(* 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,-1.8452280782328,0,8.333711750,0,-95.16896,0,1457.62,0,-25891};
Gls:={-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]]]&,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]]],Evaluate[Gs2Cs[Gls]]]];
Thread[SetDelayed[Evaluate[Gh/@Range[Length[Ghs]]],Evaluate[Gs2Cs[Ghs]]]];
Thread[SetDelayed[Evaluate[Φ/@(Range[Length[Φs]]-1)],Evaluate[Φs2Cs[Φs]]]];

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

𝒴series[A_,B_][n_]:=(A/π) (B^(2+n))(((-1))^n)  Gamma[n+3]
𝒴fseries[A_,B_,h0_,θc_,f_,Y_][n_]:=With[{x=h0(-f'[θc])/(((((θc)^2)-1))^Δ)},(((-1))^n)𝒴series[A,x B][n]+SeriesCoefficient[𝒴[A,x B][θc+θ],{θ,θc,n},Assumptions->{θ<θc,f'[θc]<0,θc>1,B>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[tc_,hc_,Δ_][Φ_][t_,h_]:=(h^(2/Δ))With[{η=tc t (((hc h))^(-1/Δ))},Φ[t (h^(-1/Δ))]-((η^2)/(8π))Log[(η^2)]]

dF2[tc_,hc_,δ_][ϕ_][t_,h_]:=((-8 (h^(2/δ)) (-2+δ) ϕ[(h^(-1/δ)) t]+t (-(((h^(2/δ)) (((h hc))^(-2/δ)) t (tc^2) δ)/π)+4 (h^(1/δ)) (-3+δ) ϕ'[(h^(-1/δ)) t]+4 t ϕ''[(h^(-1/δ)) t]))/(4 (h^2) (δ^2)))

coeff[F_[ℱ_],R_,h0_,f_,θ0_,ass_:{}][n_]:=SeriesCoefficient[-Derivative[0,2][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]

dcoeff[dF_[ℱ_],R_,h0_,f_,θ0_,ass_:{}][n_]:=SeriesCoefficient[-dF[ℱ][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,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,h0,θc,f,YY[b,λ,r]][θ],{θ,0,n},Assumptions->{θ>0,h0>0,f'[0]>0}]/.f->ff[θc,a,m]
eqm[n_,m_,r_]:=(dcoeff[dF2[ArcSinh[1],Ch,Δ][Φ],1,h0,f,1][n]==SeriesCoefficient[𝒴f[A,B,h0,θc,f,YY[b,λ,r]][θ],{θ,1,n},Assumptions->{θ>0,h0>0,f'[0]>0}])/.f->ff[θc,a,m]

Nf:=14;NY=4;

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

ysave=Table[𝒴fseries[A,B,h0,θc,f,YY[b,λ,NY]][n],{n,0,13}]/.f->ff[θc,a,Nf];

yy0=Table[(SeriesCoefficient[(((((g[x])^2)-1))^(7/4))y[g[x]],{x,0,i}]),{i,0,13}];

ds=Table[D[h0 ff[θc,a,Nf][θ]/((((θ^2)-1))^Δ),{θ,n}],{n,0,13}]/.θ->θc;

dinvs=Table[D[InverseFunction[X][x],{x,n}],{n,0,13}]/.InverseFunction[X][x]->θc/.Derivative[n_][X][θc]:>ds[[n+1]];

yy1=yy0/.{g[0]->dinvs[[1]],Derivative[n_][g][0]:>dinvs[[n+1]],Derivative[m_][y][_]:>m!ysave[[m+1]],y[x_]:>ysave[[1]]};

θcf[xc_]:=1+1/6(1+Tanh[xc])

θcf[-0.85]N

sol1=Join[FindRoot[Join[eqsl[[;;1]],eqsm[[;;1]],eqsh[[;;1]]]/.θc->θcf[xc]/.#,{{A,0.1},{h0,1.5},{b[0],0.91}}+({0,0.01RandomReal[{-1,1}]}&/@Range[3]),MaxIterations->3000],#]&@Join[{λ->0,a[0]->1,xc->-0.85},Table[a[n]->0,{n,1,14}],Table[b[n]->0,{n,1,4}]]

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

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

Show[ListLogLogPlot[sus],ListLogLogPlot[{#[[1]],(((((#[[2]])^2)-1))^(7/4))𝒴f[A,B,h0,θcf[xc],ff[θcf[xc],a,Nf],YY[b,λ,NY]][#[[2]]]/.sol1}&/@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,θcf[xc],ff[θcf[xc],a,Nf],YY[b,λ,NY]][#[[2]]]/.sol1}&/@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]]

yy2=yy1/.Join[{θc->θcf[xc/.sol1]},sol1]

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

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

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

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

Show[ListLogLogPlot[sus],ListLogLogPlot[{#[[1]],(((((#[[2]])^2)-1))^(7/4))𝒴f[A,B,h0,θcf[xc],ff[θcf[xc],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,θcf[xc],ff[θcf[xc],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]]

yy2=yy1/.Join[{θc->θcf[xc/.sol2]},sol2]

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

sol3=Join[FindRoot[Join[eqsl[[;;3]],eqsm[[;;3]],eqsh[[;;2]]]/.θc->θcf[xc]/.#,{{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],a[3]/.sol2},{a[4],a[4]/.sol2}}+({0,0.1RandomReal[{-1,1}]}&/@Range[8]),MaxIterations->500],#]&@Join[{λ->0,a[0]->1,xc->(xc/.sol2)},Table[a[n]->0,{n,5,14}],Table[b[n]->0,{n,2,4}]]

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

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

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

Show[ListLogLogPlot[sus],ListLogLogPlot[{#[[1]],(((((#[[2]])^2)-1))^(7/4))𝒴f[A,B,h0,θcf[xc],ff[θcf[xc],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,θcf[xc],ff[θcf[xc],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]]

yy2=yy1/.Evaluate[Join[{θc->θcf[xc/.sol3]},sol3]]

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

sol4=Join[FindRoot[Join[eqsl[[;;4]],eqsm[[;;4]],eqsh[[;;3]]]/.θc->θcf[xc]/.#,{{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],a[5]/.sol3},{xc,xc/.sol3},{b[2],b[2]/.sol3}}+({0,0.1RandomReal[{-1,1}]}&/@Range[11]),MaxIterations->500],#]&@Join[{λ->0,a[0]->1},Table[a[n]->0,{n,6,14}],Table[b[n]->0,{n,3,4}]]

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

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

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

Show[ListLogLogPlot[sus],ListLogLogPlot[{#[[1]],(((((#[[2]])^2)-1))^(7/4))𝒴f[A,B,h0,θcf[xc],ff[θcf[xc],a,Nf],YY[b,λ,NY]][#[[2]]]/.sol4}&/@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,θcf[xc],ff[θcf[xc],a,Nf],YY[b,λ,NY]][#[[2]]]/.sol4}&/@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,θcf[xc],ff[θcf[xc],a,Nf],YY[b,λ,NY]][#[[2]]]/.sol4}&/@func3);

i2=Table[{0,N@ℳ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]

yy2=yy1/.Evaluate[Join[{θc->θcf[xc/.sol4]},sol4]]

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

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

runtilsuccess[expr_]:=Check[ReleaseHold[expr],runtilsuccess[expr]]

FindRoot::imag="Imaginary solution."

sol5=runtilsuccess[Hold[If[0==Im[#[[1,2]]],#,Message[FindRoot::imag];#]&@Join[FindRoot[Join[eqsl[[;;7]],eqsm[[;;7]],eqsh[[;;5]]]/.θc->θcf[xc]/.#,{{A,A/.sol4},{b[0],b[0]/.sol4},{h0,h0/.sol4},{a[1],a[1]/.sol4},{a[2],a[2]/.sol4},{b[1],b[1]/.sol4},{a[3],a[3]/.sol4},{a[4],a[4]/.sol4},{a[5],a[5]/.sol4},{xc,xc/.sol4},{b[2],b[2]/.sol4},{a[6],a[6]/.sol4},{a[7],a[7]/.sol4},{a[8],a[8]/.sol4},{a[9],a[9]/.sol4},{a[10],a[10]/.sol4},{a[11],a[11]/.sol4},{b[3],b[3]/.sol4},{b[4],b[4]/.sol4}}+({0,1RandomReal[{-1,1}]}&/@Range[19]),MaxIterations->1000],#]&@Join[{λ->0,a[0]->1},Table[a[n]->0,{n,12,14}],Table[b[n]->0,{n,5,4}]]]]

Plot[ff[θcf[xc],a,Nf][θ]/.sol5,{θ,0,θcf[xc]/.sol5}]

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

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

Show[ListLogLogPlot[sus],ListLogLogPlot[{#[[1]],(((((#[[2]])^2)-1))^(7/4))𝒴f[A,B,h0,θcf[xc],ff[θcf[xc],a,Nf],YY[b,λ,NY]][#[[2]]]/.sol5}&/@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,θcf[xc],ff[θcf[xc],a,Nf],YY[b,λ,NY]][#[[2]]]/.sol5}&/@func3);

i2=Table[{0,N@ℳ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]]

yy2=yy1/.Evaluate[Join[{θc->θcf[xc/.sol4]},sol4]]

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

sol6=runtilsuccess[Hold[If[0==Im[#[[1,2]]],#,Message[FindRoot::imag];#]&@Join[FindRoot[Join[eqsl[[;;6]],eqsm[[;;5]],eqsh[[;;3]]]/.θc->(1+Exp[xc])/.#,{{A,A/.sol5},{b[0],b[0]/.sol5},{h0,h0/.sol5},{a[1],a[1]/.sol5},{a[2],a[2]/.sol5},{b[1],b[1]/.sol5},{a[3],a[3]/.sol5},{a[4],a[4]/.sol5},{a[5],a[5]/.sol5},{xc,xc/.sol5},{b[2],b[2]/.sol5},{a[6],a[6]/.sol5},{a[7],a[7]/.sol5},{b[3],b[3]/.sol5}}+({0,RandomReal[{-1,1}]}&/@Range[14]),MaxIterations->1000],#]&@Join[{λ->0,a[0]->1},Table[a[n]->0,{n,8,14}],Table[b[n]->0,{n,4,4}]]]]

Plot[ff[(1+Exp[xc]),a,Nf][θ]/.sol6,{θ,0,(1+Exp[xc])/.sol6}]

Plot[𝒴f[A,B,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}]

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

Show[ListLogLogPlot[sus],ListLogLogPlot[{#[[1]],(((((#[[2]])^2)-1))^(7/4))𝒴f[A,B,h0,(1+Exp[xc]),ff[(1+Exp[xc]),a,Nf],YY[b,λ,NY]][#[[2]]]/.sol4}&/@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,(1+Exp[xc]),ff[(1+Exp[xc]),a,Nf],YY[b,λ,NY]][#[[2]]]/.sol4}&/@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/.θc->(1+Exp[xc])/.sol4

yy2=yy1/.θc->(1+Exp[xc])/.{Derivative[m_][y][_]:>m!ysave2[[m+1]],y[x_]:>ysave2[[1]]}/.sol4

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]]