function [ux,uy,uz]=mogi(xc,yc,d,dV,xs,ys)
%%茂木モデル（Mogi 1958）による変動を計算するMATLABコード

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

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

%%参考文献
%Mogi K (1958) Relations between the Eruptions of Various Volcanoes and the Deformations of the Ground Surfaces around them. BERI, 36, 99-134. https://doi.org/10.15083/0000033924
%Delaney PT & McTigue DF (1994) Volume of magma accumulation or withdrawal estimated from surface uplift or subsidence, with application to the 1960 collapse of Kilauea volcano. Bull. Volcanol., 56, 417-424. https://doi.org/10.1007/BF00302823

k=3*dV/(4*pi);
r=((xc-xs)^2+(yc-ys)^2)^0.5;
ur=k*r/(d^2+r^2)^1.5;
uz=k*d/(d^2+r^2)^1.5;
ux=ur*(xs-xc)/(r+eps);
uy=ur*(ys-yc)/(r+eps);
end
