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

  fmn= [0.1,0.02926,0.03007,0.03058  %provided data
        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];

  m=a/z; %backslash operator to solve m
  n=b/z;    %backslash for n
 
  %------------- Building indices ---------------------------------------
  %Two loops designed to isolate the idices in the matrix that will be
  %interpolated between in the future
  mindex = [];
  b = sort(abs(fmn(:,1)-m));
  close2m = b(1:4);
  %Build indices for m values
  for j=1:4
      for i=1:8
          if close2m(j) == abs(fmn(i,1)-m)
              mindex=[mindex;i];
          end
      end
  end
  %Build indices for desired n values
  if n<1.2
      n=2;
  elseif n>1.6
      n=4;
  else
      n=int8((round(n*5)/5-1)/0.2+1);
  end
  %----------Value organization --------------------------------------
  %reference values of the previously found indices for n & m
  c=zeros(4,1);
  close2mtrue=zeros(4,1);
  for i=1:4 %organize values from the matrix
      close2mtrue(i) = fmn(mindex(i),1);
      c(i) = fmn(mindex(i),n);
  end
  
  sigma_z = q*ppval(interp1(close2mtrue,c,'cubic','pp'),m); %interpolation between the closest n & m values

end