Get["~/doc/research/first_order_singularities/mma/schofield.wl"]

SetAttributes[untilSuccess,HoldAll];
untilSuccess[e_]:=Quiet[Check[e,untilSuccess[e]]]

randomEval[eq_]:=eq/.Around[x_,y_]:>RandomVariate[NormalDistribution[x,y]]

$DefaultRemoteKernel=RemoteKernelObject["ssh://math@sident.lassp.cornell.edu"];

Match ends

RFHigh[1][4.]

Series[RFHigh[θ0][θ],{θ,I θ0,1}]

Plot[Evaluate[ReIm@RFHigh[1][θ]],{θ,-2,2}]

Generate fits

test=Integrate[x Exp[-b/x]/(((x+y))^(n+1)),{x,0,∞},Assumptions->{b>0,n∈Integers,n>2,y>0}]

ListLogPlot[(n|->Evaluate[os/(2π)/π test/.y->3/.b->2 os/π])/@Range[0,14]]

dfdξ[ff_,h_,s_][m_,x_]:=Module[{ds,dF,df},
ds=invDerivativeList[m+1][Function[θ,h[θ]/t[θ]^Δ[2]/.s]][x];
dF=NestList[Function[f,D[f,θ]],ff,m]/. θ->x;
df=NestList[D[#,θ]&,fp[θ]/t[θ]^2,m]/.  Map[Derivative[#][fp][θ]->dF[[#+1]]&,Range[0,m]]/. θ->x;
Table[Sum[df[[k+1]] BellY[j,k,ds[[;;j-k+1]]],{k,0,j}]/(j!),{j,0,m}]
]

test=dfdξ[AH RFHigh[θ0][θ]+Sum[A[i](θ^(2i)),{i,1,NFN[6]}]/.SN[6][[2]],HN[6],SN[6][[2]]][12, θc/.SN[6][[2]]]

test=dfdξ[Sum[A[i](θ^(2i)),{i,1,NFN[6]}]/.SN[6][[2]]/.SN[6][[2]],HN[6],SN[6][[2]]][11, θc/.SN[6][[2]]]

test=dfdξ[AL fLow[B,θc][-θ]/.rules[HN[6]]/.SN[6][[2]],HN[6],SN[6][[2]]][12, θc/.SN[6][[2]]]

test=dfdξ[RFReg[NFN[6]][θ]/.rules[HN[6]]/.SN[6][[2]],HN[6],SN[6][[2]]][12, θc/.SN[6][[2]]]

dfdξSing[h_,s_][m_]:=Module[{ds,dF,df},
ds=invDerivativeList[m+1][Function[θ,h[θ]/t[θ]^Δ[2]/.s]][θc/.s];
dF=Table[dfLow[AL,B][k],{k,0,m}]/. θ-> θc//.s;
df=NestList[D[#,θ]&,fp[θ]/t[θ]^2,m]/.  Map[Derivative[#][fp][θ]->dF[[#+1]]&,Range[0,m]]/. θ->θc/.s;
Table[Sum[df[[k+1]] BellY[j,k,ds[[;;j-k+1]]],{k,0,j}]/(j!),{j,0,m}]
]

test=dfdξSing[HN[6],Join[rules[HN[6]],SN[6][[2]]]][11]

HN[7][θ]

Plot[θ(1-((θ^2)/((θc)^2)))gC[0]/.θ->θc x,{x,0,1}]

ListLogPlot[ {Abs[test],(((((2 os/π))^(#-1))Gamma[#-1]os/(2π)/π)&/@Range[0,14])}]

ListLogPlot[ Abs@SeriesCoefficient[AH RFHigh[θ0][θ]/.S[7],{θ,θc/.S[7],#}]&/@Range[0,10]]

ListLogPlot[Abs[Gls]]

oursLow2=Chop@Abs@dGdξLowList[NFN[2],HN[2],Join[rules[HN[2]],SN[2][[2]]]][14]
oursLow3=Chop@Abs@dGdξLowList[NFN[3],HN[3],Join[rules[HN[3]],SN[3][[2]]]][14]
oursLow4=Chop@Abs@dGdξLowList[NFN[4],HN[4],Join[rules[HN[4]],SN[4][[2]]]][14]
oursLow5=Chop@Abs@dGdξLowList[NFN[5],HN[5],Join[rules[HN[5]],SN[5][[2]]]][14]
oursLow6=Chop@Abs@dGdξLowList[NFN[6],HN[6],Join[rules[HN[6]],SN[6][[2]]]][14]

ListLogPlot[{
Abs[Abs[Gls]],
Abs[oursLow2],
Abs[oursLow6]
},
PlotRange->All]

madPlot=ListPlot[Rest/@{{{0,0},{0,1},{14,1}},
Abs[Abs[Gls]]/(((((2 os/π))^(#-1))Gamma[#-1]os/(2π)/π)&/@Range[0,14]),
Abs[oursLow2]/(((((2 os/π))^(#-1))Gamma[#-1]os/(2π)/π)&/@Range[0,14]),
Abs[oursLow3]/(((((2 os/π))^(#-1))Gamma[#-1]os/(2π)/π)&/@Range[0,14]),
Abs[oursLow4]/(((((2 os/π))^(#-1))Gamma[#-1]os/(2π)/π)&/@Range[0,14]),
Abs[oursLow5]/(((((2 os/π))^(#-1))Gamma[#-1]os/(2π)/π)&/@Range[0,14]),
Abs[oursLow6]/(((((2 os/π))^(#-1))Gamma[#-1]os/(2π)/π)&/@Range[0,14])
},
Joined->True,PlotRange->{0.6,2},PlotLegends->{1,"Numeric","n = 2","n = 3","n = 4","n = 5","n = 6"},AxesLabel->{m,"|((SubMinus[ℱ])^((m)))| / [((Subscript[A, 0]) / π) \!\(\*SuperscriptBox[\(b\), \(m\\  - \\ 1\)]\) Γ(m - 1)]"}]

Export["/tmp/convergence.pdf",madPlot]

ListPlot[{
Abs[Abs[Gls]/(((((2 os/π))^(#-1))Gamma[#-1]os/(2π)/π)&/@Range[0,14])],
Abs[Chop@Abs@dGdξLowList[NF[2],H[2],Join[rules[H[2]],S[2]]][10]/(((((2 os/π))^(#-1))Gamma[#-1]os/(2π)/π)&/@Range[0,10])],
Abs[Chop@Abs@dGdξLowList[NF[7],H[7],Join[rules[H[7]],S[7]]][10]/(((((2 os/π))^(#-1))Gamma[#-1]os/(2π)/π)&/@Range[0,10])]
}]

ListLogPlot[{
Abs@Ghs/((0.18930^(-#))&/@Range[0,14]),
Chop[Abs@dGdξList[NFN[6],HN[6],Join[rules[HN[6]],SN[6][[2]]]][10,0]]/((0.18930^(-#))&/@Range[0,10])
}]

ListLogPlot[{
Abs[Φs],
Chop[Abs@dΦdηList[NFN[2],HN[2],Join[rules[HN[2]],SN[2][[2]]]][10,1]],
Chop[Abs@dΦdηList[NFN[6],HN[6],Join[rules[HN[6]],SN[6][[2]]]][10,1]]
}]

Show[ListPlot[θ0/θc/.SN[#][[2]]&/@Range[2,6]],Plot[2/3,{x,0,10}]]

ClearAll[θc]

ClearAll[gC]

θc:=3/2

θcW=0;

NFN[2]=2;
HN[2]=h[1];
SN[2]=annealFit[Join[resAll[NFN[2],HN[2],(10^(-5))][2],Table[A[i]((θc)^(2i)),{i,1,NFN[2]}],Table[gC[i]((θc)^(2i-1)),{i,1,1}]],{θc->1.2,θ0->0.4,AH->-2.5,A[1]->1,A[2]->0,gC[0]->1,gC[1]->0},(10^(-1))]

SN[2]=annealFit[Join[resAll[NFN[2],HN[2],(10^(-12))][2],Table[A[i]((θc)^(2i)),{i,1,NFN[2]}],Table[gC[i]((θc)^(2i-1)),{i,1,1}]],SN[2][[2]]]

NFN[3]=2;
HN[3]=h[2];
SN[3]=annealFit[Join[resAll[NFN[3],HN[3],(10^(-10))][3],Table[A[i]((θc)^(2i)),{i,1,NFN[3]}],Table[gC[i]((θc)^(2i-1)),{i,1,2}]],Join[SN[2][[2]],{gC[2]->0}],(10^(-2))]

SN[3]=annealFit[Join[resAll[NFN[3],HN[3],(10^(-12))][3],Table[A[i]((θc)^(2i)),{i,1,NFN[3]}],Table[gC[i]((θc)^(2i-1)),{i,1,2}]],SN[3][[2]]]

NFN[4]=3;
HN[4]=h[3];
SN[4]=annealFit[Join[resAll[NFN[4],HN[4],(10^(-10))][4],Table[A[i]((θc)^(2i)),{i,1,NFN[4]}],Table[gC[i]((θc)^(2i-1)),{i,1,3}]],Join[SN[3][[2]],{gC[3]->0,A[3]->0}],(10^(-3))]

SN[4]=annealFit[Join[resAll[NFN[4],HN[4],(10^(-12))][4],Table[A[i]((θc)^(2i)),{i,1,NFN[4]}],Table[gC[i]((θc)^(2i-1)),{i,1,3}]],SN[4][[2]],(10^(-3))]

NFN[5]=4;
HN[5]=h[3];
SN[5]=annealFit[Join[resAll[NFN[5],HN[5],(10^(-10))][5],Table[A[i]((θc)^(2i)),{i,1,NFN[5]}],Table[gC[i]((θc)^(2i-1)),{i,1,3}]],Join[SN[4][[2]],{A[4]->0}],(10^(-6))]

SN[5]=annealFit[Join[resAll[NFN[5],HN[5],(10^(-12))][5],Table[A[i]((θc)^(2i)),{i,1,NFN[5]}],Table[gC[i]((θc)^(2i-1)),{i,1,3}]],SN[5][[2]],(10^(-6))]

NFN[6]=5;
HN[6]=h[4];
SN[6]=annealFit[Join[resAll[NFN[6],HN[6],(10^(-10))][6],Table[A[i]((θc)^(2i)),{i,1,NFN[6]}],Table[gC[i]((θc)^(2i-1)),{i,1,4}]],Join[SN[5][[2]],{gC[4]->0,A[5]->0}],(10^(-6))]

SN[6]=annealFit[Join[resAll[NFN[6],HN[6],(10^(-12))][6],Table[A[i]((θc)^(2i)),{i,1,NFN[6]}],Table[gC[i]((θc)^(2i-1)),{i,1,4}]],SN[6][[2]],(10^(-6))]

NFN[7]=5;
HN[7]=h[5];
SN[7]=RemoteEvaluate@annealFit[Join[resAll[NFN[7],HN[7],(10^(-8))][7],Table[A[i]i!,{i,1,NFN[7]}],Table[gC[i]i!,{i,1,5}]],Join[SN[6][[2]],{A[5]->0}],(10^(-6))]

SN[7]=RemoteEvaluate@annealFit[Join[resAll[NFN[7],HN[7],(10^(-10))][7],Table[A[i]i!,{i,1,NFN[7]}],Table[gC[i]i!,{i,1,5}]],SN[7][[2]],(10^(-7))]

SN[7]=RemoteEvaluate@annealFit[Join[resAll[NFN[7],HN[7],(10^(-9))][7],Table[A[i]i!,{i,1,NFN[7]}],Table[gC[i]i!,{i,1,5}]],s7tmp[[2]],(10^(-8))]

NFN[8]=5;
HN[8]=h[7];
SN[8]=RemoteEvaluate@annealFit[Join[resAll[NFN[8],HN[8],(10^(-10))][8],Table[A[i]i!,{i,1,NFN[8]}],Table[gC[i]i!,{i,1,7}],θcW{θc-1}],Join[SN[7][[2]],{gC[7]->0,A[5]->0}],(10^(-8))]

SN[8]=RemoteEvaluate@annealFit[Join[resAll[NFN[8],HN[8],(10^(-12))][8],Table[A[i]i!,{i,1,NFN[8]}],Table[gC[i]i!,{i,1,7}],θcW{θc-1}],SN[8][[2]]]

SN[8]=annealFit[Join[resAll[NFN[8],HN[8],(10^(-15))][8],Table[A[i]i!,{i,1,NFN[8]}],Table[gC[i]i!,{i,1,0}]],SN[8][[2]]]

NFN[9]=13;
HN[9]=h[0];
SN[9]=RemoteEvaluate@annealFit[Join[res[NFN[9],HN[9],(10^(-5))][9],Table[A[i]i!,{i,1,NFN[9]}],Table[gC[i]i!,{i,1,0}],θcW{θc-1}],Join[SN[8][[2]],{A[13]->0}],(10^(-9))]

SN[9]=RemoteEvaluate@annealFit[Join[res[NFN[9],HN[9],(10^(-10))][9],Table[A[i]i!,{i,1,NFN[9]}],Table[gC[i]i!,{i,1,0}],θcW{θc-1}],SN[9][[2]]]

NFN[10]=15;
HN[10]=h[0];
SN[10]=RemoteEvaluate@annealFit[Join[res[NFN[10],HN[10],(10^(-5))][10],Table[A[i]i!,{i,1,NFN[10]}],Table[gC[i]i!,{i,1,0}],θcW{θc-1}],Join[SN[9][[2]],{A[14]->0,A[15]->0}],(10^(-10))]

SN[10]=RemoteEvaluate@annealFit[Join[res[NFN[10],HN[10],(10^(-10))][10],Table[A[i]i!,{i,1,NFN[10]}],Table[gC[i]i!,{i,1,0}],θcW{θc-1}],SN[10][[2]]]

NFN[11]=16;
HN[11]=h[0];
SN[11]=RemoteEvaluate@annealFit[Join[res[NFN[11],HN[11],(10^(-5))][11],Table[A[i]i!,{i,1,NFN[11]}],Table[gC[i]i!,{i,1,0}]],Join[SN[10][[2]],{A[16]->0}],(10^(-12))]

SN[11]=RemoteEvaluate@annealFit[Join[res[NFN[11],HN[11],(10^(-10))][11],Table[A[i]i!,{i,1,NFN[11]}],Table[gC[i]i!,{i,1,0}]],SN[11][[2]]]

test=((Φs[[;;3+1]]-#)(*/Φs[[;;ntest+1]]*))&@dΦdηList[NF[#],H[#]][3,1]//.rules[H[#]]/.S[#]&/@Range[2,7];

testN=((Φs[[;;4+1]]-#)(*/Φs[[;;ntest+1]]*))&@dΦdηList[NFN[#],HN[#]][4,1]//.rules[HN[#]]/.SN[#][[2]]&/@Range[2,7];

ListLogPlot[Evaluate[Abs[Transpose@testN]],Joined->True]

Φcomp=Flatten/@({#,Abs@((Φs[[;;4+1]]-#)(*/Φs[[;;ntest+1]]*))&@dΦdηList[NF[#],H[#]][4,1]//.rules[H[#]]/.S[#]}&/@Range[2,6])

Export["~/doc/research/first_order_singularities/paper/data/phi_comparison.dat",Map[If[Head[#]===Around,#["Value"],#]&,Φcomp,2]];

pCov=Show[ListLogPlot[Thread[{Range[2,7],#}]&/@Abs@Re@Transpose[test],Joined->True,
AxesLabel->{"Order","((Subscript[Δℱ, 0])^((n)))"},PlotLegends->LineLegend[{0,1,2,3},LegendLabel->n]],
LogLogPlot[0.1(x^(-4)),{x,1,15},PlotStyle->Black]]

N[Gls[[3]],50]

(10^(-5))//N

Export["/tmp/convergence.pdf",pCov];

Export["~/doc/research/first_order_singularities/paper/data/yl_comparison.dat",((n|->{#,n["Value"],n["Uncertainty"]})@(Abs@ruleθ0[H[#]]/.S[#]/.gC[0]->1))&/@Range[2,7]]

ListLogPlot[(Abs@ruleθ0[H[#]]/.S[#]/.gC[0]->1)&/@Range[2,7],PlotRange->{(10^(-6)),.1}]

ListPlot[Abs@(Re@ruleAH[H[#]]/.S[#]/.gC[0]->1)&/@Range[2,7]]

ListLogLogPlot[(n|->Table[{1/(2i+1),Abs@((SeriesCoefficient[H[n][θ],{θ,0,2i+1}])/(SeriesCoefficient[H[n][θ],{θ,0,2i-1}]))},{i,1,11}]/.S[n]/.gC[0]->1)/@Range[4,7],Joined->True]

ListLogLogPlot[Table[{1/(2i),Abs@((A[i])/(A[i-1]))},{i,1,11}]/.(S/@Range[4,7]),Joined->True]

Plot[Evaluate[HN[#][θ ]/.SN[#][[2]]&/@Range[2,6]],{θ,0,2}]

Plot[Evaluate[ℱ[#][θ ]//.rules[H[#]]/.S[#]&/@Range[2,6]],{θ,0,2}]

ParametricPlot[Evaluate[{-t[θ]((H[#][θ ])^(-1/Δ[])),((H[#][θ])^(-2/Δ[]))(ℱ[#][θ ]-((t[θ])^2)((Log[((H[#][θ])^2)])/(8 π Δ[])))-(((t[θ]((H[#][θ ])^(-1/Δ[]))))^2)Log[(((t[θ]((H[#][θ ])^(-1/Δ[]))))^2)](1/(8π))}//.rules[H[#]]/.S[#]&/@Range[6,7]],{θ,0.00001,1.4},AspectRatio->1/GoldenRatio,PlotRange->{{-50,50},{-2.5,0.5}},PlotPoints->400]

ParametricPlot[Evaluate[{H[#][θ]((Abs[t[θ]])^(-Δ[])),((t[θ])^(-2))ℱ[#][θ ]-((Log[((t[θ])^2)])/(8 π ))}//.rules[H[#]]/.S[#]&/@Range[2,6]],{θ,1,4},AspectRatio->1/GoldenRatio]

ParametricPlot[Evaluate[{H[#][θ]((Abs[t[θ]])^(-Δ[])),((t[θ])^(-2))ℱ[#][θ ]-((Log[((t[θ])^2)])/(8 π ))}//.rules[H[#]]/.S[#]&@7],{θ,-0.9,0.9},AspectRatio->1/GoldenRatio]

complexField=1.4;ComplexPlot[RFC[NF[7]][θ]//.rules[H[7]]/.S[7],{θ,-complexField-complexField/GoldenRatio  I,complexField+complexField/GoldenRatio I},PlotPoints->200,Exclusions->None,ColorFunction->{Automatic,"None"}]

ComplexPlot[Evaluate[RFC[0][θ]/.AH->0/.{B->1,θc->1,AL->1}],{θ,1-0.05-0.05I,1+0.05+0.05I},Exclusions->None,ColorFunction->{"Automatic",None},PlotPoints->100,WorkingPrecision->500]

complexField=0.2;ComplexPlot[RFC[NF[9]][θ]//.rules[H[9]]/.S[9],{θ,θc-complexField-complexField/GoldenRatio  I,θc+complexField+complexField/GoldenRatio I}/.S[9],PlotPoints->100,Exclusions->None,ColorFunction->{Automatic,"None"},WorkingPrecision->200]

Plot[Evaluate[Arg@RFC[0][θ+I (10^(-7))]/.AH->0/.{B->1,θc->1,AL->1}],{θ,0.9,1.1},WorkingPrecision->200,PlotPoints->50]

Stored fits

H[2]=h[1];
NF[2]=2;
ℱ[2]=RF[2];
S[2]={θc->1.148407773492004`,θ0->0.9896669889911205`,AH->-0.172823989504767`,A[1]->-0.31018352388662596`,A[2]->0.2474537923130002`,gC[0]->0.37369093055254343`,gC[1]->-0.021636313152585823`};

H[3]=h[2];
NF[3]=2;
ℱ[3]=RF[2];
S[3]={θc->1.2542120477507488`,θ0->0.6020557328641167`,AH->-0.38566364361428684`,A[1]->-0.3527514794812415`,A[2]->0.2582430860166863`,gC[0]->0.4483788209731592`,gC[1]->-0.022032295172535358`,gC[2]->0.00022200608228654115`};

H[4]=h[3];
NF[4]=3;
ℱ[4]=RF[3];
S[4]={θc->1.3164928721109121`,θ0->0.6400189996493497`,AH->-0.3563974694580203`,A[1]->-0.3550547624920263`,A[2]->0.23465947408509413`,gC[0]->0.4410742751152714`,gC[1]->-0.034817777358116885`,gC[2]->0.000678172648789985`,gC[3]->-0.00004305140578834467`,A[3]->-0.0019083731028066697`};

H[5]=h[4];
NF[5]=3;
ℱ[5]=RF[3];
S[5]={θc->1.3403205742656135`,θ0->0.6238113973493433`,AH->-0.38002950945224295`,A[1]->-0.35127522582179693`,A[2]->0.23704589676915347`,gC[0]->0.44371885415894785`,gC[1]->-0.04609943321005163`,gC[2]->-0.0007458341071947777`,gC[3]->0.00005966875622885447`,A[3]->-0.007319731639727028`,gC[4]->-4.403083529955303`*^-6};

H[6]=h[5];
NF[6]=4;
ℱ[6]=RF[4];
S[6]={θc->1.3626103817690176`,θ0->0.6462147447024515`,AH->-0.35576386594103865`,A[1]->-0.3520586281920383`,A[2]->0.23316561297622435`,gC[0]->0.43845335615925396`,gC[1]->-0.05312704168994819`,gC[2]->-0.003914782631377569`,gC[3]->-0.0004080160912692574`,A[3]->-0.006649030656179257`,gC[4]->0.000026262906640471588`,gC[5]->-1.0974538440746828`*^-6,A[4]->-0.0016899077640685814`};