(*1章3節 付図1の計算例，式(1-20) etc.  *)

(*はじめに柱先端の角度gammaを与えたときのたわみ計算を行う*)

(*gamma*)

gm=20

(* modulus : p *)

p=N[Sin[gm*Pi/2/180]]

(* 無次元荷重 Q=k*L :Nondimensional load *)

Q=EllipticK[p*p]

(* たわみのリストの準備 *)

Print["   phi       x/L         y/L "];
deflst={{0,0}};
Do[ph1=N[phi*Pi/180];
y=2*p/Q*Cos[ph1];
x=2/Q*EllipticE[ph1,p*p]-1/Q*EllipticF[ph1,p*p];
gz=2*p/EllipticK[p*p]-y;
eta=x;
deflst=Append[deflst,{N[gz],N[eta]}];
Print[N[ph1]," , ",gz," , ",eta]
,{phi,2,88,4}]
ListLinePlot[deflst,Frame->True]

(* 以上の計算に基づき 先端たわみ角 γ を20,40,60,80,100,120 と順次変えてたわみ形状を計算*)

Clear["Global`*"];
Do[gm=N[20*ig];

p=N[Sin[gm*Pi/2/180]];
Q=EllipticK[p*p];(* Q=k*L*)
Print[];
Print["gamma=",ig*20," ( P=",4*Q^2/Pi^2," P_cr )"];
Print["    phi        (x/L ,   y/L)"];deflst[ig]={{0,0}};

Do[ph1=N[phi*Pi/180];
y=2*p/Q*Cos[ph1];
x=2/Q*EllipticE[ph1,p*p]-1/Q*EllipticF[ph1,p*p];
gz=2*p/EllipticK[p*p]-y;
eta=x;
deflst[ig]=Append[N[deflst[ig],4],{N[gz,4],N[eta,4]}];
Print[PaddedForm[phi,{4,1}]," , ",PaddedForm[gz,{7,5}]," , ",PaddedForm[eta,{7,5}]]
,{phi,2,90,4}]
,{ig,1,6}];
ListLinePlot[{deflst[1],deflst[2],deflst[3],deflst[4],deflst[5],deflst[6]},Frame->True,GridLines->Automatic,AspectRatio->1,FrameLabel->{"y/L","x/L"}]