From 613e0265bd0ac2229e0640e4f78fe13e535e6a04 Mon Sep 17 00:00:00 2001
From: Jaron Kent-Dobias <jaron@kent-dobias.com>
Date: Mon, 19 Apr 2021 21:01:26 +0200
Subject: Lots of work.

---
 schofield.wl | 64 ++++++++++++++++++++++++++++++++++++------------------------
 1 file changed, 39 insertions(+), 25 deletions(-)

(limited to 'schofield.wl')

diff --git a/schofield.wl b/schofield.wl
index 1c9482c..b26b662 100644
--- a/schofield.wl
+++ b/schofield.wl
@@ -30,27 +30,36 @@ $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}]]]
+
 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}]
+*)
+h[n_][θ_] := (1 - (θ/θc)^2) Sum[gC[i]LegendreP[(2 * i + 1), θ/θc], {i, 0, n}]
 
-RFLow[B_, θc_][θ_] := (2 E^(1/(
+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_][ξ_] := (1+(ξ/ξ0)^2)^(5/6) - 1
+RFHigh[ξ0_][ξ_] := (ξ^2+ξ0^2)^(5/6)
 
-RF[n_][θ_] := AL RFLow[B, θc][θ] + AH RFHigh[θ0][θ] + Sum[A[i] θ^(2(i+1)), {i, 0, n}]
+RF[n_][θ_] := AL RFLow[B, θc][θ] + AH RFHigh[θ0][θ] + Sum[A[i] LegendreP[2i, θ], {i, 0, n}]
 
-ruleB[g_] := B - π / (2 1.3578383417066) (- g'[θc] / t[θc]^Δ[2])
+RFC[n_][θ_] := RF[n][θ] + AL I Sign[Im[θ]] ((θ-θc)Exp[-1/(B(θ-θc))]-(-θ-θc)Exp[-1/(B(-θ-θc))])
+
+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])
+ruleg0[g_] := 1 - g'[0]
 
 eqLowRHS[F_][m_] := SeriesCoefficient[F[θ], {θ, θc, m}, Assumptions -> Join[{θ < θc, θ > 1}, $Assumptions]]
 
@@ -71,6 +80,8 @@ eqMid[F_, h_][m_] := D[
     /. η -> t[θ] / h[θ]^(1 / Δ[2]),
   {θ, m} ] /. θ -> 1
 
+δ0 = 10^-16;
+
 Φs = {
   -1.197733383797993,
   -0.318810124891,
@@ -83,15 +94,15 @@ eqMid[F_, h_][m_] := D[
 }
 
 Gls = {
-  0,
-  -1.3578383417066,
-  -0.048953289720,
-  0.038863932,
-  -0.068362119,
-  0.18388370,
-  -0.6591714,
-  2.937665,
-  -15.61,
+  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,
   -679,
   5340
@@ -100,15 +111,15 @@ Gls = {
 Ghs = {
   0,
   0,
-  -1.8452280782328,
+  Around[ -1.84522807823, 10^(-11)],
   0,
-  8.333711750,
+  Around[  8.3337117508, 10^(-10)],
   0,
-  -95.16896,
+  Around[-95.16897, 10^(-5)],
   0,
-  1457.62,
+  Around[1457.62, 3 10^(-2)],
   0,
-  -25891,
+  Around[-25891, 2],
   0,
   502000
 }
@@ -119,20 +130,23 @@ 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[ΦRules, GlRules, GhRules]
+rules := Join[ΦRules, GlRules, GhRules]
 (*ξ0 := 0.18930*)
-gC[0] := 1
+(*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]
+eq[F_, g_][m_] := Flatten[Join[{ruleg0[g], ruleB[g], ruleθ0[g]},{eqLow[F, g][#](*, eqMid[F, g][#]*)} & /@ Range[0, m], eqHigh[F, g] /@ Range[0, m, 2]]] //. rules /. Around[x_, _] :> x
+
+chiSquaredLow[F_, g_][m_] := Total[(((#[[1]] /. rules)["Value"] - #[[2]])^2 / (#[[1]] /. rules)["Uncertainty"]^2)& /@ Solve[0 == (eqLow[ff, g] /@ Range[0, m]), Derivative[#][Gl][0]& /@ Range[0, m]][[1]]] /. ff -> F
+chiSquaredHigh[F_, g_][m_] := Total[(((#[[1]] /. rules)["Value"] - #[[2]])^2 / (#[[1]] /. rules)["Uncertainty"]^2)& /@ Solve[0 == (eqHigh[ff, g] /@ Range[2, m, 2]) /. ff'[0] -> 0, Derivative[#][Gh][0]& /@ Range[2, m, 2]][[1]]] /. ff -> F
 
-newSol[eqs_, oldSol_, newVars_, δ_:0, opts___] := FindRoot[
+newSol[eqs_, oldSol_, newVars_, δ_:0, γ_:0, opts___] := FindRoot[
   eqs,
   Join[
-    oldSol /. Rule -> List,
+    {#1, #2 + γ * RandomVariate[NormalDistribution[]]} & @@@ (oldSol /. Rule -> List),
     Thread[{newVars, δ * RandomVariate[NormalDistribution[], Length[newVars]]}]
   ],
-  MaxIterations -> 500,
+  MaxIterations -> 1000,
   opts
 ]
 
-- 
cgit v1.2.3-54-g00ecf