In [1]:

    

import sys
import sympy
from sympy import *
from sympy import Rational as frac
init_printing()

from memoize import memoize

vartheta, varphi = sympy.symbols('vartheta, varphi', real=True)
v, nu, delta, R_L = sympy.symbols('v, nu, delta, R_L', real=True)
Sigma_n, Sigma_ell, Sigma_lambda = sympy.symbols('Sigma_n, Sigma_ell, Sigma_lambda', real=True)
S_n, S_ell, S_lambda = sympy.symbols('S_n, S_ell, S_lambda', real=True)

nHat = Matrix([1,0,0])
lambdaHat = Matrix([0,1,0])
ellHat = Matrix([0,0,1])

NHat = Matrix([sin(vartheta)*cos(varphi), sin(vartheta)*sin(varphi), cos(vartheta)])
PHat = Matrix([sin(varphi), -cos(varphi), 0])
QHat = Matrix([cos(vartheta)*cos(varphi), cos(vartheta)*sin(varphi), -sin(vartheta)])

Sigma = Matrix([Sigma_n, Sigma_lambda, Sigma_ell])
S = Matrix([S_n, S_lambda, S_ell])

def hComplex(jComponent, kComponent):
    return sympy.trigsimp((nu/R_L)*( (jComponent.dot(PHat)*kComponent.dot(PHat) - jComponent.dot(QHat)*kComponent.dot(QHat))
                                    - sympy.I*(jComponent.dot(PHat)*kComponent.dot(QHat) + jComponent.dot(QHat)*kComponent.dot(PHat)) ) )


    

In [25]:

    

# Memoize (cache) this function for speed
@memoize
def sYlm(s, ell, m, vartheta, varphi):
    r = sympy.Symbol('r', integer=True)
    return sympy.simplify( sympy.sqrt(frac( factorial(ell+m)*factorial(ell-m),factorial(ell+s)*factorial(ell-s) )
                                      * (2*ell+1)/(4*sympy.pi) )
                           * exp(sympy.I*m*varphi)
                           * summation(binomial(ell-s,r) * binomial(ell+s,r+s-m) * (-1)**(ell-r-s-m)
                                       *sin(vartheta/2)**(2*ell-2*r-s+m) * cos(vartheta/2)**(2*r+s-m),
                                       (r,0,ell-s)) )
@memoize
def sYlmConjugate(s, ell, m, vartheta, varphi):
    r = sympy.Symbol('r', integer=True)
    return sympy.simplify( sympy.sqrt(frac( factorial(ell+m)*factorial(ell-m),factorial(ell+s)*factorial(ell-s) )
                                      * (2*ell+1)/(4*sympy.pi) )
                           * exp(-sympy.I*m*varphi)
                           * summation(binomial(ell-s,r) * binomial(ell+s,r+s-m) * (-1)**(ell-r-s-m)
                                       *sin(vartheta/2)**(2*ell-2*r-s+m) * cos(vartheta/2)**(2*r+s-m),
                                       (r,0,ell-s)) )
# Calculate these now, just so everything else will be fast
for ell in range(2,4+1):
    for m in range(-ell,ell+1):
        tmp = sYlm(-2,ell,m,vartheta,varphi)
        tmp = sYlmConjugate(-2,ell,m,vartheta,varphi)


    

In [3]:

    

h0NS = 2*v**2*sympy.simplify(hComplex(lambdaHat, lambdaHat)-hComplex(nHat, nHat))


    

In [4]:

    

h2Ss = 2*v**4*hComplex(Sigma.cross(NHat), nHat)


    

In [5]:

    

h3SO = -2*v**5*sympy.simplify((5*S_ell + 3*delta*Sigma_ell)*hComplex(lambdaHat, lambdaHat)
                              - 6*(2*S_ell + delta*Sigma_ell)*hComplex(nHat, nHat))


    

In [6]:

    

h3Sp = -2*v**5*sympy.simplify(2*hComplex(lambdaHat, nHat.cross(S + delta*Sigma))
                              + hComplex(nHat, lambdaHat.cross(9*S + 5*delta*Sigma)) )


    

In [7]:

    

h3Spv = -2*v**5*sympy.simplify(2*NHat.dot(lambdaHat)*hComplex((S + delta*Sigma).cross(NHat), nHat)
                               + NHat.dot(nHat)*hComplex((S + delta*Sigma).cross(NHat), lambdaHat) )


    

In [39]:

    

def IntegrateFunctionOnTheSphere(func):
    return integrate( integrate( sympy.expand(func*sin(vartheta)), (varphi, 0, 2*pi)), (vartheta, 0, pi))
def IntegrateFunctionAgainstsYlm(func, ell, m, s=-2):
    func = sympy.expand(func*sYlmConjugate(s,ell,m,vartheta,varphi)*sin(vartheta))
    return integrate( integrate( func, (varphi, 0, 2*pi)), (vartheta, 0, pi))


    

In [2]:

    

integrate(sin(vartheta)*cos(vartheta/2)**8, vartheta)


    

    




---------------------------------------------------------------------------
KeyboardInterrupt                         Traceback (most recent call last)
<ipython-input-2-606966c9472f> in <module>()
----> 1 integrate(sin(vartheta)*cos(vartheta/2)**8, vartheta)

