In [1]:

    

import numpy as np
import scipy.special as num


    

In [6]:

    

A=np.array([0,0,0])*1.8897259886
C=np.array([1.13858,0,0])*1.8897259886
B=np.array([0,0,0])*1.8897259886
D=np.array([1.13858,0,0])*1.8897259886
decay_A=2.263101
decay_B=2.263101
decay_C=47.105518
decay_D=47.105518
Normalisation_A=(2*decay_A/np.pi)**0.75
Normalisation_B=(2*decay_B/np.pi)**0.75
Normalisation_C=(2*decay_C/np.pi)**0.75
Normalisation_D=(2*decay_D/np.pi)**0.75
zeta=decay_A+decay_B
nu=decay_C+decay_D
xi=decay_A*decay_B/zeta
rho=zeta*nu/(zeta+nu)
P=(decay_A*A+decay_B*B)/zeta
Q=1/nu*(decay_C*C+decay_D*D)
T=rho*np.dot(P-Q,P-Q)


    

In [3]:

    

def F(T):
    if T!=0.0:
        return 0.5*num.erf(np.sqrt(T))*(np.pi/T)**0.5
    else:
        return 1


    

In [4]:

    

def K(decay1,decay2,R1,R2):
    return np.sqrt(2)*np.pi**1.25/(decay1+decay2)*np.exp(-decay1*decay2/(decay1+decay2)*np.dot(R1-R2,R1-R2))


    

In [8]:

    

ssss=(rho+nu)**(-0.5)*K(decay_A,decay_B,A,B)*K(decay_C,decay_D,C,D)*F(T)


    

In [9]:

    

ssss


    

    
Out[9]:





0.0016382887691304441

In [ ]: