function [ux,uy,uz]=fss(xc,yc,d,a,dV,xs,ys)
%%球状変動源モデル（McTigue 1987）による変動を計算するMATLABコード

%%入力
%xc,yc,d：変動源のxy座標と深さ
%a：変動源の半径
%dV：変動源の体積変化
%xs,ys：観測点のxy座標

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

%%参考文献
%McTigue DF (1987) Elastic stress and deformation near a finite spherical magma body: Resolution of the point source paradox. JGR, 92:12931-12940. https://doi.org/10.1029/JB092iB12p12931
%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;%ポアソン比
dx=xs-xc;
dy=ys-yc;
[ux,uy,uz]=disp(dx,dy,d,dV,a,nu);
end

function [ux,uy,uz]=disp(dx,dy,d,dV,a,nu)
K=dV/(pi*a^3*(1+(a/d)^2));
r=(dx^2+dy^2)^0.5;
ur=(1-nu)*(K*a^3)*(r/(r^2+d^2)^1.5)*(1-(a/d)^3*((1+nu)/(2*(7-5*nu)) ...
  -(15*(2-nu))/(4*(7-5*nu))*d^2/(r^2+d^2)));
ux=ur*dx/(r+eps);
uy=ur*dy/(r+eps);
end
uz=(1-nu)*(K*a^3)*(d/(r^2+d^2)^1.5)*(1-(a/d)^3*((1+nu)/(2*(7-5*nu)) ...
  -(15*(2-nu))/(4*(7-5*nu))*d^2/(d^2+r^2)));
end