function [ux,uy,uz,uxx,uxy,uyx,uyy,uzx,uzy]...
  =tensile(xc,yc,d,L,W,dip,strike,opening,xs,ys)
%入力
%xc,yc,d：割れ目上端の中心のxy座標と深さ
%L,W：割れ目の長さと幅
%dip：開口割れ目の傾斜角（水平から時計回り）
%strike：開口割れ目の走向（北から時計回り）
%opening：開口量
%xs,ys：観測点のxy座標
%出力
%ux,uy,uz：変位のx,y,z成分
%uxx,uxy：uxのx,y微分量
%uyx,uyy：uyのx,y微分量
%uzx,uzy：uzのx,y微分量
nu=0.25;%ポアソン比
cd=cosd(dip);
sd=sind(dip);
cs=cosd(strike);
ss=sind(strike);
k=1-2*nu;
u=opening/(2*pi);
x=ss.*(xs-xc)+cs.*(ys-yc);
y=-cs.*(xs-xc)+ss.*(ys-yc);
p=y.*cd+d.*sd;
q=y.*sd-d.*cd;
[ux1,uy1,uz1,uxx1,uxy1,uyx1,uyy1,uzx1,uzy1]=ui(x+L/2,p+W,cd,sd,k,q);
[ux2,uy2,uz2,uxx2,uxy2,uyx2,uyy2,uzx2,uzy2]=ui(x+L/2,p,cd,sd,k,q);
[ux3,uy3,uz3,uxx3,uxy3,uyx3,uyy3,uzx3,uzy3]=ui(x-L/2,p+W,cd,sd,k,q);
[ux4,uy4,uz4,uxx4,uxy4,uyx4,uyy4,uzx4,uzy4]=ui(x-L/2,p,cd,sd,k,q);
ux0=u.*(ux1-ux2-ux3+ux4);
uy0=u.*(uy1-uy2-uy3+uy4);
uz0=u.*(uz1-uz2-uz3+uz4);
uxx0=-u.*(uxx1-uxx2-uxx3+uxx4);
uxy0=-u.*(uxy1-uxy2-uxy3+uxy4);
uyx0=-u.*(uyx1-uyx2-uyx3+uyx4);
uyy0=-u.*(uyy1-uyy2-uyy3+uyy4);
uzx0=-u.*(uzx1-uzx2-uzx3+uzx4);
uzy0=-u.*(uzy1-uzy2-uzy3+uzy4);
ux=ux0.*ss-uy0.*cs;
uy=ux0.*cs+uy0.*ss;
uz=uz0;
uxx=(uxx0-uyy0)/2*(ss.^2-cs.^2)+(uxx0+uyy0)/2+(uxy0+uyx0).*cs.*ss;
uxy=-(uxx0-uyy0).*cs.*ss+(uxy0-uyx0)/2+(uxy0+uyx0)/2*(ss.^2-cs.^2);
uyx=(uxx0-uyy0).*cs.*ss+(uxy0-uyx0)/2-(uxy0+uyx0)/2*(ss.^2-cs.^2);
uyy=-(uxx0-uyy0)/2*(ss.^2-cs.^2)+(uxx0+uyy0)/2-(uxy0+uyx0).*cs.*ss;
uzx=uzx0.*ss-uzy0.*cs;
uzy=uzx0.*cs+uzy0.*ss;
end

function [I1,I2,I3,I4,I5]=In(xi,eta,cd,sd,k,q)
ytil=eta.*cd+q.*sd;
dtil=eta.*sd-q.*cd;
R=sqrt(xi.^2+eta.^2+q.^2);
X=sqrt(xi.^2+q.^2);
if cd==0
  I5=-k.*xi.*sd/(R+dtil);
  I4=-k.*q/(R+dtil);
  if R+eta==0
    I3=k/2*(eta/(R+dtil)+ytil.*q/(R+dtil).^2+log(R-eta));
  else
    I3=k/2*(eta/(R+dtil)+ytil.*q/(R+dtil).^2-log(R+eta));
  end
  I1=-k/2*xi.*q/(R+dtil).^2;
else
  if xi==0
    I5=0;
  else
    I5=k.*2/cd.*atan((eta.*(X+q.*cd)+X.*(R+X).*sd)/(xi.*(R+X).*cd));
  end
  if R+eta==0
    I4=k/cd.*(log(R+dtil)+sd.*log(R-eta));
    I3=k.*(1/cd.*ytil/(R+dtil)+log(R-eta))+sd/cd.*I4;
  else
    I4=k/cd.*(log(R+dtil)-sd.*log(R+eta));
    I3=k.*(1/cd.*ytil/(R+dtil)-log(R+eta))+sd/cd.*I4;
  end
  I1=k.*(-1/cd.*xi/(R+dtil))-sd/cd.*I5;
end
if R+eta==0
  I2=k.*log(R-eta)-I3;
else
  I2=-k.*log(R+eta)-I3;
end
end

