diff options
Diffstat (limited to 'schofield.wl')
| -rw-r--r-- | schofield.wl | 24 |
1 files changed, 12 insertions, 12 deletions
diff --git a/schofield.wl b/schofield.wl index ecf61ff..358d5a4 100644 --- a/schofield.wl +++ b/schofield.wl @@ -41,8 +41,7 @@ h[n_][θ_] := (1 - (θ/θc)^2) Sum[gC[i] θ^(2*i+1), {i, 0, n}] fLow[B_, θc_][θ_] := (θc Exp[1/(B θc)] ExpIntegralEi[-1/(B θc)] + (θ - θc) Exp[-1/(B (θ - θc))] ExpIntegralEi[1/(B (θ - θc))]) / π RFLow[B_, θc_][θ_] := fLow[B, θc][θ] + fLow[B, θc][-θ] -RFHigh[ξ0_][ξ_] := (ξ^2+ξ0^2)^(5/6) - (ξ0^2)^(5/6) - +RFHigh[ξ0_][ξ_] := (-I ξ+ ξ0)^(5/6) + (I ξ+ ξ0)^(5/6) - 2 ξ0^(5/6) RFReg[n_][θ_] := AL fLow[B, θc][-θ] + AH RFHigh[θ0][θ] + Sum[A[i] θ^(2 i), {i, 1, n}] RF[n_][θ_] := RFReg[n][θ] + AL fLow[B, θc][θ] dfLow[AL_, B_][m_] := Piecewise[{{AL m! Gamma[m - 1] B^(m - 1) / π, m > 1}, {AL θc Exp[1/(B θc)] ExpIntegralEi[-1/(B θc)] / π, m == 0}, {0, True}}] @@ -97,7 +96,7 @@ dΦdηList[n_, h_, params_:{}][m_, θp_] := Module[ ruleB[g_] := B -> (2 * os / π) * (- g'[θc] / t[θc]^Δ[2]) ξYL[g_] := g[I θ0]/(-t[I θ0])^Δ[2]/I -AYL[g_] := AH / ((g[I θ0]/ I)^(2 / Δ[2]) * (-η[g]'[I θ0] / (2 θ0 I))^(5/6)) +AYL[g_] := AH / ((g[I θ0]/ I)^(2 / Δ[2]) * (-η[g]'[I θ0] / (I))^(5/6)) ruleθ0[g_] := ξYL[g] - Around[0.18930, 0.00005] ruleAL[g_] := AL -> Exp[Δ[2] t[θc]^(Δ[2] - 1) t'[θc] / (2 os / π g'[θc]) - t[θc]^Δ[2] g''[θc] / (4 os / π g'[θc]^2)] t[θc]^(1/8) os / (2 π) * g'[θc] ruleAH[g_] := AYL[g] + Around[1.37, 0.02] @@ -137,12 +136,12 @@ eqMid[F_, h_][m_] := D[ ]], *) Around[0.110886196683, 2.0 10^-12], - Around[0.01642689465, 1.0 10^-11], + Around[0.01642689465, 2.0 10^-11], Around[-2.639978 10^-4, 1.0 10^-10], Around[-5.140526 10^-4, 1.0 10^-10], Around[2.08865 10^-4, 1.0 10^-9], Around[-4.4819 10^-5, 1.0 10^-9], - Around[3.194 10^-7, 1.0 10^-9], + Around[3.16 10^-7, 1.0 10^-9], Around[4.31 10^-6, 0.01 10^-6], Around[-1.99 10^-6, 0.01 10^-6] } @@ -153,15 +152,15 @@ Gls = { −1.000960328725262189480934955172097320572505951770117 Sqrt[2]/((2 )^(-7/8) (2^(3/16)/os)^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[ 0.18388370, 1.0 10^(-8)], + Around[-0.6591714, 1.0 10^(-7)], Around[ 2.937665, 3.0 10^(-6)], Around[-15.61, 1.0 10^(-2)], Around[ 96.76, 1.0 10^(-2)], Around[-6.79 10^2, 1.0], - Around[5.34 10^3, 10.], + Around[ 5.34 10^3, 10.], Around[-4.66 10^4, 0.01 10^4], - Around[4.46 10^5, 0.01 10^5], + Around[ 4.46 10^5, 0.01 10^5], Around[-4.66 10^6, 0.01 10^6] } @@ -172,7 +171,7 @@ Ghs = { 0, Around[ 8.333711750, 5.0 10^(-9)], 0, - Around[-95.16897, 3.0 10^(-5)], + Around[-95.16896, 1.0 10^(-5)], 0, Around[1457.62, 3.0 10^(-2)], 0, @@ -189,7 +188,7 @@ dRule[sym_][f_, i_] := Derivative[i[[1]] - 1][sym][0] -> f (i[[1]] - 1)! GlRules = MapIndexed[dRule[Gl], Gls]; GhRules = MapIndexed[dRule[Gh], Ghs]; -rules[g_] := Join[ΦRules, GlRules, GhRules, {ruleAL[g], ruleB[g], gC[0] -> 1}] +rules[g_] := Join[ΦRules, GlRules, GhRules, {ruleAL[g], ruleB[g]}] eqAround[n_, g_][m_, p_, q_] := Flatten[Join[{ruleAH[g], ruleθ0[g]}, eqLow[n, g][#] & /@ Range[0, m],eqMid[RF[n], g][#] & /@ Range[0, p], eqHigh[n, g] /@ Range[2, q, 2]]] //. rules[g] @@ -206,6 +205,7 @@ resMid[n_, g_, δ_][m_] := formResiduals[Φs[[;;m+1]], dΦdηList[n, g][m, 1], Πres[F_, g_, δ_][m_] := Join[resLow[F, g, δ][m], resHigh[F, g, δ][m]] //. rules[g] chiSquared[F_, g_, δ_][m_] := Total[res[F, g, δ][m]^2] resAll[F_, g_, δ_][m_] := Join[resLow[F, g, δ][m], resHigh[F, g, δ][m], (*resMid[F, g, δ][m],*) {ruleθ0[g] / 0.00005, ruleAH[g] / 0.02} /. Around[x_, _] :> x] //. rules[g] +resHighAll[F_, g_, δ_][m_] := Join[resLow[F, g, δ][m][[{1}]], resHigh[F, g, δ][m], (*resMid[F, g, δ][m],*) {ruleθ0[g] / 0.00005, ruleAH[g] / 0.02} /. Around[x_, _] :> x] //. rules[g] newSol[eqs_, oldSol_, newVars_, δ_:0, γ_:0, opts___] := FindRoot[ eqs, @@ -296,7 +296,7 @@ levenbergMarquardtHelper[Δ_, rf_, Jf_, β0_, λ0_ : 1, ν_ : 1.5, ε_ : 10^-15] {newC, β} ] -annealFit[r_, β0_, γ_ : 0, λ0_ : 1, ν_ : 2, μ_ : 3, ε_: 10^-14, h_ : 0.1, α_ : 0.1] := +annealFit[r_, β0_, γ_ : 0, λ0_ : 1, ν_ : 2, μ_ : 3, ε_: 10^-15, h_ : 0.1, α_ : 0.1] := Module[ { rf, Jf, n = Length[β0], |