(*3.1 節, 図3.3作成*)

(*式(3.5)より y_max/l を求める．Gospodnetic論文に基づく計算*)

(*荷重の大きさPl^2/EI => q2 は，式（3．7）で評価する*)

Clear["Global`*"]

Print["psi[deg]   y_max/l　　(Pl^2/EI) "]
lst1={};lst2={};p2=(1/Sqrt[2])^2;
Do[
psi0=Pi/40*i;phi=ArcCos[Sqrt[Sin[psi0]]];
(*psi0=0 のときには，不定形が現れるので，それを回避する*)
If [i<=0 ,{psi0=0.0001;phi=ArcCos[Sqrt[Sin[psi0]]]}];
fn=N[Re[(-EllipticK[p2]+EllipticF[phi,p2]
+2*EllipticE[p2]
-2*EllipticE[phi,p2])]];
ymaxl=(Sqrt[2]*Sin[psi0]*Cos[phi]-Cos[psi0]*fn)/(Sqrt[2]*Cos[psi0]*Cos[phi]+Sin[psi0]*fn);
q2=Cos[psi0]*(Sqrt[2]*Cos[psi0]*Cos[phi]+Sin[psi0]*fn)^2;
psi0deg=N[psi0*180/Pi];
lst1=Append[lst1,{N[psi0deg],ymaxl}];
lst2=Append[lst2,{N[psi0deg],q2}];
Print[N[psi0deg]," , ",PaddedForm[ymaxl,5]," , ",PaddedForm[q2,5]],{i,0,20}]

ListLinePlot[{lst1,lst2},Frame->True,GridLines->Automatic,Joined->True,FrameLabel->{"Psi","(Subscript[y, max])/l, (PL^2)/EI"}]

(*(Pl^2/EI)の最大値付近 と y_max/l *)

Print["psi[deg] psi[rad]  (Pl^2/EI)     y_max/l"]
Do[
psi0=Pi/180*i;
phi=ArcCos[Sqrt[Sin[psi0]]];
fn=N[Re[(-EllipticK[p2]+EllipticF[phi,p2]
+2*EllipticE[p2]
-2*EllipticE[phi,p2])]];
ymaxl=(Sqrt[2]*Sin[psi0]*Cos[phi]-Cos[psi0]*fn)/(Sqrt[2]*Cos[psi0]*Cos[phi]+Sin[psi0]*fn);
q2=Cos[psi0]*(Sqrt[2]*Cos[psi0]*Cos[phi]+Sin[psi0]*fn)^2;
Print[i, " , ",N[psi0]," , ",PaddedForm[q2,9]," , ", PaddedForm[ymaxl,7]],{i,38.2,38.4,0.01}]

(*Psi0=38.3 deg 付近でy_max/l が最大値*)