From bdfbf32258d85d5cccae59c843cfe07aade07b8a Mon Sep 17 00:00:00 2001
From: Jaron Kent-Dobias <jaron@kent-dobias.com>
Date: Mon, 25 Oct 2021 09:48:56 +0200
Subject: Work.

---
 schofield.py | 34 ++++++++++++++++++++++++++++++++++
 1 file changed, 34 insertions(+)
 create mode 100644 schofield.py

(limited to 'schofield.py')

diff --git a/schofield.py b/schofield.py
new file mode 100644
index 0000000..96fe68a
--- /dev/null
+++ b/schofield.py
@@ -0,0 +1,34 @@
+from sympy import *
+from mpmath import glaisher, pi
+
+β = Rational(1/8)
+δ = 15
+α = 0
+σ = Rational(5, 6)
+Δ = β * δ
+
+s = 2**Rational(1, 12) * exp(-β) * glaisher**Rational(3, 2)
+
+def f(B, θc, θ):
+    return (θc * Ei(-1 / (B * θc)) * exp(1 / (B * θc)) + (θ - θc) * Ei(1/(B * (θ - θc))) * exp(-1/(B * (θ - θc)))) / pi
+
+def FH(θ0, θ):
+    return (θ**2 + θ0**2)**σ - (θ0**2)**σ
+
+def FR(AL, AH, As, B, θc, θ0, θ):
+    return AL * f(B, θc, -θ) + AH * FH(θ0, θ) + sum([A * θ**(2 * i + 2) for i,A in enumerate(As)])
+
+def F(AL, AH, As, B, θc, θ0, θ):
+    return FR(AL, AH, As, B, θc, θ0, θ) + AL * f(B, θc, θ)
+
+def dfc(AL, B, m):
+    if m == 0:
+        return AL * θc * Ei(-1 / (B * θc)) * exp(1 / (B * θc)) / pi
+    elif m == 1:
+        return 0
+    else:
+        return AL * factorial(m) * factorial(m - 2) * B**(m-1) / pi
+
+def dFc(AL, AH, As, B, θc, θ0, m):
+    θ = Symbol['θ']
+    return dfc(AL, B, m) + diff(FR(AL, AH, As, B, θc, θ0, θ), θ).subs(θ, θc) 
-- 
cgit v1.2.3-54-g00ecf