(*p.79, Eq3-17の計算*)

Clear["Global`*"]; 
Do[
lambda=10*ilamb*Pi/180;
(*If[ilamb==0,lambda=0.01];*)
 p=Sin[Pi/4+lambda/2];
Print["lambda= ",lambda*180/Pi];Print[" psi0    Pl^2/EI"];
lstpdata[ilamb]={};
Do[
n=ArcSin[(Cos[(psi0-lambda)/2]-Sin[(psi0-lambda)/2])/(Cos[lambda/2]+Sin[lambda/2])];
Lphi=EllipticF[n,p*p]-EllipticK[p*p]+2*EllipticE[p*p]-2*EllipticE[n,p*p];
Ltht=Cos[psi0-lambda]*(Sin[psi0-lambda]*Lphi+2*p*Cos[psi0-lambda]*Cos[n])^2;
lstpdata[ilamb]=Append[lstpdata[ilamb],{psi0*180/Pi,Ltht}];
Print[" ",psi0*180/Pi," , ",Ltht]
,{psi0,0,Pi/2,0.05*Pi/2}],{ilamb,0,8,1}];
ListLinePlot[{lstpdata[0],lstpdata[1],lstpdata[2],lstpdata[3],lstpdata[4],lstpdata[5],lstpdata[6],lstpdata[7],lstpdata[8]},Frame->True,GridLines->Automatic,FrameLabel->{" (Subscript[ψ, 0])","(Pl^2)/(EI)"}]