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}
  }
]

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_][θ_] := (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][θ] + AH RFHigh[θ0][θ] + Sum[A[i] LegendreP[2i, θ], {i, 0, n}]

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]]

eqLowLHS[h_][m_] :=D[
  t[θ]^2 (Gl[h[θ] t[θ]^-Δ[2]] + Log[t[θ]^2]/(8 π)),
  {θ, m} ] / m! /. θ -> θc

eqLow[F_, h_][m_] := eqLowRHS[F][m] - eqLowLHS[h][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[F_, h_][m_] := eqHighRHS[F][m] - eqHighLHS[h][m]

eqMid[F_, h_][m_] := D[
  F[θ] - t[θ]^2 Log[h[θ]^2]/(8 Δ[2]π) - h[θ]^((2-α[2])/Δ[2]) Φ[η]
    /. η -> t[θ] / h[θ]^(1 / Δ[2]),
  {θ, m} ] /. θ -> 1

δ0 = 10^-16;

Φ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
}

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,
  -679,
  5340
}

Ghs = {
  0,
  0,
  Around[ -1.84522807823, 10^(-11)],
  0,
  Around[  8.3337117508, 10^(-10)],
  0,
  Around[-95.16897, 10^(-5)],
  0,
  Around[1457.62, 3 10^(-2)],
  0,
  Around[-25891, 2],
  0,
  502000
}

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];

rules := Join[ΦRules, GlRules, GhRules]
(*ξ0 := 0.18930*)
(*gC[0] := 1*)
tC[0] := 1

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, γ_:0, opts___] := FindRoot[
  eqs,
  Join[
    {#1, #2 + γ * RandomVariate[NormalDistribution[]]} & @@@ (oldSol /. Rule -> List),
    Thread[{newVars, δ * RandomVariate[NormalDistribution[], Length[newVars]]}]
  ],
  MaxIterations -> 1000,
  opts
]

EndPackage[]