/Users/boyle/.continuum/anaconda/lib/python2.7/site-packages/sympy/utilities/decorator.pyc in threaded_func(expr, *args, **kwargs)
     33                                       func(expr.rhs, *args, **kwargs))
     34             else:
---> 35                 return func(expr, *args, **kwargs)
     36 
     37     return threaded_func

/Users/boyle/.continuum/anaconda/lib/python2.7/site-packages/sympy/integrals/integrals.pyc in integrate(*args, **kwargs)
   1287     if isinstance(integral, Integral):
   1288         return integral.doit(deep=False, meijerg=meijerg, conds=conds,
-> 1289                              risch=risch, manual=manual)
   1290     else:
   1291         return integral

/Users/boyle/.continuum/anaconda/lib/python2.7/site-packages/sympy/integrals/integrals.pyc in doit(self, **hints)
    553                     function, xab[0],
    554                     meijerg=meijerg1, risch=risch, manual=manual,
--> 555                     conds=conds)
    556                 if antideriv is None and meijerg1 is True:
    557                     ret = try_meijerg(function, xab)

/Users/boyle/.continuum/anaconda/lib/python2.7/site-packages/sympy/integrals/integrals.pyc in _eval_integral(self, f, x, meijerg, risch, manual, conds)
    917                 try:
    918                     if conds == 'piecewise':
--> 919                         h = heurisch_wrapper(g, x, hints=[])
    920                     else:
    921                         h = heurisch(g, x, hints=[])

/Users/boyle/.continuum/anaconda/lib/python2.7/site-packages/sympy/integrals/heurisch.pyc in heurisch_wrapper(f, x, rewrite, hints, mappings, retries, degree_offset, unnecessary_permutations)
    126 
    127     res = heurisch(f, x, rewrite, hints, mappings, retries, degree_offset,
--> 128                    unnecessary_permutations)
    129     if not isinstance(res, Basic):
    130         return res

/Users/boyle/.continuum/anaconda/lib/python2.7/site-packages/sympy/integrals/heurisch.pyc in heurisch(f, x, rewrite, hints, mappings, retries, degree_offset, unnecessary_permutations)
    564 
    565     if not (F.atoms(Symbol) - set(V)):
--> 566         solution = _integrate('Q')
    567 
    568         if solution is None:

/Users/boyle/.continuum/anaconda/lib/python2.7/site-packages/sympy/integrals/heurisch.pyc in _integrate(field)
    519         # revise this code.
    520         candidate = poly_part/poly_denom + Add(*log_part)
--> 521         h = F - _derivation(candidate) / denom
    522         raw_numer = h.as_numer_denom()[0]
    523 

/Users/boyle/.continuum/anaconda/lib/python2.7/site-packages/sympy/integrals/heurisch.pyc in _derivation(h)
    382     numers = [ cancel(denom*g) for g in diffs ]
    383     def _derivation(h):
--> 384         return Add(*[ d * h.diff(v) for d, v in zip(numers, V) ])
    385 
    386     def _deflation(p):

/Users/boyle/.continuum/anaconda/lib/python2.7/site-packages/sympy/core/expr.pyc in diff(self, *symbols, **assumptions)
   2773         new_symbols = list(map(sympify, symbols))  # e.g. x, 2, y, z
   2774         assumptions.setdefault("evaluate", True)
-> 2775         return Derivative(self, *new_symbols, **assumptions)
   2776 
   2777     ###########################################################################

/Users/boyle/.continuum/anaconda/lib/python2.7/site-packages/sympy/core/function.pyc in __new__(cls, expr, *variables, **assumptions)
   1100                     old_v = v
   1101                     v = new_v
-> 1102                 obj = expr._eval_derivative(v)
   1103                 nderivs += 1
   1104                 if not is_symbol:

/Users/boyle/.continuum/anaconda/lib/python2.7/site-packages/sympy/core/add.pyc in _eval_derivative(self, s)
    352 
    353     def _eval_derivative(self, s):
--> 354         return self.func(*[f.diff(s) for f in self.args])
    355 
    356     def _eval_nseries(self, x, n, logx):

/Users/boyle/.continuum/anaconda/lib/python2.7/site-packages/sympy/core/expr.pyc in diff(self, *symbols, **assumptions)
   2773         new_symbols = list(map(sympify, symbols))  # e.g. x, 2, y, z
   2774         assumptions.setdefault("evaluate", True)
-> 2775         return Derivative(self, *new_symbols, **assumptions)
   2776 
   2777     ###########################################################################

/Users/boyle/.continuum/anaconda/lib/python2.7/site-packages/sympy/core/function.pyc in __new__(cls, expr, *variables, **assumptions)
   1050         if evaluate:
   1051             symbol_set = set(sc[0] for sc in variable_count if sc[0].is_Symbol)
-> 1052             if symbol_set.difference(expr.free_symbols):
   1053                 return S.Zero
   1054 

/Users/boyle/.continuum/anaconda/lib/python2.7/site-packages/sympy/core/basic.pyc in free_symbols(self)
    558         method."""
    559         union = set.union
--> 560         return reduce(union, [arg.free_symbols for arg in self.args], set())
    561 
    562     @property

/Users/boyle/.continuum/anaconda/lib/python2.7/site-packages/sympy/core/symbol.pyc in free_symbols(self)
    134         return False
    135 
--> 136     @property
    137     def free_symbols(self):
    138         return set([self])

KeyboardInterrupt: 

In [ ]:

    




    

In [ ]: