1 files changed, 176 insertions, 0 deletions
diff --git a/schofield-mult.wl b/schofield-mult.wl
new file mode 100644
index 0000000..330e451
--- /dev/null
+++ b/schofield-mult.wl
@@ -0,0 +1,176 @@
+
+BeginPackage["Schofield`"]
+
+$Assumptions = {θc > 0, θc > 1, gC[_] ∈ Reals, B > 0, γ > 0, ξ0 > 0}
+
+β[D_:2] := Piecewise[
+  {
+    {1/8, D == 2},
+    {0.326419, D == 3},
+    {1/2, D == 4},
+    {β, True}
+  }
+]
+
+δ[D_:2] := Piecewise[
+  {
+    {15, D == 2},
+    {4.78984, D == 3},
+    {3, D == 4},
+    {δ, True}
+  }
+]
+
+α[D_:2] := Piecewise[
+  {
+    {0, D == 2},
+    {0.11008, D == 3},
+    {0, D == 4},
+    {α, True}
+  }
+]
+
+Δ[D_:2] := β[D] δ[D]
+
+OverBar[s] := 1.357838341706595496
+
+t[θ_] := ((θ)^2 - 1)
+h[n_][θ_] := (1 - (θ/θc)^2) Sum[gC[i]LegendreP[(2 * i + 1), θ/θc], {i, 0, n}]
+η[g_][θ_] := t[θ] / (g[θ] / I)^(1 / Δ[2])
+
+RFLow[B_, θc_][θ_] := (1/\[Pi])(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_][ξ_] := (ξ^2+ξ0^2)^(5/6)
+
+RF[n_][θ_] := AL RFLow[B, θc][θ] RFHigh[θ0][θ] + Sum[A[i] LegendreP[(2 i), θ/θc] , {i, 1, n}]
+RFReg[n_][θ_] := AL RFHigh[θ0][θ] (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))]) + Sum[A[i] LegendreP[(2 i), θ/θc], {i, 1, n}]
+dRFc[n_][m_] := AL m! Sum[Piecewise[{{ Gamma[j - 1] B^(j - 1) / π, j>1}, {0, True}}] (D[RFHigh[θ0][θ], {θ, m - j}] / (m - j)! /. θ -> θc), {j, 0, m}] + 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
+ruleAL[g_] := AL RFHigh[θ0][θc] + t[θc]^2 OverBar[s] / (2 π) * (- g'[θc] / t[θc]^Δ[2])
+ruleAH[g_] := AL Re[RFLow[B, θc][θ0 I]]+ 1.37 * (g[I θ0]/ I)^(2 / Δ[2]) * (-η[g]'[I θ0] / (2 θ0 I))^(5/6)
+
+eqLowRHSReg[n_][m_] := dRFc[n][m]
+
+eqLowLHS[h_][m_] :=D[
+  t[θ]^2 (Gl[h[θ] t[θ]^-Δ[2]] + Log[t[θ]^2]/(8 π)),
+  {θ, m} ] /. θ -> θc
+
+eqLow[n_, h_][m_] := (eqLowRHSReg[n][m] - eqLowLHS[h][m]) / m!
+
+eqHighRHS[F_][m_] := D[F[θ], {θ, m} ] /. θ -> 0
+
+eqHighLHS[h_][m_] := D[(-t[θ])^2 (Gh[h[θ] (-t[θ])^-Δ[2]] + Log[(-t[θ])^2]/(8 π)), {θ, m} ] /. θ -> 0
+
+eqHigh[n_, h_][m_] := (eqHighRHS[RF[n]][m] - eqHighLHS[h][m]) / m!
+
+eqMid[F_, h_][m_] := D[
+  F[θ] - t[θ]^2 Log[h[θ]^2]/(8 Δ[2]π) - h[θ]^((2-α[2])/Δ[2]) Φ[η]
+    /. η -> t[θ] / h[θ]^(1 / Δ[2]),
+  {θ, m} ] / m! /. θ -> 1
+
+δ0 = 10^(-14);
+
+Φs = {
+  -1.197733383797993,
+  -0.318810124891,
+  0.110886196683,
+  0.01642689465,
+  -2.639978 10^-4,
+  -5.140526 10^-4,
+  2.08856 10^-4,
+  -4.4819 10^-5,
+  3.16 10^-7,
+  4.31 10^-6,
+  -1.99 10^-6
+}
+
+Gls = {
+  Around[0, δ0],
+  Around[-OverBar[s], δ0],
+  Around[−0.048953289720, 2 10^(-12)],
+  Around[ 0.0388639290, 1 10^(-10)],
+  Around[-0.068362121, 1 10^(-9)],
+  Around[ 0.18388371, 1 10^(-8)],
+  Around[-0.659170, 1 10^(-6)],
+  Around[ 2.937665, 3 10^(-6)],
+  Around[-15.61, 10^(-2)],
+  96.76,
+  -6.79 10^2,
+  5.34 10^3,
+  -4.66 10^4,
+  4.46 10^5,
+  -4.66 10^6
+}
+
+Ghs = {
+  Around[0, δ0],
+  Around[0, δ0],
+  Around[ -1.84522807823, 10^(-11)],
+  Around[0, δ0],
+  Around[  8.3337117508, 10^(-10)],
+  Around[0, δ0],
+  Around[-95.16897, 10^(-5)],
+  Around[0, δ0],
+  Around[1457.62, 3 10^(-2)],
+  0,
+  Around[-2.5891 10^4, 2],
+  0,
+  5.02 10^5,
+  0,
+  -1.04 10^7
+}
+
+dRule[sym_][f_, i_] := Derivative[i[[1]] - 1][sym][0] -> f (i[[1]] - 1)!
+
+ΦRules = MapIndexed[dRule[Φ], Φs];
+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*)
+
+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[2, 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])]
+chiSquaredHigh[n_, g_][m_] := Total[(((#[[1]] /. rules)["Value"] - #[[2]])^2 / (#[[1]] /. rules)["Uncertainty"]^2)& /@ ({Ghs[[#+1]], dFdξHigh[n, g][#] / #!} & /@ Range[0, m])]
+chiSquared[F_, g_][m_] := chiSquaredLow[F, g][m] + chiSquaredHigh[F, g][m] + ruleB[g]^2 / δ0^2 + ruleθ0[g]^2 / 0.00005^2
+
+newSol[eqs_, oldSol_, newVars_, δ_:0, γ_:0, opts___] := FindRoot[
+  eqs,
+  Join[
+    {#1, #2 + γ * RandomVariate[NormalDistribution[]]} & @@@ (oldSol /. Rule -> List),
+    Thread[{newVars, δ * RandomVariate[NormalDistribution[], Length[newVars]]}]
+  ],
+  MaxIterations -> 50000,
+  opts
+]
+
+EndPackage[]
+