change opticalpotential parameter to be class, add alpha

This commit is contained in:
Ryan@Home 2025-03-13 19:16:35 -04:00
parent 0af018f0d0
commit 623fc99c27
4 changed files with 209 additions and 165 deletions

View File

@ -2,7 +2,8 @@
import time
import matplotlib.pyplot as plt
from distortedWave import WoodsSaxonPot, CoulombPotential, SpinOrbit_Pot, WS_SurfacePot, DistortedWave
from solveSE import WoodsSaxonPot, CoulombPotential, SpinOrbit_Pot, WS_SurfacePot
from distortedWave import DistortedWave
from dwba_zr import DWBA_ZR
# haha = DWBA_ZR("16O", "d", "p", "17O", "1/2+", "1s1/2", 0.0, 10)

View File

@ -66,29 +66,29 @@ class DWBA_ZR:
print("====================== Incoming wave function ")
if A_a == 2 and Z_a == 1:
op.AnCai(A_A, Z_A, self.ELab)
pot = op.AnCai(A_A, Z_A, self.ELab)
if A_a == 1 and Z_a == 1:
op.Koning(A_A, Z_A, self.ELab , Z_a)
pot = op.Koning(A_A, Z_A, self.ELab , Z_a)
self.dwI = self.reactDigest.dwI
self.dwI.spin_A = self.spin_A
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.AddPotential(WoodsSaxonPot( -pot.v, pot.r0, pot.a), False)
self.dwI.AddPotential(WoodsSaxonPot(-1j*pot.vi, pot.ri0, pot.ai), False)
self.dwI.AddPotential(WS_SurfacePot(-1j*pot.vsi, pot.rsi0, pot.asi), False)
self.dwI.AddPotential(SpinOrbit_Pot( -pot.vso, pot.rso0, pot.aso), False)
self.dwI.AddPotential(SpinOrbit_Pot(-1j*pot.vsoi, pot.rsoi0, pot.asoi), False)
self.dwI.AddPotential(CoulombPotential( pot.rc0), False)
self.dwI.PrintPotentials()
self.maxL1 = self.dwI.maxL
print("====================== Outgoing wave function ")
if A_b == 1 and Z_b == 1:
op.Koning(A_B, Z_B, self.Eout, Z_b)
pot = op.Koning(A_B, Z_B, self.Eout, Z_b)
if A_b == 2 and Z_b == 1:
op.AnCai(A_B, Z_B, self.Eout)
pot = op.AnCai(A_B, Z_B, self.Eout)
self.maxL2 = self.maxL1 + self.l
@ -97,12 +97,12 @@ class DWBA_ZR:
self.dwO.maxL = self.maxL2
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.AddPotential(WoodsSaxonPot( -pot.v, pot.r0, pot.a), False)
self.dwO.AddPotential(WoodsSaxonPot(-1j*pot.vi, pot.ri0, pot.ai), False)
self.dwO.AddPotential(WS_SurfacePot(-1j*pot.vsi, pot.rsi0, pot.asi), False)
self.dwO.AddPotential(SpinOrbit_Pot( -pot.vso, pot.rso0, pot.aso), False)
self.dwO.AddPotential(SpinOrbit_Pot(-1j*pot.vsoi, pot.rsoi0, pot.asoi), False)
self.dwO.AddPotential(CoulombPotential( pot.rc0), False)
self.dwO.PrintPotentials()

View File

