1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
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[]
|
