finsihed DWBA zero range, need optimization
This commit is contained in:
Ryan@Home
parent
3b67db12e2
commit
337adb2446
|
|
@ -1,474 +1,53 @@
|
||||||
#!/usr/bin/env python3
|
#!/usr/bin/env python3
|
||||||
|
|
||||||
from boundState import BoundState
|
|
||||||
from solveSE import WoodsSaxonPot, CoulombPotential, SpinOrbit_Pot, WS_SurfacePot
|
|
||||||
import matplotlib.pyplot as plt
|
|
||||||
|
|
||||||
# boundState = BoundState(16, 8, 1, 0, 1, 0, 0.5, -3.273)
|
|
||||||
# boundState = BoundState(16, 8, 1, 0, 0, 2, 2.5, -4.14)
|
|
||||||
# boundState.SetPotential(1.10, 0.65, -6, 1.25, 0.65, 1.30)
|
|
||||||
# boundState.FindPotentialDepth(-75, -50, 0.1)
|
|
||||||
# # boundState.PrintWF()
|
|
||||||
# boundState.PlotBoundState()
|
|
||||||
|
|
||||||
# exit()
|
|
||||||
|
|
||||||
from distortedWave import DistortedWave
|
|
||||||
|
|
||||||
# dw = DistortedWave("60Ni", "p", 30)
|
|
||||||
# dw.ClearPotential()
|
|
||||||
# dw.AddPotential(WoodsSaxonPot(-47.937-2.853j, 1.20, 0.669), False)
|
|
||||||
# dw.AddPotential(WS_SurfacePot(-6.878j, 1.28, 0.550), False)
|
|
||||||
# dw.AddPotential(SpinOrbit_Pot(-5.250 + 0.162j, 1.02, 0.590), False)
|
|
||||||
# dw.AddPotential(CoulombPotential(1.258), False)
|
|
||||||
|
|
||||||
# dw = DistortedWave("60Ni", "d", 60)
|
|
||||||
# dw.PrintInput()
|
|
||||||
# dw.ClearPotential()
|
|
||||||
# dw.AddPotential(WoodsSaxonPot(-81.919, 1.15, 0.768), False)
|
|
||||||
# dw.AddPotential(WoodsSaxonPot(-4.836j, 1.33, 0.464), False)
|
|
||||||
# dw.AddPotential(WS_SurfacePot(-8.994j, 1.373, 0.774), False)
|
|
||||||
# dw.AddPotential(SpinOrbit_Pot(-3.557, 0.972, 1.011), False)
|
|
||||||
# dw.AddPotential(CoulombPotential(1.303), False)
|
|
||||||
|
|
||||||
|
|
||||||
# dw.CalScatteringMatrix()
|
|
||||||
# dw.PrintScatteringMatrix()
|
|
||||||
|
|
||||||
# dw.PlotDCSUnpolarized(180, 1)
|
|
||||||
|
|
||||||
# exit()
|
|
||||||
|
|
||||||
# for i in range(1, 19):
|
|
||||||
# theta = 10*i
|
|
||||||
# # ruth = dw.RutherFord(theta)
|
|
||||||
# # coulAmp = dw.CoulombScatterintAmp(theta)
|
|
||||||
# dw.CalLegendre(theta)
|
|
||||||
# nuAmp1 = dw.NuclearScatteringAmp(-0.5, 0.5, 14)
|
|
||||||
# nuAmp2 = dw.NuclearScatteringAmp(0.5, -0.5, 14)
|
|
||||||
# # dsc = dw.DCSUnpolarized(theta, 14)
|
|
||||||
# # print(f"{theta:3.0f}, {nuAmp1:15.5f}, {nuAmp2:15.5f}, {dsc:10.6f}, {ruth:10.6f}")
|
|
||||||
# print(f"{theta:3.0f}, {nuAmp1:15.5f}, {nuAmp2:15.5f}")
|
|
||||||
|
|
||||||
|
|
||||||
import sys, os
|
|
||||||
import re
|
|
||||||
import numpy as np
|
|
||||||
from scipy.integrate import simpson
|
|
||||||
import matplotlib.pyplot as plt
|
|
||||||
|
|
||||||
sys.path.append(os.path.join(os.path.dirname(__file__), '../Cleopatra'))
|
|
||||||
from IAEANuclearData import IsotopeClass
|
|
||||||
|
|
||||||
from clebschGordan import clebsch_gordan, quantum_factorial, obeys_triangle_rule
|
|
||||||
from sympy.physics.quantum.cg import wigner_9j
|
|
||||||
|
|
||||||
# Woods-Saxon
|
|
||||||
v = 0
|
|
||||||
r0 = 0
|
|
||||||
a = 0
|
|
||||||
vi = 0
|
|
||||||
ri0 = 0
|
|
||||||
ai = 0
|
|
||||||
# Woods-Saxon Surface
|
|
||||||
vsi = 0
|
|
||||||
rsi0 = 0
|
|
||||||
asi = 0
|
|
||||||
# Spin-orbit
|
|
||||||
vso = 0
|
|
||||||
rso0 = 0
|
|
||||||
aso = 0
|
|
||||||
vsoi = 0
|
|
||||||
rsoi0 = 0
|
|
||||||
asoi = 0
|
|
||||||
# Coulomb
|
|
||||||
rc0 = 0
|
|
||||||
|
|
||||||
def AnCai(A : int, Z : int, E : float):
|
|
||||||
global v, r0, a, vi, ri0, ai, vsi, rsi0, asi, vso, rso0, aso, vsoi, rsoi0, asoi, rc0
|
|
||||||
|
|
||||||
A3 = A**(1./3.)
|
|
||||||
v = 91.85 - 0.249*E + 0.000116*pow(E,2) + 0.642 * Z / A3
|
|
||||||
r0 = 1.152 - 0.00776 / A3
|
|
||||||
a = 0.719 + 0.0126 * A3
|
|
||||||
|
|
||||||
vi = 1.104 + 0.0622 * E
|
|
||||||
ri0 = 1.305 + 0.0997 / A3
|
|
||||||
ai = 0.855 - 0.1 * A3
|
|
||||||
|
|
||||||
vsi = 10.83 - 0.0306 * E
|
|
||||||
rsi0 = 1.334 + 0.152 / A3
|
|
||||||
asi = 0.531 + 0.062 * A3
|
|
||||||
|
|
||||||
vso = 3.557
|
|
||||||
rso0 = 0.972
|
|
||||||
aso = 1.011
|
|
||||||
|
|
||||||
vsoi = 0.0
|
|
||||||
rsoi0 = 0.0
|
|
||||||
asoi = 0.0
|
|
||||||
|
|
||||||
rc0 = 1.303
|
|
||||||
|
|
||||||
def Koning(A : int, Z : int, E : float, Zproj : float):
|
|
||||||
global v, r0, a, vi, ri0, ai, vsi, rsi0, asi, vso, rso0, aso, vsoi, rsoi0, asoi, rc0
|
|
||||||
|
|
||||||
N = A-Z
|
|
||||||
A3 = A**(1./3.)
|
|
||||||
|
|
||||||
vp1 = 59.3 + 21.*(N-Z)/A - 0.024*A
|
|
||||||
vn1 = 59.3 - 21.*(N-Z)/A - 0.024*A
|
|
||||||
|
|
||||||
vp2 = 0.007067 + 0.00000423*A
|
|
||||||
vn2 = 0.007228 - 0.00000148*A
|
|
||||||
|
|
||||||
vp3 = 0.00001729 + 0.00000001136 * A
|
|
||||||
vn3 = 0.00001994 - 0.00000002 * A
|
|
||||||
|
|
||||||
vp4 = 7e-9 # = vn4
|
|
||||||
vn4 = vp4
|
|
||||||
|
|
||||||
wp1 = 14.667 + 0.009629*A
|
|
||||||
wn1 = 12.195 + 0.0167*A
|
|
||||||
|
|
||||||
wp2 = 73.55 + 0.0795*A # = wn2
|
|
||||||
wn2 = wp2
|
|
||||||
|
|
||||||
dp1 = 16 + 16.*(N-Z)/A
|
|
||||||
dn1 = 16 - 16.*(N-Z)/A
|
|
||||||
|
|
||||||
dp2 = 0.018 + 0.003802/(1 + np.exp((A-156.)/8)) # = dn2
|
|
||||||
dn2 = dp2
|
|
||||||
|
|
||||||
dp3 = 11.5 # = dn3
|
|
||||||
dn3 = dp3
|
|
||||||
|
|
||||||
vso1 = 5.922 + 0.003 * A
|
|
||||||
vso2 = 0.004
|
|
||||||
|
|
||||||
wso1 = -3.1
|
|
||||||
wso2 = 160
|
|
||||||
|
|
||||||
epf = -8.4075 + 0.01378 *A
|
|
||||||
enf = -11.2814 + 0.02646 *A
|
|
||||||
|
|
||||||
rc = 1.198 + 0.697/pow(A3,2) + 12.995/pow(A3,5)
|
|
||||||
vc = 1.73/rc * Z / A3
|
|
||||||
|
|
||||||
v = vp1*(1 - vp2*(E-epf) + vp3*pow(E-epf,2) - vp4*pow(E-epf,3)) + vc * vp1 * (vp2 - 2*vp3*(E-epf) + 3*vp4*pow(E-epf,2))
|
|
||||||
#neutron
|
|
||||||
if Zproj == 0 :
|
|
||||||
v = vn1*(1 - vn2*(E-enf) + vn3*pow(E-enf,2) - vn4*pow(E-enf,3))
|
|
||||||
|
|
||||||
r0 = 1.3039 - 0.4054 / A3
|
|
||||||
a = 0.6778 - 0.000148 * A
|
|
||||||
|
|
||||||
vi = wp1 * pow(E-epf,2)/(pow(E-epf,2) + pow(wp2,2))
|
|
||||||
if Zproj == 0 :
|
|
||||||
vi = wn1 * pow(E-enf,2)/(pow(E-enf,2) + pow(wn2,2))
|
|
||||||
|
|
||||||
ri0 = 1.3039 - 0.4054 / A3
|
|
||||||
ai = 0.6778 - 0.000148 * A
|
|
||||||
|
|
||||||
vsi = dp1 * pow(E-epf,2)/(pow(E-epf,2)+pow(dp3,2)) * np.exp(-dp2*(E-epf))
|
|
||||||
if Zproj == 0 :
|
|
||||||
vsi = dn1 * pow(E-enf,2)/(pow(E-enf,2)+pow(dn3,2)) * np.exp(-dn2*(E-enf))
|
|
||||||
|
|
||||||
rsi0 = 1.3424 - 0.01585 * A3
|
|
||||||
asi = 0.5187 + 0.0005205 * A
|
|
||||||
if Zproj == 0:
|
|
||||||
asi = 0.5446 - 0.0001656 * A
|
|
||||||
|
|
||||||
vso = vso1 * np.exp(-vso2 * (E-epf))
|
|
||||||
if Zproj == 0:
|
|
||||||
vso = vso1 * np.exp(-vso2 * (E-enf))
|
|
||||||
|
|
||||||
rso0 = 1.1854 - 0.647/A3
|
|
||||||
aso = 0.59
|
|
||||||
|
|
||||||
vsoi = wso1 * pow(E-epf,2)/(pow(E-epf,2)+pow(wso2,2))
|
|
||||||
if Zproj == 0 :
|
|
||||||
vsoi = wso1 * pow(E-enf,2)/(pow(E-enf,2)+pow(wso2,2))
|
|
||||||
|
|
||||||
rsoi0 = 1.1854 - 0.647/A3
|
|
||||||
asoi = 0.59
|
|
||||||
|
|
||||||
rc0 = rc
|
|
||||||
|
|
||||||
def ConvertLSym(LSym :str) -> int:
|
|
||||||
if LSym == "s" :
|
|
||||||
return 0
|
|
||||||
elif LSym == "p" :
|
|
||||||
return 1
|
|
||||||
elif LSym == "d" :
|
|
||||||
return 2
|
|
||||||
elif LSym == "f" :
|
|
||||||
return 3
|
|
||||||
elif LSym == "g" :
|
|
||||||
return 4
|
|
||||||
elif LSym == "h" :
|
|
||||||
return 5
|
|
||||||
elif LSym == "i" :
|
|
||||||
return 6
|
|
||||||
elif LSym == "j" :
|
|
||||||
return 7
|
|
||||||
elif LSym == "k" :
|
|
||||||
return 8
|
|
||||||
else :
|
|
||||||
return -1
|
|
||||||
|
|
||||||
#==========================================
|
|
||||||
|
|
||||||
nu_A = "16O"
|
|
||||||
nu_a = "d"
|
|
||||||
nu_b = "p"
|
|
||||||
nu_B = "17O"
|
|
||||||
ELabPreU = 10 # MeV/u
|
|
||||||
Ex = 0.87
|
|
||||||
J_B = "1/2+"
|
|
||||||
orbital = "1s1/2"
|
|
||||||
|
|
||||||
import time
|
import time
|
||||||
start_time = time.time() # Start the timer
|
from dwba_zr import DWBA_ZR
|
||||||
|
|
||||||
iso = IsotopeClass()
|
haha = DWBA_ZR("16O", "d", "p", "17O", "1/2+", "1s1/2", 0.0, 10)
|
||||||
|
|
||||||
A_A, Z_A = iso.GetAZ(nu_A)
|
# haha.boundState.PlotBoundState()
|
||||||
A_a, Z_a = iso.GetAZ(nu_a)
|
|
||||||
A_b, Z_b = iso.GetAZ(nu_b)
|
|
||||||
A_B, Z_B = iso.GetAZ(nu_B)
|
|
||||||
|
|
||||||
A_x = abs(A_a - A_b)
|
haha.PrintRadialIntegral()
|
||||||
Z_x = abs(Z_a - Z_b)
|
# haha.PlotRadialIntegral()
|
||||||
|
|
||||||
mass_A = iso.GetMassFromSym(nu_A)
|
# haha.PlotDistortedWave(True, 2, 1, 20)
|
||||||
mass_a = iso.GetMassFromSym(nu_a)
|
# haha.PlotDistortedWave(False, 2, 1.5, 20)
|
||||||
mass_b = iso.GetMassFromSym(nu_b)
|
|
||||||
mass_B = iso.GetMassFromSym(nu_B)
|
|
||||||
|
|
||||||
mass_x = iso.GetMassFromAZ( A_x, Z_x)
|
j = 1/2
|
||||||
|
sa = 1
|
||||||
|
sb = 1/2
|
||||||
|
|
||||||
if A_A < A_B : # (d,p)
|
# hehe = haha.Gamma(0, 1, 0, 0.5, 0, 1, 0.5)
|
||||||
A_c = A_A
|
# print(hehe)
|
||||||
Z_c = Z_A
|
|
||||||
BindingEnergy = mass_B - mass_A - mass_x + Ex
|
|
||||||
else: #(p,d)
|
|
||||||
A_c = A_B
|
|
||||||
Z_c = Z_B
|
|
||||||
BindingEnergy = mass_A - mass_B - mass_x
|
|
||||||
|
|
||||||
sym_A = iso.GetSymbol(A_A, Z_A)
|
# haha.PreCalLegendreP(10)
|
||||||
sym_B = iso.GetSymbol(A_B, Z_B)
|
# print(haha.legendrePArray)
|
||||||
|
|
||||||
spin_A_str = iso.GetJpi(A_A, Z_A)
|
# jaja = haha.Beta(-1, 1, -0.5) * haha.ffactor
|
||||||
spin_B_str = J_B
|
# print(jaja)
|
||||||
|
|
||||||
spin_A = float(eval(re.sub(r'[+-]', '', spin_A_str)))
|
# lala = haha.AngDist(10)
|
||||||
spin_B = float(eval(re.sub(r'[+-]', '', J_B)))
|
# print(lala)
|
||||||
|
|
||||||
if A_a == 2 and Z_a == 1:
|
angList = []
|
||||||
spin_a = 1.0
|
xsec = []
|
||||||
spin_b = 0.5
|
|
||||||
else:
|
|
||||||
spin_a = 0.5
|
|
||||||
spin_b = 1.0
|
|
||||||
|
|
||||||
s = 0.5 # spin of x, neutron or proton
|
start_time = time.time()
|
||||||
|
for i in range(0, 181, 5):
|
||||||
|
angList.append(i)
|
||||||
|
kaka = haha.AngDist(i)
|
||||||
|
xsec.append(kaka)
|
||||||
|
print(i, kaka)
|
||||||
|
|
||||||
#=================== digest orbital
|
stop_time = time.time()
|
||||||
match = re.search(r'[a-zA-Z]', orbital) # Find first letter
|
print(f"Total time {(stop_time - start_time) :.2f} sec")
|
||||||
if match:
|
|
||||||
index = match.start() # Get position of the first letter
|
|
||||||
|
|
||||||
node = int(orbital[:index])
|
|
||||||
l_sym = orbital[index:index+1]
|
|
||||||
j_sym = orbital[index+1:]
|
|
||||||
j = eval(j_sym)
|
|
||||||
l = ConvertLSym(l_sym)
|
|
||||||
|
|
||||||
#==== check the angular conservasion
|
import matplotlib.pyplot as plt
|
||||||
passJ = False
|
|
||||||
if obeys_triangle_rule(spin_A, spin_B, j):
|
|
||||||
passJ = True
|
|
||||||
else:
|
|
||||||
print(f"the orbital spin-J ({j}) does not consver J({nu_A}) + J({nu_B}) = {spin_A} + {spin_B}.")
|
|
||||||
|
|
||||||
passS = False
|
plt.plot(angList, xsec)
|
||||||
if obeys_triangle_rule(spin_a, spin_b, s):
|
plt.xlim(-1, 181)
|
||||||
passS = True
|
plt.yscale("log")
|
||||||
else:
|
plt.grid()
|
||||||
print(f"the orbital spin-s ({s})does not consver S({nu_a}) + S({nu_b}) = {spin_a} + {spin_b}.")
|
|
||||||
|
|
||||||
passl = False
|
|
||||||
if obeys_triangle_rule(j, s, l):
|
|
||||||
passl = True
|
|
||||||
else:
|
|
||||||
print(f"the orbital spin-l ({l})does not consver J({j}) + J({s}).")
|
|
||||||
|
|
||||||
if passJ == False or passS == False or passl == False :
|
|
||||||
print("Fail angular momentum conservation.")
|
|
||||||
exit()
|
|
||||||
|
|
||||||
reactionStr = f"{nu_A}({spin_A_str})({nu_a},{nu_b}){nu_B}({Ex:.3f}|{spin_B_str}, {orbital}) @ {ELabPreU:.1f} MeV/u"
|
|
||||||
|
|
||||||
print("==================================================")
|
|
||||||
print(reactionStr)
|
|
||||||
|
|
||||||
Q_value = mass_A + mass_a - mass_b - mass_B - Ex
|
|
||||||
print(f"Transfer Orbtial : {orbital}")
|
|
||||||
print(f"Q-value : {Q_value:10.6f} MeV")
|
|
||||||
print(f"Binding : {BindingEnergy:10.6f} MeV")
|
|
||||||
|
|
||||||
#=================== find the maximum L for partial wave
|
|
||||||
mass_I = mass_A * mass_a / (mass_A + mass_a) # reduced mass of incoming channel
|
|
||||||
hbarc = 197.3269788 # MeV.fm
|
|
||||||
k_I = np.sqrt(2*mass_I * A_a * ELabPreU)/hbarc # wave number of incoming channel
|
|
||||||
touching_Radius = 1.25*(A_A**(1./3) + A_a**(1./3)) + 10 # add 10 fm
|
|
||||||
maxL = int(touching_Radius * k_I) # maximum partial wave
|
|
||||||
print(f"max L : {maxL}")
|
|
||||||
|
|
||||||
#================== Bound state
|
|
||||||
print("====================== Bound state ")
|
|
||||||
boundState = BoundState(A_c, Z_c, A_x, Z_x, node, l, j, BindingEnergy)
|
|
||||||
boundState.SetPotential(1.10, 0.65, -6, 1.25, 0.65, 1.30)
|
|
||||||
boundState.FindPotentialDepth(-70, -55, 0.1)
|
|
||||||
# # boundState.PrintWF()
|
|
||||||
# boundState.PlotBoundState()
|
|
||||||
|
|
||||||
# exit()
|
|
||||||
|
|
||||||
#================== incoming wave function
|
|
||||||
print("====================== Incoming wave function ")
|
|
||||||
AnCai(A_A, Z_A, A_a * ELabPreU)
|
|
||||||
|
|
||||||
dwI = DistortedWave(nu_A, nu_a, ELabPreU * A_a)
|
|
||||||
dwI.maxL = maxL
|
|
||||||
dwI.PrintInput()
|
|
||||||
dwI.ClearPotential()
|
|
||||||
dwI.AddPotential(WoodsSaxonPot(-v, r0, a), False)
|
|
||||||
dwI.AddPotential(WoodsSaxonPot(-1j*vi, ri0, ai), False)
|
|
||||||
dwI.AddPotential(WS_SurfacePot(-1j*vsi, rsi0, asi), False)
|
|
||||||
dwI.AddPotential(SpinOrbit_Pot(-vso , rso0, aso), False)
|
|
||||||
dwI.AddPotential(SpinOrbit_Pot(- 1j* vsoi, rsoi0, asoi), False)
|
|
||||||
dwI.AddPotential(CoulombPotential(rc0), False)
|
|
||||||
dwI.PrintPotentials()
|
|
||||||
|
|
||||||
sm_I, wfu_I = dwI.CalScatteringMatrix()
|
|
||||||
|
|
||||||
dwI.PrintScatteringMatrix()
|
|
||||||
|
|
||||||
# dwI.PlotDistortedWave(1, 1, 20)
|
|
||||||
# dwI.PlotScatteringMatrix()
|
|
||||||
|
|
||||||
#================= outgoing wave function
|
|
||||||
print("====================== Outgoing wave function ")
|
|
||||||
Koning(A_B, Z_B, A_a*ELabPreU + Q_value - Ex, Z_b)
|
|
||||||
|
|
||||||
dwO = DistortedWave(nu_B, nu_b, ELabPreU * A_a + Q_value - Ex)
|
|
||||||
dwO.maxL = maxL
|
|
||||||
dwO.ClearPotential()
|
|
||||||
dwO.AddPotential(WoodsSaxonPot(-v, r0, a), False)
|
|
||||||
dwO.AddPotential(WoodsSaxonPot(-1j*vi, ri0, ai), False)
|
|
||||||
dwO.AddPotential(WS_SurfacePot(-1j*vsi, rsi0, asi), False)
|
|
||||||
dwO.AddPotential(SpinOrbit_Pot(-vso , rso0, aso), False)
|
|
||||||
dwO.AddPotential(SpinOrbit_Pot(- 1j* vsoi, rsoi0, asoi), False)
|
|
||||||
dwO.AddPotential(CoulombPotential(rc0), False)
|
|
||||||
dwO.PrintPotentials()
|
|
||||||
|
|
||||||
sm_O, wfu_O = dwO.CalScatteringMatrix()
|
|
||||||
|
|
||||||
dwO.PrintScatteringMatrix()
|
|
||||||
|
|
||||||
# dwO.PlotDistortedWave(1, 1.5, 20)
|
|
||||||
|
|
||||||
end_time = time.time() # End the timer
|
|
||||||
print(f"Time used {(end_time - start_time) * 1000:.2f} milliseconds")
|
|
||||||
|
|
||||||
#=================== Calculate radial integral
|
|
||||||
print("====================== Calculating Radial integrals")
|
|
||||||
|
|
||||||
def FormatSpin(spin : float) -> str:
|
|
||||||
if int(2*spin) % 2 == 0 :
|
|
||||||
return f"{int(spin):+d}"
|
|
||||||
else:
|
|
||||||
return f"{int(2*spin):+d}/2"
|
|
||||||
|
|
||||||
spin_a = dwI.spin_a
|
|
||||||
spin_b = dwO.spin_a
|
|
||||||
|
|
||||||
radialInt = np.zeros((maxL+1, int(2*spin_a+1), int(2*spin_b+1)), dtype=complex)
|
|
||||||
|
|
||||||
bs = boundState.GetBoundStateWF()
|
|
||||||
|
|
||||||
for L in range(0, maxL+1):
|
|
||||||
for index1 in range(0, len(wfu_I[L])):
|
|
||||||
wf1 = wfu_I[L][index1]
|
|
||||||
for index2 in range(0, len(wfu_O[L])):
|
|
||||||
wf2 = wfu_O[L][index2]
|
|
||||||
# if L == 0 and index1 == 2 and index2 == 1 :
|
|
||||||
# for i in range(0, len(bs)):
|
|
||||||
# if i%50 == 0 :
|
|
||||||
# print(bs[i], wf1[i], wf2[i], bs[i]* wf1[i]* wf2[i])
|
|
||||||
pf1 = np.exp(1j*dwI.CoulombPhaseShift(L))
|
|
||||||
pf2 = np.exp(1j*dwI.CoulombPhaseShift(L))
|
|
||||||
integral = simpson (bs*wf1*wf2, dx=boundState.dr)
|
|
||||||
radialInt[L][index1][index2] = integral * pf1 * pf2
|
|
||||||
|
|
||||||
#print radial integral
|
|
||||||
for index1 in range(0, int(2*spin_a) + 1):
|
|
||||||
for index2 in range(0, int(2*spin_b) + 1):
|
|
||||||
print(f"======================= J1 = L{FormatSpin(index1-spin_a)}, J2 = L{FormatSpin(index2-spin_b)}")
|
|
||||||
for L in range(0, maxL+1):
|
|
||||||
J1 = L + index1 - spin_a
|
|
||||||
J2 = L + index2 - spin_b
|
|
||||||
print(f"{L:2d}, {J1:4.1f}, {J2:4.1f}, {np.real(radialInt[L][index1][index2]):12.4e} + {np.imag(radialInt[L][index1][index2]):12.4e}I")
|
|
||||||
|
|
||||||
# Plot radial integral
|
|
||||||
fig, axes = plt.subplots(int(2*spin_b+1), int(2*spin_a+1), figsize=(6*int(2*spin_a+1), 4*int(2*spin_b+1)))
|
|
||||||
|
|
||||||
for index2 in range(0, int(2*spin_b) + 1):
|
|
||||||
for index1 in range(0, int(2*spin_a) + 1):
|
|
||||||
haha = []
|
|
||||||
l_list = []
|
|
||||||
for L in range(0, maxL+1):
|
|
||||||
J1 = L + index1 - spin_a
|
|
||||||
J2 = L + index2 - spin_b
|
|
||||||
if J1 < 0 or J2 < 0 :
|
|
||||||
continue
|
|
||||||
l_list.append(L)
|
|
||||||
haha.append(radialInt[L][index1][index2])
|
|
||||||
axes[index2, index1].plot(l_list, np.real(haha), label="Real", marker='o')
|
|
||||||
axes[index2, index1].plot(l_list, np.imag(haha), label="Imag", marker='x')
|
|
||||||
axes[index2, index1].legend()
|
|
||||||
axes[index2, index1].set_xlabel('L')
|
|
||||||
axes[index2, index1].set_ylabel('Value')
|
|
||||||
axes[index2, index1].set_title(f'Radial Int. vs L for Spin J1 = L{FormatSpin(index1-spin_a)}, J2 = L{FormatSpin(index2-spin_b)}.')
|
|
||||||
axes[index2, index1].set_xlim(-1, maxL+1)
|
|
||||||
axes[index2, index1].grid()
|
|
||||||
|
|
||||||
plt.tight_layout()
|
|
||||||
plt.show(block=False)
|
plt.show(block=False)
|
||||||
|
|
||||||
input("Press Enter to continue...")
|
input("Press Enter to continue...")
|
||||||
|
|
||||||
|
|
||||||
def Gamma(L1, J1, L2, J2, m, ma, mb):
|
|
||||||
if int(L1 + L2 + l)%2 != 0:
|
|
||||||
return 0
|
|
||||||
else:
|
|
||||||
fact0 = wigner_9j(j, l, s, J1, L1, spin_a, J2, L2, spin_b)
|
|
||||||
if fact0 == 0:
|
|
||||||
return 0
|
|
||||||
else:
|
|
||||||
fact1 = pow(-1, m) * np.power(1j, L1-L2-l) * (2*L2+1) * np.sqrt((2*l+1)*(2*s+1)*(2*L1+1)*(2*J2+1))
|
|
||||||
fact2 = np.sqrt( quantum_factorial(L2-m) / quantum_factorial(L2+m) )
|
|
||||||
fact3 = clebsch_gordan(J2, mb-m, j, m-mb+ma, J1, ma)
|
|
||||||
fact4 = clebsch_gordan(L1, 0, spin_a, ma, J1, ma)
|
|
||||||
fact5 = clebsch_gordan(L2, -m, spin_b, mb, J2, mb-m)
|
|
||||||
fact6 = clebsch_gordan(L1, 0, l, 0, L2, 0)
|
|
||||||
return fact0 * fact1 * fact2 * fact3 * fact4 * fact5 * fact6
|
|
||||||
|
|
||||||
|
|
||||||
def Beta(m, ma, mb, theta_deg):
|
|
||||||
return 0
|
|
||||||
|
|
||||||
|
|
||||||
|
|
|
||||||
|
|
@ -39,6 +39,7 @@ class BoundState(SolvingSE):
|
||||||
|
|
||||||
def FindPotentialDepth(self, Vmin, Vmax, Vstep=1, isPathWhittaker = True):
|
def FindPotentialDepth(self, Vmin, Vmax, Vstep=1, isPathWhittaker = True):
|
||||||
start_time = time.time() # Start the timer
|
start_time = time.time() # Start the timer
|
||||||
|
print(f"Potential Depth Search from {Vmin:.2f} to {Vmax:.2f} MeV, step {Vstep:.2f}")
|
||||||
V0List = np.arange(Vmin, Vmax, Vstep)
|
V0List = np.arange(Vmin, Vmax, Vstep)
|
||||||
lastSolU = []
|
lastSolU = []
|
||||||
minLastSolU = 0
|
minLastSolU = 0
|
||||||
|
|
|
||||||
|
|
@ -53,16 +53,18 @@ class DistortedWave(SolvingSE):
|
||||||
self.ScatMatrix = []
|
self.ScatMatrix = []
|
||||||
self.distortedWaveU = []
|
self.distortedWaveU = []
|
||||||
|
|
||||||
|
tempZeroList = np.zeros(len(self.rpos), dtype=np.complex128)
|
||||||
|
|
||||||
for L in range(0, maxL+1):
|
for L in range(0, maxL+1):
|
||||||
sigma = self.CoulombPhaseShift()
|
# sigma = self.CoulombPhaseShift()
|
||||||
|
|
||||||
temp_ScatMatrix = []
|
temp_ScatMatrix = []
|
||||||
temp_distortedWaveU = []
|
temp_distortedWaveU = []
|
||||||
|
|
||||||
for J in np.arange(L-self.S, L + self.S+1, 1):
|
for J in np.arange(L-self.S, L + self.S+1, 1):
|
||||||
if J < 0:
|
if J < 0 or (L==0 and J != self.S):
|
||||||
temp_ScatMatrix.append(0)
|
temp_ScatMatrix.append(0)
|
||||||
temp_distortedWaveU.append(0)
|
temp_distortedWaveU.append(tempZeroList)
|
||||||
continue
|
continue
|
||||||
|
|
||||||
self.SetLJ(L, J)
|
self.SetLJ(L, J)
|
||||||
|
|
|
||||||
372
Raphael/dwba_zr.py
Executable file
372
Raphael/dwba_zr.py
Executable file
|
|
@ -0,0 +1,372 @@
|
||||||
|
#!/usr/bin/env python3
|
||||||
|
import sys, os
|
||||||
|
import re
|
||||||
|
import numpy as np
|
||||||
|
from scipy.integrate import simpson
|
||||||
|
from scipy.interpolate import interp1d
|
||||||
|
import matplotlib.pyplot as plt
|
||||||
|
import time
|
||||||
|
from sympy import S
|
||||||
|
from sympy.physics.quantum.cg import wigner_9j
|
||||||
|
|
||||||
|
sys.path.append(os.path.join(os.path.dirname(__file__), '../Cleopatra'))
|
||||||
|
from IAEANuclearData import IsotopeClass
|
||||||
|
|
||||||
|
from assLegendreP import associated_legendre_array
|
||||||
|
from clebschGordan import clebsch_gordan, quantum_factorial, obeys_triangle_rule
|
||||||
|
from boundState import BoundState
|
||||||
|
from solveSE import WoodsSaxonPot, CoulombPotential, SpinOrbit_Pot, WS_SurfacePot
|
||||||
|
from distortedWave import DistortedWave
|
||||||
|
import opticalPotentials as op
|
||||||
|
|
||||||
|
class DWBA_ZR:
|
||||||
|
def __init__(self, nu_A:str, nu_a:str, nu_b:str, nu_B:str, JB:str, orbital:str, ExB:float, ELabPerU:float):
|
||||||
|
start_time = time.time()
|
||||||
|
iso = IsotopeClass()
|
||||||
|
|
||||||
|
A_A, Z_A = iso.GetAZ(nu_A)
|
||||||
|
A_a, Z_a = iso.GetAZ(nu_a)
|
||||||
|
A_b, Z_b = iso.GetAZ(nu_b)
|
||||||
|
A_B, Z_B = iso.GetAZ(nu_B)
|
||||||
|
|
||||||
|
self.ELab = A_a * ELabPerU
|
||||||
|
|
||||||
|
mass_A = iso.GetMassFromSym(nu_A)
|
||||||
|
mass_a = iso.GetMassFromSym(nu_a)
|
||||||
|
mass_b = iso.GetMassFromSym(nu_b)
|
||||||
|
mass_B = iso.GetMassFromSym(nu_B)
|
||||||
|
self.ExB = ExB
|
||||||
|
|
||||||
|
# sym_A = iso.GetSymbol(A_A, Z_A)
|
||||||
|
# sym_B = iso.GetSymbol(A_B, Z_B)
|
||||||
|
|
||||||
|
spin_A_str = iso.GetJpi(A_A, Z_A)
|
||||||
|
self.spin_A = float(eval(re.sub(r'[+-]', '', spin_A_str)))
|
||||||
|
self.spin_B = float(eval(re.sub(r'[+-]', '', JB)))
|
||||||
|
|
||||||
|
if A_a == 2 and Z_a == 1:
|
||||||
|
self.spin_a = 1.0
|
||||||
|
self.spin_b = 0.5
|
||||||
|
else:
|
||||||
|
self.spin_a = 0.5
|
||||||
|
self.spin_b = 1.0
|
||||||
|
|
||||||
|
#====== transfering nucleon
|
||||||
|
self.s = 1/2 # spin of x, neutron or proton
|
||||||
|
A_x = abs(A_a - A_b)
|
||||||
|
Z_x = abs(Z_a - Z_b)
|
||||||
|
|
||||||
|
mass_x = iso.GetMassFromAZ( A_x, Z_x)
|
||||||
|
|
||||||
|
#======== core
|
||||||
|
if A_A < A_B : # (d,p)
|
||||||
|
A_c = A_A
|
||||||
|
Z_c = Z_A
|
||||||
|
BindingEnergy = mass_B - mass_A - mass_x + self.ExB
|
||||||
|
else: #(p,d)
|
||||||
|
A_c = A_B
|
||||||
|
Z_c = Z_B
|
||||||
|
BindingEnergy = mass_A - mass_B - self.ExB - mass_x
|
||||||
|
|
||||||
|
#=================== digest orbital
|
||||||
|
match = re.search(r'[a-zA-Z]', orbital) # Find first letter
|
||||||
|
if match:
|
||||||
|
index = match.start() # Get position of the first letter
|
||||||
|
|
||||||
|
node = int(orbital[:index])
|
||||||
|
l_sym = orbital[index:index+1]
|
||||||
|
j_sym = orbital[index+1:]
|
||||||
|
self.j = eval(j_sym)
|
||||||
|
self.l = op.ConvertLSym(l_sym)
|
||||||
|
|
||||||
|
passJ = False
|
||||||
|
if obeys_triangle_rule(self.spin_A, self.spin_B, self.j):
|
||||||
|
passJ = True
|
||||||
|
else:
|
||||||
|
print(f"the orbital spin-J ({self.j}) does not consver J({nu_A}) + J({nu_B}) = {self.spin_A} + {self.spin_B}.")
|
||||||
|
|
||||||
|
passS = False
|
||||||
|
if obeys_triangle_rule(self.spin_a, self.spin_b, self.s):
|
||||||
|
passS = True
|
||||||
|
else:
|
||||||
|
print(f"the orbital spin-s ({self.s}) does not consver S({nu_a}) + S({nu_b}) = {self.spin_a} + {self.spin_b}.")
|
||||||
|
|
||||||
|
passL = False
|
||||||
|
if obeys_triangle_rule(self.j, self.s, self.l):
|
||||||
|
passL = True
|
||||||
|
else:
|
||||||
|
print(f"the orbital spin-L ({self.l}) does not consver J({self.j}) + S({self.s}).")
|
||||||
|
|
||||||
|
self.isSpinBalanced = passJ * passS * passL
|
||||||
|
if self.isSpinBalanced == False :
|
||||||
|
print("Fail angular momentum conservation.")
|
||||||
|
return
|
||||||
|
else:
|
||||||
|
print("All Spin are balance.")
|
||||||
|
|
||||||
|
self.reactionStr = f"{nu_A}({spin_A_str})({nu_a},{nu_b}){nu_B}({ExB:.3f}|{JB}, {orbital}) @ {ELabPerU:.1f} MeV/u"
|
||||||
|
print("==================================================")
|
||||||
|
print(self.reactionStr)
|
||||||
|
|
||||||
|
self.Q_value = mass_A + mass_a - mass_b - mass_B - ExB
|
||||||
|
print(f"Transfer Orbtial : {orbital}")
|
||||||
|
print(f"Q-value : {self.Q_value:10.6f} MeV")
|
||||||
|
print(f"Binding : {BindingEnergy:10.6f} MeV")
|
||||||
|
|
||||||
|
print("====================== Bound state ")
|
||||||
|
self.boundState = BoundState(A_c, Z_c, A_x, Z_x, node, self.l, self.j, BindingEnergy)
|
||||||
|
self.boundState.SetPotential(1.10, 0.65, -6, 1.25, 0.65, 1.30)
|
||||||
|
self.boundState.FindPotentialDepth(-70, -45, 0.5)
|
||||||
|
|
||||||
|
print("====================== Incoming wave function ")
|
||||||
|
op.AnCai(A_A, Z_A, self.ELab)
|
||||||
|
|
||||||
|
self.maxL = 15
|
||||||
|
|
||||||
|
self.dwI = DistortedWave(nu_A, nu_a, self.ELab)
|
||||||
|
self.dwI.maxL = self.maxL
|
||||||
|
self.dwI.PrintInput()
|
||||||
|
self.dwI.ClearPotential()
|
||||||
|
self.dwI.AddPotential(WoodsSaxonPot( -op.v, op.r0, op.a), False)
|
||||||
|
self.dwI.AddPotential(WoodsSaxonPot(-1j*op.vi, op.ri0, op.ai), False)
|
||||||
|
self.dwI.AddPotential(WS_SurfacePot(-1j*op.vsi, op.rsi0, op.asi), False)
|
||||||
|
self.dwI.AddPotential(SpinOrbit_Pot( -op.vso, op.rso0, op.aso), False)
|
||||||
|
self.dwI.AddPotential(SpinOrbit_Pot(-1j*op.vsoi, op.rsoi0, op.asoi), False)
|
||||||
|
self.dwI.AddPotential(CoulombPotential( op.rc0), False)
|
||||||
|
self.dwI.PrintPotentials()
|
||||||
|
|
||||||
|
self.mass_I = self.dwI.mu # reduced mass of incoming channel
|
||||||
|
self.k_I = self.dwI.k # wave number of incoming channel
|
||||||
|
self.maxL = self.dwI.maxL
|
||||||
|
|
||||||
|
sm_I, wfu_I = self.dwI.CalScatteringMatrix()
|
||||||
|
|
||||||
|
self.dwI.PrintScatteringMatrix()
|
||||||
|
|
||||||
|
Ecm_I = self.dwI.Ecm
|
||||||
|
Ecm_O = Ecm_I + self.Q_value
|
||||||
|
Eout = ((Ecm_O + mass_b + mass_B + self.ExB)**2 - (mass_b + mass_B + ExB)**2)/2/mass_B
|
||||||
|
|
||||||
|
print("====================== Outgoing wave function ")
|
||||||
|
op.Koning(A_B, Z_B, self.ELab + self.Q_value - ExB, Z_b)
|
||||||
|
|
||||||
|
self.dwO = DistortedWave(nu_B, nu_b, Eout)
|
||||||
|
self.dwO.spin_A = self.spin_B
|
||||||
|
self.dwO.maxL = self.maxL
|
||||||
|
self.dwO.PrintInput()
|
||||||
|
self.dwO.ClearPotential()
|
||||||
|
self.dwO.AddPotential(WoodsSaxonPot( -op.v, op.r0, op.a), False)
|
||||||
|
self.dwO.AddPotential(WoodsSaxonPot(-1j*op.vi, op.ri0, op.ai), False)
|
||||||
|
self.dwO.AddPotential(WS_SurfacePot(-1j*op.vsi, op.rsi0, op.asi), False)
|
||||||
|
self.dwO.AddPotential(SpinOrbit_Pot( -op.vso, op.rso0, op.aso), False)
|
||||||
|
self.dwO.AddPotential(SpinOrbit_Pot(-1j*op.vsoi, op.rsoi0, op.asoi), False)
|
||||||
|
self.dwO.AddPotential(CoulombPotential( op.rc0), False)
|
||||||
|
|
||||||
|
self.dwO.PrintPotentials()
|
||||||
|
|
||||||
|
sm_O, wfu_O_temp = self.dwO.CalScatteringMatrix()
|
||||||
|
|
||||||
|
#============ rescale the outgoing wave
|
||||||
|
rpos_O_temp = self.dwO.rpos * A_B/A_A
|
||||||
|
self.rpos_O = []
|
||||||
|
rpos_O_filled = False
|
||||||
|
self.wfu_O = []
|
||||||
|
for L in range(0, self.maxL+1):
|
||||||
|
temp_wfu_L = []
|
||||||
|
for k in range(0, int(2*self.spin_b)+1):
|
||||||
|
wfu_O_inter_real = interp1d(rpos_O_temp, np.real(wfu_O_temp[L][k]), kind='cubic')
|
||||||
|
wfu_O_inter_imag = interp1d(rpos_O_temp, np.imag(wfu_O_temp[L][k]), kind='cubic')
|
||||||
|
temp_wfu = []
|
||||||
|
for r in self.dwI.rpos:
|
||||||
|
if r > 20 :
|
||||||
|
break
|
||||||
|
if rpos_O_filled == False:
|
||||||
|
self.rpos_O.append(r)
|
||||||
|
temp_wfu.append(wfu_O_inter_real(r) + 1j * wfu_O_inter_imag(r))
|
||||||
|
rpos_O_filled = True
|
||||||
|
temp_wfu_L.append(temp_wfu)
|
||||||
|
self.wfu_O.append(temp_wfu_L)
|
||||||
|
|
||||||
|
self.dwO.PrintScatteringMatrix()
|
||||||
|
|
||||||
|
print("====================== Calculating Radial integrals")
|
||||||
|
|
||||||
|
self.radialInt = np.zeros((self.maxL+1, int(2*self.spin_a+1), int(2*self.l+1), int(2*self.spin_b+1)), dtype=complex)
|
||||||
|
|
||||||
|
bs = self.boundState.GetBoundStateWF()
|
||||||
|
|
||||||
|
for L1 in range(0, self.maxL+1):
|
||||||
|
for index1 in range(0, len(wfu_I[L1])):
|
||||||
|
wf1 = wfu_I[L1][index1]
|
||||||
|
for L2 in np.arange(abs(L1 - self.l), L1 + self.l + 1, 1):
|
||||||
|
for index2 in range(0, len(self.wfu_O[int(L2)])):
|
||||||
|
wf2 = self.wfu_O[int(L2)][index2]
|
||||||
|
pf1 = np.exp(1j*self.dwI.CoulombPhaseShift(L1))
|
||||||
|
pf2 = np.exp(1j*self.dwO.CoulombPhaseShift(L2))
|
||||||
|
integral = simpson (bs[:200]*wf1[:200]*wf2[:200], dx=self.boundState.dr)
|
||||||
|
indexL2 = int(L2 - abs(L1-self.l))
|
||||||
|
self.radialInt[L1][index1][indexL2][index2] = integral * pf1 * pf2
|
||||||
|
|
||||||
|
mass_I = self.dwI.mu
|
||||||
|
k_I = self.dwI.k
|
||||||
|
mass_O = self.dwO.mu # reduced mass of outgoing channel
|
||||||
|
k_O = self.dwO.k # wave number of outgoing channel
|
||||||
|
|
||||||
|
D0 = 1.55e+4 # for (d,p)
|
||||||
|
|
||||||
|
self.massFactor = A_B/A_A
|
||||||
|
self.ffactor = np.sqrt(4*np.pi)/k_I /k_O * A_B/A_A
|
||||||
|
|
||||||
|
self.xsecScalingfactor = A_B / A_A *D0 * mass_I * mass_O / np.pi / self.dwI.hbarc**4 / k_I**3 / k_O * (2*self.spin_B + 1) / (2*self.spin_A+1) / (2*self.spin_a +1)
|
||||||
|
|
||||||
|
stop_time = time.time()
|
||||||
|
print(f"Total time {(stop_time - start_time) * 1000:.2f} msec")
|
||||||
|
|
||||||
|
#========== end of contructor
|
||||||
|
|
||||||
|
def FormatSpin(self, spin : float) -> str:
|
||||||
|
if int(2*spin) % 2 == 0 :
|
||||||
|
return f"{int(spin):+d}"
|
||||||
|
else:
|
||||||
|
return f"{int(2*spin):+d}/2"
|
||||||
|
|
||||||
|
|
||||||
|
def PrintRadialIntegral(self):
|
||||||
|
for index1 in range(0, int(2*self.spin_a) + 1):
|
||||||
|
for index2 in range(0, int(2*self.spin_b) + 1):
|
||||||
|
print(f"======================= J1 = L{self.FormatSpin(index1-self.spin_a)}, J2 = L{self.FormatSpin(index2-self.spin_b)}")
|
||||||
|
for L1 in range(0, self.maxL+1):
|
||||||
|
print("{", end="")
|
||||||
|
for L2 in np.arange(abs(L1 - self.l), L1 + self.l + 1):
|
||||||
|
J1 = L1 + index1 - self.spin_a
|
||||||
|
J2 = int(L2) + index2 - self.spin_b
|
||||||
|
indexL2 = int(L2 - abs(L1-self.l))
|
||||||
|
print(f"{{{L1:2d}, {J1:4.1f}, {int(L2):2d}, {J2:4.1f}, {np.real(self.radialInt[L1][index1][indexL2][index2]):12.4e} + {np.imag(self.radialInt[L1][index1][indexL2][index2]):12.4e}I}},", end="")
|
||||||
|
print("},")
|
||||||
|
print("=========================== end of Radial Integrals.")
|
||||||
|
|
||||||
|
def PlotRadialIntegral(self):
|
||||||
|
spin_b = self.spin_b
|
||||||
|
spin_a = self.spin_a
|
||||||
|
l = int(self.l)
|
||||||
|
maxL = self.maxL
|
||||||
|
|
||||||
|
fig, axes = plt.subplots(int(2*spin_b+1)*int(2*l+1), int(2*spin_a+1), figsize=(6*int(2*spin_a+1), 4*int(2*spin_b+1)*int(2*l+1)))
|
||||||
|
|
||||||
|
for index2 in range(0, int(2*spin_b) + 1):
|
||||||
|
for index1 in range(0, int(2*spin_a) + 1):
|
||||||
|
l_list = []
|
||||||
|
for indexL2 in range(0, int(2*l) + 1):
|
||||||
|
haha = []
|
||||||
|
for L1 in range(0, maxL+1):
|
||||||
|
J1 = L1 + index1 - spin_a
|
||||||
|
L2 = int(abs(L1-l) + indexL2)
|
||||||
|
J2 = L2 + index2 - spin_b
|
||||||
|
if J1 < 0 or J2 < 0 :
|
||||||
|
continue
|
||||||
|
l_list.append(L1)
|
||||||
|
haha.append(self.radialInt[L1][index1][indexL2][index2])
|
||||||
|
axes[int(2*l+1)*index2 + indexL2, index1].plot(l_list, np.real(haha), label="Real", marker='o')
|
||||||
|
axes[int(2*l+1)*index2 + indexL2, index1].plot(l_list, np.imag(haha), label="Imag", marker='x')
|
||||||
|
axes[int(2*l+1)*index2 + indexL2, index1].legend()
|
||||||
|
axes[int(2*l+1)*index2 + indexL2, index1].set_xlabel('L1')
|
||||||
|
axes[int(2*l+1)*index2 + indexL2, index1].set_ylabel('Value')
|
||||||
|
axes[int(2*l+1)*index2 + indexL2, index1].set_title(f'Radial Int. vs L for Spin J1 = L{self.FormatSpin(index1-spin_a)}, L2 = L1{indexL2-l:+d}, J2 = L{self.FormatSpin(index2-spin_b)}.')
|
||||||
|
axes[int(2*l+1)*index2 + indexL2, index1].set_xlim(-1, maxL+1)
|
||||||
|
axes[int(2*l+1)*index2 + indexL2, index1].grid()
|
||||||
|
|
||||||
|
plt.tight_layout()
|
||||||
|
plt.show(block=False)
|
||||||
|
input("Press Enter to continue...")
|
||||||
|
|
||||||
|
def PlotScatteringMatrix(self, isIncoming):
|
||||||
|
if isIncoming :
|
||||||
|
self.dwI.PlotScatteringMatrix()
|
||||||
|
else:
|
||||||
|
self.dwO.PlotScatteringMatrix()
|
||||||
|
|
||||||
|
def PlotDistortedWave(self, isIncoming, L, J, maxR = None):
|
||||||
|
if isIncoming:
|
||||||
|
self.dwI.PlotDistortedWave(L, J, maxR)
|
||||||
|
else:
|
||||||
|
plt.plot(self.rpos_O, np.real(self.wfu_O[L][int(J-L + self.spin_b)]), label="Real")
|
||||||
|
plt.plot(self.rpos_O, np.imag(self.wfu_O[L][int(J-L + self.spin_b)]), label="Imag")
|
||||||
|
plt.title(f"Radial wave function for L={L} and J={J}")
|
||||||
|
if maxR != None :
|
||||||
|
plt.xlim(-1, maxR)
|
||||||
|
plt.legend()
|
||||||
|
plt.grid()
|
||||||
|
plt.show(block=False)
|
||||||
|
input("Press Enter to continue...")
|
||||||
|
|
||||||
|
def Gamma(self, L1:int, J1, L2:int, J2, m, ma, mb):
|
||||||
|
if int(L1 + L2 + self.l)%2 != 0: #check if the sum of L1 + L2 + l is even
|
||||||
|
return 0
|
||||||
|
else:
|
||||||
|
fact0 = wigner_9j(S(2*self.j)/2, S(2*self.l)/2, S(2*self.s)/2, S(2*J1)/2, S(2*L1)/2, S(2*self.spin_a)/2, S(2*J2)/2, S(2*L2)/2, S(2*self.spin_b)/2).evalf()
|
||||||
|
if fact0 == 0:
|
||||||
|
return 0
|
||||||
|
else:
|
||||||
|
fact1 = pow(-1, m) * np.power(1j, L1-L2-self.l) * (2*L2+1) * np.sqrt((2*self.l+1)*(2*self.s+1)*(2*L1+1)*(2*J2+1))
|
||||||
|
fact2 = np.sqrt( quantum_factorial(L2-m) / quantum_factorial(L2+m) )
|
||||||
|
fact3 = clebsch_gordan(J2, mb-m, self.j, m-mb+ma, J1, ma)
|
||||||
|
fact4 = clebsch_gordan(L1, 0, self.spin_a, ma, J1, ma)
|
||||||
|
fact5 = clebsch_gordan(L2, -m, self.spin_b, mb, J2, mb-m)
|
||||||
|
fact6 = clebsch_gordan(L1, 0, int(self.l), 0, L2, 0)
|
||||||
|
return fact0 * fact1 * fact2 * fact3 * fact4 * fact5 * fact6
|
||||||
|
|
||||||
|
def Beta(self, m:int, ma, mb):
|
||||||
|
result = 0
|
||||||
|
for L1 in np.arange(0, self.maxL+1):
|
||||||
|
for J1 in np.arange(abs(L1 - self.spin_a), L1 + self.spin_a + 1, 1):
|
||||||
|
for L2 in np.arange(abs(L1 - self.l), L1 + self.l + 1, 1):
|
||||||
|
for J2 in np.arange(abs(L2 - self.spin_b), L2 + self.spin_b + 1, 1):
|
||||||
|
|
||||||
|
if not obeys_triangle_rule(J1, self.j, J2):
|
||||||
|
continue
|
||||||
|
if not(abs(m) <= L2):
|
||||||
|
continue
|
||||||
|
if int(L1 + L2 + self.l) % 2 != 0:
|
||||||
|
continue
|
||||||
|
|
||||||
|
index1 = int(J1 - L1 + self.spin_a)
|
||||||
|
index2 = int(J2 - L2 + self.spin_b)
|
||||||
|
indexL2 = int(L1 - abs(L1 - self.l))
|
||||||
|
|
||||||
|
gg = self.Gamma(L1, J1, L2, J2, m, ma, mb)
|
||||||
|
if gg == 0:
|
||||||
|
continue
|
||||||
|
lp = self.legendrePArray[int(L2)][int(abs(m))]
|
||||||
|
if m < 0 :
|
||||||
|
lp *= (-1)**m * quantum_factorial(int(L2)+m)/ quantum_factorial(int(L2)-m)
|
||||||
|
ri = self.radialInt[int(L1)][index1][indexL2][index2]
|
||||||
|
# print(f"{L1:2d}, {J1:4.1f}({index1:d}), {L2:2d}({indexL2:d}), {J2:4.1f}({index2:d}), {gg:10.6f}, {ri *self.ffactor :.10f}, {lp:10.6f}")
|
||||||
|
|
||||||
|
result += gg * lp * ri
|
||||||
|
|
||||||
|
return result
|
||||||
|
|
||||||
|
def PreCalLegendreP(self, theta_deg:float, maxL:int = None, maxM:int = None):
|
||||||
|
if maxL is None:
|
||||||
|
maxL = self.maxL
|
||||||
|
if maxM is None:
|
||||||
|
maxM = int(self.j + self.spin_b + self.spin_a)
|
||||||
|
self.legendrePArray = associated_legendre_array(maxL, maxM, theta_deg)
|
||||||
|
|
||||||
|
def AngDist(self, theta_deg):
|
||||||
|
xsec = 0
|
||||||
|
self.PreCalLegendreP(theta_deg)
|
||||||
|
for ma in np.arange(-self.spin_a, self.spin_a + 1, 1):
|
||||||
|
for mb in np.arange(-self.spin_b, self.spin_b + 1, 1):
|
||||||
|
for m in np.arange(-self.j + mb - ma, self.j + mb -ma + 1, 1):
|
||||||
|
haha = self.Beta(m, ma, mb)
|
||||||
|
xsec += np.abs(haha)**2
|
||||||
|
|
||||||
|
return xsec * self.xsecScalingfactor * 10 # factor 10 for fm^2 = 10 mb
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
159
Raphael/opticalPotentials.py
Executable file
159
Raphael/opticalPotentials.py
Executable file
|
|
@ -0,0 +1,159 @@
|
||||||
|
#!/usr/bin/env python3
|
||||||
|
|
||||||
|
#need to use import
|
||||||
|
|
||||||
|
import numpy as np
|
||||||
|
|
||||||
|
# Woods-Saxon
|
||||||
|
v = 0
|
||||||
|
r0 = 0
|
||||||
|
a = 0
|
||||||
|
vi = 0
|
||||||
|
ri0 = 0
|
||||||
|
ai = 0
|
||||||
|
# Woods-Saxon Surface
|
||||||
|
vsi = 0
|
||||||
|
rsi0 = 0
|
||||||
|
asi = 0
|
||||||
|
# Spin-orbit
|
||||||
|
vso = 0
|
||||||
|
rso0 = 0
|
||||||
|
aso = 0
|
||||||
|
vsoi = 0
|
||||||
|
rsoi0 = 0
|
||||||
|
asoi = 0
|
||||||
|
# Coulomb
|
||||||
|
rc0 = 0
|
||||||
|
|
||||||
|
def AnCai(A : int, Z : int, E : float):
|
||||||
|
global v, r0, a, vi, ri0, ai, vsi, rsi0, asi, vso, rso0, aso, vsoi, rsoi0, asoi, rc0
|
||||||
|
|
||||||
|
A3 = A**(1./3.)
|
||||||
|
v = 91.85 - 0.249*E + 0.000116*pow(E,2) + 0.642 * Z / A3
|
||||||
|
r0 = 1.152 - 0.00776 / A3
|
||||||
|
a = 0.719 + 0.0126 * A3
|
||||||
|
|
||||||
|
vi = 1.104 + 0.0622 * E
|
||||||
|
ri0 = 1.305 + 0.0997 / A3
|
||||||
|
ai = 0.855 - 0.1 * A3
|
||||||
|
|
||||||
|
vsi = 10.83 - 0.0306 * E
|
||||||
|
rsi0 = 1.334 + 0.152 / A3
|
||||||
|
asi = 0.531 + 0.062 * A3
|
||||||
|
|
||||||
|
vso = 3.557
|
||||||
|
rso0 = 0.972
|
||||||
|
aso = 1.011
|
||||||
|
|
||||||
|
vsoi = 0.0
|
||||||
|
rsoi0 = 0.0
|
||||||
|
asoi = 0.0
|
||||||
|
|
||||||
|
rc0 = 1.303
|
||||||
|
|
||||||
|
def Koning(A : int, Z : int, E : float, Zproj : float):
|
||||||
|
global v, r0, a, vi, ri0, ai, vsi, rsi0, asi, vso, rso0, aso, vsoi, rsoi0, asoi, rc0
|
||||||
|
|
||||||
|
N = A-Z
|
||||||
|
A3 = A**(1./3.)
|
||||||
|
|
||||||
|
vp1 = 59.3 + 21.*(N-Z)/A - 0.024*A
|
||||||
|
vn1 = 59.3 - 21.*(N-Z)/A - 0.024*A
|
||||||
|
|
||||||
|
vp2 = 0.007067 + 0.00000423*A
|
||||||
|
vn2 = 0.007228 - 0.00000148*A
|
||||||
|
|
||||||
|
vp3 = 0.00001729 + 0.00000001136 * A
|
||||||
|
vn3 = 0.00001994 - 0.00000002 * A
|
||||||
|
|
||||||
|
vp4 = 7e-9 # = vn4
|
||||||
|
vn4 = vp4
|
||||||
|
|
||||||
|
wp1 = 14.667 + 0.009629*A
|
||||||
|
wn1 = 12.195 + 0.0167*A
|
||||||
|
|
||||||
|
wp2 = 73.55 + 0.0795*A # = wn2
|
||||||
|
wn2 = wp2
|
||||||
|
|
||||||
|
dp1 = 16 + 16.*(N-Z)/A
|
||||||
|
dn1 = 16 - 16.*(N-Z)/A
|
||||||
|
|
||||||
|
dp2 = 0.018 + 0.003802/(1 + np.exp((A-156.)/8)) # = dn2
|
||||||
|
dn2 = dp2
|
||||||
|
|
||||||
|
dp3 = 11.5 # = dn3
|
||||||
|
dn3 = dp3
|
||||||
|
|
||||||
|
vso1 = 5.922 + 0.003 * A
|
||||||
|
vso2 = 0.004
|
||||||
|
|
||||||
|
wso1 = -3.1
|
||||||
|
wso2 = 160
|
||||||
|
|
||||||
|
epf = -8.4075 + 0.01378 *A
|
||||||
|
enf = -11.2814 + 0.02646 *A
|
||||||
|
|
||||||
|
rc = 1.198 + 0.697/pow(A3,2) + 12.995/pow(A3,5)
|
||||||
|
vc = 1.73/rc * Z / A3
|
||||||
|
|
||||||
|
v = vp1*(1 - vp2*(E-epf) + vp3*pow(E-epf,2) - vp4*pow(E-epf,3)) + vc * vp1 * (vp2 - 2*vp3*(E-epf) + 3*vp4*pow(E-epf,2))
|
||||||
|
#neutron
|
||||||
|
if Zproj == 0 :
|
||||||
|
v = vn1*(1 - vn2*(E-enf) + vn3*pow(E-enf,2) - vn4*pow(E-enf,3))
|
||||||
|
|
||||||
|
r0 = 1.3039 - 0.4054 / A3
|
||||||
|
a = 0.6778 - 0.000148 * A
|
||||||
|
|
||||||
|
vi = wp1 * pow(E-epf,2)/(pow(E-epf,2) + pow(wp2,2))
|
||||||
|
if Zproj == 0 :
|
||||||
|
vi = wn1 * pow(E-enf,2)/(pow(E-enf,2) + pow(wn2,2))
|
||||||
|
|
||||||
|
ri0 = 1.3039 - 0.4054 / A3
|
||||||
|
ai = 0.6778 - 0.000148 * A
|
||||||
|
|
||||||
|
vsi = dp1 * pow(E-epf,2)/(pow(E-epf,2)+pow(dp3,2)) * np.exp(-dp2*(E-epf))
|
||||||
|
if Zproj == 0 :
|
||||||
|
vsi = dn1 * pow(E-enf,2)/(pow(E-enf,2)+pow(dn3,2)) * np.exp(-dn2*(E-enf))
|
||||||
|
|
||||||
|
rsi0 = 1.3424 - 0.01585 * A3
|
||||||
|
asi = 0.5187 + 0.0005205 * A
|
||||||
|
if Zproj == 0:
|
||||||
|
asi = 0.5446 - 0.0001656 * A
|
||||||
|
|
||||||
|
vso = vso1 * np.exp(-vso2 * (E-epf))
|
||||||
|
if Zproj == 0:
|
||||||
|
vso = vso1 * np.exp(-vso2 * (E-enf))
|
||||||
|
|
||||||
|
rso0 = 1.1854 - 0.647/A3
|
||||||
|
aso = 0.59
|
||||||
|
|
||||||
|
vsoi = wso1 * pow(E-epf,2)/(pow(E-epf,2)+pow(wso2,2))
|
||||||
|
if Zproj == 0 :
|
||||||
|
vsoi = wso1 * pow(E-enf,2)/(pow(E-enf,2)+pow(wso2,2))
|
||||||
|
|
||||||
|
rsoi0 = 1.1854 - 0.647/A3
|
||||||
|
asoi = 0.59
|
||||||
|
|
||||||
|
rc0 = rc
|
||||||
|
|
||||||
|
def ConvertLSym(LSym :str) -> int:
|
||||||
|
if LSym == "s" :
|
||||||
|
return 0
|
||||||
|
elif LSym == "p" :
|
||||||
|
return 1
|
||||||
|
elif LSym == "d" :
|
||||||
|
return 2
|
||||||
|
elif LSym == "f" :
|
||||||
|
return 3
|
||||||
|
elif LSym == "g" :
|
||||||
|
return 4
|
||||||
|
elif LSym == "h" :
|
||||||
|
return 5
|
||||||
|
elif LSym == "i" :
|
||||||
|
return 6
|
||||||
|
elif LSym == "j" :
|
||||||
|
return 7
|
||||||
|
elif LSym == "k" :
|
||||||
|
return 8
|
||||||
|
else :
|
||||||
|
return -1
|
||||||
|
|
@ -197,7 +197,7 @@ class SolvingSE:
|
||||||
self.eta = 0
|
self.eta = 0
|
||||||
else:
|
else:
|
||||||
self.eta = self.Z * self.ee * self.k /2 /self.Ecm
|
self.eta = self.Z * self.ee * self.k /2 /self.Ecm
|
||||||
self.maxL = int(self.k * (1.4 * (self.A_A**(1/3) + self.A_a**(1/3)) + 5))
|
self.maxL = int(self.k * (1.4 * (self.A_A**(1/3) + self.A_a**(1/3)) + 10))
|
||||||
|
|
||||||
# def SetA_ExSpin(self, ExA, sA ):
|
# def SetA_ExSpin(self, ExA, sA ):
|
||||||
# self.ExA = ExA
|
# self.ExA = ExA
|
||||||
|
|
|
||||||
Loading…
Reference in New Issue
Block a user