From 11fab12f99701c8207e509a86f3e898f7411880d Mon Sep 17 00:00:00 2001
From: Jaron Kent-Dobias <jaron@kent-dobias.com>
Date: Thu, 2 Sep 2021 18:22:57 +0200
Subject: Fixed many errors.

---
 schofield.wl | 18 ++++++++++--------
 1 file changed, 10 insertions(+), 8 deletions(-)

(limited to 'schofield.wl')

diff --git a/schofield.wl b/schofield.wl
index f9bae8d..83cf6ac 100644
--- a/schofield.wl
+++ b/schofield.wl
@@ -36,6 +36,7 @@ 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))] + 
@@ -67,7 +68,8 @@ dFdηList[n_, h_][m_][tt_] := Module[{ff, hh}, NestList[ddη[hh], h[θ]^(-2 / Δ
 
 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 + t[θc]^2 OverBar[s] / (2 π) * (- g'[θc] / t[θc]^Δ[2])
+ruleAL[g_] := AL - Exp[Δ[2] t[θc]^(Δ[2] - 1) t'[θc] / (2 OverBar[s] / π g'[θc]) - t[θc]^Δ[2] g''[θc] / (4 OverBar[s] / π g'[θc]^2)] t[θc]^(1/8) OverBar[s] / (2 π) * g'[θc]
+ruleAH[g_] := AH + 1.37 * (g[I θ0]/ I)^(2 / Δ[2]) * (-η[g]'[I θ0] / (2 θ0 I))^(5/6)
 
 eqLowRHSReg[n_][m_] := dRFc[n][m]
 
@@ -91,8 +93,8 @@ eqMid[F_, h_][m_] := D[
 δ0 = 10^(-14);
 
 Φs = {
-  -1.197733383797993,
-  -0.318810124891,
+  -1.19773338379799339,
+  -0.31881012489061,
   0.110886196683,
   0.01642689465,
   -2.639978 10^-4,
@@ -106,8 +108,8 @@ eqMid[F_, h_][m_] := D[
 
 Gls = {
   Around[0, δ0],
-  Around[-OverBar[s], δ0],
-  Around[−0.048953289720, 2 10^(-12)],
+  -OverBar[s],
+  −0.0489532897203,
   Around[ 0.0388639290, 1 10^(-10)],
   Around[-0.068362121, 1 10^(-9)],
   Around[ 0.18388371, 1 10^(-8)],
@@ -125,7 +127,7 @@ Gls = {
 Ghs = {
   Around[0, δ0],
   Around[0, δ0],
-  Around[ -1.84522807823, 10^(-11)],
+  -1.845228078232838,
   Around[0, δ0],
   Around[  8.3337117508, 10^(-10)],
   Around[0, δ0],
@@ -153,7 +155,7 @@ rules := Join[ΦRules, GlRules, GhRules]
 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[0, q, 2]]] //. rules /. Around[x_, _] :> x
+eq[n_, g_][m_, p_, q_] := Flatten[Join[{ruleB[g], ruleθ0[g], g'[0] - 1, ruleAH[g], ruleAL[g]}, 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])]
@@ -166,7 +168,7 @@ newSol[eqs_, oldSol_, newVars_, δ_:0, γ_:0, opts___] := FindRoot[
     {#1, #2 + γ * RandomVariate[NormalDistribution[]]} & @@@ (oldSol /. Rule -> List),
     Thread[{newVars, δ * RandomVariate[NormalDistribution[], Length[newVars]]}]
   ],
-  MaxIterations -> 50000,
+  MaxIterations -> 500,
   opts
 ]
 
-- 
cgit v1.2.3-54-g00ecf