Much work.

1 files changed, 16 insertions, 13 deletions

diff --git a/schofield.wl b/schofield.wl
index 7f6bed8..c4ebc32 100644
--- a/schofield.wl
+++ b/schofield.wl
@@ -30,18 +30,12 @@ $Assumptions = {θc > 0, θc > 1, gC[_] ∈ Reals, B > 0, γ > 0, ξ0 > 0}
   }
 ]
 
-OverBar[s] := 1.357838341706595496
-
 Δ[D_:2] := β[D] δ[D]
-t[θ_] := ((θ)^2 - 1)
 
-(*
-hBasis = Orthogonalize[x^# & /@ Range[0, 20], Function[{f, g}, Integrate[f g (1 - x^2)^2, {x, -1, 1}]]]
+OverBar[s] := 1.357838341706595496
 
-h[n_][θ_] := (1 - (θ/θc)^2) Sum[gC[i] hBasis[[2 * i + 2]], {i, 0, n}] /. x -> θ / θc
+t[θ_] := ((θ)^2 - 1)
 h[n_][θ_] := (1 - (θ/θc)^2) Sum[gC[i]LegendreP[(2 * i + 1), θ/θc], {i, 0, n}]
-*)
-h[n_][θ_] := (1 - (θ/θc)^2) Sum[gC[i]θ^(2 * i + 1), {i, 0, n}]
 
 RFLow[B_, θc_][θ_] := (1/\[Pi])(2 E^(1/(
    B \[Theta]c)) \[Theta]c ExpIntegralEi[-(1/(B \[Theta]c))] + 
@@ -52,19 +46,24 @@ ExpIntegralEi[1/(B \[Theta] - B \[Theta]c)] -
      B \[Theta] + B \[Theta]c))])
 RFHigh[ξ0_][ξ_] := (ξ^2+ξ0^2)^(5/6)
 
-RF[n_][θ_] := AL RFLow[B, θc][θ] + AH RFHigh[θ0][θ] + Sum[A[i] θ^(2 i), {i, 0, n}]
+RF[n_][θ_] := AL RFLow[B, θc][θ] + AH RFHigh[θ0][θ] + Sum[A[i] LegendreP[(2 i), θ/θc] , {i, 0, n}]
 RFReg[n_][θ_] := AL (1/\[Pi])(2 E^(1/(
    B \[Theta]c)) \[Theta]c ExpIntegralEi[-(1/(B \[Theta]c))] - 
   E^(1/(B \[Theta] + 
     B \[Theta]c)) (\[Theta] + \[Theta]c) ExpIntegralEi[-(1/(
-     B \[Theta] + B \[Theta]c))]) + AH RFHigh[θ0][θ] + Sum[A[i] LegendreP[2i, θ], {i, 0, n}]
+     B \[Theta] + B \[Theta]c))]) + AH RFHigh[θ0][θ] + Sum[A[i] LegendreP[(2 i), θ/θc], {i, 0, n}]
 dRFc[n_][m_] := Piecewise[{{AL m! Gamma[m - 1] B^(m - 1) / π, m>1}, {0, True}}] + D[RFReg[n][θ], {θ, m}] /. θ -> θc
 
 RFC[n_][θ_] := RF[n][θ] + AL I Sign[Im[θ]] ((θ-θc)Exp[-1/(B(θ-θc))]-(-θ-θc)Exp[-1/(B(-θ-θc))])
 
 ddξ[h_][f_] := D[f, θ] / D[h[θ] / RealAbs[t[θ]]^Δ[2], θ]
