boussinesq_spline.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
        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;
  n=b/z;
%The boussinesq_lookup.m code is copied 
%to a file called boussinesq_spline.m 
%using a cubic spline to interpolate in 
%two dimensions, both m and n, that returns sigma_z.
  mv=fmn(:,1);
  n1=fmn(:,2);
  n2=fmn(:,3);
  n3=fmn(:,4);
  
  if n<1.3
      n=1.2;
      nv=n1;
  elseif 1.3<n && n<1.5
      n=1.4;
      nv=n2;
  elseif n>=1.5
      n=1.6;
      nv=n3;
  end 
  
  i=round(m,1);
  
  for i = i:0.7
    if m<=i && m>i+0.1
        i=i;
    end
    x=[i-0.1 i i+0.1 i+0.2];
    m1=i;
  end 
end

nf0=nv((m1+10)-1);
nf1=nv(m1*10);
nf2=nv((m1*10)+1);
nf3=nv((m1*10)+2);
nf=[nf0 nf1 nf2 nf3];
p=polyfit(x,nf,3); 
v=polyval(p,m);
sigma_z=q*v;
end