From d843c9d2d5c25ecc10760b1fc834348b1290921b Mon Sep 17 00:00:00 2001 From: Jaron Kent-Dobias Date: Tue, 17 Aug 2021 19:10:01 +0200 Subject: Lots of work. --- schofield-3d.wl | 189 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 189 insertions(+) create mode 100644 schofield-3d.wl (limited to 'schofield-3d.wl') diff --git a/schofield-3d.wl b/schofield-3d.wl new file mode 100644 index 0000000..87fbe9c --- /dev/null +++ b/schofield-3d.wl @@ -0,0 +1,189 @@ + +BeginPackage["Schofield`"] + +$Assumptions = {θc > 0, θc > 1, gC[_] ∈ Reals, B > 0, γ > 0, ξ0 > 0} + +β[D_:3] := Piecewise[ + { + {1/8, D == 2}, + {0.326419, D == 3}, + {1/2, D == 4}, + {β, True} + } +] + +δ[D_:3] := Piecewise[ + { + {15, D == 2}, + {4.78984, D == 3}, + {3, D == 4}, + {δ, True} + } +] + +α[D_:3] := Piecewise[ + { + {0, D == 2}, + {0.11008, D == 3}, + {0, D == 4}, + {α, True} + } +] + +Δ[D_:3] := β[D] δ[D] + +t[θ_] := ((θ)^2 - 1) +h[n_][θ_] := (1 - (θ/θc)^2) Sum[gC[i]LegendreP[(2 * i + 1), θ/θc], {i, 0, n}] + +RFLow[B_, θc_][θ_] := -(1/(12 \[Pi]))(( + E^(-(1/(B^2 (\[Theta] - \[Theta]c)^2))) + ExpIntegralE[7/6, -(1/(B^2 (\[Theta] - \[Theta]c)^2))] Gamma[1/ + 6])/(\[Theta] - \[Theta]c)^2 + (-(( + 2 E^(-(1/(B^2 \[Theta]c^2))) + ExpIntegralE[7/6, -(1/(B^2 \[Theta]c^2))])/\[Theta]c^2) + ( + E^(-(1/(B^2 (\[Theta] + \[Theta]c)^2))) + ExpIntegralE[7/ + 6, -(1/(B^2 (\[Theta] + \[Theta]c)^2))])/(\[Theta] + \ + \[Theta]c)^2) Gamma[1/6] + ( + 4 B E^(-(1/(B^2 (\[Theta] - \[Theta]c)^2))) + ExpIntegralE[5/3, -(1/(B^2 (\[Theta] - \[Theta]c)^2))] Gamma[2/ + 3])/(\[Theta] - \[Theta]c) + (( + 8 B E^(-(1/(B^2 \[Theta]c^2))) + ExpIntegralE[5/3, -(1/(B^2 \[Theta]c^2))])/\[Theta]c - ( + 4 B E^(-(1/(B^2 (\[Theta] + \[Theta]c)^2))) + ExpIntegralE[5/ + 3, -(1/(B^2 (\[Theta] + \[Theta]c)^2))])/(\[Theta] + \ + \[Theta]c)) Gamma[2/3]) +RFHigh[ξ0_][ξ_] := (ξ^2+ξ0^2)^(1+0.085) + +RF[n_][θ_] := AL RFLow[B, θc][θ] + AH RFHigh[θ0][θ] + Sum[A[i] LegendreP[(2 i), θ/θc] , {i, 0, n}] +RFReg[n_][θ_] := -AL (1/(12 \[Pi]))(( + (-(( + 2 E^(-(1/(B^2 \[Theta]c^2))) + ExpIntegralE[7/6, -(1/(B^2 \[Theta]c^2))])/\[Theta]c^2) + ( + E^(-(1/(B^2 (\[Theta] + \[Theta]c)^2))) + ExpIntegralE[7/ + 6, -(1/(B^2 (\[Theta] + \[Theta]c)^2))])/(\[Theta] + \ + \[Theta]c)^2) Gamma[1/6] + \ + (( + 8 B E^(-(1/(B^2 \[Theta]c^2))) + ExpIntegralE[5/3, -(1/(B^2 \[Theta]c^2))])/\[Theta]c - ( + 4 B E^(-(1/(B^2 (\[Theta] + \[Theta]c)^2))) + ExpIntegralE[5/ + 3, -(1/(B^2 (\[Theta] + \[Theta]c)^2))])/(\[Theta] + \ + \[Theta]c)) Gamma[2/3]) + AH RFHigh[θ0][θ] + Sum[A[i] LegendreP[(2 i), θ/θc], {i, 0, n}] +dRFc[n_][m_] := AL m! Gamma[7/6 + m/2] B^(m + 2) / (2 π) + 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[θ]]^Δ[3], θ] +ddη[h_][f_] := D[f, θ] / D[t[θ] / h[θ]^(1 / Δ[3]), θ] +dFdξLow[n_, h_][m_] := Module[{ff, hh}, Nest[ddξ[hh], ff[θ] / t[θ]^(3 * ν[3]), 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[θ]^(3 * ν[3]), 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[θ]^(-3 * ν[3] / Δ[3]) ff[θ], 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[θ]^(3 * ν[3]), 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[θ]^(3 * ν[3]), 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[θ]^(-3 * ν[3] / Δ[3]) (ff[θ]), m] /. θ -> tt /. Map[Derivative[#][ff][tt] -> Derivative[#][RF[n]][tt] &, Range[0, m]] /. hh -> h] + +ruleB[g_] := B - 0.0464597 +ruleθ0[g_] := Simplify[g[I θ0]/(-t[I θ0])^Δ[3]/I] - +ruleAL[g_] := AL + t[θc]^2 OverBar[s] / (2 π) * (- g'[θc] / t[θc]^Δ[2]) + +eqLowRHSReg[n_][m_] := dRFc[n][m] + +eqLowLHS[h_][m_] :=D[ + t[θ]^2 (Gl[h[θ] t[θ]^-Δ[2]] + Log[t[θ]^2]/(8 π)), + {θ, m} ] /. θ -> θc + +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[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} ] / m! /. θ -> 1 + +δ0 = 10^(-14); + +Φ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, + 3.16 10^-7, + 4.31 10^-6, + -1.99 10^-6 +} + +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 = { + Around[0, δ0], + Around[0, δ0], + Around[ -1.84522807823, 10^(-11)], + 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[-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]; + +ClearAll[gC] +rules := Join[ΦRules, GlRules, GhRules] +(*ξ0 := 0.18930*) +(*gC[0] := 1*) +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 + + (* *) +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, + Join[ + {#1, #2 + γ * RandomVariate[NormalDistribution[]]} & @@@ (oldSol /. Rule -> List), + Thread[{newVars, δ * RandomVariate[NormalDistribution[], Length[newVars]]}] + ], + MaxIterations -> 50000, + opts +] + +EndPackage[] + -- cgit v1.2.3-54-g00ecf