(* Constants fixed for the 3D simple cubic-lattice Ising model *)
δ=4.78984;
β=0.326419;
ν=0.629971;
γ=1.237075;
α=0.11008;
Δ=β δ;
Tc=1/0.22165455;
ℳ0=1.6919045; (*Talapov_1996_The*)
𝒮0=10.083/(Tc^(2ν));(*Caselle_2007_The*)
B=(1/2)ℳ0(Sqrt[(3/π)]) ((((1/(𝒮0))))^(3/2));

(2.08^(3/2))

Exp[0.4]

10.083/(Tc^(2  ν))

𝒮0

data=Dataset[With[
{
info=ToExpression/@Rest[StringSplit[FileBaseName[#],"_"]],
dat=Import[#]
},
<|
"D"->info[[1]],
"L"->info[[2]],
"T"->info[[3]],
"H"->info[[4]],
"Me"->Around@@dat[[1]],
"Mm"->Around@@dat[[2]],
"Mt"->Around@@dat[[3]],
"Ma"->Around@@dat[[4]],
"χe"->Around@@dat[[5]],
"χm"->Around@@dat[[6]],
"χt"->Around@@dat[[7]],
"χa"->Around@@dat[[8]]
|>]&/@FileNames["~/doc/research/first_order_singularities/data/3d_ising/metastable_3_*_m.dat.processed"]];

t[Tc_][T_]:=Abs[((T-Tc)/T)]
sus=data[All,{#H/#T ((t[Tc][#T])^(-Δ)),#χe["Value"]#T ((t[Tc][#T])^γ)/((#L)^(#D))}&];
mag=data[All,{#H/#T ((t[Tc][#T])^(-Δ)),#T #Me["Value"] ((t[Tc][#T])^(-β))/((#L)^(#D))}&];

(* Coefficients from Butera and Pernici *)
Ghs=(#((#)!&/@(Range[Length[#]])))&@-{0,0.5608,0,-0.1608,0,0.146,0,-0.187,0,0.26,0,-0.19,0,0.026};
Gls={-ℳ0,-0.99051/(2 1), Tc 6.015/ ℳ0(((Ghs[[2]]/4.78))^2)/(2 3),Tc Ghs[[4]]/9.2/(3 4)};
Φs={-4.268216450407312`};

Ghs

Gls

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

Integrate[Exp[-1/X]Exp[-1/(X^2)](X^(7/3))

𝒴[A_,B_][X_]:=With[{Y=B X},((2 A (B^2))/π)((1/(54 (Y^8)))(-6 (Y^2) (Y (-6+13 (Y^2)) Gamma[(2/3)]+3 (3-8 (Y^2)) Gamma[(7/6)])+(E^(-(1/(Y^2)))) (18-69 (Y^2)+35 (Y^4)) (2 Y ExpIntegralE[(5/3),-(1/(Y^2))] Gamma[(2/3)]-3 ExpIntegralE[(7/6),-(1/(Y^2))] Gamma[(7/6)])))]
𝒴f[A_,B_,h0_,θc_,f_,Y_][θ_]:=With[{x=h0(-f'[θc])/((Abs[((θc)^2)-1])^Δ)},𝒴[A,B][x(θc-θ)]+𝒴[A,B][x(θc+θ)]+Y[θ]]

ClearAll[𝒴num]

𝒴dat=Table[{10^x,Re@𝒴[1,1][SetPrecision[(10^x),60]]},{x,-7,7,1/10}];

𝒴num=Interpolation[𝒴dat];

𝒴fnum[A_,B_,h0_,θc_,f_,Y_][θ_]:=With[{x=h0(-f'[θc])/((Abs[((θc)^2)-1])^Δ)},A (((B))^2)𝒴num[x B(θc-θ)]+A (((B))^2)𝒴num[x B(θc+θ)]+Y[θ]]

Show[ListLogLogPlot[{#1 ,#2}&@@@sus],
LogLogPlot[-2Gls[[2]]𝒴num[B x]/𝒴series[1,1][0],{x,0.001,10}]
]

𝒴series[A_,B_][n_]:=(A/π) (B^(2+n))(((-1))^n)   Gamma[(13/6)+(n/2)]
𝒴fseries[A_,B_,h0_,θc_,f_,Y_][n_]:=With[{x=h0(-f'[θc])/((Abs[((θc)^2)-1])^Δ)},(((-1))^n)𝒴series[A,x B][n]/(x^2)+SeriesCoefficient[𝒴[A,B][x(θ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-α)/Δ))Φ[t (h^(-1/Δ))]

dF2[tc_,hc_,ΔΔ_][ϕ_][t_,h_]:=((1.88992` (h^(-2+((1.88992`)/(ΔΔ)))) (-1+((1.88992`)/(ΔΔ))) ϕ[(h^(-1/ΔΔ)) t])/(ΔΔ))-((3.77984` (h^(-2+((0.88992`)/(ΔΔ)))) t ϕ'[(h^(-1/ΔΔ)) t])/((ΔΔ)^2))+(h^(1.88992`/ΔΔ)) (-(((h^(-2-(1/(ΔΔ)))) t (-1-(1/(ΔΔ))) ϕ'[(h^(-1/ΔΔ)) t])/(ΔΔ))+(((h^(-2-(2/(ΔΔ)))) (t^2) ϕ''[(h^(-1/ΔΔ)) t])/((ΔΔ)^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:=0;NY=0;

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

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

yy0=Table[(SeriesCoefficient[(((((g[x])^2)-1))^γ)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])

Gh[2]/Gl[2]

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

sol1=Thread[Rule[sol1[[All,1]],SetPrecision[Re[sol1[[All,2]]],20]]]

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

θcf[xc]/.sol1

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

func[[1]]

susData={#[[1]],Re[(((((#[[2]])^2)-1))^γ) 𝒴fnum[A,B,h0,θcf[xc],ff[θcf[xc],a,Nf],YY[b,λ,NY]][#[[2]]]/.sol1]}&/@func;

Show[ListLogLogPlot[sus],ListLogLogPlot[susData,Joined->True,PlotStyle->Green,PlotRange->All],LogLogPlot[{
𝒴[-π Gl[2]/((B^2))/Gamma[13/3],B][X]
},{X,0.001,1000},PlotStyle->Red]]

Export["~/doc/research/thesis/data/essential/3d_scaling_1.dat",susData];

i1=Table[{(#[[i,1]]+#[[i+1,1]])/2,(#[[i+1,1]]-#[[i,1]])(#[[i,2]]+#[[i+1,2]])/2},{i,1,Length[#]-1}]&@susData;

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

i2[[-10]]

Show[ListPlot[mag,PlotRange->All],ListPlot[i2,Joined->True,PlotStyle->Red,PlotRange->All]]

Export["~/doc/research/thesis/data/essential/3d_scaling_m.dat",i2]

θcf[xc/.sol1]

sol1

SetPrecision[yy1[[8]],50]/.Join[{θc->θcf[xc/.sol1]},sol1[[;;6]]]

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

Gls

ListLogPlot[{Abs[𝒴series[-2Gl[2]/𝒴series[1,B][0],B][#]&/@(Range[13]-1)],Abs[-(#(#+1)&/@Range[3])Gls[[2;;]]],Abs[yy2]}]

Export[
"~/doc/research/thesis/data/essential/3d_coefficients_1.dat"
,Thread[{Abs[-(#(#+1)&/@Range[13])PadRight[Gls[[2;;]],13]],Abs[𝒴series[-2Gl[2]/𝒴series[1,B][0],B][#]&/@(Range[13]-1)],Abs[yy2][[;;13]]}]
]