diff options
Diffstat (limited to 'IsingScalingFunction.wl')
| -rw-r--r-- | IsingScalingFunction.wl | 292 |
1 files changed, 292 insertions, 0 deletions
diff --git a/IsingScalingFunction.wl b/IsingScalingFunction.wl new file mode 100644 index 0000000..81d5e47 --- /dev/null +++ b/IsingScalingFunction.wl @@ -0,0 +1,292 @@ +BeginPackage["IsingScalingFunction`"] + +InverseCoordinates::usage = "Numerically convert Schofield coordinates to t and h." + +g::usage = "g[θ0, gs][θ] gives the Schofield coordinate transformation defined in (14)." + +ut::usage = "ut[θ] gives the scaling field u_t as a function of Schofield coordinates." + +uh::usage = "uh[θ0, gs][θ] gives the scaling field u_h as a function of Schofield coordinates." + +η::usage = "η[θ0, gs][θ] gives the invariant scaling combination η." + +ξ::usage = "ξ[θ0, gs][θ] gives the invariant scaling combination ξ." + +ReScriptF::usage = "ReScriptF[θ0, θYL, B, C0, CYL, Gs][θ] gives the free energy scaling function defined in (19)." + +ScriptF::usage = "ScriptF[θ0, θYL, B, C0, CYL, Gs][θ] gives the free energy scaling function defined in (35)." + +DScriptFPlusMinusDξθ0List::usage = + "DScriptFPlusMinusDξθ0List computes the first m derivatives of the scaling function F_- evaluated at θ_0." + +DScriptFPlusMinusDξList::usage = + "DScriptFPlusMinusDξList computes the first m derivatives of the scaling function F_+/-." + +DScriptF0DηList::usage = + "DScriptF0DηList computes the first m derivatives of the scaling function F_0." + +DScriptF0Dη::usage = + "DScriptF0Dη computes the mth derivative of the scaling function F_0." + +DScriptMCasDξList::usage = "Computes the first m derivatives of the scaling function M given by Caselle et al." + +uf::usage = "uf computes the singular free energy u_f." + +DufDut::usage = + "DufDut computes derivatives of the singular free energy u_f with respect to the temperature-like scaling field u_t." + +DufDuh::usage = + "DufDuh computes derivatives of the singular free energy u_f with respect to the temperature-like scaling field u_h." + +ruleB::usage = "Fixes B given other data as in (38)." + +ruleC0::usage = "Fixes C0 given other data as in (39)." + +Data::usage = "Data[n] gives data from the fit to nth order from Table II." + +PrepareArgument::usage = "Converts scaling function data into appropriate argument to function interfaces." + +θ0Cas::usage = "θ0 from Caselle et al." + +gsCas::usage = "g function coefficients from Caselle et al." + +Φs::usage = "List of numeric coefficients for the scaling function F_0" + +Gls::usage = "List of numeric coefficients for the scaling function F_-" + +Ghs::usage = "List of numeric coefficients for the scaling function F_+" + +Begin["Private`"] + +β := 1/8 + +δ := 15 + +Δ := β δ + +OverlineS := 2^(1/12) Exp[-1/8] Glaisher^(3/2) + +Φs := { + -Gamma[1/3]Gamma[1/5]Gamma[7/15]/(2 π Gamma[2/3]Gamma[4/5]Gamma[8/15])(4 π^2 Gamma[13/16]^2 Gamma[3/4]/(Gamma[3/16]^2 Gamma[1/4]))^(8/15), + -0.31881012489061, + Around[0.110886196683, 2.0 10^-12], + 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.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] +} + +Gls := { + 0, + -OverlineS, + −1.000960328725262189480934955172097320572505951770117 Sqrt[2]/((2 )^(-7/8) (2^(3/16)/OverlineS)^2)/2/(12 \[Pi]), + Around[ 0.038863932, 3.0 10^(-9)], + Around[−0.068362119, 2.0 10^(-9)], + 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[-4.66 10^4, 0.01 10^4], + Around[ 4.46 10^5, 0.01 10^5], + Around[-4.66 10^6, 0.01 10^6] +} + +Ghs := { + 0, + 0, + -1.000815260440212647119476363047210236937534925597789 Sqrt[2]/((2 )^(-7/8) (2^(3/16)/OverlineS)^2)/2, + 0, + Around[ 8.333711750, 5.0 10^(-9)], + 0, + Around[-95.16896, 1.0 10^(-5)], + 0, + Around[1457.62, 3.0 10^(-2)], + 0, + Around[-2.5891 10^4, 2.0], + 0, + Around[5.02 10^5, 0.01 10^5], + 0, + Around[-1.04 10^7, 0.01 10^7] +} + +Data[2] = Rationalize[#, 10^-20] & /@ <| + "θ0" -> 1.148407773492004`, + "θYL" -> 0.9896669889911205`, + "CYL" -> -0.172823989504767`, + "Gs" -> {-0.31018352388662596`, 0.2474537923130002`}, + "gs" -> {0.37369093055254343`, -0.021636313152585823`} +|> + +Data[3] = Rationalize[#, 10^-20] & /@ <| + "θ0" -> 1.2542120477507488`, + "θYL" -> 0.6020557328641167`, + "CYL" -> -0.38566364361428684`, + "Gs" -> {-0.3527514794812415`, 0.2582430860166863`}, + "gs" -> {0.4483788209731592`, -0.022032295172535358`, 0.00022200608228654115`} +|> + +Data[4] = Rationalize[#, 10^-20] & /@ <| + "θ0" -> 1.3164928721109121`, + "θYL" -> 0.6400189996493497`, + "CYL" -> -0.3563974694580203`, + "Gs" -> {-0.3550547624920263`, 0.23465947408509413`, -0.0019083731028066697`}, + "gs" -> {0.4410742751152714`, -0.034817777358116885`, 0.000678172648789985`, -0.00004305140578834467`} +|> + +Data[5] = Rationalize[#, 10^-20] & /@ <| + "θ0" -> 1.3403205742656135`, + "θYL" -> 0.6238113973493433`, + "CYL" -> -0.38002950945224295`, + "Gs" -> {-0.35127522582179693`, 0.23704589676915347`, -0.007319731639727028`}, + "gs" -> {0.44371885415894785`, -0.04609943321005163`, -0.0007458341071947777`, 0.00005966875622885447`, -4.403083529955303`*^-6} +|> + +Data[6] = Rationalize[#, 10^-20] & /@ <| + "θ0" -> 1.3626103817690176`, + "θYL" -> 0.6462147447024515`, + "CYL" -> -0.35576386594103865`, + "Gs" -> {-0.3520586281920383`, 0.23316561297622435`, -0.006649030656179257`, -0.0016899077640685814`}, + "gs" -> {0.43845335615925396`, -0.05312704168994819`, -0.003914782631377569`, -0.0004080160912692574`, 0.000026262906640471588`, -1.0974538440746828`*^-6} +|> + +PrepareArgument[data_] := With[ + { + θ0 = data["θ0"], + gs = data["gs"] + }, + { + θ0, + data["θYL"], + ruleB[θ0, gs], + ruleC0[θ0, gs], + data["CYL"], + data["Gs"], + gs + } +] + +t[θ_] := θ^2 - 1 + +g[θ0_, gs_][θ_] := (1 - (θ/θ0)^2) Total[MapIndexed[Function[{gi, i}, gi θ^(2*i[[1]]-1)], gs]] + +ut[R_, θ_] := R t[θ] + +uh[θ0_, gs_][R_, θ_] := R^Δ g[θ0, gs][θ] + +InverseCoordinates[\[Theta]0_, gs_, wp_:20][tn_, hn_] := + ({Exp[logR], \[Theta]0 Tanh[x]} /. + FindRoot[{ + Rationalize[tn, 10^-30] == ut[Exp[logR], \[Theta]0 Tanh[x]], + Rationalize[hn, 10^-30] == uh[\[Theta]0, gs][Exp[logR], \[Theta]0 Tanh[x]] + }, {{logR, 2}, {x, Sign[hn]/2}}, WorkingPrecision -> wp]) /; + NumericQ[tn] && NumericQ[hn] + +η[θ0_, gs_][θ_] := t[θ] / RealAbs[g[θ0, gs][θ]]^(1 / Δ) + +ξ[θ0_, gs_][θ_] := g[θ0, gs][θ] / RealAbs[t[θ]]^Δ + +ScriptR[θc_, B_][θ_] := (θc Exp[1/(B θc)] ExpIntegralEi[-1/(B θc)] + (θ - θc) Exp[-1/(B (θ - θc))] ExpIntegralEi[1/(B (θ - θc))]) / π + +ReScriptF0[C0_, θc_, B_][θ_] := C0 (ScriptR[θc, B][θ] + ScriptR[θc, B][-θ]) + +ScriptFYL[θYL_, CYL_][θ_] := CYL ((-I θ + θYL)^(5/6) + (I θ + θYL)^(5/6) - 2 θYL^(5/6)) + +ReScriptFRegular[θ0_, θYL_, B_, C0_, CYL_, Gs_][θ_] := C0 ScriptR[θ0, B][-θ] + ScriptFYL[θYL, CYL][θ] + Total[MapIndexed[Function[{G, i}, G θ^(2*i[[1]])], Gs]] + +ReScriptF[θ0_, θYL_, B_, C0_, CYL_, Gs_][θ_] := ReScriptFRegular[θ0, θYL, B, C0, CYL, Gs][θ] + C0 ScriptR[θ0, B][θ] + +DReScriptFIrregular[θ0_, B_, C0_][m_] := Piecewise[{{C0 m! Gamma[m - 1] B^(m - 1) / π, m > 1}, {C0 θ0 Exp[1/(B θ0)] ExpIntegralEi[-1/(B θ0)] / π, m == 0}}, 0] + +ScriptF[θ0_, θYL_, B_, C0_, CYL_, Gs_][θ_] := ReScriptF[θ0, θYL, B, C0, CYL, Gs][θ] + C0 I Sign[Im[θ]] ((θ-θ0)Exp[-1/(B(θ-θ0))]-(-θ-θ0)Exp[-1/(B(-θ-θ0))]) + +ScriptFPlusMinus[ScriptF_][θ_] := ScriptF[θ] / t[θ]^2 - 1/(8 \[Pi]) Log[t[θ]^2] + +ScriptF0[θ0_, gs_][ScriptF_][θ_] := RealAbs[g[θ0, gs][θ]]^(-2 / Δ) ScriptF[θ] - η[θ0, gs][θ]^2 Log[g[θ0, gs][θ]^2] / (8 π Δ) + +uf[params__][R_, θ_] := R^2 ReScriptF[params][θ] + t[θ]^2 R^2 / (8 π) Log[R^2] + +EfficientDerivativeList[n_][f_][x_] := Module[ + {xp}, NestList[Function[g, D[g, xp]], f[xp], n] /. xp -> x +] + +InverseDerivativeList[n_][f_][x_] := Module[ + {xp, dfs, fp, Pns}, + dfs = Rest[EfficientDerivativeList[n][f][x]]; + Pns = FoldList[Function[{Pm, m}, + fp'[xp] D[Pm, xp] - (2 m - 1) fp''[xp] Pm], 1, Range[n - 1]] /. + Derivative[m_][fp][xp] :> dfs[[m]]; + MapIndexed[{Pn, i} \[Function] Pn/dfs[[1]]^(2 i[[1]] - 1), Pns] +] + +CompositeFunctionDerivativeList[G_, F_, X_, FSupp_:(0&)][m_, θ_] := Module[ + { ds, dF, df, fp }, + ds = InverseDerivativeList[m+1][X][θ]; + dF = EfficientDerivativeList[m][F][θ] + FSupp /@ Range[0, m]; + df = EfficientDerivativeList[m][G[fp]][θ] /. + Map[Derivative[#][fp][θ] -> dF[[# + 1]] &, Range[0, m]]; + Table[Sum[df[[k+1]] BellY[j, k, ds[[;; j - k + 1]]], {k, 0, j}]/(j!), {j, 0, m}] +] + +DScriptFPlusMinusDξθ0List[θ0_, θYL_, B_, C0_, CYL_, Gs_, gs_][m_] := CompositeFunctionDerivativeList[ + ScriptFPlusMinus, ReScriptFRegular[θ0, θYL, B, C0, CYL, Gs], + ξ[θ0, gs], DReScriptFIrregular[θ0, B, C0] + ][m, θ0] + +DScriptFPlusMinusDξList[θ0_, θYL_, B_, C0_, CYL_, Gs_, gs_][m_, θ_] := CompositeFunctionDerivativeList[ + ScriptFPlusMinus, ReScriptF[θ0, θYL, B, C0, CYL, Gs], ξ[θ0, gs] + ][m, θ] + +DScriptF0DηList[θ0_, θYL_, B_, C0_, CYL_, Gs_, gs_][m_, θ_] := CompositeFunctionDerivativeList[ + ScriptF0[θ0, gs], ReScriptF[θ0, θYL, B, C0, CYL, Gs], η[θ0, gs] + ][m, θ] + +DScriptFPlusMinusDξθ0[params__][m_] := Last[DScriptFPlusMinusDξθ0List[params][m]] + +DScriptFPlusMinusDξ[params__][m_, θ_] := Last[DScriptFPlusMinusDξList[params][m, θ]] + +DScriptF0Dη[params__][m_, θ_] := Last[DScriptF0DηList[params][m, θ]] + +DufDut[θ0_, θYL_, B_, C0_, CYL_, Gs_, gs_][m_][R_, θ_] := m! RealAbs[uh[θ0, gs][R, θ]]^(2 / Δ - m / Δ) DScriptF0Dη[θ0, θYL, B, C0, CYL, Gs, gs][m, θ] + Log[uh[θ0, gs][R, θ]^2] / (8 π Δ) Derivative[m][Function[utp, utp^2]][ut[R, θ]] + +DufDuh[θ0_, θYL_, B_, C0_, CYL_, Gs_, gs_][m_][R_, θ_] := m! RealAbs[ut[R, θ]]^(2-m Δ) DScriptFPlusMinusDξ[θ0, θYL, B, C0, CYL, Gs, gs][m, θ] + Derivative[m][1&][θ] ut[R, θ]^2 / (8 π) Log[ut[R, θ]^2] + +ruleB[θ0_, gs_] := (2 * OverlineS / π) * (- g[θ0, gs]'[θ0] / t[θ0]^Δ) + +ruleC0[θ0_, gs_] := Exp[Δ t[θ0]^(Δ - 1) t'[θ0] / (2 OverlineS / π g[θ0, gs]'[θ0]) - t[θ0]^Δ g[θ0, gs]''[θ0] / (4 OverlineS / π g[θ0, gs]'[θ0]^2)] t[θ0]^(1/8) OverlineS / (2 π) * g[θ0, gs]'[θ0] + +θ0Cas := Sqrt[1.16951]; + +h0Cas := a b ρ /. { + a -> c2^2/c4, + b -> (-c4/c2^3)^(1/2), + ρ -> 2.00881 + } /. { + c2->Ghs[[3]] 2!, + c4->Ghs[[5]]["Value"] 4! + } + +gsCas := h0Cas * { + 1, + -0.222389, + -0.043547, + -0.014809, + -0.007168 + } + +m0Cas := -Ghs[[3]]2! h0Cas + +DScriptMCasDξList[m_, θ_] := CompositeFunctionDerivativeList[ + Identity, Function[θp, m0Cas * θp / RealAbs[θp^2 - 1]^β], ξ[θ0Cas, gsCas] + ][m, θ] + +End[] + +EndPackage[] + |