+ddη[h_][f_] := D[f, θ] / D[t[θ] / h[θ]^(1 / Δ[2]), θ]
 dFdξLow[n_, h_][m_] := Module[{ff, hh}, Nest[ddξ[hh], ff[θ] / t[θ]^2 - Log[t[θ]^2] / (8 π), m] /. θ -> θc /. Map[Derivative[#][ff][θc] -> dRFc[n][#] &, Range[0, m]] /. hh -> h]
 dFdξHigh[n_, h_][m_] := Module[{ff, hh}, Nest[ddξ[hh], ff[θ] / t[θ]^2 - Log[t[θ]^2] / (8 π), m] /. θ -> 0 /. Map[Derivative[#][ff][0] -> eqHighRHS[RF[n]][#] &, Range[0, m]] /. hh -> h]
+dFdη[n_, h_][m_][tt_] := Module[{ff, hh}, Nest[ddη[hh], h[θ]^(-2 / Δ[]) (ff[θ] - t[θ]^2 Log[hh[θ]^2] / (8 π Δ[])), m] /. θ -> tt /. Map[Derivative[#][ff][tt] -> Derivative[#][RF[n]][tt] &, Range[0, m]] /. hh -> h]
+dFdξLowList[n_, h_][m_] := Module[{ff, hh}, NestList[ddξ[hh], ff[θ] / t[θ]^2 - Log[t[θ]^2] / (8 π), m] /. θ -> θc /. Map[Derivative[#][ff][θc] -> dRFc[n][#] &, Range[0, m]] /. Map[Derivative[#][hh][θc] -> Derivative[#][h][θc] &, Range[0, m]]]
+dFdξHighList[n_, h_][m_] := Module[{ff, hh}, NestList[ddξ[hh], ff[θ] / t[θ]^2 - Log[t[θ]^2] / (8 π), m] /. θ -> 0 /. Map[Derivative[#][ff][0] -> eqHighRHS[RF[n]][#] &, Range[0, m]] /. hh -> h]
+dFdηList[n_, h_][m_][tt_] := Module[{ff, hh}, NestList[ddη[hh], h[θ]^(-2 / Δ[2]) (ff[θ] - t[θ]^2 Log[hh[θ]^2] / (8 π Δ[2])), m] /. θ -> tt /. Map[Derivative[#][ff][tt] -> Derivative[#][RF[n]][tt] &, Range[0, m]] /. hh -> h]
 
 ruleB[g_] := B - (2 * OverBar[s] / π) * (- g'[θc] / t[θc]^Δ[2])
 ruleθ0[g_] := Simplify[g[I θ0]/(-t[I θ0])^Δ[2]/I] - 0.18930
@@ -99,7 +98,10 @@ eqMid[F_, h_][m_] := D[
   -2.639978 10^-4,
   -5.140526 10^-4,
   2.08856 10^-4,
-  -4.4819 10^-5
+  -4.4819 10^-5,
+  3.16 10^-7,
+  4.31 10^-6,
+  -1.99 10^-6
 }
 
 Gls = {
@@ -139,13 +141,14 @@ dRule[sym_][f_, i_] := Derivative[i[[1]] - 1][sym][0] -> f (i[[1]] - 1)!
 GlRules = MapIndexed[dRule[Gl], Gls];
 GhRules = MapIndexed[dRule[Gh], Ghs];
 
+ClearAll[gC]
 rules := Join[ΦRules, GlRules, GhRules]
 (*ξ0 := 0.18930*)
 (*gC[0] := 1*)
 tC[0] := 1
-gC[0] := 1
+(*gC[0] := 1*)
 
-eq[n_, g_][m_] := Flatten[Join[{ruleB[g], ruleθ0[g]},{eqLow[n, g][#](*, eqMid[F, g][#]*)} & /@ Range[0, m], eqHigh[n, g] /@ Range[0, m, 2]]] //. rules /. Around[x_, _] :> x
+eq[n_, g_][m_, p_, q_] := Flatten[Join[{ruleB[g], ruleθ0[g], g'[0] - 1}, eqLow[n, g][#] & /@ Range[0, m],eqMid[RF[n], g][#] & /@ Range[0, p], eqHigh[n, g] /@ Range[0, q, 2]]] //. rules /. Around[x_, _] :> x
 
  (* *)
 chiSquaredLow[n_, g_][m_] := Total[(((#[[1]] /. rules)["Value"] - #[[2]])^2 / (#[[1]] /. rules)["Uncertainty"]^2)& /@ ({Gls[[#+1]], dFdξLow[n, g][#] / #!} & /@ Range[0, m])]