function [J1,J2,J3,J4]=Jn(xi,eta,cd,sd,k,q)
ytil=eta.*cd+q.*sd;
dtil=eta.*sd-q.*cd;
R=sqrt(xi.^2+eta.^2+q.^2);
X=sqrt(xi.^2+q.^2);
[K1,K2,K3]=Kn(xi,eta,cd,sd,k,q);
if cd==0
  J1=k/2*q/(R+dtil).^2*(2*xi.^2/(R.*(R+dtil))-1);
  J2=k/2*xi.*sd/(R+dtil).^2*(2*q.^2/(R.*(R+dtil))-1);
else
  J1=k/cd.*(xi.^2/(R.*(R+dtil).^2)-1/(R+dtil))-sd/cd.*K3;
  J2=k/cd.*(xi.*ytil/(R.*(R+dtil).^2))-sd/cd.*K1;
end
if R+eta==0
  J3=-J2;
  J4=-k.*cd/R-J1;
else
  J3=k.*(-xi/(R.*(R+eta)))-J2;
  J4=k.*(-cd/R-q.*sd/(R.*(R+eta)))-J1;
end
end

function [K1,K2,K3]=Kn(xi,eta,cd,sd,k,q)
ytil=eta.*cd+q.*sd;
dtil=eta.*sd-q.*cd;
R=sqrt(xi.^2+eta.^2+q.^2);
X=sqrt(xi.^2+q.^2);
if cd==0
  K1=k.*xi.*q/(R.*(R+dtil).^2);
  K3=k.*sd/(R+dtil).*(xi.^2/(R.*(R+dtil))-1);
else
  if R+eta==0
    K1=k.*xi/cd/(R.*(R+dtil));
    K3=-k/cd.*ytil/(R.*(R+dtil));
  else
    K1=k.*xi/cd.*(1/(R.*(R+dtil))-sd/(R.*(R+eta)));
    K3=k/cd.*(q/(R.*(R+eta))-ytil/(R.*(R+dtil)));
  end
end
if R+eta==0
  K2=-k.*sd/R-K3;
else
  K2=k.*(-sd/R+q.*cd/(R.*(R+eta)))-K3;
end
end

function [uxi,uyi,uzi,uxxi,uxyi,uyxi,uyyi,uzxi,uzyi]=ui(xi,eta,cd,sd,k,q)
ytil=eta.*cd+q.*sd;
dtil=eta.*sd-q.*cd;
R=sqrt(xi.^2+eta.^2+q.^2);
X=sqrt(xi.^2+q.^2);
[I1,I2,I3,I4,I5]=In(xi,eta,cd,sd,k,q);
[J1,J2,J3,J4]=Jn(xi,eta,cd,sd,k,q);
[K1,K2,K3]=Kn(xi,eta,cd,sd,k,q);
Axi=(2*R+xi)/(R.^3*(R+xi).^2);
if R+eta==0
  Aeta=0;
else
  Aeta=(2*R+eta)/(R.^3*(R+eta).^2);
end
if R+eta==0
  uxi=-I3.*sd.^2;
else
  uxi=q.^2/(R.*(R+eta))-I3.*sd.^2;
end
if q==0
  if R+eta==0
    uyi=-dtil.*q/(R.*(R+xi))-I1.*sd.^2;
  else
    uyi=-dtil.*q/(R.*(R+xi))-sd.*xi.*q/(R.*(R+eta))-I1.*sd.^2;
  end
else
  if R+eta==0
    uyi=-dtil.*q/(R.*(R+xi))+sd.*atan(xi.*eta/(q.*R))-I1.*sd.^2;
  else
    uyi=-dtil.*q/(R.*(R+xi))-sd.*(xi.*q/(R.*(R+eta))-atan(xi.*eta/(q.*R)))-I1.*sd.^2;
  end
end
if q==0
  if R+eta==0
    uzi=ytil.*q/(R.*(R+xi))-I5.*sd.^2;
  else
    uzi=ytil.*q/(R.*(R+xi))+cd.*xi.*q/(R.*(R+eta))-I5.*sd.^2;
  end
else
  if R+eta==0
    uzi=ytil.*q/(R.*(R+xi))-cd.*atan(xi.*eta/(q.*R))-I5.*sd.^2;
  else
    uzi=ytil.*q/(R.*(R+xi))+cd.*(xi.*q/(R.*(R+eta))-atan(xi.*eta/(q.*R)))-I5.*sd.^2;
  end
end
uxxi=xi.*q.^2*Aeta+J3.*sd.^2;
uxyi=-dtil.*q/R.^3-xi.^2*q.*Aeta.*sd+J1.*sd.^2;
uyxi=q.^2/R.^3*cd+q.^3*Aeta.*sd+J1.*sd.^2;
uyyi=(ytil.*cd-dtil.*sd).*q.^2*Axi-2*q.*cd.*sd/(R.*(R+xi))-(xi.*q.^2*Aeta-J2).*sd.^2;
uzxi=q.^2/R.^3*sd-q.^3*Aeta.*cd+K3.*sd.^2;
uzyi=(ytil.*sd+dtil.*cd).*q.^2*Axi+xi.*q.^2*Aeta.*sd.*cd-(2*q/(R.*(R+xi))-K1).*sd.^2;
end
