function [ux,uy,uz]=pss(xc,yc,d,a,b,dip,strike,dP,xs,ys)
%%長球状変動源モデル（Yang et al. 1988）による変動を計算するMATLABコード

%%入力
%xc,yc,d：変動源の中心のxy座標と深さ
%a,b：変動源の長半径と短半径
%dip：変動源の傾斜角（水平から時計回り）
%strike：変動源の走向（北から時計回り）
%dP：変動源の内部圧力変化
%xs,ys：観測点のxy座標

%%出力
%ux,uy,uz：変位のx,y,z成分

%%参考文献
%Yang X, Davis PM, Dieterich JH (1988) Deformation from inflation of a dipping finite prolate spheroid in an elastic half-space as a model for volcanic stressing. JGR, 93:4249-4257. https://doi.org/10.1029/JB093iB05p04249
%Battaglia M, Cervelli PF, Murray JR (2013) Modeling crustal deformation near active faults and volcanic centers: a catalog of deformation models and modeling approaches. Techniques and Methods 13-B1. https://doi.org/10.3133/tm13B1

nu=0.25;%ポアソン比
E=1e9;%ヤング率
zs=0;
mu=E/(2*(1+nu));
cs=cosd(strike);
ss=sind(strike);
xp=cs*(xs-xc)-ss*(ys-yc);
yp=ss*(xs-xc)+cs*(ys-yc);
c=(a^2-b^2)^0.5;
[u11,u21,u31]=primitive(mu,nu,zs,xp,yp,d,dP,dip,a,b,c,c);
[u12,u22,u32]=primitive(mu,nu,zs,xp,yp,d,dP,dip,a,b,c,-c);
upx=-u11-u12;
upy=-u21-u22;
upz=u31+u32;
ux=cs*upx+ss*upy;
uy=-ss*upx+cs*upy;
uz=upz;
end
      
function [u1,u2,u3]=primitive(mu,nu,zs,xp,yp,d,dP,dip,a,b,c,xi)
cd=cosd(dip);
sd=sind(dip);
xi2=xi*cd;
xi3=xi*sd;
x1=xp;
x2=yp;
x3=zs-d;
x3b=zs+d;
y1=x1;
y2=x2-xi2;
y3=x3-xi3;
y3b=x3b+xi3;
r2=x2*sd-x3*cd;
q2=x2*sd+x3b*cd;
r3=x2*cd+x3*sd;
q3=-x2*cd+x3b*sd;
r3b=r3-xi;
q3b=q3+xi;
rr1=(y1^2+y2^2+y3^2)^0.5;
rr2=(y1^2+y2^2+y3b^2)^0.5;
c0=d/(sd+eps);
beta=(q2*cd+(1+sd)*(rr2+q3b))/(cd*y1+eps);
dr3b=rr1+r3b;
dq3b=rr2+q3b;
dy3b=rr2+y3b;
lr3b=log(rr1+r3b);
lq3b=log(rr2+q3b);
ly3b=log(rr2+y3b);
[P1,a1,b1]=general(mu,nu,dP,a,b,c);
aa1s=a1/(rr1*dr3b)+b1*(lr3b+(r3+xi)/dr3b);
aa1sb=-a1/(rr2*dq3b)-b1*(lq3b+(q3-xi)/dq3b);
aa1=xi/rr1+lr3b;
aa1b=xi/rr2-lq3b;
aa2=rr1-r3*lr3b;
aa2b=rr2-q3*lq3b;
aa3=xi*r3b/rr1+rr1;
aa3b=xi*q3b/rr2-rr2;
bbs=(a1/rr1+2*b1*aa2)+(3-4*nu)*(a1/rr2+2*b1*aa2b);
bb=xi*(xi+c0)/rr2-aa2b-c0*lq3b;
ff1s=0;
ff2s=0;
ff1=0;
ff2=0;
f1=xi*y1/dy3b+(3/(cd+eps)^2)*(y1*sd*ly3b-y1*lq3b+2*q2*atan2(beta))+2*y1*lq3b-4*x3b*atan2(beta)/cd;
f2=xi*y2/dy3b+(3/(cd+eps)^2)*(q2*sd*lq3b-q2*ly3b+2*y1*sd*atan2(beta)+cd*(rr2-y3b))-2*cd*aa2b+(2/(cd+eps))*(x3b*ly3b-q3*lq3b);
f3=(1/(cd+eps))*(q2*lq3b-q2*sd*ly3b+2*y1*atan2(beta))+2*sd*aa2b+q3*ly3b-xi;
c1=(a*b^2/(xi+eps)^3)/(16*mu*(1-nu));
c2=2*c1*P1;
uu11=c1*(aa1s*y1+(3-4*nu)*aa1sb*y1+ff1s*y1);
uu12=c1*(sd*(aa1s*r2+(3-4*nu)*aa1sb*q2+ff1s*q2)+cd*(bbs-ff2s));
uu13=c1*(-cd*(aa1s*r2+(3-4*nu)*aa1sb*q2-ff1s*q2)+sd*(bbs+ff2s));
uu21=c2*((aa1*y1+(3-4*nu)*aa1b*y1+ff1*y1)-4*(1-nu)*(1-2*nu)*f1);
uu22=c2*(sd*(aa1*r2+(3-4*nu)*aa1b*q2+ff1*q2)-4*(1-nu)*(1-2*nu)*f2+4*(1-nu)*cd*(aa2+aa2b)+cd*(aa3-(3-4*nu)*aa3b-ff2));
uu23=c2*(cd*(-aa1*r2+(3-4*nu)*aa1b*q2+ff1*q2)+4*(1-nu)*(1-2*nu)*f3+4*(1-nu)*sd*(aa2+aa2b)+sd*(aa3+(3-4*nu)*aa3b+ff2-2*(3-4*nu)*bb));
u1=uu11+uu21;
u2=uu12+uu22;
u3=uu13+uu23;
end
      
function [P1,a1,b1]=general(mu,nu,dP,a,b,c)
lambda=2*mu*nu/(1-2*nu);
acratio=(a-c)/(a+c);
coef1=2*pi*a*b^2;
den1=8*pi*(1-nu);
qq=3/den1;
rr=(1-2*nu)/den1;
ia=-coef1*(2/(a*c^2)+log(acratio)/c^2);
iaa=-coef1*(2/(3*a^3*c^2)+2/(a*c^4)+log(acratio)/c^5);
a11=2*rr*(ia-4*pi);
a12=-2*rr*(ia+4*pi);
a21=qq*a^2*iaa+rr*ia-1;
a22=-qq*a^2*iaa-ia*(2*rr-qq);
den2=3*lambda+2*mu;
num2=3*a22-a12;
den3=a11*a22-a12*a21;
num3=a11-3*a21;      
P1=dP*(2*mu/den2)*(num2-num3)/den3;
P2=dP*(1/den2)*(num2*lambda+2*(lambda+mu)*num3)/den3;
a1=-2*b^2*P1;
b1=3*(b/c)^2*P1+2*(1-2*nu)*P2;
end