@ -4,137 +4,180 @@
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
class OpticalPotential:
# 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
def Print(self):
print(f"Wood-Saxon (Re) : {self.v:7.3f}, {self.r0:5.3f}, {self.a:5.3f}")
print(f" (Im) : {self.vi:7.3f}, {self.ri0:5.3f}, {self.ai:5.3f}")
print(f"WS Surface (Im) : {self.vsi:7.3f}, {self.rsi0:5.3f}, {self.asi:5.3f}")
print(f"Spin-Orbit (Re) : {self.vso:7.3f}, {self.rso0:5.3f}, {self.aso:5.3f}")
print(f" (Im) : {self.vsoi:7.3f}, {self.rsoi0:5.3f}, {self.asoi:5.3f}")
print(f"Coulomb (Re) : {'':7s}, {self.rc0:5.3f}")
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
class SuAndHan(OpticalPotential):
def __init__(self, A : int, Z :int , E : float):
N = A - Z
A3 = A**(1./3.)
vsi = 10.83 - 0.0306 * E
rsi0 = 1.334 + 0.152 / A3
asi = 0.531 + 0.062 * A3
vsiCOND = 27.5816 - 0.0797 * E + 48.0*(N-Z)/A
viCOND = -4.0174 + 0.1409 * E
self.v = 175.0881 - 0.6236 * E + 0.0006*E*E + 30.*(N-Z)/A - 0.236 * Z/A3
self.r0 = 1.3421
self.a = 0.6578
vso = 3.557
rso0 = 0.972
aso = 1.011
self.vi = viCOND
if viCOND < 0 :
self.vi = 0.0
self.ri0 = 1.4259
self.ai = 0.6578
vsoi = 0.0
rsoi0 = 0.0
asoi = 0.0
self.vsi = vsiCOND
if vsiCOND < 0 :
self.vsi = 0.0
self.rsi0 = 1.2928
self.asi = 0.6359
rc0 = 1.303
self.vso = 0.0
self.rso0 = 1.2686
self.aso = 0.85
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
self.vsoi = 0.0
self.rsoi0 = 0.0
self.asoi = 0.0
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
self.rc0 = 1.350
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
class AnCai(OpticalPotential):
def __init__(self, A : int, Z :int , E : float):
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
A3 = A**(1./3.)
self.v = 91.85 - 0.249*E + 0.000116*pow(E,2) + 0.642 * Z / A3
self.r0 = 1.152 - 0.00776 / A3
self.a = 0.719 + 0.0126 * A3
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))
self.vi = 1.104 + 0.0622 * E
self.ri0 = 1.305 + 0.0997 / A3
self.ai = 0.855 - 0.1 * A3
rsi0 = 1.3424 - 0.01585 * A3
asi = 0.5187 + 0.0005205 * A
if Zproj == 0:
asi = 0.5446 - 0.0001656 * A
self.vsi = 10.83 - 0.0306 * E
self.rsi0 = 1.334 + 0.152 / A3
self.asi = 0.531 + 0.062 * A3
vso = vso1 * np.exp(-vso2 * (E-epf))
if Zproj == 0:
vso = vso1 * np.exp(-vso2 * (E-enf))
self.vso = 3.557
self.rso0 = 0.972
self.aso = 1.011
rso0 = 1.1854 - 0.647/A3
aso = 0.59
self.vsoi = 0.0
self.rsoi0 = 0.0
self.asoi = 0.0
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))
self.rc0 = 1.303
rsoi0 = 1.1854 - 0.647/A3
asoi = 0.59
class Koning(OpticalPotential):
def __init__(self, A : int, Z :int , E : float, Zproj : float):
rc0 = rc
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
self.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 :
self.v = vn1*(1 - vn2*(E-enf) + vn3*pow(E-enf,2) - vn4*pow(E-enf,3))
self.r0 = 1.3039 - 0.4054 / A3
self.a = 0.6778 - 0.000148 * A
self.vi = wp1 * pow(E-epf,2)/(pow(E-epf,2) + pow(wp2,2))
if Zproj == 0 :
self.vi = wn1 * pow(E-enf,2)/(pow(E-enf,2) + pow(wn2,2))
self.ri0 = 1.3039 - 0.4054 / A3
self.ai = 0.6778 - 0.000148 * A
self.vsi = dp1 * pow(E-epf,2)/(pow(E-epf,2)+pow(dp3,2)) * np.exp(-dp2*(E-epf))
if Zproj == 0 :
self.vsi = dn1 * pow(E-enf,2)/(pow(E-enf,2)+pow(dn3,2)) * np.exp(-dn2*(E-enf))
self.rsi0 = 1.3424 - 0.01585 * A3
self.asi = 0.5187 + 0.0005205 * A
if Zproj == 0:
self.asi = 0.5446 - 0.0001656 * A
self.vso = vso1 * np.exp(-vso2 * (E-epf))
if Zproj == 0:
self.vso = vso1 * np.exp(-vso2 * (E-enf))
self.rso0 = 1.1854 - 0.647/A3
self.aso = 0.59
self.vsoi = wso1 * pow(E-epf,2)/(pow(E-epf,2)+pow(wso2,2))
if Zproj == 0 :
self.vsoi = wso1 * pow(E-enf,2)/(pow(E-enf,2)+pow(wso2,2))
self.rsoi0 = 1.1854 - 0.647/A3
self.asoi = 0.59
def ConvertLSym(LSym :str) -> int:
if LSym == "s" :

