Δ:=15/8
γ:=7/4
β:=1/8

B:=(Tc^2)ℳ0/(16 π)
ℳ0:=((((2^(5/2))ArcSinh[1]))^β)
Tc:=2/Log[1+Sqrt[2]]

Ch:=0.838677624411
C0:=0.0255369719
Gs:={-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,-5.140526,2.08856,-4.4819};

t[Tc_][T_]:=Abs[((T-Tc)/T)]

data=Import["~/doc/research/sethna/doc/essential-ising_doc/data/data_square-4096.dat"];

assData=Thread[Rule[{"L","T","H","E","δE","M","δM","C","δC","χ","δχ"},#]]&/@data;

htΔvχtγ={"H"/"T"(((t[Tc]["T"]))^(-Δ)),"T""χ" (((t[Tc]["T"]))^γ)}/.#&/@assData;

G2C[G_,n_]:=(((2ArcSinh[1]))^2)(((Ch (((2ArcSinh[1]))^(-Δ))))^n)G

Gs2Cs[Gs_]:=MapIndexed[G2C[#1,#2[[1]]]&,Gs]

Φ2C[Φ_,n_]:=(((2ArcSinh[1]))^2)(((Ch (((2ArcSinh[1]))^(-15/8))))^(16/15-8/15n))Φ

Φ2f[Φ_,n_]:=((1/15)-((8 n)/15)) ((16/15)-((8 n)/15))Φ2C[Φ,n]

Φs2fs[Φs_]:=MapIndexed[Φ2f[#1,#2[[1]]-1]&,Φs]

𝒴[A_,B_][X_]:=With[{Y=B X},((A (B^2))/(π (Y^3))) (Y(Y-1)-Exp[1/Y]ExpIntegralEi[-1/Y])]

f[c_,λ_,n_][X_]:=((c (X^n))/(1+(((λ X))^(n+1))))

F[λ_,{Cs___}][X_]:=Total@MapIndexed[f[#1,λ,#2[[1]]-1][X]&,{Cs}]

t𝒴[A_,B_,λ_,λ2_,{Cs___}][X_]:=(((1+B X))^(1/15))With[{Y=X},(𝒴[A,B][Y]+F[λ,{Cs}][B Y])]

coeff[A_,λ_,λ2_,{cs___}][m_]:=SeriesCoefficient[-1/((2+m)(1+m))t𝒴[A,B,λ,λ2,{cs}][X],{X,0,m},Assumptions->{X>0}];

tail[A_,λ_,λ2_,{cs___}][m_]:=SeriesCoefficient[-t𝒴[A,B,λ,λ2,{cs}][X],{X,∞,6/15+m 8/15},Assumptions->{X>0}]

coeq[gs_,Φs_]:=Join[MapThread[Equal,{MapIndexed[G2C[#1,1+#2[[1]]]&,gs],coeff[A,λ,λ2,c/@(Range[Length[gs]+Length[Φs]-2]+1)][#]&/@(Range[Length[gs]]-1)}],MapThread[Equal,{tail[A,λ,λ2,c/@(Range[Length[gs]+Length[Φs]-2]+1)][#]&/@Range[Length[Φs]],MapIndexed[Φ2f[#1,#2[[1]]-1]&,Φs]}]]

coeffSolve[gs_,Φs_]:=Solve[coeq[gs,Φs],Join[{A,λ},c/@(Range[Length[gs]+Length[Φs]-2]+1)],Reals]

n1:=6
n2:=1

sol=coeffSolve[Take[Gs,{2,2+n1}],Take[Φs,n2]];

resol=Select[sol,Re@#[[2,2]]>0&&AllTrue[#[[All,2]],Im[#]==0&]&]

Show[LogLogPlot[{𝒴[π/2 C0/(B^2),B][Y],𝒴[π/2 C0/(B^2),B][Y](((1+(((0.0705991/(((π/2 C0/(B^2))/(2 (2^(11/16)) (π^2) ((ArcSinh[1])^(15/8)))))))^15)Y))^(1/15)),t𝒴[A,B,λ,λ2,c/@(Range[n1+n2-1]+1)][Y]/.resol[[1]],0.0705991 (Y^(-14/15))},{Y,0.001,10000},PlotRange->{0.00001,.05}],ListLogLogPlot[htΔvχtγ]]

χfitData={"H"/"T" ((t[Tc]["T"])^(-Δ)),"T" "χ" ((t[Tc]["T"])^γ)}/.#&/@assData;

χfitWeights=1/(((Tc "δχ"((t[Tc]["T"])^γ)))^2)/.#&/@assData;

nn:=1

fit=NonlinearModelFit[χfitData,t𝒴[π/2 (C0-c[2])/(B^2),B,λ,λ2,c/@(Range[nn]+1)][X],Join[{λ,λ2},Thread[{c/@(Range[nn]+1),ConstantArray[0,nn]}]],X,PrecisionGoal->70,MaxIterations->1000,Method->Automatic,Weights->χfitWeights];

SeriesCoefficient[Evaluate[t𝒴[A,B,λ,λ2,c/@(Range[n1+n2-1]+1)][X]/.sol[[1]]],{X,0,0},Assumptions->{X>0}]

fit=NonlinearModelFit[χfitData,{Evaluate[t𝒴[A,B,λ,λ2,c/@(Range[n1+n2-1]+1)][X]/.sol[[1]]]},{λ},X,PrecisionGoal->90,MaxIterations->1000,Method->Automatic,Weights->χfitWeights];

fit["ParameterTable"]

pχlog=Show[ListPlot[Log10@htΔvχtγ,PlotRange->{{-3,6},{-5,-1}}],Plot[{Log10@t𝒴[π/2 C0/((B^2)),B,0,λ2,{}][10^X],Log10@(fit[10^X])},{X,-3,4},PlotStyle->{Black,{Black,Dashed}}]]

fit["BestFitParameters"]

ListLogPlot[{Join[{ℳ0},Abs@coeff[π/2 C0/(B^2),1,0,{}][#]&/@(Range[15]-1)],Abs@Gs2Cs[Gs],Join[{ℳ0},(x|->Abs@coeff[A,λ,λ2,c/@(Range[n1+n2-1]+1)][#]/.x&/@(Range[15]-1))@sol[[1]]],Join[{ℳ0},Abs@SeriesCoefficient[t𝒴[A,B,λ,λ2,c/@(Range[n1+n2-1]+1)][X]/((#+2)(#+1))/.sol[[1]]/.fit["BestFitParameters"],{X,0,#},Assumptions->{X>0}]&/@(Range[15]-1)]}]

ratioData=With[{X="H"/"T" ((t[Tc]["T"])^(-Δ))},{X,"T" "χ" ((t[Tc]["T"])^γ)/𝒴[π/2 C0/((B^2)),B][X]}]/.#&/@assData;

ratioSmall[X_]:=-2Rest@Gs2Cs[Gs].((X^(#))&/@(Range[Length[Gs]-1]-1))/𝒴[π C0/(2(B^2)),B][X]

ratioLarge[X_]:=-Φs2fs[Φs].((X^(-(6/15+8/15#)))&/@Range[Length[Φs]])/𝒴[π C0/(2(B^2)),B][X]

ratioLarge2[X_]:=-Φs2fs[Φs].((X^(-(#)))&/@Range[Length[Φs]])/𝒴[π C0/(2(B^2)),B][X]

Show[ListLogLogPlot[ratioData,PlotRange->{{0.001,100000},{0.6,1.6}}],LogLogPlot[{ratioSmall[X],ratioLarge[X],ratioLarge2[X],(((1+(((-Φs2fs[Φs][[1]]/(C0/(2B))))^15)X))^(1/15)),-Φs2fs[Φs][[1]]/(C0/(2B))(X^(1/15)),fit[X]/𝒴[π C0/(2(B^2)),B][X],t𝒴[A,B,λ,λ2,c/@(Range[n1+n2-1]+1)][X]/𝒴[π C0/(2(B^2)),B][X]/.resol[[1]]},{X,0.001,100000}]]

D[Sum[a[n](x^n),{n,0,8}],{x,2}]

1/((2+#)(1+#))&/@(Range[10]-1)