Much work.

1 files changed, 56 insertions, 23 deletions

diff --git a/schofield.wl b/schofield.wl
index 646881d..8ac7ffc 100644
--- a/schofield.wl
+++ b/schofield.wl
@@ -1,7 +1,7 @@
 
 BeginPackage["Schofield`"]
 
-$Assumptions = {θc > 0, θc > 1, gC[_] ∈ Reals, B > 0, γ > 0}
+$Assumptions = {θc > 0, θc > 1, gC[_] ∈ Reals, B > 0, γ > 0, ξ0 > 0}
 
 β[D_:2] := Piecewise[
   {
@@ -31,28 +31,43 @@ $Assumptions = {θc > 0, θc > 1, gC[_] ∈ Reals, B > 0, γ > 0}
 ]
 
 Δ[D_:2] := β[D] δ[D]
-t[θ_] := (θ)^2 - 1
-h[n_][θ_] := (1 - (θ/θc)^2) Sum[gC[i] θ^(2i+1), {i, 0, n}]
-
-RFf[y_] := π^-1 y Exp[1/y] ExpIntegralEi[-1/y]
-RFnonsingular[n_][θ_] := A RFf[B(θc+θ)] + Sum[FC[i] θ^(2i) , {i, 0, n}]
-RF[n_][θ_] := A RFf[B(θc-θ)] + RFnonsingular[n][θ]
+t[θ_] := ((θ)^2 - 1)
+(*
+hBasis = Orthogonalize[x^# & /@ Range[0, 20], Function[{f, g}, Integrate[f g (1 - x^2)^2, {x, -1, 1}]]]
+h[n_][θ_] := (1 - (θ/θc)^2) Sum[gC[i] hBasis[[2 * i + 2]], {i, 0, n}] /. x -> θ / θc
+*)
+h[n_][θ_] := (1 - (θ/θc)^2) Sum[gC[i]θ^(2 * i + 1), {i, 0, n}]
+
+RFLow[B_, θc_][θ_] := (2 E^(1/(
+   B \[Theta]c)) \[Theta]c ExpIntegralEi[-(1/(B \[Theta]c))] + 
+  E^(1/(B (-\[Theta] + \[Theta]c))) (\[Theta] - \[Theta]c) \
+ExpIntegralEi[1/(B \[Theta] - B \[Theta]c)] - 
+  E^(1/(B \[Theta] + 
+    B \[Theta]c)) (\[Theta] + \[Theta]c) ExpIntegralEi[-(1/(
+     B \[Theta] + B \[Theta]c))])
+RFHigh[ξ0_][ξ_] := (1+(ξ/ξ0)^2)^(5/6) - 1
+
+RF[θ_] := AL RFLow[B, θc][θ] + AH RFHigh[θ0][θ] + AM θ^2
+
+ruleB[g_] := B - π / (2 1.3578383417066) (- g'[θc] / t[θc]^Δ[2])
+ruleθ0[g_] := Simplify[g[I θ0]/(-t[I θ0])^Δ[2]/I] - 0.18930
+
+eqLowRHS[F_][m_] := SeriesCoefficient[F[θ], {θ, θc, m}, Assumptions -> Join[{θ < θc, θ > 1}, $Assumptions]]
+
+eqLowLHS[h_][m_] :=D[
+  t[θ]^2 (Gl[h[θ] t[θ]^-Δ[2]] + Log[t[θ]^2]/(8 π)),
+  {θ, m} ] / m! /. θ -> θc
 
-RFfCoeff[m_] := 0 /; m < 2
-RFfCoeff[m_] := Gamma[m - 1] / π /; m > 1
+eqLow[F_, h_][m_] := eqLowRHS[F][m] - eqLowLHS[h][m]
 
-ruleB[g_] := π / (2 1.3578383417066) (- g'[θc] / t[θc]^Δ[2])
+eqHighRHS[F_][m_] := D[F[θ], {θ, m} ] /. θ -> 0
 
-eqLow[n_, h_][m_] := -A B^m RFfCoeff[m] + D[
-  RFnonsingular[n][θ] - t[θ]^2 (Gl[h[θ] t[θ]^-Δ[2]] + Log[t[θ]^2]/(8 π)),
-  {θ, m} ] / m! /. θ -> θc
+eqHighLHS[h_][m_] := D[(-t[θ])^2 (Gh[h[θ] (-t[θ])^-Δ[2]] + Log[(-t[θ])^2]/(8 π)), {θ, m} ] /. θ -> 0
 
-eqHigh[n_, h_][m_] := D[
-  RF[n][θ] - (-t[θ])^2 (Gh[h[θ] (-t[θ])^-Δ[2]] + Log[(-t[θ])^2]/(8 π)),
-  {θ, m} ] /. θ -> 0
+eqHigh[F_, h_][m_] := eqHighRHS[F][m] - eqHighLHS[h][m]
 
-eqMid[n_, h_][m_] := D[
-  RF[n][θ] - t[θ]^2 Log[h[θ]^2]/(8 Δ[2]π) - h[θ]^((2-α[2])/Δ[2]) Φ[η]
+eqMid[F_, h_][m_] := D[
+  F[θ] - t[θ]^2 Log[h[θ]^2]/(8 Δ[2]π) - h[θ]^((2-α[2])/Δ[2]) Φ[η]
     /. η -> t[θ] / h[θ]^(1 / Δ[2]),
   {θ, m} ] /. θ -> 1
 
@@ -76,7 +91,10 @@ Gls = {
   0.18388370,
   -0.6591714,
   2.937665,
-  -15.61
+  -15.61,
+  96.76,
+  -679,
+  5340
 }
 
 Ghs = {
@@ -90,7 +108,9 @@ Ghs = {
   0,
   1457.62,
   0,
-  -25891
+  -25891,
+  0,
+  502000
 }
 
 dRule[sym_][f_, i_] := Derivative[i[[1]] - 1][sym][0] -> f (i[[1]] - 1)!
@@ -99,9 +119,22 @@ dRule[sym_][f_, i_] := Derivative[i[[1]] - 1][sym][0] -> f (i[[1]] - 1)!
 GlRules = MapIndexed[dRule[Gl], Gls];
 GhRules = MapIndexed[dRule[Gh], Ghs];
 
-rules[g_] := Join[{B -> ruleB[g]}, ΦRules, GlRules, GhRules]
-
-eq[F_, g_][m_] := Flatten[Join[{eqLow[F, g][#], eqMid[F, g][#]} & /@ Range[0, m], eqHigh[F, g] /@ Range[0, m, 2]]] /. rules[g]
+rules[g_] := Join[ΦRules, GlRules, GhRules]
+(*ξ0 := 0.18930*)
+gC[0] := 1
+tC[0] := 1
+
+eq[F_, g_][m_] := Flatten[Join[{ruleB[g], ruleθ0[g]},{eqLow[F, g][#](*, eqMid[F, g][#]*)} & /@ Range[0, m], eqHigh[F, g] /@ Range[2, m, 2]]] //. rules[g]
+
+newSol[eqs_, oldSol_, newVars_, δ_:0, opts___] := FindRoot[
+  eqs,
+  Join[
+    oldSol /. Rule -> List,
+    Thread[{newVars, δ * RandomVariate[NormalDistribution[], Length[newVars]]}]
+  ],
+  MaxIterations -> 10000,
+  opts
+]
 
 EndPackage[]