In [23]:
%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
import numpy.matlib
from pylab import *
import scipy.special
import cmath
In [24]:
# Set the numerical parameters
RMin = 0.0
RMax = 12.0
RMatch = RMax
lOrbital = 1
Step = 0.01
Dim = (RMax-RMin)/Step
# Define your physical parameters
hc = 197.32705
beta = 0.0478450
Rws = 1.2 * 10.0**(0.3333) # This assumes that 10Be is independent from that other neutron
aws = 0.65
V0 = -61.1
E = 0.1
# Set your initial conditions
y0 = 0
yprime0 = 1
# Here we set up the harmonic oscillator potential
def f(R,v):
F = np.matlib.zeros((2,2))
F[0,1] = 1
F[1,0] = lOrbital*(lOrbital+1)/(R*R) + beta * ( V0/(1+exp( (R-Rws)/aws )) - E);
return F*v
In [25]:
# Define your Hankel functions in terms of the cylindrical Bessel functions
def F(l,r):
return (math.pi*r/2)**(0.5)*scipy.special.jv(l,r)
def G(l,r):
return -(math.pi*r/2)**(0.5)*scipy.special.yv(l,r)
def Hplus(l,r):
return G(l,r) + 1j*F(l,r)
def Hminus(l,r):
return G(l,r) - 1j*F(l,r)
# And also define their derivatives
def Fprime(l,r):
return 0.5*math.pi/2*(math.pi*r/2)**(-0.5)*scipy.special.jv(l,r) + (math.pi*r/2)**(0.5)*scipy.special.jvp(l, r, 1)
def Gprime(l,r):
return -0.5*math.pi/2*(math.pi*r/2)**(-0.5)*scipy.special.yv(l,r) - (math.pi*r/2)**(0.5)*scipy.special.yvp(l, r, 1)
def Hplusp(l,r):
return Gprime(l,r) + 1j*Fprime(l,r)
def Hminusp(l,r):
return Gprime(l,r) - 1j*Fprime(l,r)
In [26]:
# Runge-Kutta
# Start the algorithm with your initial conditions
y = np.matlib.zeros((2,int(Dim)))
y[0,0] = y0
y[1,0] = yprime0
r = np.linspace(RMin,RMax,Dim)
for i in xrange(1,int(Dim)):
# print r[i]
k1 = f(r[i],y[:,i-1])
k2 = f(r[i] + Step/2, y[:,i-1] + Step/2*k1)
k3 = f(r[i] + Step/2, y[:,i-1] + Step/2*k2)
k4 = f(r[i] + Step, y[:,i-1] + Step*k3)
y[:,i] = y[:,i-1] + Step/6 * (k1 + 2*k2 + 2*k3 + k4)
# print k1, k2, k3, k4
# print y
In [27]:
# Match the solution for r < RMatch (which you get from Runge-Kutta) to the asymptotic solution for r > RMatch
mat = np.zeros((2,2), dtype=complex)
mat[0,0] = y[0,int(Dim-1)]
mat[1,0] = y[1,int(Dim-1)]
mat[0,1] = Hplus(lOrbital, RMatch)
mat[1,1] = Hplusp(lOrbital, RMatch)
vec = np.zeros((2,1), dtype=complex)
vec[0,0] = Hminus(lOrbital, RMatch)
vec[1,0] = Hminusp(lOrbital, RMatch)
params = np.dot(np.linalg.inv(mat),vec)
cl = params[0,0] # This is the scaling factor to account for your derivative guess at R=0, such that u'(0)=cl*y'(0)
Sl = params[1,0] # This is your S-matrix value
print 'The S-matrix element is $S_L$ = ', Sl # See eqns 6.1.7 and 6.1.8 in Filomena's book for details
print 'which has magnitude', abs(Sl), 'and corresponds to a phase $\delta_L = $', 1/(2*1j)*cmath.log(Sl)
In [32]:
plt.plot(r, (cl*y.A[0,:]).real, 'b-')
plt.plot(r+RMatch, (Hminus(lOrbital,r+RMatch)-Sl*Hplus(lOrbital,r+RMatch)).real, 'b--')
plt.xlabel(r'$r\ (fm)$')
plt.ylabel(r'$u_L(r)$')
plt.title(r'Scattering wavefunction for L = %s (real part)' %(lOrbital))
#plt.savefig('realpart.pdf')
plt.savefig('realpart.png')
plt.show()
In [29]:
plt.plot(r, (cl*y.A[0,:]).imag, 'b-')
plt.plot(r+RMatch, (Hminus(lOrbital,r+RMatch)-Sl*Hplus(lOrbital,r+RMatch)).imag, 'b--')
plt.xlabel(r'$r\ (fm)$')
plt.ylabel(r'$u_L(r)$')
plt.title(r'Scattering wavefunction for L = %s (imaginary part)' %(lOrbital))
#plt.savefig('imagpart.pdf')
#plt.savefig('imagpart.png')
plt.show()
In [30]:
plt.plot(r, abs(cl*y.A[0,:]), 'b-')
plt.plot(r+RMatch, abs(Hminus(lOrbital,r+RMatch)-Sl*Hplus(lOrbital,r+RMatch)), 'b--')
plt.xlabel(r'$r\ (fm)$')
plt.ylabel(r'$u_L(r)$')
plt.title(r'Scattering wavefunction for L = %s (magnitude)' %(lOrbital))
#plt.savefig('magnitude.pdf')
#plt.savefig('magnitude.png')
plt.show()
In [31]:
print cl*y.A[0,int(Dim-1)]
print Hminus(lOrbital,RMatch)-Sl*Hplus(lOrbital,RMatch)
print cl*y.A[1,int(Dim-1)]
print Hminusp(lOrbital,RMatch)-Sl*Hplusp(lOrbital,RMatch)