(* Constants fixed for the 2D square-lattice Ising model *)
Δ:=15/8
γ:=7/4
β:=1/8
ℳ0:=((((2^(5/2))ArcSinh[1]))^β)
Tc:=2/Log[1+Sqrt[2]]
B:=(Tc^2)ℳ0/(16 π)

data=Import["~/doc/research/sethna/doc/essential-ising_doc/data/data_square-4096.dat"];
t[Tc_][T_]:=Abs[((T-Tc)/T)]
assData=Thread[Rule[{"L","T","H","E","δE","M","δM","C","δC","χ","δχ"},#]]&/@data;
sus={"H"/"T"(((t[Tc]["T"]))^(-Δ)),"T""χ" (((t[Tc]["T"]))^γ)}/.#&/@assData;
mag={"H"/"T"(((t[Tc]["T"]))^(-Δ)),"M" (((t[Tc]["T"]))^(-β))}/.#&/@assData;

(* Constants from Mangazeev et al *)
Ch:=0.838677624411
Ghs={0,0,-1.8452280782328,0,8.333711750,0,-95.16896,0,1457.62,0,-25891};
Gls:={0,-1.3578383417066,-0.048953289720,0.038863932,-0.068362119,0.18388370,-0.6591714,2.937665,-15.61};
Φs:={-1.197733383797993,-0.318810124891,0.110886196683,0.01642689465,-2.639978 (10^(-4)),-5.140526 (10^(-4)),2.08856 (10^(-4)),-4.4819 (10^(-5))};

(* Parametric variables from Caselle et al *)
hCas[θ_]:=(θ−(θ^3)/1.16951)(1−0.222389(θ^2)−0.043547(θ^4)−0.014809(θ^6)−0.007168(θ^8))
h0Cas:=0.940647
θcCas:=1.08144

(* Translation between the variables of Mangazeev et al and our own *)
G2C[G_,n_]:=(((2ArcSinh[1]))^2)(((Ch (((2ArcSinh[1]))^(-Δ))))^n)G
Gs2Cs[Gs_]:=MapIndexed[G2C[#1,#2[[1]]-1]&,Gs]
Φ2C[Φ_,n_]:=(Ch^(2/Δ))(((-2ArcSinh[1](Ch^(-1/Δ))))^n)Φ
Φs2Cs[Gs_]:=MapIndexed[Φ2C[#1,#2[[1]]-1]&,Gs]

(* Set the rules for expansion coefficients *)
ClearAll[Gl];ClearAll[Gh];ClearAll[Φ];
Thread[SetDelayed[Evaluate[Gl/@(Range[Length[Gls]]-1)],Evaluate[Gs2Cs[Gls]]]];
Thread[SetDelayed[Evaluate[Gh/@(Range[Length[Ghs]]-1)],Evaluate[Gs2Cs[Ghs]]]];
Thread[SetDelayed[Evaluate[Φ/@(Range[Length[Φs]]-1)],Evaluate[Φs2Cs[Φs]]]];

ℱ[A_,B_,C_][X_]:=With[{Y=B X},-A Y(C-(1/π)Exp[1/Y]ExpIntegralEi[-1/Y])]
ℱf[A_,B_,C_,h0_,θc_,f_,Y_][θ_]:=With[{x=h0(-f'[θc])/(((((θc)^2)-1))^Δ)},ℱ[A,x B,C][θc-θ]+ℱ[A,x B,C][θc+θ]-2ℱ[A,x B,C][θc]+Y[θ]]

ℱseries[A_,B_,C_][n_]:=(A/π) (B^n)(((-1))^(n+1))  Gamma[n-1]/;n≥2
ℱseries[A_,B_,C_][1]:=-A B C
ℱseries[A_,B_,C_][0]:=0
ℱfseries[A_,B_,C_,h0_,θc_,f_,Y_][n_]:=With[{x=h0(-f'[θc])/(((((θc)^2)-1))^Δ)},(((-1))^n)ℱseries[A,x B,C][n]+SeriesCoefficient[ℱ[A,x B,C][θc+θ],{θ,θc,n},Assumptions->{θ<θc,f'[θc]<0,θc>1,B>0,C>0}]+SeriesCoefficient[Y[θ],{θ,θc,n}]]

ff[θc_,a_,N_][θ_]:=(1-((θ^2)/((θc)^2)))θ Sum[a[n](θ^(2n)),{n,0,N}]
YY[a_,λ_,N_][θ_]:=Sum[a[n](θ^(2n))Exp[-λ (θ^2)],{n,0,N}]

F1[G_][t_,h_]:=((Abs[t])^2)G[h ((Abs[t])^(-Δ))]
F2[Φ_][t_,h_]:=(h^(2/Δ))Φ[t (h^(-1/Δ))]

coeff[F_[ℱ_],R_,h0_,f_,θ0_,ass_:{}][n_]:=SeriesCoefficient[F[ℱ][R(1-(θ^2)),h0 (R^Δ)f[θ]],{θ,θ0,n},Assumptions->Join[{R>0,θ>0,θc>1,θ<θc},ass]]/.{f[0]->0,Derivative[m_][f][0]:>0/;EvenQ[m],f[θc]->0}/.Derivative[m_][ℱ][0]:>m!ℱ[m]

eql[n_,m_,r_]:=(coeff[F1[Gl],1,h0,f,θc,{θ>1}][n]==ℱfseries[A,B,m0,h0,θc,f,YY[b,λ,r]][n])/.f->ff[θc,a,m]
eqh[n_,m_,r_]:=coeff[F1[Gh],1,h0,f,0,{θ<1}][n]==SeriesCoefficient[ℱf[A,B,m0,h0,θc,f,YY[b,λ,r]][θ],{θ,0,n},Assumptions->{θ>0,h0>0,f'[0]>0,m0>0}]/.f->ff[θc,a,m]
eqm[n_,m_,r_]:=(coeff[F2[Φ],1,h0,f,1][n]==SeriesCoefficient[ℱf[A,B,m0,h0,θc,f,YY[b,λ,r]][θ],{θ,1,n},Assumptions->{θ>0,h0>0,f'[0]>0,m0>0}])/.f->ff[θc,a,m]

Nf:=14;NY=4;

eqsl=eql[#,Nf,NY]&/@(Range[8]-1);
eqsh=eqh[#,Nf,NY]&/@(2(Range[5]-1));
eqsm=eqm[#,Nf,NY]&/@(Range[8]-1);

sol1=Join[FindRoot[Join[eqsl[[;;3]],eqsm[[;;3]],eqsh[[2;;2]]]/.θc->1+Exp[xc]/.#,{{A,10},{h0,1},{m0,0},{b[1],0},{xc,-2},{a[1],0},{b[2],0}}+({0,5RandomReal[{-1,1}]}&/@Range[7]),MaxIterations->1000],#]&@Join[{a[0]->1,λ->0,b[0]->0},Table[a[n]->0,{n,2,14}],Table[b[n]->0,{n,3,4}]]

Plot[ℱf[A,B,m0,h0,1+Exp[xc],ff[1+Exp[xc],a,Nf],YY[b,λ,NY]][θ]/.sol1,{θ,-(1+Exp[xc])-.1/.sol1,1+Exp[xc]+.1/.sol1},GridLines->{{{1+Exp[xc]/.sol1,{Red,Thick}},{-(1+Exp[xc])/.sol1,{Red,Thick}}},None}]

func=Sort@(Reverse/@Table[{θ,h0 ff[1+Exp[xc],a,Nf][θ]/((Abs[1-(θ^2)])^Δ)/.sol1},{θ,1.0001,1+Exp[xc]-0.0001/.sol1,0.0001}]);

ffunc={#[[1]],-(((((#[[2]])^2)-1))^(-2))ℱf[A,B,m0,h0,1+Exp[xc],ff[1+Exp[xc],a,Nf],YY[b,λ,NY]][#[[2]]]/.sol1}&/@func;

Show[Plot[-ℱ[-π Gl[2]/((B^2)),B,ℳ0/B/(-π Gl[2]/((B^2)))][X],{X,0.01,1300}],ListPlot[ffunc]]

sol2=Join[FindRoot[Join[eqsl[[;;3]],eqsm[[;;3]]]/.θc->1+Exp[xc]/.#,{{A,A/.sol1},{xc,xc/.sol1},{m0,m0/.sol1},{a[1],a[1]/.sol1},{h0,h0/.sol1},{b[1],b[1]/.sol1}}+({0,0.1RandomReal[{-1,1}]}&/@Range[6]),MaxIterations->1000],#]&@Join[{a[0]->1,λ->0,b[0]->0},Table[a[n]->0,{n,2,14}],Table[b[n]->0,{n,2,4}]]

Plot[ℱf[A,B,m0,h0,1+Exp[xc],ff[1+Exp[xc],a,Nf],YY[b,λ,NY]][θ]/.sol2,{θ,-(1+Exp[xc])-.1/.sol2,1+Exp[xc]+.1/.sol2},GridLines->{{{1+Exp[xc]/.sol2,{Red,Thick}},{-(1+Exp[xc])/.sol2,{Red,Thick}}},None}]

func=Sort@(Reverse/@Table[{θ,h0 ff[1+Exp[xc],a,Nf][θ]/((Abs[1-(θ^2)])^Δ)/.sol2},{θ,1.0001,1+Exp[xc]-0.0001/.sol2,0.0001}]);

ffunc={#[[1]],-(((((#[[2]])^2)-1))^(-2))ℱf[A,B,m0,h0,1+Exp[xc],ff[1+Exp[xc],a,Nf],YY[b,λ,NY]][#[[2]]]/.sol2}&/@func;

Show[Plot[-ℱ[-π Gl[2]/((B^2)),B,ℳ0/B/(-π Gl[2]/((B^2)))][X],{X,0.01,1300}],ListPlot[ffunc]]

sol3=Join[FindRoot[Join[eqsl[[;;4]],eqsm[[;;4]],eqsh[[;;2]]]/.θc->1+Exp[xc]/.#,{{A,A/.sol2},{b[0],b[0]/.sol2},{h0,h0/.sol2},{m0,m0/.sol2},{b[1],b[1]/.sol2},{a[1],a[1]/.sol2},{xc,xc/.sol2},{a[2],a[2]/.sol2},{a[3],a[3]/.sol2},{b[2],b[2]/.sol2}}+({0,0.1RandomReal[{-1,1}]}&/@Range[10]),MaxIterations->1000],#]&@Join[{a[0]->1,λ->1},Table[a[n]->0,{n,4,14}],Table[b[n]->0,{n,3,4}]]

Plot[ℱf[A,B,m0,h0,1+Exp[xc],ff[1+Exp[xc],a,Nf],YY[b,λ,NY]][θ]/.sol3,{θ,-(1+Exp[xc])-1/.sol3,1+Exp[xc]+1/.sol3},GridLines->{{{1+Exp[xc]/.sol3,{Red,Thick}},{-(1+Exp[xc])/.sol3,{Red,Thick}}},None}]

func=Sort@(Reverse/@Table[{θ,h0 ff[1+Exp[xc],a,Nf][θ]/((Abs[1-(θ^2)])^Δ)/.sol3},{θ,1.0001,1+Exp[xc]-0.0001/.sol3,0.0001}]);

ffunc={#[[1]],-(((((#[[2]])^2)-1))^(-2))ℱf[A,B,m0,h0,1+Exp[xc],ff[1+Exp[xc],a,Nf],YY[b,λ,NY]][#[[2]]]/.sol3}&/@func;

Show[Plot[-ℱ[-π Gl[2]/((B^2)),B,ℳ0/B/(-π Gl[2]/((B^2)))][X],{X,0.01,1300}],ListPlot[ffunc]]

sol4=Join[FindRoot[Join[eqsl[[;;5]],eqsm[[;;5]],eqsh[[;;3]]]/.θc->1+Exp[xc]/.#,{{A,A/.sol3},{b[0],b[0]/.sol3},{h0,h0/.sol3},{m0,m0/.sol3},{b[1],b[1]/.sol3},{a[1],a[1]/.sol3},{xc,xc/.sol3},{a[2],a[2]/.sol3},{a[3],a[3]/.sol3},{b[2],b[2]/.sol3},{a[4],a[4]/.sol3},{a[5],a[5]/.sol3},{b[3],b[3]/.sol3}}+({0,0.1RandomReal[{-1,1}]}&/@Range[13]),MaxIterations->500],#]&@Join[{λ->1,a[0]->1},Table[a[n]->0,{n,6,14}],Table[b[n]->0,{n,4,4}]]

Plot[ℱf[A,B,m0,h0,1+Exp[xc],ff[1+Exp[xc],a,Nf],YY[b,λ,NY]][θ]/.sol4,{θ,-(1+Exp[xc])-1/.sol4,1+Exp[xc]+1/.sol4},GridLines->{{{1+Exp[xc]/.sol4,{Red,Thick}},{-(1+Exp[xc])/.sol4,{Red,Thick}}},None}]

func=Sort@(Reverse/@Table[{θ,h0 ff[1+Exp[xc],a,Nf][θ]/((Abs[1-(θ^2)])^Δ)/.sol4},{θ,1.0001,1+Exp[xc]-0.0001/.sol4,0.0001}]);

ffunc={#[[1]],-(((((#[[2]])^2)-1))^(-2))ℱf[A,B,m0,h0,1+Exp[xc],ff[1+Exp[xc],a,Nf],YY[b,λ,NY]][#[[2]]]/.sol4}&/@func;

Show[Plot[-ℱ[-π Gl[2]/((B^2)),B,ℳ0/B/(-π Gl[2]/((B^2)))][X],{X,0.01,1300}],ListPlot[ffunc]]

sol5=Join[FindRoot[Join[eqsl[[;;4]],eqsm[[;;4]],eqsh[[;;3]]]/.θc->1+Exp[xc]/.#,{{A,A/.sol4},{b[0],b[0]/.sol4},{xc,xc/.sol4},{h0,h0/.sol4},{m0,m0/.sol4},{a[1],a[1]/.sol4},{b[1],b[1]/.sol4},{a[2],a[2]/.sol4},{a[3],a[3]/.sol4},{a[4],a[4]/.sol4},{b[2],b[2]/.sol4}}+({0,0.1RandomReal[{-1,1}]}&/@Range[11]),MaxIterations->100],#]&@Join[{λ->0,a[0]->1},Table[a[n]->0,{n,5,14}],Table[b[n]->0,{n,3,4}]]

Plot[ℱf[A,B,m0,h0,1+Exp[xc],ff[1+Exp[xc],a,Nf],YY[b,λ,NY]][θ]/.sol5,{θ,-(1+Exp[xc])-.1/.sol5,1+Exp[xc]+.1/.sol5},GridLines->{{{1+Exp[xc]/.sol5,{Red,Thick}},{-(1+Exp[xc])/.sol5,{Red,Thick}}},None}]

sol6=Join[FindRoot[Join[eqsl[[;;5]],eqsm[[;;5]],eqsh[[;;3]]]/.θc->1+Exp[xc]/.#,{{A,A/.sol5},{b[0],b[0]/.sol5},{xc,xc/.sol5},{h0,h0/.sol5},{m0,m0/.sol5},{a[1],a[1]/.sol5},{b[1],b[1]/.sol5},{a[2],a[2]/.sol5},{a[3],a[3]/.sol5},{a[4],a[4]/.sol5},{b[2],b[2]/.sol5},{a[5],a[5]/.sol5},{a[6],a[6]/.sol5}}+({0,0.1RandomReal[{-1,1}]}&/@Range[13]),MaxIterations->100],#]&@Join[{λ->0,a[0]->1},Table[a[n]->0,{n,7,14}],Table[b[n]->0,{n,3,4}]]

Plot[ℱf[A,B,m0,h0,1+Exp[xc],ff[1+Exp[xc],a,Nf],YY[b,λ,NY]][θ]/.sol6,{θ,-(1+Exp[xc])-.1/.sol6,1+Exp[xc]+.1/.sol6},GridLines->{{{1+Exp[xc]/.sol6,{Red,Thick}},{-(1+Exp[xc])/.sol6,{Red,Thick}}},None}]

Plot[ff[1+Exp[xc],a,Nf][θ]/.sol3,{θ,0,2}]

func=Sort@(Reverse/@Table[{θ,h0 ff[1+Exp[xc],a,Nf][θ]/((Abs[1-(θ^2)])^Δ)/.sol1},{θ,1.0001,1+Exp[xc]-0.0001/.sol1,0.0001}]);

ffunc={#[[1]],(((((#[[2]])^2)-1))^(-2))ℱf[A,B,m0,h0,1+Exp[xc],ff[1+Exp[xc],a,Nf],YY[b,λ,NY]][#[[2]]]/.sol1}&/@func;

ListPlot[ffunc]

dFfunc=Table[{(ffunc[[i+1,1]]+ffunc[[i,1]])/2,-(ffunc[[i+1,2]]-ffunc[[i,2]])/(ffunc[[i+1,1]]-ffunc[[i,1]])},{i,1,Length[ffunc]-1}];

ListPlot[{Log[#[[1]]],#[[2]]}&/@#&/@{mag,dFfunc}]

Show[LogLogPlot[-ℱ[-π Gl[2]/((B^2)),B,ℳ0][X],{X,0.001,1000},PlotStyle->Red],ListLogLogPlot[ffunc,Joined->True,PlotStyle->Green]]

sol7=Join[FindRoot[Join[eqsl[[;;5]],eqsm[[;;5]],eqsh[[;;4]]]/.θc->1+Exp[xc]/.#,{{A,A/.sol6},{b[0],b[0]/.sol6},{xc,xc/.sol6},{h0,h0/.sol6},{m0,m0/.sol6},{a[1],a[1]/.sol6},{b[1],b[1]/.sol6},{a[2],a[2]/.sol6},{a[3],a[3]/.sol6},{a[4],a[4]/.sol6},{b[2],b[2]/.sol6},{a[5],a[5]/.sol6},{a[6],a[6]/.sol6},{b[3],a[3]/.sol6}}+({0,0.1RandomReal[{-1,1}]}&/@Range[14]),MaxIterations->100],#]&@Join[{λ->0,a[0]->1},Table[a[n]->0,{n,7,14}],Table[b[n]->0,{n,4,4}]]

Plot[ℱf[A,B,m0,h0,1+Exp[xc],ff[1+Exp[xc],a,Nf],YY[b,λ,NY]][θ]/.sol6,{θ,-(1+Exp[xc])-.1/.sol6,1+Exp[xc]+.1/.sol6},GridLines->{{{1+Exp[xc]/.sol6,{Red,Thick}},{-(1+Exp[xc])/.sol6,{Red,Thick}}},None}]

Plot[ℱf[A,B,m0,h0,θc,ff[θc,a,Nf],YY[b,λ,NY]][θ]/.sol2,{θ,-θc-1/.sol2,θc+1/.sol2},GridLines->{{{θc/.sol2,{Red,Thick}},{-θc/.sol2,{Red,Thick}}},None}]

Plot[𝒴f[A,B,h0,θc,ff[θc,a,Nf],YY[b,λ,NY]][θ]/.sol2,{θ,-θc-1/.sol2,θc+1/.sol2},GridLines->{{{θc/.sol2,{Red,Thick}},{-θc/.sol2,{Red,Thick}}},None}]

func=Sort@(Reverse/@Table[{θ,h0 ff[θc,a,Nf][θ]/((Abs[1-(θ^2)])^Δ)/.sol2},{θ,1.0001,θc-0.0001/.sol2,0.0001}]);

Show[ListLogLogPlot[sus],ListLogLogPlot[{#[[1]],(((((#[[2]])^2)-1))^(7/4))𝒴f[A,B,h0,θc,ff[θc,a,Nf],YY[b,λ,NY]][#[[2]]]/.sol2}&/@func,Joined->True,PlotStyle->Green],LogLogPlot[𝒴[-π Gl[2]/((B^2)),B][X],{X,0.001,1000},PlotStyle->Red]]

i1=Table[{(#[[i,1]]+#[[i+1,1]])/2,(#[[i+1,1]]-#[[i,1]])(#[[i,2]]+#[[i+1,2]])/2},{i,1,Length[#]-1}]&@({#[[1]],(((((#[[2]])^2)-1))^(7/4))𝒴f[A,B,h0,θc,ff[θc,a,Nf],YY[b,λ,NY]][#[[2]]]/.sol2}&/@func);

i2=Table[{0,ℳ0}+{#[[i,1]],Total[Take[#[[All,2]],i-1]]},{i,1,Length[#]}]&@i1;

Show[ListPlot[i2,Joined->True,PlotRange->{{0,1300},{1.22,1.75}},PlotStyle->Red],ListPlot[mag,PlotRange->All]]

ysave2=ysave/.sol2

yy2=yy1/.{Derivative[m_][y][_]:>m!ysave2[[m+1]],y[x_]:>ysave2[[1]]}/.sol2

ListLogPlot[{Abs[-(#(#+1)&/@Range[7])Gs2Cs[Gls][[2;;]]],Abs[yy2],Abs[𝒴series[-π Gl[2]/((B^2)),B][#]&/@(Range[7]-1)]}]

sol3=Join[FindRoot[Join[eqsl[[;;3]],eqsm[[;;3]],eqsh[[;;2]]]/.#,{{A,A/.sol2},{b[0],b[0]/.sol2},{h0,h0/.sol2},{a[1],a[1]/.sol2},{a[2],a[2]/.sol2},{b[1],b[1]/.sol2},{a[3],0},{a[4],0}}+({0,0.5RandomReal[{-1,1}]}&/@Range[8]),MaxIterations->5000],#]&@Join[{θc->(θc/.sol1),λ->0,a[0]->1},Table[a[n]->0,{n,5,14}],Table[b[n]->0,{n,2,4}]]

Plot[ff[θc,a,Nf][θ]/.sol3,{θ,0,θc/.sol3}]

Plot[𝒴f[A,B,h0,θc,ff[θc,a,Nf],YY[b,λ,NY]][θ]/.sol3,{θ,-θc-.1/.sol3,θc+.1/.sol3},GridLines->{{{θc/.sol3,{Red,Thick}},{-θc/.sol3,{Red,Thick}}},None}]

func3=Sort@(Reverse/@Table[{θ,h0 ff[θc,a,Nf][θ]/((Abs[1-(θ^2)])^Δ)/.sol3},{θ,1.00001,θc-0.00001/.sol3,0.00001}]);

Show[ListLogLogPlot[sus],ListLogLogPlot[{#[[1]],(((((#[[2]])^2)-1))^(7/4))𝒴f[A,B,h0,θc,ff[θc,a,Nf],YY[b,λ,NY]][#[[2]]]/.sol3}&/@func3,Joined->True,PlotStyle->Green],LogLogPlot[𝒴[-π Gl[2]/((B^2)),B][X],{X,0.001,1000},PlotStyle->Red]]

Export["~/doc/research/sethna/doc/essential-ising_doc/figs/fig-sus_scaling-func-mod.dat",{#[[1]],(((((#[[2]])^2)-1))^(7/4))𝒴f[A,B,h0,θc,ff[θc,a,Nf],YY[b,λ,NY]][#[[2]]]/.sol3}&/@func3]

i1=Table[{(#[[i,1]]+#[[i+1,1]])/2,(#[[i+1,1]]-#[[i,1]])(#[[i,2]]+#[[i+1,2]])/2},{i,1,Length[#]-1}]&@({#[[1]],(((((#[[2]])^2)-1))^(7/4))𝒴f[A,B,h0,θc,ff[θc,a,Nf],YY[b,λ,NY]][#[[2]]]/.sol3}&/@func3);

i2=Table[{0,ℳ0}+{#[[i,1]],Total[Take[#[[All,2]],i-1]]},{i,1,Length[#]}]&@i1;

Show[ListPlot[i2,Joined->True,PlotRange->{{0,1300},{1.22,1.75}},PlotStyle->Red],ListPlot[mag,PlotRange->All]]

Export["~/doc/research/sethna/doc/essential-ising_doc/figs/fig-mag_scaling-func-mod.dat",i2]

N[ℳ0]

ysave2=ysave/.sol3

yy2=yy1/.{Derivative[m_][y][_]:>m!ysave2[[m+1]],y[x_]:>ysave2[[1]]}/.sol3

ListLogPlot[{Abs[-(#(#+1)&/@Range[7])Gs2Cs[Gls][[2;;]]],Abs[yy2],Abs[𝒴series[-π Gl[2]/((B^2)),B][#]&/@(Range[7]-1)]}]

Export["~/doc/research/sethna/doc/essential-ising_doc/figs/fig-series-data.dat",Transpose@{Range[7]-1,Abs[𝒴series[-π Gl[2]/((B^2)),B][#]&/@(Range[7]-1)],Abs[yy2],Abs[-(#(#+1)&/@Range[7])Gs2Cs[Gls][[2;;]]]}]

sol4=Join[FindRoot[Join[eqsl[[;;4]],eqsm[[;;4]],eqsh[[;;2]]]/.#,{{A,A/.sol3},{b[0],b[0]/.sol3},{h0,h0/.sol3},{a[1],a[1]/.sol3},{a[2],a[2]/.sol3},{b[1],b[1]/.sol3},{a[3],a[3]/.sol3},{a[4],a[4]/.sol3},{a[5],0},{θc,1.055}}+({0,0.01RandomReal[{-1,1}]}&/@Range[10]),MaxIterations->2000],#]&@Join[{λ->0,a[0]->1},Table[a[n]->0,{n,6,14}],Table[b[n]->0,{n,2,4}]]

Plot[ff[θc,a,Nf][θ]/.sol3,{θ,0,θc/.sol3}]

Plot[𝒴f[A,B,h0,θc,ff[θc,a,Nf],YY[b,λ,NY]][θ]/.sol3,{θ,-θc-.1/.sol3,θc+.1/.sol3},GridLines->{{{θc/.sol3,{Red,Thick}},{-θc/.sol3,{Red,Thick}}},None}]

func3=Sort@(Reverse/@Table[{θ,h0 ff[θc,a,Nf][θ]/((Abs[1-(θ^2)])^Δ)/.sol3},{θ,1.0001,θc-0.0001/.sol3,0.0001}]);

Show[ListLogLogPlot[sus],ListLogLogPlot[{#[[1]],(((((#[[2]])^2)-1))^(7/4))𝒴f[A,B,h0,θc,ff[θc,a,Nf],YY[b,λ,NY]][#[[2]]]/.sol3}&/@func3,Joined->True,PlotStyle->Green],LogLogPlot[𝒴[-π Gl[2]/((B^2)),B][X],{X,0.001,1000},PlotStyle->Red]]

i1=Table[{(#[[i,1]]+#[[i+1,1]])/2,(#[[i+1,1]]-#[[i,1]])(#[[i,2]]+#[[i+1,2]])/2},{i,1,Length[#]-1}]&@({#[[1]],(((((#[[2]])^2)-1))^(7/4))𝒴f[A,B,h0,θc,ff[θc,a,Nf],YY[b,λ,NY]][#[[2]]]/.sol3}&/@func3);

i2=Table[{0,ℳ0}+{#[[i,1]],Total[Take[#[[All,2]],i-1]]},{i,1,Length[#]}]&@i1;

Show[ListPlot[i2,Joined->True,PlotRange->{{0,1300},{1.22,1.75}},PlotStyle->Red],ListPlot[mag,PlotRange->All]]

ysave2=ysave/.sol3

yy2=yy1/.{Derivative[m_][y][_]:>m!ysave2[[m+1]],y[x_]:>ysave2[[1]]}/.sol3

ListLogPlot[{Abs[-(#(#+1)&/@Range[7])Gs2Cs[Gls][[2;;]]],Abs[yy2],Abs[𝒴series[-π Gl[2]/((B^2)),B][#]&/@(Range[7]-1)]}]

sol3=Join[FindRoot[Join[eqsl[[;;3]],eqsm[[;;2]]]/.#,{{A,A/.sol2},{b[0],b[0]/.sol2},{h0,h0/.sol2},{a[1],a[1]/.sol2},{a[2],0}}+({0,0.1RandomReal[{-1,1}]}&/@Range[5]),MaxIterations->10000],#]&@Join[{θc->(θc/.sol2),a[0]->1,λ->0},Table[a[n]->0,{n,3,14}],Table[b[n]->0,{n,1,4}]]

Plot[𝒴f[A,B,h0,θc,ff[θc,a,Nf],YY[b,λ,NY]][θ]/.sol3,{θ,-θc-.1/.sol3,θc+.1/.sol3},GridLines->{{{θc/.sol3,{Red,Thick}},{-θc/.sol3,{Red,Thick}}},None}]

func3=Sort@(Reverse/@Table[{θ,h0 ff[θc,a,Nf][θ]/((Abs[1-(θ^2)])^Δ)/.sol3},{θ,1.0001,θc-0.0001/.sol3,0.0001}]);

Plot[ff[θc,a,Nf][θ]/.sol3,{θ,0,θc/.sol3}]

Show[ListLogLogPlot[sus],ListLogLogPlot[{#[[1]],(((((#[[2]])^2)-1))^(7/4))𝒴f[A,B,h0,θc,ff[θc,a,Nf],YY[b,λ,NY]][#[[2]]]/.sol3}&/@func3,Joined->True,PlotStyle->Green],LogLogPlot[𝒴[-π Gl[2]/((B^2)),B][X],{X,0.001,1000},PlotStyle->Red]]

sol4=Join[FindRoot[Join[eqsl[[;;3]],eqsm[[;;3]]]/.#,{{A,A/.sol3},{b[0],b[0]/.sol3},{h0,h0/.sol3},{a[1],a[1]/.sol3},{a[2],a[2]/.sol3},{a[3],0}}+({0,0.1RandomReal[{-1,1}]}&/@Range[6]),MaxIterations->100000],#]&@Join[{θc->(θc/.sol3),a[0]->1,λ->0},Table[a[n]->0,{n,4,14}],Table[b[n]->0,{n,1,4}]]

sol4=Join[FindRoot[eqset1[[{1,2,3,4,5,6,7}]]/.#,Evaluate[{{A,A/.sol3},{b[0],b[0]/.sol3},{a[1],a[1]/.sol3},{a[3],a[3]/.sol3},{b[2],b[2]/.sol3},{a[5],0.},{b[4],0}}+({0,RandomReal[{-1,1}]}&/@Range[7])],MaxIterations->5000],#]&@{h0->1,a[7]->0,a[9]->0,b[8]->0,b[6]->0,θc->1.05 }

Plot[ff[θc,a,Nf][θ]/.sol4,{θ,0,θc/.sol4}]

Plot[𝒴f[A,B,h0,θc,ff[θc,a,Nf],YY[b,NY]][θ]/.sol4,{θ,-θc-1/.sol4,θc+1/.sol4},GridLines->{{{θc/.sol4,{Red,Thick}},{-θc/.sol4,{Red,Thick}}},None}]

func4=Sort@(Reverse/@Table[{θ,ff[θc,a,Nf][θ]/((Abs[1-(θ^2)])^Δ)/.sol4},{θ,1.0001,θc/.sol4,0.0001}]);

Show[ListLogLogPlot[sus],ListLogLogPlot[{#[[1]],(((((#[[2]])^2)-1))^(7/4))𝒴f[A,B,h0,θc,ff[θc,a,Nf],YY[b,NY]][#[[2]]]/.sol4}&/@func4,Joined->True,PlotStyle->Green],LogLogPlot[𝒴[-π Gl[2]/((B^2)),B][X],{X,0.001,1000},PlotStyle->Red]]

i1[[3]]

sol5=Join[FindRoot[eqset1[[{1,2,3,4,6,5}]]/.#,{{A,A/.sol4},{b[0],b[0]/.sol4},{a[5],a[5]/.sol4},{a[3],a[3]/.sol4},{a[7],a[7]/.sol4},{h0,0.15}},MaxIterations->20000],#]&@{a[9]->0,a[1]->1,b[2]->0,b[4]->0,b[6]->0,b[8]->0,θc->θcCas}

Plot[ff[θc,a,Nf][θ]/.sol5,{θ,0,θc/.sol4}]

Plot[𝒴f[A,B,h0,θc,ff[θc,a,Nf],YY[b,NY]][θ]/.sol5,{θ,-θc-.1/.sol5,θc+.1/.sol5},GridLines->{{{θc/.sol5,{Red,Thick}},{-θc/.sol5,{Red,Thick}}},None}]

func5=Sort@(Reverse/@Table[{θ,ff[θc,a,Nf][θ]/((Abs[1-(θ^2)])^Δ)/.sol5},{θ,1.0001,θc/.sol5,0.0001}]);

Show[ListLogLogPlot[sus],ListLogLogPlot[{#[[1]],(((((#[[2]])^2)-1))^(7/4))𝒴f[A,B,h0,θc,ff[θc,a,Nf],YY[b,NY]][#[[2]]]/.sol5}&/@func5,Joined->True,PlotStyle->Green],LogLogPlot[𝒴[-π Gl[2]/((B^2)),B][X],{X,0.001,1000},PlotStyle->Red]]