boussinesq_lookup.m

function sigma_z=boussinesq_lookup(q,a,b,z)
  % function that determines stress under corner of an a by b rectangular platform
  % z-meters below the platform. The calculated solutions are in the fmn data
  % m=fmn(:,1)
  % in column 2, fmn(:,2), n=1.2
  % in column 3, fmn(:,2), n=1.4
  % in column 4, fmn(:,2), n=1.6
  
% the data table used to interpolate
  fmn= [0.1,0.02926,0.03007,0.03058
        0.2,0.05733,0.05894,0.05994
        0.3,0.08323,0.08561,0.08709
        0.4,0.10631,0.10941,0.11135
        0.5,0.12626,0.13003,0.13241
        0.6,0.14309,0.14749,0.15027
        0.7,0.15703,0.16199,0.16515
        0.8,0.16843,0.17389,0.17739];
 
  %fins the m and n values
  m=a/z;
  n=b/z;
 
  %the nth order plus 1
  n1=4;

  %find the values near the m and n values
  k=find(fmn(:,1)>=m);
  j=find(fmn(:,1)<=m);
  k=min(k);
  j=max(j);
  

  %sets the four rows the table will be interpolated from 
  f(1)=fmn(j-1,1);
  f(2)=fmn(j,1);
  f(3)=fmn(k,1);
  f(4)=fmn(k+1,1);
  

  %a series of if statements to determine what column to interpolate from
  %for the n value
  if n<1.2 
      c=2;
  end
  if n>1.2 && n<1.3
      c=2;
  end
  if n>=1.3 && n<1.5
      c=3;
  else 
      c=4;
  end
  
  %this if statement is for if an exact value in the table is selected it
  %calculates the stress, and ends the function
  if j==k
      sigma_z=q*fmn(j,c);
        return  
  end
  
  b = zeros(n1,n1);
% assign dependent variables to the first column of b.

b(:,1) = fmn(j-1:k+1,c); % the (:) ensures that y is a column vector.

%uses parts of the Newtint.m function to calculate the interpolated value
for j = 2:n1
  for i = 1:n1-j+1
    b(i,j) = (b(i+1,j-1)-b(i,j-1))/(f(i+j-1)-f(i));
  end
end

% use the finite divided differences to interpolate
  xt = 1;
yint = b(1,1);
for j = 1:n1-1
  xt = xt*(m-f(j));
  yint = yint+b(1,j+1)*xt;
end
  sigma_z=yint*q;
  
end