View File

@ -85,87 +85,87 @@ with open(fileOutName, "w") as file:
file.write(f"{0.1:+08.4f}{15:+08.4f}\n")
#===== Block 5
if A_a == 2 :
op.AnCai(A_A, Z_A, A_a*ELab)
pot = op.AnCai(A_A, Z_A, A_a*ELab)
else:
op.Koning(A_A, Z_A, A_a*ELab, Z_a)
pot = op.Koning(A_A, Z_A, A_a*ELab, Z_a)
file.write(f"{A_a*ELab:+08.4f}")
file.write(f"{A_a:+08.4f}")
file.write(f"{Z_a:+08.4f}")
file.write(f"{A_A:+08.4f}")
file.write(f"{Z_A:+08.4f}")
file.write(f"{op.rc0:+08.4f}")
file.write(f"{pot.rc0:+08.4f}")
file.write(f"{"":8s}")
file.write(f"{"":8s}")
file.write(f"{2*spin_a:+08.4f}\n")
# Woods-Saxon
file.write(f"{1:+08.4f}")
file.write(f"{-op.v:+08.4f}") # real
file.write(f"{op.r0:+08.4f}") #
file.write(f"{op.a:+08.4f}") #
file.write(f"{-pot.v:+08.4f}") # real
file.write(f"{pot.r0:+08.4f}") #
file.write(f"{pot.a:+08.4f}") #
file.write(f"{"":8s}") # spin-orbit skipped
file.write(f"{-op.vi:+08.4f}") # imag
file.write(f"{op.ri0:+08.4f}") #
file.write(f"{op.ai:+08.4f}\n") #
file.write(f"{-pot.vi:+08.4f}") # imag
file.write(f"{pot.ri0:+08.4f}") #
file.write(f"{pot.ai:+08.4f}\n") #
# Woods-Saxon surface
file.write(f"{2:+08.4f}")
file.write(f"{"":8s}") # real
file.write(f"{"":8s}") #
file.write(f"{"":8s}") #
file.write(f"{"":8s}") # spin-orbit skipped
file.write(f"{4*op.vsi:+08.4f}") # imag
file.write(f"{op.rsi0:+08.4f}") #
file.write(f"{op.asi:+08.4f}\n") #
file.write(f"{4*pot.vsi:+08.4f}") # imag
file.write(f"{pot.rsi0:+08.4f}") #
file.write(f"{pot.asi:+08.4f}\n") #
# Spin-Orbit
file.write(f"{-4:+08.4f}")
file.write(f"{-4*op.vso:+08.4f}") # real
file.write(f"{op.rso0:+08.4f}") #
file.write(f"{op.aso:+08.4f}") #
file.write(f"{-4*pot.vso:+08.4f}") # real
file.write(f"{pot.rso0:+08.4f}") #
file.write(f"{pot.aso:+08.4f}") #
file.write(f"{"":8s}") # spin-orbit skipped
file.write(f"{-4*op.vsoi:+08.4f}") # imag
file.write(f"{op.rsoi0:+08.4f}") #
file.write(f"{op.asoi:+08.4f}\n") #
file.write(f"{-4*pot.vsoi:+08.4f}") # imag
file.write(f"{pot.rsoi0:+08.4f}") #
file.write(f"{pot.asoi:+08.4f}\n") #
#===== Block 6
if A_a == 2 :
op.Koning(A_B, Z_B, A_a*ELab + Q_value - Ex, Z_b)
pot = op.Koning(A_B, Z_B, A_a*ELab + Q_value - Ex, Z_b)
else:
op.AnCai(A_B, Z_B, A_a*ELab + Q_value - Ex)
pot = op.AnCai(A_B, Z_B, A_a*ELab + Q_value - Ex)
file.write(f"{Q_value:+08.4f}")
file.write(f"{A_b:+08.4f}")
file.write(f"{Z_b:+08.4f}")
file.write(f"{A_B:+08.4f}")
file.write(f"{Z_B:+08.4f}")
file.write(f"{op.rc0:+08.4f}")
file.write(f"{pot.rc0:+08.4f}")
file.write(f"{"":8s}")
file.write(f"{"":8s}")
file.write(f"{2*spin_b:+08.4f}\n")
# Woods-Saxon
file.write(f"{1:+08.4f}")
file.write(f"{-op.v:+08.4f}") # real
file.write(f"{op.r0:+08.4f}") #
file.write(f"{op.a:+08.4f}") #
file.write(f"{-pot.v:+08.4f}") # real
file.write(f"{pot.r0:+08.4f}") #
file.write(f"{pot.a:+08.4f}") #
file.write(f"{"":8s}") # spin-orbit skipped
file.write(f"{-op.vi:+08.4f}") # imag
file.write(f"{op.ri0:+08.4f}") #
file.write(f"{op.ai:+08.4f}\n") #
file.write(f"{-pot.vi:+08.4f}") # imag
file.write(f"{pot.ri0:+08.4f}") #
file.write(f"{pot.ai:+08.4f}\n") #
# Woods-Saxon surface
file.write(f"{2:+08.4f}")
file.write(f"{"":8s}") # real
file.write(f"{"":8s}") #
file.write(f"{"":8s}") #
file.write(f"{"":8s}") # spin-orbit skipped
file.write(f"{4*op.vsi:+08.4f}") # imag
file.write(f"{op.rsi0:+08.4f}") #
file.write(f"{op.asi:+08.4f}\n") #
file.write(f"{4*pot.vsi:+08.4f}") # imag
file.write(f"{pot.rsi0:+08.4f}") #
file.write(f"{pot.asi:+08.4f}\n") #
# Spin-Orbit
file.write(f"{-4:+08.4f}")
file.write(f"{-4*op.vso:+08.4f}") # real
file.write(f"{op.rso0:+08.4f}") #
file.write(f"{op.aso:+08.4f}") #
file.write(f"{-4*pot.vso:+08.4f}") # real
file.write(f"{pot.rso0:+08.4f}") #
file.write(f"{pot.aso:+08.4f}") #
file.write(f"{"":8s}") # spin-orbit skipped
file.write(f"{-4*op.vsoi:+08.4f}") # imag
file.write(f"{op.rsoi0:+08.4f}") #
file.write(f"{op.asoi:+08.4f}\n") #
file.write(f"{-4*pot.vsoi:+08.4f}") # imag
file.write(f"{pot.rsoi0:+08.4f}") #
file.write(f"{pot.asoi:+08.4f}\n") #
#====== bound state
file.write(f"{BindingEnergy:+08.4f}")
file.write(f"{A_x:+08.4f}")