Merge branch 'master' of git:research/first_order_singularities/paper

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[]
+