from sympy import Rational, Symbol, Ei, diff from mpmath import glaisher, pi, exp, factorial β = 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 FA(As, θ): return sum([A * θ**(2 * i + 2) for i, A in enumerate(As)]) def FR(AL, AH, As, B, θc, θ0, θ): return AL * f(B, θc, -θ) + AH * FH(θ0, θ) + FA(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, θc, 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)