function main 

%*********************參數說明 
%L為樁長，r為計算點與樁軸線的水平距離，z為計算點距離地面深度。都取正值 
%f1,f2,f3依次是點荷載、均布荷載以及三角形荷載下的Geddes應力解函數 
%Ip,Ir,It依次是對應三個荷載下的應力系數，負表示拉，正表示壓 
%參考文獻:J.D.Geddes Stresses in Foundation Soils Due to Vertical Subsurface Loading, 1966 
%計算結果于Geddes的附加應力表是一致的 
clc 
clear 
L=10; 
r=2; 
v=0.5; 
z=3; 

disp(’Ip, It, Ir’) 

[Ip,Ir,It]=Geddes(L,r,v,z) 

disp(’***************’) 

function [Ip,Ir,It]=Geddes(L,r,v,z) 

m=z/L; 
n=r/L; 
disp(’m and n’) 
[m n] 

if n==0 
disp(’計算點位于樁軸線上，采用n=0.002進行修正，并模擬刺入沉降!’) 
n=0.002; 
end 

Ip=f1(m,n,v); 
Ir=f2(m,n,v); 
It=f3(m,n,v); 

function Ip_sub=f1(m,n,v) %點荷載 

F=sqrt(m^2+n^2); 
A=sqrt(n^2+(m-1)^2); 
B=sqrt(n^2+(m+1)^2); 
temp0=1/(8*pi*(1-v)); 
temp1=(1-2*v)*(m-1)/A^3-(1-2*v)*(m-1)/B^3+3*(m-1)^3/A^5; 
temp2=(3*(3-4*v)*m*(m+1)^2-3*(m+1)*(5*m-1))/B^5; 
temp3=30*m*(m+1)^3/B^7; 
Ip_sub=-temp0*(temp1+temp2+temp3); 

function Ir_sub=f2(m,n,v) %矩形荷載 
F=sqrt(m^2+n^2); 
A=sqrt(n^2+(m-1)^2); 
B=sqrt(n^2+(m+1)^2); 

temp0=1/(8*pi*(1-v)) 
temp1=-2*(2-v)/A; 
temp2=(2*(2-v)+2*(1-2*v)*(1+m)*m/n^2)/B; 
temp3=-2*(1-2*v)*m^2/n^2/F+n^2/A^3; 
temp4=(4*m^2-4*(1+v)*(m/n)^2*m^2)/F^3; 
temp5=(4*m*(1+v)*(m+1)*(m/n+1/n)^2-(4*m^2+n^2))/B^3; 
temp6=6*m^2*((m^4-n^4)/n^2)/F^5; 
temp7=6*m*(m*n^2-(m+1)^5/n^2)/B^5; 
Ir_sub=temp0*(temp1+temp2+temp3+temp4+temp5+temp6+temp7); 

function It_sub=f3(m,n,v) %三角形荷載 
F=sqrt(m^2+n^2); 
A=sqrt(n^2+(m-1)^2); 
B=sqrt(n^2+(m+1)^2); 
temp0=1/(4*pi*(1-v)); 
temp1=-2*(2-v)/A; 
temp2=(2*(2-v)*(4*m+1)-2*(1-2*v)*(1+m)*m^2/n^2)/B; 
temp3=(2*(1-2*v)*m^3/n^2-8*(2-v)*m)/F+(m*n^2+(m-1)^3)/A^3; 
temp4=(4*v*m*n^2+4*m^3-15*n^2*m-2*(5+2*v)*(m/n)^2*(m+1)^3+(m+1)^3)/B^3; 
temp5=(2*(7-2*v)*m*n^2-6*m^3+2*(5+2*v)*(m/n)^2*m^3)/F^3; 
temp6=(6*m*n^2*(n^2-m^2)+12*(m/n)^2*(m+1)^5)/B^5; 
temp7=-(12*(m/n)^2*m^5+6*m*n^2*(n^2-m^2))/F^5; 
temp100=((A+m-1)/(F+m))*((B+m+1)/(F+m)); 
temp8=-2*(2-v)*log(temp100); 
It_sub=temp0*(temp1+temp2+temp3+temp4+temp5+temp6+temp7+temp8);