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(:,3), n=1.4
% in column 4, fmn(:,4), 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;

if (n>1.5)
    n=1.6;
elseif (1.3<=n)&&(n<=1.5)
    n=1.4;
elseif 1.3<=n
    n=1.2;
end

f=ones(1,4);
x=ones(1,4);
%for loop to find the values in the fmn matrix that should be used.
for i=1:4
[~,p]=min(abs(m-fmn(:,1)));
M =fmn(p,1);
fmn(p,1)=0;
if n==1.2
   t=fmn(p,2);
elseif n==1.4
    t=fmn(p,3);
elseif n==1.6
   t=fmn(p,4);
end

x(i)= M;
f(i)=t;
end

%coefficients for the polynomial are found with the b values.
b1=f(1);
b2=(f(2)-f(1))/(x(2)-x(1));
b3=(((f(3)-f(2))/(x(3)-x(2)))-((f(2)-f(1))/(x(2)-x(1)))/(x(3)-x(1)));
b4=(((f(4)-f(3))/(x(4)-x(3)))-((f(3)-f(2))/(x(3)-x(2)))-((f(2)-f(1))/(x(2)-x(1))))/(x(4)-x(1));
%the coefficients are  plugged into the third order polynoomial
%interpolation.
f3=b1+(b2*(m-x(1)))+(b3*(m-x(1))*(m-x(2)))+(b4*(m-x(1))*(m-x(2))*(m-x(3)));
%Stress in the vertical direction is equal to force*the third order
%polynomial.
sigma_z=q*f3;

end