1 files changed, 129 insertions, 43 deletions
diff --git a/schofield.wl b/schofield.wl
index 83cf6ac..1c750a9 100644
--- a/schofield.wl
+++ b/schofield.wl
@@ -1,7 +1,7 @@
 
 BeginPackage["Schofield`"]
 
-$Assumptions = {θc > 0, θc > 1, gC[_] ∈ Reals, B > 0, γ > 0, ξ0 > 0}
+$Assumptions = {θc > 0, θc > 1, gC[_] ∈ Reals, B > 0}
 
 β[D_:2] := Piecewise[
   {
@@ -32,10 +32,10 @@ $Assumptions = {θc > 0, θc > 1, gC[_] ∈ Reals, B > 0, γ > 0, ξ0 > 0}
 
 Δ[D_:2] := β[D] δ[D]
 
-OverBar[s] := 1.357838341706595496
+OverBar[s] := 2^(1/12) Exp[-1/8] Glaisher^(3/2)
 
-t[θ_] := ((θ)^2 - 1)
-h[n_][θ_] := (1 - (θ/θc)^2) Sum[gC[i]LegendreP[(2 * i + 1), θ/θc], {i, 0, n}]
+t[θ_] := ((θ/1)^2 - 1)
+h[n_][θ_] := (1 - (θ/θc)^2) Sum[gC[i] θ^(2*i+1), {i, 0, n}]
 η[g_][θ_] := t[θ] / (g[θ] / I)^(1 / Δ[2])
 
 RFLow[B_, θc_][θ_] := (1/\[Pi])(2 E^(1/(
@@ -45,31 +45,63 @@ 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)
+RFHigh[ξ0_][ξ_] := (ξ^2+ξ0^2)^(5/6) - (ξ0^2)^(5/6)
 
-RF[n_][θ_] := AL RFLow[B, θc][θ] + AH RFHigh[θ0][θ] + Sum[A[i] LegendreP[(2 i), θ/θc] , {i, 0, n}]
+RF[n_][θ_] := AL RFLow[B, θc][θ] + AH RFHigh[θ0][θ] + Sum[A[i] θ^(2 i), {i, 1, 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[(2 i), θ/θc], {i, 0, n}]
+     B \[Theta] + B \[Theta]c))]) + AH RFHigh[θ0][θ] + Sum[A[i] θ^(2 i), {i, 1, 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]), θ]
+
+invDerivativeList[n_][f_][x_] := Module[
+    {xp, dfs, fp, Pns},
+  dfs = Rest[NestList[D[#, xp] &, f[xp], n]] /. xp -> x;
+  Pns = FoldList[{Pm, m} |-> 
+        fp'[xp] D[Pm, xp] - (2 m - 1) fp''[xp] Pm, 1, Range[n - 1]] /. 
+    Derivative[m_][fp][xp] :> dfs[[m]];
+  MapIndexed[{Pn, i} |-> Pn/dfs[[1]]^(2 i[[1]] - 1), Pns]
+  ]
+
+dFdξLowList[n_, h_][m_] := Module[
+  { ds, dF, df },
+  ds = invDerivativeList[m+1][Function[θ, h[θ] / t[θ]^Δ[2]]][θc];
+  dF = NestList[Function[f, D[f, θ]], RFReg[n][θ], m] + Table[Piecewise[{{AL k! Gamma[k - 1] B^(k - 1)/\[Pi], k > 1}, {0, True}}], {k, 0, m}] /. θ -> θc;
+  df = NestList[D[#, \[Theta]] &, 
+      fp[\[Theta]]/t[\[Theta]]^2 - 1/(8 \[Pi]) Log[t[\[Theta]]^2], 
+      m] /. 
+   Map[Derivative[#][fp][\[Theta]] -> dF[[# + 1]] &, 
+    Range[0, m]] /. θ -> θc;
+  Table[Sum[df[[k+1]] BellY[j, k, ds[[;; j - k + 1]]], {k, 0, j}]/(j!), {j, 0, m}]
+]
+
+dFdξHighList[n_, h_][m_] := Module[
+  { ds, dF, df },
+  ds = invDerivativeList[m+1][Function[θ, h[θ] / (-t[θ])^Δ[2]]][0];
+  dF = NestList[Function[f, D[f, θ]], RF[n][θ], m] /. θ -> 0;
+  df = NestList[D[#, \[Theta]] &, 
+      fp[\[Theta]]/t[\[Theta]]^2 - 1/(8 \[Pi]) Log[t[\[Theta]]^2], 
+      m] /. 
+   Map[Derivative[#][fp][\[Theta]] -> dF[[# + 1]] &, 
+    Range[0, m]] /. θ -> 0;
+  Table[Sum[df[[k+1]] BellY[j, k, ds[[;; j - k + 1]]], {k, 0, j}]/(j!), {j, 0, m}]
+]
+
 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 - 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)
+ruleB[g_] := B -> (2 * OverBar[s] / π) * (- g'[θc] / t[θc]^Δ[2])
+ruleθ0[g_] := Around[0.18930, 0.00005] - Simplify[g[I θ0]/(-t[I θ0])^Δ[2]/I]
+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 / ((g[I θ0]/ I)^(2 / Δ[2]) * (-η[g]'[I θ0] / (2 θ0 I))^(5/6)) + Around[1.37, 0.02]
 
 eqLowRHSReg[n_][m_] := dRFc[n][m]
 
@@ -107,16 +139,16 @@ eqMid[F_, h_][m_] := D[
 }
 
 Gls = {
-  Around[0, δ0],
+  0,
   -OverBar[s],
-  −0.0489532897203,
-  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,
+  −1.000960328725262189480934955172097320572505951770117 Sqrt[2]/((2 )^(-7/8) (2^(3/16)/OverBar[s])^2)/2/(12 \[Pi]),
+  Around[ 0.038863932, 3.0 10^(-9)],
+  Around[−0.068362119, 2.0 10^(-9)],
+  Around[ 0.18388371,  1.0 10^(-8)],
+  Around[-0.659170,    1.0 10^(-6)],
+  Around[ 2.937665,    3.0 10^(-6)],
+  Around[-15.61,       1.0 10^(-2)],
+  Around[ 96.76,       1.0 10^(-2)],
   -6.79 10^2,
   5.34 10^3,
   -4.66 10^4,
@@ -125,17 +157,17 @@ Gls = {
 }
 
 Ghs = {
-  Around[0, δ0],
-  Around[0, δ0],
-  -1.845228078232838,
-  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,
+  -1.000815260440212647119476363047210236937534925597789 Sqrt[2]/((2 )^(-7/8) (2^(3/16)/OverBar[s])^2)/2,
+  0,
+  Around[  8.333711750, 5.0 10^(-9)],
+  0,
+  Around[-95.16897,     3.0 10^(-5)],
+  0,
+  Around[1457.62,       3.0 10^(-2)],
+  0,
+  Around[-2.5891 10^4,  2.0],
   0,
   5.02 10^5,
   0,
@@ -148,19 +180,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];
 
-ClearAll[gC]
-rules := Join[ΦRules, GlRules, GhRules]
-(*ξ0 := 0.18930*)
-(*gC[0] := 1*)
-tC[0] := 1
-(*gC[0] := 1*)
+ClearAll[rules]
+rules[g_] := Join[ΦRules, GlRules, GhRules, {ruleAL[g], ruleB[g], gC[0]->1}]
+
+eq[n_, g_][m_, p_, q_] := Flatten[Join[{ruleθ0[g], ruleAH[g], g'[0] θc - 1}, eqLow[n, g][#] & /@ Range[0, m],eqMid[RF[n], g][#] & /@ Range[0, p], eqHigh[n, g] /@ Range[2, q, 2]]] //. rules[g] /. 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
+eqAround[n_, g_][m_, p_, q_] := Flatten[Join[{ruleθ0[g], ruleAH[g]}, eqLow[n, g][#] & /@ Range[0, m],eqMid[RF[n], g][#] & /@ Range[0, p], eqHigh[n, g] /@ Range[2, q, 2]]] //. rules[g]
 
- (* *)
-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
+formResiduals[data_, functions_, δ_:0] := If[Head[#1]===Around,
+  (#1["Value"] - #2) / Max[#1["Uncertainty"], δ],
+  (#1 - #2) / δ] & @@@ Thread[{data, functions}]
+
+resLow[n_, g_, δ_][m_] := formResiduals[Gls[[;;m+1]], dFdξLowList[n, g][m], δ]
+
+resHigh[n_, g_, δ_][m_] := Rest[formResiduals[Ghs[[;;m+1]], dFdξHighList[n, g][m], δ][[;;;;2]]]
+
+res[F_, g_, δ_][m_] := Join[resLow[F, g, δ][m], resHigh[F, g, δ][m], {ruleθ0[g] / 0.00005, ruleAH[g] / 0.02} /. Around[x_, _] :> x]
+chiSquared[F_, g_, δ_][m_] := Total[res[F, g, δ][m]^2]
 
 newSol[eqs_, oldSol_, newVars_, δ_:0, γ_:0, opts___] := FindRoot[
   eqs,
@@ -172,5 +208,55 @@ newSol[eqs_, oldSol_, newVars_, δ_:0, γ_:0, opts___] := FindRoot[
   opts
 ]
 
+levenburgMarquardt[r_, \[Beta]0_, \[Lambda]0_ : 1, \[Nu]_ : 2] := 
+ Module[
+    {
+     n, \[Beta], new\[Beta], rc, J, I, oldJ, newJ, oldr, newr, M, 
+   g, \[Delta], newcost, \[Lambda], oldcost, temp
+   },
+  n = Length[\[Beta]0]; \[Lambda] = \[Lambda]0; \[Beta] = \[Beta]0;
+  PrintTemporary["Compiling the Jacobian..."];
+  rc = Compile[{{x, _Real, 1}}, 
+      Evaluate[
+       r /. Thread[Rule[First /@ \[Beta], Part[x, #] & /@ Range[n]]]]];
+  J = Compile[{{x, _Real, 1}}, 
+      Evaluate[
+       D[r, {First /@ \[Beta]}] /. 
+      Thread[Rule[First /@ \[Beta], Part[x, #] & /@ Range[n]]]]];
+  PrintTemporary["Beginning the algorithm."];
+  oldJ = J[\[Beta][[All, 2]]];
+  oldr = rc[\[Beta][[All, 2]]];
+  oldcost = Re[Total[oldr^2]];
+  new\[Beta] = \[Beta];
+  g = Re[Transpose[oldJ] . oldr];
+  M = Re[Transpose[oldJ] . oldJ];
+  \[Delta] = LinearSolve[M + \[Lambda] DiagonalMatrix[Diagonal[M]], g];
+  PrintTemporary[Dynamic[oldcost]]
+  While[Norm[\[Delta]] > 10^-15,
+   \[Delta] = 
+    LinearSolve[M + \[Lambda]/\[Nu] DiagonalMatrix[Diagonal[M]], g];
+   new\[Beta][[All, 2]] -= \[Delta];
+   newr = rc[new\[Beta][[All, 2]]];
+   newcost = Re[Total[newr^2]];
+   While[newcost > oldcost,
+      \[Delta] = 
+     LinearSolve[M + \[Lambda] DiagonalMatrix[Diagonal[M]], g];
+    new\[Beta] = \[Beta];
+    new\[Beta][[All, 2]] -= \[Delta];
+    newr = rc[new\[Beta][[All, 2]]];
+    newcost = Re[Total[newr^2]];
+    \[Lambda] *= \[Nu];
+    ];
+   \[Lambda] /= \[Nu];
+   oldcost = newcost;
+   oldr = newr;
+   \[Beta] = new\[Beta];
+   oldJ = J[\[Beta][[All, 2]]];
+   g = Re[Transpose[oldJ] . oldr];
+   M = Re[Transpose[oldJ] . oldJ];
+   ];
+  {newcost, \[Beta]}
+ ]
+
 EndPackage[]