(*1章3節， h_maxの誘導*)

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

(*EllipticK[p*p]]<-->E[p](Fay)とすればよい*)

(*第１種完全楕円積分EllipticK[p*p]，第１種不完全楕円積分EllipticF[]*)

Clear["Global`*"]

(*はじめに，K(p)-B(p)p^2/(1-p^2)の符号反転箇所を調べ，p の存在範囲を探る．*)

Do[
p2=p*p;
B=NIntegrate[(Cos[phi])^2/(1-p*p*(Sin[phi])^2)^0.5,{phi,0,Pi/2}];
fn=EllipticK[p2]-B*p*p/(1-p*p);
Print["p= ",p," , fn=",fn]
,{p,0.5,0.95,0.05}]

(*これで，p=0.8, 0.85 付近に根が存在することがわかった．次は，FindRoot を用いた求根．*)

ans1=FindRoot[EllipticK[p*p]-p*p/(1-p*p)*Integrate[(Cos[phi])^2/(1-p*p*(Sin[phi])^2)^0.5,{phi,0,Pi/2}],{p,0.80}]

fn/.ans1

(*水平の最大たわみ hmax/L*)

hmax=(2*p/EllipticK[p*p])/.ans1

(* 荷重 P= 4P_cr*K^2(p*p)/(Pi)^2 *)

P=4*(((EllipticK[p*p])^2)/Pi^2)/.ans1