diff options
Diffstat (limited to 'schofield.wl')
| -rw-r--r-- | schofield.wl | 46 |
1 files changed, 29 insertions, 17 deletions
diff --git a/schofield.wl b/schofield.wl index b26b662..7f6bed8 100644 --- a/schofield.wl +++ b/schofield.wl @@ -39,9 +39,9 @@ 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}] +*) +h[n_][θ_] := (1 - (θ/θc)^2) Sum[gC[i]θ^(2 * i + 1), {i, 0, n}] RFLow[B_, θc_][θ_] := (1/\[Pi])(2 E^(1/( B \[Theta]c)) \[Theta]c ExpIntegralEi[-(1/(B \[Theta]c))] + @@ -52,35 +52,44 @@ ExpIntegralEi[1/(B \[Theta] - B \[Theta]c)] - 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}] +RF[n_][θ_] := AL RFLow[B, θc][θ] + AH RFHigh[θ0][θ] + Sum[A[i] θ^(2 i), {i, 0, 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[2i, θ], {i, 0, 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], θ] +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] + 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]] +eqLowRHSReg[n_][m_] := dRFc[n][m] eqLowLHS[h_][m_] :=D[ t[θ]^2 (Gl[h[θ] t[θ]^-Δ[2]] + Log[t[θ]^2]/(8 π)), - {θ, m} ] / m! /. θ -> θc + {θ, m} ] /. θ -> θc -eqLow[F_, h_][m_] := eqLowRHS[F][m] - eqLowLHS[h][m] +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[F_, h_][m_] := eqHighRHS[F][m] - eqHighLHS[h][m] +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} ] /. θ -> 1 -δ0 = 10^-16; +δ0 = 10^-14; Φs = { -1.197733383797993, @@ -109,14 +118,14 @@ Gls = { } Ghs = { - 0, - 0, + Around[0, δ0], + Around[0, δ0], Around[ -1.84522807823, 10^(-11)], - 0, + Around[0, δ0], Around[ 8.3337117508, 10^(-10)], - 0, + Around[0, δ0], Around[-95.16897, 10^(-5)], - 0, + Around[0, δ0], Around[1457.62, 3 10^(-2)], 0, Around[-25891, 2], @@ -134,11 +143,14 @@ rules := Join[ΦRules, GlRules, GhRules] (*ξ0 := 0.18930*) (*gC[0] := 1*) tC[0] := 1 +gC[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 +eq[n_, g_][m_] := Flatten[Join[{ruleB[g], ruleθ0[g]},{eqLow[n, g][#](*, eqMid[F, g][#]*)} & /@ Range[0, m], eqHigh[n, 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 + (* *) +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, |