def transmutagen_cram_error(degree, t0, prec=200):
from ..partialfrac import (thetas_alphas, thetas_alphas_to_expr_complex,
customre, t)
from ..cram import get_CRAM_from_cache
from sympy import re, exp, nsolve, diff
expr = get_CRAM_from_cache(degree, prec)
thetas, alphas, alpha0 = thetas_alphas(expr, prec)
part_frac = thetas_alphas_to_expr_complex(thetas, alphas, alpha0)
part_frac = part_frac.replace(customre, re)
E = part_frac - exp(-t)
E = E.evalf(20)
return E.subs(t, nsolve(diff(E, t), t0, prec=prec))
python类exp()的实例源码
def can_do(ap, bq, numerical=True, div=1, lowerplane=False):
from sympy import exp_polar, exp
r = hyperexpand(hyper(ap, bq, z))
if r.has(hyper):
return False
if not numerical:
return True
repl = {}
for n, a in enumerate(r.free_symbols - set([z])):
repl[a] = randcplx(n)/div
[a, b, c, d] = [2, -1, 3, 1]
if lowerplane:
[a, b, c, d] = [2, -2, 3, -1]
return tn(
hyper(ap, bq, z).subs(repl),
r.replace(exp_polar, exp).subs(repl),
z, a=a, b=b, c=c, d=d)
def test_hyperexpand_bases():
assert hyperexpand(hyper([2], [a], z)) == \
a + z**(-a + 1)*(-a**2 + 3*a + z*(a - 1) - 2)*exp(z)* \
lowergamma(a - 1, z) - 1
# TODO [a+1, a-S.Half], [2*a]
assert hyperexpand(hyper([1, 2], [3], z)) == -2/z - 2*log(-z + 1)/z**2
assert hyperexpand(hyper([S.Half, 2], [S(3)/2], z)) == \
-1/(2*z - 2) + atanh(sqrt(z))/sqrt(z)/2
assert hyperexpand(hyper([S(1)/2, S(1)/2], [S(5)/2], z)) == \
(-3*z + 3)/4/(z*sqrt(-z + 1)) \
+ (6*z - 3)*asin(sqrt(z))/(4*z**(S(3)/2))
assert hyperexpand(hyper([1, 2], [S(3)/2], z)) == -1/(2*z - 2) \
- asin(sqrt(z))/(sqrt(z)*(2*z - 2)*sqrt(-z + 1))
assert hyperexpand(hyper([-S.Half - 1, 1, 2], [S.Half, 3], z)) == \
sqrt(z)*(6*z/7 - S(6)/5)*atanh(sqrt(z)) \
+ (-30*z**2 + 32*z - 6)/35/z - 6*log(-z + 1)/(35*z**2)
assert hyperexpand(hyper([1 + S.Half, 1, 1], [2, 2], z)) == \
-4*log(sqrt(-z + 1)/2 + S(1)/2)/z
# TODO hyperexpand(hyper([a], [2*a + 1], z))
# TODO [S.Half, a], [S(3)/2, a+1]
assert hyperexpand(hyper([2], [b, 1], z)) == \
z**(-b/2 + S(1)/2)*besseli(b - 1, 2*sqrt(z))*gamma(b) \
+ z**(-b/2 + 1)*besseli(b, 2*sqrt(z))*gamma(b)
# TODO [a], [a - S.Half, 2*a]
def test_quantum_fourier():
assert QFT(0, 3).decompose() == \
SwapGate(0, 2)*HadamardGate(0)*CGate((0,), PhaseGate(1)) * \
HadamardGate(1)*CGate((0,), TGate(2))*CGate((1,), PhaseGate(2)) * \
HadamardGate(2)
assert IQFT(0, 3).decompose() == \
HadamardGate(2)*CGate((1,), RkGate(2, -2))*CGate((0,), RkGate(2, -3)) * \
HadamardGate(1)*CGate((0,), RkGate(1, -2))*HadamardGate(0)*SwapGate(0, 2)
assert represent(QFT(0, 3), nqubits=3) == \
Matrix([[exp(2*pi*I/8)**(i*j % 8)/sqrt(8) for i in range(8)] for j in range(8)])
assert QFT(0, 4).decompose() # non-trivial decomposition
assert qapply(QFT(0, 3).decompose()*Qubit(0, 0, 0)).expand() == qapply(
HadamardGate(0)*HadamardGate(1)*HadamardGate(2)*Qubit(0, 0, 0)
).expand()
def test_wavefunction():
a = 1/Z
R = {
(1, 0): 2*sqrt(1/a**3) * exp(-r/a),
(2, 0): sqrt(1/(2*a**3)) * exp(-r/(2*a)) * (1 - r/(2*a)),
(2, 1): S(1)/2 * sqrt(1/(6*a**3)) * exp(-r/(2*a)) * r/a,
(3, 0): S(2)/3 * sqrt(1/(3*a**3)) * exp(-r/(3*a)) *
(1 - 2*r/(3*a) + S(2)/27 * (r/a)**2),
(3, 1): S(4)/27 * sqrt(2/(3*a**3)) * exp(-r/(3*a)) *
(1 - r/(6*a)) * r/a,
(3, 2): S(2)/81 * sqrt(2/(15*a**3)) * exp(-r/(3*a)) * (r/a)**2,
(4, 0): S(1)/4 * sqrt(1/a**3) * exp(-r/(4*a)) *
(1 - 3*r/(4*a) + S(1)/8 * (r/a)**2 - S(1)/192 * (r/a)**3),
(4, 1): S(1)/16 * sqrt(5/(3*a**3)) * exp(-r/(4*a)) *
(1 - r/(4*a) + S(1)/80 * (r/a)**2) * (r/a),
(4, 2): S(1)/64 * sqrt(1/(5*a**3)) * exp(-r/(4*a)) *
(1 - r/(12*a)) * (r/a)**2,
(4, 3): S(1)/768 * sqrt(1/(35*a**3)) * exp(-r/(4*a)) * (r/a)**3,
}
for n, l in R:
assert simplify(R_nl(n, l, r, Z) - R[(n, l)]) == 0
def test_Znm():
# http://en.wikipedia.org/wiki/Solid_harmonics#List_of_lowest_functions
th, ph = Symbol("theta", real=True), Symbol("phi", real=True)
from sympy.abc import n,m
assert Znm(0, 0, th, ph) == Ynm(0, 0, th, ph)
assert Znm(1, -1, th, ph) == (-sqrt(2)*I*(Ynm(1, 1, th, ph)
- exp(-2*I*ph)*Ynm(1, 1, th, ph))/2)
assert Znm(1, 0, th, ph) == Ynm(1, 0, th, ph)
assert Znm(1, 1, th, ph) == (sqrt(2)*(Ynm(1, 1, th, ph)
+ exp(-2*I*ph)*Ynm(1, 1, th, ph))/2)
assert Znm(0, 0, th, ph).expand(func=True) == 1/(2*sqrt(pi))
assert Znm(1, -1, th, ph).expand(func=True) == (sqrt(3)*I*sqrt(-cos(th)**2 + 1)*exp(I*ph)/(4*sqrt(pi))
- sqrt(3)*I*sqrt(-cos(th)**2 + 1)*exp(-I*ph)/(4*sqrt(pi)))
assert Znm(1, 0, th, ph).expand(func=True) == sqrt(3)*cos(th)/(2*sqrt(pi))
assert Znm(1, 1, th, ph).expand(func=True) == (-sqrt(3)*sqrt(-cos(th)**2 + 1)*exp(I*ph)/(4*sqrt(pi))
- sqrt(3)*sqrt(-cos(th)**2 + 1)*exp(-I*ph)/(4*sqrt(pi)))
assert Znm(2, -1, th, ph).expand(func=True) == (sqrt(15)*I*sqrt(-cos(th)**2 + 1)*exp(I*ph)*cos(th)/(4*sqrt(pi))
- sqrt(15)*I*sqrt(-cos(th)**2 + 1)*exp(-I*ph)*cos(th)/(4*sqrt(pi)))
assert Znm(2, 0, th, ph).expand(func=True) == 3*sqrt(5)*cos(th)**2/(4*sqrt(pi)) - sqrt(5)/(4*sqrt(pi))
assert Znm(2, 1, th, ph).expand(func=True) == (-sqrt(15)*sqrt(-cos(th)**2 + 1)*exp(I*ph)*cos(th)/(4*sqrt(pi))
- sqrt(15)*sqrt(-cos(th)**2 + 1)*exp(-I*ph)*cos(th)/(4*sqrt(pi)))
def eval(cls, nu, z):
from sympy import (unpolarify, expand_mul, uppergamma, exp, gamma,
factorial)
nu2 = unpolarify(nu)
if nu != nu2:
return expint(nu2, z)
if nu.is_Integer and nu <= 0 or (not nu.is_Integer and (2*nu).is_Integer):
return unpolarify(expand_mul(z**(nu - 1)*uppergamma(1 - nu, z)))
# Extract branching information. This can be deduced from what is
# explained in lowergamma.eval().
z, n = z.extract_branch_factor()
if n == 0:
return
if nu.is_integer:
if (nu > 0) != True:
return
return expint(nu, z) \
- 2*pi*I*n*(-1)**(nu - 1)/factorial(nu - 1)*unpolarify(z)**(nu - 1)
else:
return (exp(2*I*pi*nu*n) - 1)*z**(nu - 1)*gamma(1 - nu) + expint(nu, z)
def _getunbranched(cls, ar):
from sympy import exp_polar, log, polar_lift
if ar.is_Mul:
args = ar.args
else:
args = [ar]
unbranched = 0
for a in args:
if not a.is_polar:
unbranched += arg(a)
elif a.func is exp_polar:
unbranched += a.exp.as_real_imag()[1]
elif a.is_Pow:
re, im = a.exp.as_real_imag()
unbranched += re*unbranched_argument(
a.base) + im*log(abs(a.base))
elif a.func is polar_lift:
unbranched += arg(a.args[0])
else:
return None
return unbranched
def _mul_args(f):
"""
Return a list ``L`` such that Mul(*L) == f.
If f is not a Mul or Pow, L=[f].
If f=g**n for an integer n, L=[g]*n.
If f is a Mul, L comes from applying _mul_args to all factors of f.
"""
args = Mul.make_args(f)
gs = []
for g in args:
if g.is_Pow and g.exp.is_Integer:
n = g.exp
base = g.base
if n < 0:
n = -n
base = 1/base
gs += [base]*n
else:
gs.append(g)
return gs
def test_as_integral():
from sympy import Function, Integral
f = Function('f')
assert mellin_transform(f(x), x, s).rewrite('Integral') == \
Integral(x**(s - 1)*f(x), (x, 0, oo))
assert fourier_transform(f(x), x, s).rewrite('Integral') == \
Integral(f(x)*exp(-2*I*pi*s*x), (x, -oo, oo))
assert laplace_transform(f(x), x, s).rewrite('Integral') == \
Integral(f(x)*exp(-s*x), (x, 0, oo))
assert str(inverse_mellin_transform(f(s), s, x, (a, b)).rewrite('Integral')) \
== "Integral(x**(-s)*f(s), (s, _c - oo*I, _c + oo*I))"
assert str(inverse_laplace_transform(f(s), s, x).rewrite('Integral')) == \
"Integral(f(s)*exp(s*x), (s, _c - oo*I, _c + oo*I))"
assert inverse_fourier_transform(f(s), s, x).rewrite('Integral') == \
Integral(f(s)*exp(2*I*pi*s*x), (s, -oo, oo))
# NOTE this is stuck in risch because meijerint cannot handle it
def test_recursive():
from sympy import symbols, exp_polar, expand
a, b, c = symbols('a b c', positive=True)
r = exp(-(x - a)**2)*exp(-(x - b)**2)
e = integrate(r, (x, 0, oo), meijerg=True)
assert simplify(e.expand()) == (
sqrt(2)*sqrt(pi)*(
(erf(sqrt(2)*(a + b)/2) + 1)*exp(-a**2/2 + a*b - b**2/2))/4)
e = integrate(exp(-(x - a)**2)*exp(-(x - b)**2)*exp(c*x), (x, 0, oo), meijerg=True)
assert simplify(e) == (
sqrt(2)*sqrt(pi)*(erf(sqrt(2)*(2*a + 2*b + c)/4) + 1)*exp(-a**2 - b**2
+ (2*a + 2*b + c)**2/8)/4)
assert simplify(integrate(exp(-(x - a - b - c)**2), (x, 0, oo), meijerg=True)) == \
sqrt(pi)/2*(1 + erf(a + b + c))
assert simplify(integrate(exp(-(x + a + b + c)**2), (x, 0, oo), meijerg=True)) == \
sqrt(pi)/2*(1 - erf(a + b + c))
def main():
x = Symbol("x")
a = Symbol("a")
h = Symbol("h")
show( limit(sqrt(x**2 - 5*x + 6) - x, x, oo), -Rational(5)/2 )
show( limit(x*(sqrt(x**2 + 1) - x), x, oo), Rational(1)/2 )
show( limit(x - sqrt3(x**3 - 1), x, oo), Rational(0) )
show( limit(log(1 + exp(x))/x, x, -oo), Rational(0) )
show( limit(log(1 + exp(x))/x, x, oo), Rational(1) )
show( limit(sin(3*x)/x, x, 0), Rational(3) )
show( limit(sin(5*x)/sin(2*x), x, 0), Rational(5)/2 )
show( limit(((x - 1)/(x + 1))**x, x, oo), exp(-2))
def real(self, nested_scope=None):
"""Return the correspond floating point number."""
op = self.children[0].name
expr = self.children[1]
dispatch = {
'sin': sympy.sin,
'cos': sympy.cos,
'tan': sympy.tan,
'asin': sympy.asin,
'acos': sympy.acos,
'atan': sympy.atan,
'exp': sympy.exp,
'ln': sympy.log,
'sqrt': sympy.sqrt
}
if op in dispatch:
arg = expr.real(nested_scope)
return dispatch[op](arg)
else:
raise NodeException("internal error: undefined external")
def sym(self, nested_scope=None):
"""Return the corresponding symbolic expression."""
op = self.children[0].name
expr = self.children[1]
dispatch = {
'sin': sympy.sin,
'cos': sympy.cos,
'tan': sympy.tan,
'asin': sympy.asin,
'acos': sympy.acos,
'atan': sympy.atan,
'exp': sympy.exp,
'ln': sympy.log,
'sqrt': sympy.sqrt
}
if op in dispatch:
arg = expr.sym(nested_scope)
return dispatch[op](arg)
else:
raise NodeException("internal error: undefined external")
def _print_Pow(self, expr):
PREC = sympy.printing.precedence.precedence(expr)
if expr.exp == 0:
return "1.0"
elif expr.exp == 0.5:
return 'sqrt(%s)' % self._print(expr.base)
elif expr.exp == 1:
return self.parenthesize(expr.base,PREC)
elif expr.exp.is_constant and int(expr.exp) == expr.exp:
base = self.parenthesize(expr.base,PREC)
if expr.exp > 0:
return "(" + "*".join(repeat(base,int(expr.exp))) + ")"
else:
return "(1.0/" + "/".join(repeat(base,int(expr.exp))) + ")"
return 'pow(%s, %s)' % (self._print(expr.base),
self._print(expr.exp))
def __init__(self, space, nuclearHamiltonian, overRideDT=None):
self.myHamiltonian = nuclearHamiltonian
self.mySpace = space
if overRideDT is None:
self.dt = self.mySpace.dt
else:
self.dt = overRideDT
#set up the kinetic propagator,
kineticSurface = self.myHamiltonian.myKineticOperator.surface
kineticSurface = np.exp(-1.0j *self.dt * kineticSurface / self.mySpace.hbar)
self.myKineticOperator = momentumOperator(self.mySpace, kineticSurface)
#set up the potential propagator
potentialSurface = self.myHamiltonian.myPotentialOperator.surface
potentialSurface = np.exp(-1.0j * self.dt * potentialSurface / self.mySpace.hbar)
self.myPotentialOperator = positionOperator(self.mySpace, potentialSurface)
def __init__(self, space, nuclearHamiltonian, overRideDT=None):
self.myHamiltonian = nuclearHamiltonian
self.mySpace = space
if overRideDT is None:
self.dt = self.mySpace.dt
else:
self.dt = overRideDT
#set up the kinetic propagator,
kineticSurface = self.myHamiltonian.myKineticOperator.surface
kineticSurface = np.exp(-1.0j * .5 *self.dt * kineticSurface / self.mySpace.hbar)
self.myKineticOperator = momentumOperator(self.mySpace, kineticSurface)
#set up the potential propagator
potentialSurface = self.myHamiltonian.myPotentialOperator.surface
potentialSurface = np.exp(-1.0j * self.dt * potentialSurface / self.mySpace.hbar)
self.myPotentialOperator = positionOperator(self.mySpace, potentialSurface)
def __init__(self, space, nuclearHamiltonian, overRideDT=None):
#raise Exception("DO NOT USE HIGH-ACCURACY PROPAGATOR; FUNCTIONALITY NOT CODED YET")
self.myHamiltonian = nuclearHamiltonian
self.mySpace = space
if overRideDT is None:
self.dt = self.mySpace.dt
else:
self.dt = overRideDT
gamma = 1.0 / (2.0 - 2.0**(1.0/3.0))
lam = -1.0j * self.dt / self.mySpace.hbar
#set up the kinetic propagator,
kineticSurfaceLambda = self.myHamiltonian.myKineticOperator.surface * lam * .5
kineticSurface1 = np.exp( gamma * kineticSurfaceLambda )
kineticSurface2 = np.exp( (1.0 - gamma) * kineticSurfaceLambda )
self.myKineticOperator1 = momentumOperator(self.mySpace, kineticSurface1)
self.myKineticOperator2 = momentumOperator(self.mySpace, kineticSurface2)
#set up the potential propagator
potentialSurfaceLambda = self.myHamiltonian.myPotentialOperator.surface * lam
potentialSurface1 = np.exp( gamma * potentialSurfaceLambda )
potentialSurface2 = np.exp( (1.0 - 2.0 *gamma) * potentialSurfaceLambda )
self.myPotentialOperator1 = positionOperator(self.mySpace, potentialSurface1)
self.myPotentialOperator2 = positionOperator(self.mySpace, potentialSurface2)
def potentialFunction(self, x):
naturalPotential = self.De * (1 - np.exp(-self.a * (x - self.center)))**2 + self.energyOffset
imaginaryPotential = 0
try:
imaginaryPotentialZeros = np.zeros(x.shape)
except:
if x > self.startOfAbsorbingPotential:
imaginaryPotential = -1.0j *self.strength * np.cosh( (np.real(x) - self.mySpace.xMax)**2 / (self.mySpace.Dx*30)**2)**(-2)
else:
imaginaryPotential = 0
return naturalPotential + imaginaryPotential
if self.absorbingPotential:
imaginaryPotential = -1.0j *self.strength * np.cosh( (np.real(x) - self.mySpace.xMax)**2 / (self.mySpace.Dx*30)**2)**(-2)
imaginaryPotential = np.where(x > self.startOfAbsorbingPotential, imaginaryPotential, imaginaryPotentialZeros)
return naturalPotential + imaginaryPotential
def can_do(ap, bq, numerical=True, div=1, lowerplane=False):
from sympy import exp_polar, exp
r = hyperexpand(hyper(ap, bq, z))
if r.has(hyper):
return False
if not numerical:
return True
repl = {}
randsyms = r.free_symbols - {z}
while randsyms:
# Only randomly generated parameters are checked.
for n, a in enumerate(randsyms):
repl[a] = randcplx(n)/div
if not any([b.is_Integer and b <= 0 for b in Tuple(*bq).subs(repl)]):
break
[a, b, c, d] = [2, -1, 3, 1]
if lowerplane:
[a, b, c, d] = [2, -2, 3, -1]
return tn(
hyper(ap, bq, z).subs(repl),
r.replace(exp_polar, exp).subs(repl),
z, a=a, b=b, c=c, d=d)
def test_hyperexpand_bases():
assert hyperexpand(hyper([2], [a], z)) == \
a + z**(-a + 1)*(-a**2 + 3*a + z*(a - 1) - 2)*exp(z)* \
lowergamma(a - 1, z) - 1
# TODO [a+1, a-S.Half], [2*a]
assert hyperexpand(hyper([1, 2], [3], z)) == -2/z - 2*log(-z + 1)/z**2
assert hyperexpand(hyper([S.Half, 2], [S(3)/2], z)) == \
-1/(2*z - 2) + atanh(sqrt(z))/sqrt(z)/2
assert hyperexpand(hyper([S(1)/2, S(1)/2], [S(5)/2], z)) == \
(-3*z + 3)/4/(z*sqrt(-z + 1)) \
+ (6*z - 3)*asin(sqrt(z))/(4*z**(S(3)/2))
assert hyperexpand(hyper([1, 2], [S(3)/2], z)) == -1/(2*z - 2) \
- asin(sqrt(z))/(sqrt(z)*(2*z - 2)*sqrt(-z + 1))
assert hyperexpand(hyper([-S.Half - 1, 1, 2], [S.Half, 3], z)) == \
sqrt(z)*(6*z/7 - S(6)/5)*atanh(sqrt(z)) \
+ (-30*z**2 + 32*z - 6)/35/z - 6*log(-z + 1)/(35*z**2)
assert hyperexpand(hyper([1 + S.Half, 1, 1], [2, 2], z)) == \
-4*log(sqrt(-z + 1)/2 + S(1)/2)/z
# TODO hyperexpand(hyper([a], [2*a + 1], z))
# TODO [S.Half, a], [S(3)/2, a+1]
assert hyperexpand(hyper([2], [b, 1], z)) == \
z**(-b/2 + S(1)/2)*besseli(b - 1, 2*sqrt(z))*gamma(b) \
+ z**(-b/2 + 1)*besseli(b, 2*sqrt(z))*gamma(b)
# TODO [a], [a - S.Half, 2*a]
def test_meijerg_lookup():
from sympy import uppergamma, Si, Ci
assert hyperexpand(meijerg([a], [], [b, a], [], z)) == \
z**b*exp(z)*gamma(-a + b + 1)*uppergamma(a - b, z)
assert hyperexpand(meijerg([0], [], [0, 0], [], z)) == \
exp(z)*uppergamma(0, z)
assert can_do_meijer([a], [], [b, a + 1], [])
assert can_do_meijer([a], [], [b + 2, a], [])
assert can_do_meijer([a], [], [b - 2, a], [])
assert hyperexpand(meijerg([a], [], [a, a, a - S(1)/2], [], z)) == \
-sqrt(pi)*z**(a - S(1)/2)*(2*cos(2*sqrt(z))*(Si(2*sqrt(z)) - pi/2)
- 2*sin(2*sqrt(z))*Ci(2*sqrt(z))) == \
hyperexpand(meijerg([a], [], [a, a - S(1)/2, a], [], z)) == \
hyperexpand(meijerg([a], [], [a - S(1)/2, a, a], [], z))
assert can_do_meijer([a - 1], [], [a + 2, a - S(3)/2, a + 1], [])
def _eval_power(self, expt):
"""
b is I = sqrt(-1)
e is symbolic object but not equal to 0, 1
I**r -> (-1)**(r/2) -> exp(r/2*Pi*I) -> sin(Pi*r/2) + cos(Pi*r/2)*I, r is decimal
I**0 mod 4 -> 1
I**1 mod 4 -> I
I**2 mod 4 -> -1
I**3 mod 4 -> -I
"""
if isinstance(expt, Number):
if isinstance(expt, Integer):
expt = expt.p % 4
if expt == 0:
return S.One
if expt == 1:
return S.ImaginaryUnit
if expt == 2:
return -S.One
return -S.ImaginaryUnit
return (S.NegativeOne)**(expt*S.Half)
return
def test_quantum_fourier():
assert QFT(0, 3).decompose() == \
SwapGate(0, 2)*HadamardGate(0)*CGate((0,), PhaseGate(1)) * \
HadamardGate(1)*CGate((0,), TGate(2))*CGate((1,), PhaseGate(2)) * \
HadamardGate(2)
assert IQFT(0, 3).decompose() == \
HadamardGate(2)*CGate((1,), RkGate(2, -2))*CGate((0,), RkGate(2, -3)) * \
HadamardGate(1)*CGate((0,), RkGate(1, -2))*HadamardGate(0)*SwapGate(0, 2)
assert represent(QFT(0, 3), nqubits=3) == \
Matrix([[exp(2*pi*I/8)**(i*j % 8)/sqrt(8) for i in range(8)] for j in range(8)])
assert QFT(0, 4).decompose() # non-trivial decomposition
assert qapply(QFT(0, 3).decompose()*Qubit(0, 0, 0)).expand() == qapply(
HadamardGate(0)*HadamardGate(1)*HadamardGate(2)*Qubit(0, 0, 0)
).expand()
def test_wavefunction():
a = 1/Z
R = {
(1, 0): 2*sqrt(1/a**3) * exp(-r/a),
(2, 0): sqrt(1/(2*a**3)) * exp(-r/(2*a)) * (1 - r/(2*a)),
(2, 1): S(1)/2 * sqrt(1/(6*a**3)) * exp(-r/(2*a)) * r/a,
(3, 0): S(2)/3 * sqrt(1/(3*a**3)) * exp(-r/(3*a)) *
(1 - 2*r/(3*a) + S(2)/27 * (r/a)**2),
(3, 1): S(4)/27 * sqrt(2/(3*a**3)) * exp(-r/(3*a)) *
(1 - r/(6*a)) * r/a,
(3, 2): S(2)/81 * sqrt(2/(15*a**3)) * exp(-r/(3*a)) * (r/a)**2,
(4, 0): S(1)/4 * sqrt(1/a**3) * exp(-r/(4*a)) *
(1 - 3*r/(4*a) + S(1)/8 * (r/a)**2 - S(1)/192 * (r/a)**3),
(4, 1): S(1)/16 * sqrt(5/(3*a**3)) * exp(-r/(4*a)) *
(1 - r/(4*a) + S(1)/80 * (r/a)**2) * (r/a),
(4, 2): S(1)/64 * sqrt(1/(5*a**3)) * exp(-r/(4*a)) *
(1 - r/(12*a)) * (r/a)**2,
(4, 3): S(1)/768 * sqrt(1/(35*a**3)) * exp(-r/(4*a)) * (r/a)**3,
}
for n, l in R:
assert simplify(R_nl(n, l, r, Z) - R[(n, l)]) == 0
def chapman(z, Nm, Hm, H_O, exp=NP.exp):
"""
Return Chapman function electron density at height *z*, maximum
electron density at the F-peak *Nm*, height at the maximum *Hm*,
and the scale height of atomic oxygen *H_O*. The exponential
function can be overridden with *exp*.
"""
return Nm * exp((1 - ((z - Hm) / H_O) - exp(-(z - Hm) / H_O)) / 2)
def chapman_sym(z, Nm, Hm, H_O):
"""
Return symbolic (i.e., :module:`sympy`) Chapman electron density
profile function.
"""
return chapman(z, Nm, Hm, H_O, exp=SYM.exp)
def paper_cram_error(degree, t0, prec=200):
from ..partialfrac import t, customre
from sympy import re, exp, nsolve, diff
paper_part_frac = get_paper_part_frac(degree).replace(customre, re)
E = paper_part_frac - exp(-t)
return E.subs(t, nsolve(diff(E, t), t0, prec=prec))
def _plot(args):
n = args.degree
from sympy import exp
from .. import plot_in_terminal
rat_func = create_expression(n)
print(rat_func)
plot_in_terminal(rat_func - exp(-t), (0, 100), prec=20)
def main():
setup_matplotlib_rc()
expr = get_CRAM_from_cache(degree, prec)
c = 0.6*degree
# Get the translated approximation on [-1, 1]. This is similar logic from CRAM_exp().
n, d = map(Poly, fraction(expr))
inv = -c*(t + 1)/(t - 1)
p, q = map(lambda i: Poly(i, t), fraction(inv))
n, d = n.transform(p, q), d.transform(p, q)
rat_func = n/d.TC()/(d/d.TC())
rat_func = rat_func.evalf(prec)
plt.clf()
fig, ax = plt.subplots()
fig.set_size_inches(4, 4)
plot_in_terminal(expr - exp(-t), (0, 100), prec=prec, points=points,
axes=ax)
# Hide grid lines
ax.grid(False)
# Hide axes ticks
ax.set_xticks([])
ax.set_yticks([])
plt.axis('off')
plt.tight_layout()
plt.savefig('logo.png', transparent=True, bbox_inches='tight', pad_inches=0)
