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

Definitions

F1=RF[1];
h1=h[2];
eq1=eq[F1,h1][4];

sol1=newSol[eq1,{θc->1+RandomReal[],θ0->0.2,B->RandomReal[]},{AH,AL,A[0],A[1],gC[1],gC[2]},1]

GraphicsRow[{
Plot[h1[θ]/.sol1,{θ,-θc/.sol1,θc/.sol1}],
Plot[F1[θ]/.sol1,{θ,-θc/.sol1,θc/.sol1}]
}]

GraphicsRow[{
ListLogPlot[Evaluate[Transpose[-{eqLowRHS[F1][#],eqLowLHS[h1][#]}&/@Range[6]]/.rules[h1]/.sol1],PlotRange->All],
ListLogPlot[Chop@Evaluate[Transpose[Abs@{eqHighRHS[F1][#],eqHighLHS[h1][#]}&/@Range[11]]/.rules[h1]/.sol1],PlotRange->All]
}]

F2=RF[2];
h2=h[4];
eq2=eq[F2,h2][6];

sol2=newSol[eq2,sol1,{A[2],gC[3],gC[4]},.1]

GraphicsRow[{
Plot[h2[θ]/.sol2,{θ,0,θc/.sol2}],
Plot[F2[θ]/.sol2,{θ,0,θc/.sol2}]
}]

GraphicsRow[{
ListLogPlot[Evaluate[Transpose[-{eqLowRHS[F2][#],eqLowLHS[h2][#]}&/@Range[8]]/.rules[h2]/.sol2],PlotRange->All],
ListLogPlot[Chop@Evaluate[Transpose[-{eqHighRHS[F2][#],eqHighLHS[h2][#]}&/@Range[11]]/.rules[h2]/.sol2],PlotRange->All]
}]

h3=h[5];
F3=RF[2];
eq3=eq[F3,h3][7];

sol3=newSol[eq3,sol2,{gC[5]},0.01]

GraphicsRow[{
Plot[h3[θ]/.sol3,{θ,0,θc/.sol3}],
Plot[F3[θ]/.sol3,{θ,0,θc/.sol3}]
}]

GraphicsRow[{
ListLogPlot[Evaluate[Transpose[-{eqLowRHS[F3][#],eqLowLHS[h3][#]}&/@Range[8]]/.rules[h3]/.sol3],PlotRange->All],
ListLogPlot[Evaluate[Transpose[-{eqHighRHS[F3][#],eqHighLHS[h3][#]}&/@Range[11]]/.rules[h3]/.sol3],PlotRange->All]
}]

h4=h[6];
F4=RF[3];
eq4=eq[F4,h4][8];

sol4=newSol[eq4,sol3,{gC[6],A[3]},0.01]

GraphicsRow[{
Plot[h4[θ]/.sol4,{θ,0,θc/.sol4}],
Plot[F4[θ]/.sol4,{θ,0,θc/.sol4}]
}]

GraphicsRow[{
ListLogPlot[Evaluate[Transpose[-{eqLowRHS[F4][#],eqLowLHS[t4,h4][#]}&/@Range[11]]/.rules[t4,h4]/.sol4],PlotRange->All],
ListLogPlot[Evaluate[Transpose[-{eqHighRHS[F4][#],eqHighLHS[t4,h4][#]}&/@Range[11]]/.rules[t4,h4]/.sol4],PlotRange->All]
}]

h5=h[11];
F5=RF;
eq5=eq[F5,h5][8];

sol5=newSol[eq5,sol4,{gC[9],gC[10],gC[11]},0.0001]

GraphicsRow[{
Plot[h4[θ]/.sol4,{θ,0,θc/.sol4}],
Plot[F4[θ]/.B->ruleB[t4,h4]/.sol4,{θ,0,θc/.sol4}]
}]

h5=h[6];
F5=RF[6];
eq5=eq[F5,h5][10];

sol5=newSol[eq5,sol4,{gC[6],FC[6]},0.000001]

GraphicsRow[{
Plot[h5[θ]/.sol5,{θ,0,θc/.sol5}],
Plot[F5[θ]/.B->ruleB[h5]/.sol5,{θ,0,θc/.sol5}]
}]

ListLogPlot[Abs[gC/@Range[10]]/.sol5]

h6=h[6];
F6=RF[7];
eq6=eq[7,h6][10];

sol6=newSol[eq6,sol5,{FC[7],FC[6]},(10^(-4))]

GraphicsRow[{
Plot[h4[θ]/.sol4,{θ,0,θc/.sol4}],
Plot[F4[θ]/.B->ruleB[h4]/.sol4,{θ,0,θc/.sol4}]
}]

ListLogPlot[Evaluate[-Transpose[{eqLowRHS[3,h3][#],eqLowLHS[3,h3][#]}&/@Range[12]]/.rules[h3]/.sol3],PlotRange->All]

Limit[D[RF[1][θ],{θ,2}],θ->θc]//Expand

eqLowRHS[1,h1][1]

ParametricPlot[Evaluate[(R|->{R t[θ],(R^(Δ[2]))h2[θ]})/@{0.25,0.375,0.5,0.75,1}/.sol2/.θc->θc1],{θ,-θc1,θc1}]

Φη2=Table[{t[θ]((h2[θ])^(-1/Δ[2])),((h2[θ])^(-2/Δ[2]))(F2[θ]-(((t[θ])^2)/(8π))Log[((t[θ])^2)])}/.B->ruleB[h2]/.sol2,{θ,0.0001,θc-0.0001/.sol2,0.0001}];

ListPlot[Φη2,Joined->True,PlotRange->{{-20,20},{-2,0}}]

Eθ2=Table[{θ,-(Transpose[Inverse[D[{R t[θp],(R^(Δ[2]))h2[θp]},{{R,θp}}]/.θp->θ]].(D[(R^2)(F2[θp]+((t[θp])^2)Log[(R^2)]/(8π)),{{R,θp}}]/.θp->θ))[[1]]/(((R ))^1)}/.B->ruleB[h2]/.sol2/.θc->θc1/.R->1,{θ,0,θc1-0.001,0.001}];

EMθ2=Table[{θ,-(Transpose[Inverse[D[{R t[θp],(R^(Δ[2]))h2[θp]},{{R,θp}}]/.θp->θ]].(D[(R^2)(F2[θp]+((t[θp])^2)Log[(R^2)]/(8π)),{{R,θp}}]/.θp->θ))[[2]],-(Transpose[Inverse[D[{R t[θp],(R^(Δ[2]))h2[θp]},{{R,θp}}]/.θp->θ]].(D[(R^2)(F2[θp]+((t[θp])^2)Log[(R^2)]/(8π)),{{R,θp}}]/.θp->θ))[[1]]}/.B->ruleB[h2]/.sol2/.θc->θc1/.R->1,{θ,-θc1,θc1,0.01}];

Graphics3D[{Black,Sphere[{0,0,0},0.1]}]

Show[
ListPointPlot3D[Flatten[Table[(({R (t[#1]+#3),(R^(Δ[2]))h2[#1]+(R^(1/8))#2,(R^(Δ[2]))h2[#1]-(R^(1/8))#2}/.sol2/.θc->θc1)&@@@EMθ2),{R,0,5,0.1}],1],AxesLabel->{(SubPlus[T]),(SubPlus[H]),(SubMinus[H])}]
]

ListPlot3D[Flatten[(({R t[#1]-R #3,(R^(Δ[2]))h2[#1]-(R^(Δ[2]))#2,(R^(Δ[2]))h2[#1]+(R^(Δ[2]))#2}/.sol2/.θc->θc1)&@@@EMθ2)/.R->Range[0.01,1,0.01],{1}]]

ListPlot3D[{,PlotRange->All]

Mξ2=Table[{h2[θ]((t[θ])^(-Δ[2])),-(Transpose[Inverse[D[{R t[θp],(R^(Δ[2]))h2[θp]},{{R,θp}}]/.θp->θ]].(D[(R^2)(F2[θp]+((t[θp])^2)Log[(R^2)]/(8π)),{{R,θp}}]/.θp->θ))[[2]]/(((R t[θ]))^(1/8))}/.B->ruleB[h2]/.sol2/.θc->θc1/.R->1,{θ,1+0.001,θc1-0.001,0.001}];

ListPlot[Mξ2,Joined->True]

χξ2=Table[{h2[θ]((t[θ])^(-Δ[2])),-(Transpose[Inverse[D[{R t[θp],(R^(Δ[2]))h2[θp]},{{R,θp}}]]].D[(Transpose[Inverse[D[{R t[θp],(R^(Δ[2]))h2[θp]},{{R,θp}}]]].D[(R^2)(F2[θp]+((t[θp])^2)Log[(R^2)]/(8π)),{{R,θp}}])[[2]],{{R,θp}}])[[2]](((R t[θp]))^(7/4))}/.θp->θ/.B->ruleB[h2]/.sol2/.θc->θc1/.R->1,{θ,1+0.001,θc1-0.001,0.001}];

ListPlot[χξ2,Joined->True]

h3=h[2];
F3=RF[2];
eq3=eq[2,h3][4];

sol3=newSol[eq3,sol2,{gC[1],gC[2],FC[2]}]

GraphicsRow[{
Plot[h3[θ]/.sol3,{θ,0,θc/.sol3}],
Plot[F3[θ]/.B->ruleB[h3]/.sol3,{θ,0,θc/.sol3}]
}]

h4=h[3];
F4=RF[2];
eq4=eq[2,h4][5];

sol4=newSol[eq4,sol3,{gC[3]},1]

GraphicsRow[{
Plot[h4[θ]/.sol4,{θ,0,θc/.sol4}],
Plot[F4[θ]/.B->ruleB[h4]/.sol4,{θ,0,θc/.sol4}]
}]

h5=h[4];
F5=RF[3];
eq5=eq[3,h5][6];

sol5=FindRoot[eq5,{{θc,θc/.sol4},{A,A/.sol4},{gC[0],gC[0]/.sol4},{FC[0],FC[0]/.sol4},{FC[1],FC[1]/.sol4},{FC[2],FC[2]/.sol4},{gC[1],gC[1]/.sol4},{gC[2],gC[2]/.sol4},{gC[3],gC[3]/.sol4},{gC[4],0.1RandomReal[{-1,1}]},{FC[3],0.1RandomReal[{-1,1}]}},MaxIterations->50000]

GraphicsRow[{
Plot[h5[θ]/.sol5,{θ,0,θc/.sol5}],
Plot[F5[θ]/.B->ruleB[h5]/.sol5,{θ,0,θc/.sol5}]
}]

h6=h[5];
F6=RF[3];
eq6=eq[3,h6][7];

sol6=FindRoot[eq6,{{θc,θc/.sol5},{A,A/.sol5},{gC[0],gC[0]/.sol5},{FC[0],FC[0]/.sol5},{FC[1],FC[1]/.sol5},{FC[2],FC[2]/.sol5},{gC[1],gC[1]/.sol5},{gC[2],gC[2]/.sol5},{gC[3],gC[3]/.sol5},{gC[4],gC[4]/.sol5},{FC[3],FC[3]/.sol5},{gC[5],0.1 RandomReal[{-1,1}]}},MaxIterations->50000]

GraphicsRow[{
Plot[h6[θ]/.sol6,{θ,0,θc/.sol6}],
Plot[F6[θ]/.B->ruleB[h6]/.sol6,{θ,0,θc/.sol6}]
}]

h7=h[6];
F7=RF[4];
eq7=eq[4,h7][8];

sol7=FindRoot[eq7,{{θc,θc/.sol6},{A,A/.sol6},{gC[0],gC[0]/.sol6},{FC[0],FC[0]/.sol6},{FC[1],FC[1]/.sol6},{FC[2],FC[2]/.sol6},{gC[1],gC[1]/.sol6},{gC[2],gC[2]/.sol6},{gC[3],gC[3]/.sol6},{gC[4],gC[4]/.sol6},{FC[3],FC[3]/.sol6},{gC[5],gC[5]/.sol6},{gC[6],0.01 RandomReal[{-1,1}]},{FC[4],0.01 RandomReal[{-1,1}]}},MaxIterations->50000]

GraphicsRow[{
Plot[h7[θ]/.sol7,{θ,0,θc/.sol7}],
Plot[F7[θ]/.B->ruleB[h7]/.sol7,{θ,0,θc/.sol7}]
}]

h8=h[7];
F8=RF[4];
eq8=eq[4,h8][9];

sol8=FindRoot[eq8,{{θc,θc/.sol7},{A,A/.sol7},{gC[0],gC[0]/.sol7},{FC[0],FC[0]/.sol7},{FC[1],FC[1]/.sol7},{FC[2],FC[2]/.sol7},{gC[1],gC[1]/.sol7},{gC[2],gC[2]/.sol7},{gC[3],gC[3]/.sol7},{gC[4],gC[4]/.sol7},{FC[3],FC[3]/.sol7},{gC[5],gC[5]/.sol7},{gC[6],gC[6]/.sol7},{FC[4],FC[4]/.sol7},{gC[7],0.01RandomReal[{-1,1}]}},MaxIterations->50000]

GraphicsRow[{
Plot[h8[θ]/.sol8,{θ,0,θc/.sol8}],
Plot[F8[θ]/.B->ruleB[h8]/.sol8,{θ,0,θc/.sol8}]
}]

h9=h[8];
F9=RF[5];
eq9=eq[5,h9][10];

sol9=FindRoot[eq9,{{θc,θc/.sol8},{A,A/.sol8},{gC[0],gC[0]/.sol8},{FC[0],FC[0]/.sol8},{FC[1],FC[1]/.sol8},{FC[2],FC[2]/.sol8},{gC[1],gC[1]/.sol8},{gC[2],gC[2]/.sol8},{gC[3],gC[3]/.sol8},{gC[4],gC[4]/.sol8},{FC[3],FC[3]/.sol8},{gC[5],gC[5]/.sol8},{gC[6],gC[6]/.sol8},{FC[4],FC[4]/.sol8},{gC[7],gC[7]/.sol8},{gC[8],0.001 RandomReal[{-1,1}]},{FC[5],0.001 RandomReal[{-1,1}]}},MaxIterations->50000,Method->{"Newton", "StepControl" -> "TrustRegion"}]

GraphicsRow[{
Plot[h9[θ]/.sol9,{θ,0,θc/.sol9}],
Plot[F9[θ]/.B->ruleB[h8]/.sol9,{θ,0,θc/.sol9}]
}]

ListLogPlot[Abs[gC/@Range[9]]/.sol9]

Φη3=Table[{t[θ]((h3[θ])^(-1/Δ[2])),((h3[θ])^(-2/Δ[2]))(F3[θ]-(((t[θ])^2)/(8π))Log[((t[θ])^2)])}/.B->ruleB[h3]/.sol3/.θc->θc1,{θ,0.0001,θc1-0.0001,0.0001}];

Mξ3=Table[{h3[θ]((t[θ])^(-Δ[2])),-(Transpose[Inverse[D[{R t[θp],(R^(Δ[2]))h3[θp]},{{R,θp}}]/.θp->θ]].(D[(R^2)(F3[θp]+((t[θp])^2)Log[(R^2)]/(8π)),{{R,θp}}]/.θp->θ))[[2]]/(((R t[θ]))^(1/8))}/.B->ruleB[h3]/.sol3/.θc->θc1/.R->1,{θ,1+0.001,θc1-0.001,0.001}];

χξ3=Table[{h3[θ]((t[θ])^(-Δ[2])),-(Transpose[Inverse[D[{R t[θp],(R^(Δ[2]))h3[θp]},{{R,θp}}]]].D[(Transpose[Inverse[D[{R t[θp],(R^(Δ[2]))h3[θp]},{{R,θp}}]]].D[(R^2)(F3[θp]+((t[θp])^2)Log[(R^2)]/(8π)),{{R,θp}}])[[2]],{{R,θp}}])[[2]](((R t[θp]))^(7/4))}/.θp->θ/.B->ruleB[h3]/.sol3/.θc->θc1/.R->1,{θ,1+0.001,θc1-0.001,0.001}];

χ3ξ3=Table[{h3[θ]((t[θ])^(-Δ[2])),(Transpose[Inverse[D[{R t[θp],(R^(Δ[2]))h3[θp]},{{R,θp}}]]].D[(Transpose[Inverse[D[{R t[θp],(R^(Δ[2]))h3[θp]},{{R,θp}}]]].D[(Transpose[Inverse[D[{R t[θp],(R^(Δ[2]))h3[θp]},{{R,θp}}]]].D[(R^2)(F3[θp]+((t[θp])^2)Log[(R^2)]/(8π)),{{R,θp}}])[[2]],{{R,θp}}])[[2]],{{R,θp}}])[[2]](((R t[θp]))^(7/4+Δ[2]))}/.θp->θ/.B->ruleB[h3]/.sol3/.θc->θc1/.R->1,{θ,1+0.001,θc1-0.001,0.001}];

ListPlot[{Φη2,Φη3},Joined->True,PlotRange->{{-20,20},{-2,0}}]

ListPlot[{Mξ2,Mξ3},Joined->True]

ListLogLogPlot[{χξ2,χξ3},Joined->True]

h4=h[7];
F4=RF[3];
eq4=eq[3,h4][4];

sol4=FindRoot[eq4/.θc->θc1,{{A,A/.sol3},{gC[0],gC[0]/.sol3},{FC[0],FC[0]/.sol3},{FC[1],FC[1]/.sol3},{FC[2],FC[2]/.sol3},{gC[1],gC[1]/.sol3},{gC[2],gC[2]/.sol3},{gC[3],gC[3]/.sol3},{gC[4],gC[4]/.sol3},{gC[5],gC[5]/.sol3},{gC[6],0.1RandomReal[]},{gC[7],0.1RandomReal[]},{FC[3],0.1RandomReal[]}},MaxIterations->500]

GraphicsRow[{
Plot[h4[θ]/.sol4/.θc->θc1,{θ,0,θc1}],
Plot[F4[θ]/.B->ruleB[h4]/.sol4/.θc->θc1,{θ,0,θc1}]
}]

ListLogPlot[Abs[gC/@Range[9]/.sol4]]

h5=h[8];
F5=RF[4];
eq5=eq[4,h5][5];

sol5=FindRoot[eq5/.θc->θc1,{{A,A/.sol4},{gC[0],gC[0]/.sol4},{FC[0],FC[0]/.sol4},{FC[1],FC[1]/.sol4},{FC[2],FC[2]/.sol4},{gC[1],gC[1]/.sol4},{gC[2],gC[2]/.sol4},{gC[3],gC[3]/.sol4},{gC[4],gC[4]/.sol4},{gC[5],gC[5]/.sol4},{gC[6],gC[6]/.sol4},{gC[7],gC[7]/.sol4},{FC[3],FC[3]/.sol4},{gC[8],0.1RandomReal[]},{FC[4],0.1RandomReal[]}},MaxIterations->5000]

GraphicsRow[{
Plot[h5[θ]/.sol5/.θc->θc1,{θ,0,θc1}],
Plot[F5[θ]/.B->ruleB[h5]/.sol5/.θc->θc1,{θ,0,θc1}]
}]

ListLogPlot[Abs[gC/@Range[9]/.sol5]]

h6=h[10];
F6=RF[5];
eq6=eq[5,h6][6];

sol6=FindRoot[eq6/.θc->θc1,{{A,A/.sol5},{gC[0],gC[0]/.sol5},{FC[0],FC[0]/.sol5},{FC[1],FC[1]/.sol5},{FC[2],FC[2]/.sol5},{gC[1],gC[1]/.sol5},{gC[2],gC[2]/.sol5},{gC[3],gC[3]/.sol5},{gC[4],gC[4]/.sol5},{gC[5],gC[5]/.sol5},{gC[6],gC[6]/.sol5},{gC[7],gC[7]/.sol5},{FC[3],FC[3]/.sol5},{gC[8],gC[8]/.sol5},{FC[4],FC[4]/.sol5},{gC[9],(10^(-3))RandomReal[{-1,1}]},{gC[10],(10^(-3))RandomReal[{-1,1}]},{FC[5],(10^(-3))RandomReal[{-1,1}]}},MaxIterations->1500]

GraphicsRow[{
Plot[h6[θ]/.sol6/.θc->θc1,{θ,0,θc1}],
Plot[F6[θ]/.B->ruleB[h6]/.sol6/.θc->θc1,{θ,0,θc1}]
}]

ListLogPlot[Abs[gC/@Range[10]/.sol6]]

h7=h[12];
F7=RF[5];
eq7=eq[5,h7][7];

sol7=FindRoot[eq7/.θc->θc1,{{A,A/.sol6},{gC[0],gC[0]/.sol6},{FC[0],FC[0]/.sol6},{FC[1],FC[1]/.sol6},{FC[2],FC[2]/.sol6},{gC[1],gC[1]/.sol6},{gC[2],gC[2]/.sol6},{gC[3],gC[3]/.sol6},{gC[4],gC[4]/.sol6},{gC[5],gC[5]/.sol6},{gC[6],gC[6]/.sol6},{gC[7],gC[7]/.sol6},{FC[3],FC[3]/.sol6},{gC[8],gC[8]/.sol6},{FC[4],FC[4]/.sol6},{gC[9],gC[9]/.sol6},{gC[10],gC[10]/.sol6},{FC[5],FC[5]/.sol6},{gC[11],0.1RandomReal[{-1,1}]},{gC[12],0.1RandomReal[{-1,1}]}},MaxIterations->1500]

GraphicsRow[{
Plot[h7[θ]/.sol7/.θc->θc1,{θ,0,θc1}],
Plot[F7[θ]/.B->ruleB[h7]/.sol7/.θc->θc1,{θ,0,θc1}]
}]

ListLogPlot[Abs[gC/@Range[14]/.sol7]]

Caselle comparison

hCaselle[θ_]:=h0 θ(1-((θ^2)/1.16951))(1-0.222389 (θ^2)-0.043547 (θ^4)-0.014809 (θ^6)-0.007168 (θ^8))

m0Sol=Solve[-Gls[[2]]==m0 θ/((((θ^2)-1))^(β[2])) /.θ->Sqrt[1.16951],m0][[1]]

(Transpose[Inverse[D[{R t[θp],(R^(Δ[2]))hCaselle[θp]},{{R,θp}}]/.θp->θ]].(D[m0 (R^(β[2]))θp,{{R,θp}}]/.θp->θ))[[2]]/(R^(-7/4))/.θ->Sqrt[1.16951]/.m0Sol

h0Sol=Solve[-2 Gls[[3]]==(Transpose[Inverse[D[{R t[θp],(R^(Δ[2]))hCaselle[θp]},{{R,θp}}]/.θp->θ]].(D[m0 (R^(β[2]))θp,{{R,θp}}]/.θp->θ))[[2]]/(((R t[θ]))^(-7/4))/.θ->Sqrt[1.16951]/.m0Sol,h0][[1]]

MCaselle=Table[{hCaselle[θ]((t[θ])^(-Δ[2])),m0 θ/((((θ^2)-1))^(β[2]))}/.m0Sol/.h0Sol,{θ,1.0001,Sqrt[1.16951],0.0001}];

ListPlot[{Mξ3,MCaselle},Joined->True]

χCaselle=Table[{hCaselle[θ]((t[θ])^(-Δ[2])),(Transpose[Inverse[D[{R t[θp],(R^(Δ[2]))hCaselle[θp]},{{R,θp}}]/.θp->θ]].(D[m0 (R^(β[2]))θp,{{R,θp}}]/.θp->θ))[[2]]/(((R t[θ]))^(-7/4))}/.m0Sol/.h0Sol,{θ,1.0001,Sqrt[1.16951],0.0001}];

ListPlot[{χξ3,χCaselle},Joined->True]

χ3Caselle=Table[{hCaselle[θ]((t[θ])^(-Δ[2])),-(Transpose[Inverse[D[{R t[θp],(R^(Δ[2]))hCaselle[θp]},{{R,θp}}]]].D[(Transpose[Inverse[D[{R t[θp],(R^(Δ[2]))hCaselle[θp]},{{R,θp}}]/.θp->θ]].(D[m0 (R^(β[2]))θp,{{R,θp}}]))[[2]],{{R,θp}}])[[2]](((R t[θ]))^(7/4+Δ[2]))/.θp->θ}/.m0Sol/.h0Sol,{θ,1.0001,Sqrt[1.16951],0.0001}];

ListLogLogPlot[{χ3ξ3,χ3Caselle},Joined->True]