(*2章2節の例題計算*)

(*Mathematica の楕円積分の引数の定義がFayの本と違っている*)

(*EllipticK[p*p]]<-->E[p](Fay)とすればよい.また，Mathematica の楕円積分の引数は，たとえば，EllipticF[phi,p2] となっていることに注意（角度が先で，p2があと）*)

Clear["Global`*"];
L=100;P=1;EI=1000;
q=N[Sqrt[P*L^2/EI]];q

(*式（2.16）により q を与えたときのpを求める．*)

fn=q-EllipticK[p*p]+EllipticF[ArcSin[1/Sqrt[2]/p],p*p];
pans=FindRoot[fn==0,{p,0.8}][[1]][[2]]

(*複素数根が求められるが，虚数部が小さいので無視する*)

pans1=Re[pans]

(*このpを用いてphi=ArcSin(1/(Sqrt[2]p))よりphiを求める*)

phi1=ArcSin[1/Sqrt[2]/pans1]

(*式（2.18）により，はり先端の垂直たわみδ[in.]を求める*)

sdelta=L/q*(EllipticK[pans1*pans1]-EllipticF[phi1,pans1*pans1]-2*EllipticE[pans1*pans1]+2*EllipticE[phi1,pans1*pans1])

(*式（2.20）より，はり先端の水平たわみΔ [in.]を求める*)

Ldelta=L-Sqrt[2*EI*(2*pans1^2-1)/P]

(*Sin[psi_0]=2*p*p-1より，先端たわみ角 psi_0 を求める*)

psi0=ArcSin[2*pans1*pans1-1]

(*degree 計算*)

psi0*180/3.141595