(*片持ちはり大たわみ計算 Fig2-3, Fig2-4 *)

(*楕円積分を利用*)

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

EllipticK[0.25]

(*楕円関数K を数値積分によって評価した結果*)

Do[
p=N[i]/10;p2=p*p;
ans1=EllipticK[p2]-NIntegrate[1/(1-p*p*(Sin[th])^2)^(1/2),{th,0,Pi/2}];
Print[p," , ",ans1],{i,1,9}]

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

EllipticK[0.25]

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

(*なお，Mathematica の楕円積分の引数は，たとえば，EllipticF[phi,p2] となっていることに注意（角度が先で，pがあと）*)

Print["psi0     p       q^2      delta/L   (L-Delta)/L)    sigma0" ];
Do[
psi0=N[psi00*Pi/180];p=((1+Sin[psi0])/2)^(1/2);p2=p*p;
phi1=ArcSin[1/Sqrt[2]/p];
q=Re[EllipticK[p2]-EllipticF[phi1,p2]];
deltaL=Re[1/q*(EllipticK[p2]-EllipticF[phi1,p2]
-2*EllipticE[p2]
+2*EllipticE[phi1,p2])];
LdeltaL=Sqrt[2]/q*(Sin[psi0])^(1/2);
sigma0=Sqrt[2]*q*(Sin[psi0])^(1/2);
Print[psi00," , ",p," , ",q*q," , ",deltaL," , ",LdeltaL," , ",sigma0],{psi00,0,85,2.5}]

(*楕円関数を利用しないで直接，数値積分する．=> 数値積分結果は楕円関数解とほぼ等しい．*)

Print[" psi0    q^2   delta/L  (L-delta)/L"];
tabdeltaL={};tabLdL={};tabLinear={{0,0},{1,3}};
Do[
psi0=N[psideg*Pi/180];
If[psi0==0,
    q=0;deltaL=0;LdeltaL=1;,
q=1/Sqrt[2]*NIntegrate[1/(Sin[psi0]-Sin[psi])^(1/2),
{psi,0,psi0}];
deltaL=1/(Sqrt[2]*q)*NIntegrate[Sin[psi]/(Sin[psi0]-Sin[psi])^(1/2),
{psi,0,psi0}];
LdeltaL=Sqrt[2]/q*Sqrt[Sin[psi0]]];
Print[psideg," , ",q*q," , ",deltaL," , ",LdeltaL];
(*結果のテーブルへの追加*)
tabdeltaL=Append[tabdeltaL,{deltaL,q*q}];
tabLdL=Append[tabLdL,{LdeltaL,q*q}];
,{psideg,0,85,2.5}]
(*繰り返し終了*)
(*結果のプロット*)
ListPlot[{tabdeltaL,tabLdL,tabLinear},Frame->True,GridLines->Automatic,Joined->True,PlotStyle->{Thick},FrameLabel->{ "Deflections",Style["(k^2)(l^2)",FontSize->20] },LabelStyle->FontSize->15]

TtabdeltaL=Table[{tabdeltaL[[i]][[2]],tabdeltaL[[i]][[1]]},{i,Length[tabdeltaL]}];
TtabLdL=Table[{tabLdL[[i]][[2]],tabLdL[[i]][[1]]},{i,Length[tabLdL]}];TtabLinear=Table[{tabLinear[[i]][[2]],tabLinear[[i]][[1]]},{i,Length[tabLinear]}];ListPlot[{TtabdeltaL,TtabLdL,TtabLinear},Frame->True,GridLines->Automatic,Joined->True,PlotStyle->{Thick},FrameLabel->{Style["Load (k^2)(l^2)",FontSize->15],"Deflections"} ,LabelStyle->FontSize->15]

(*変形座標(x,y)の計算*)

Print[" psi0     x        y"];
grtable{};
Do[
tabxy={};
psi0=N[psideg*Pi/180];
If[psi0==0,
 q=0;,
q=1/Sqrt[2]*NIntegrate[1/(Sin[psi0]-Sin[psi])^(1/2),
{psi,0,psi0}]];

Do[
xx=1/(Sqrt[2]*q)*NIntegrate[Cos[psi]/(Sin[psi0]-Sin[psi])^(1/2),
{psi,0,psi1}];
yy=-1/(Sqrt[2]*q)*NIntegrate[Sin[psi]/(Sin[psi0]-Sin[psi])^(1/2),
{psi,0,psi1}];
(*結果のテーブルへの追加*)
tabxy=Append[tabxy,{xx,yy}];
,{psi1,0,psi0,psi0/20}];

graph[psideg]=ListPlot[tabxy,Frame->True,GridLines->Automatic,Joined->True,PlotStyle->{Thick},FrameLabel->{"x/L","y/L"},PlotLabel->"psi0="<>ToString[psideg]];
(*grtable=Append[grtable,graph[psideg]];*)
,{psideg,10,85,10}]
(*繰り返し終了*)
(*結果のプロット*)

Show[{graph[60],graph[50],graph[40],graph[30],graph[20],graph[10]},AspectRatio->1,PlotRange->{{0,1},{0,-0.8}},PlotLabel->"psi=10 to 60"]

g2=Animate[Show[graph[i],AspectRatio->1,PlotRange->{{0,1},{0,-0.8}}],{i,10,85,10}]

Export["fig2-3-animation.avi",g2]

(*無次元曲げ応力sigma_B の計算, 長さ方向分布計算*)

Print[" s/L    Dimensionless Stress"];
grtableM{};
Do[
tabxyM={};
psi0=N[psideg*Pi/180];(*Degree -> Rad*)
q=1/Sqrt[2]*NIntegrate[1/(Sin[psi0]-Sin[psi])^(1/2),
{psi,0,psi0}];
Print["q=", q, " , ","(psi_0=",psideg,")"];
Do[
ss=1/(Sqrt[2]*q)*NIntegrate[1/(Sin[psi0]-Sin[psi])^(1/2),
{psi,0,psi1}];
sigmaB=2^(1/2)*q*(Sin[psi0]-Sin[psi1])^(1/2);
(*結果のテーブルへの追加*)
tabxyM=Append[tabxyM,{ss,sigmaB}];Print[ss," , ",sigmaB];
,{psi1,0,psi0,psi0/10}];

graphM[psideg]=ListPlot[tabxyM,Frame->True,GridLines->Automatic,Joined->True,PlotStyle->{Thick},FrameLabel->{"s/L","sigma_B"},PlotLabel->"psi0="<>ToString[psideg]];
(*grtable=Append[grtable,graph[psideg]];*)
,{psideg,20,85,20}]
(*繰り返し終了*)
(*結果のプロット*)
Show[{graphM[60],graphM[40],graphM[20]},AspectRatio->1,PlotRange->Automatic]