From c911e604eb7b66ee0616a693bcf97472fc936410 Mon Sep 17 00:00:00 2001 From: "Ryan@SOLARIS_testStation" Date: Tue, 25 Feb 2025 14:35:38 -0500 Subject: [PATCH] seperate reactionData, add symbolic link of IAEANuclearData --- Raphael/IAEANuclearData.py | 1 + Raphael/README.md | 1 + Raphael/distortedWave.py | 4 +- Raphael/dwba_zr.py | 124 +----------------------------------- Raphael/reactionData.py | 126 +++++++++++++++++++++++++++++++++++++ Raphael/solveSE.py | 2 +- 6 files changed, 133 insertions(+), 125 deletions(-) create mode 120000 Raphael/IAEANuclearData.py create mode 100755 Raphael/reactionData.py diff --git a/Raphael/IAEANuclearData.py b/Raphael/IAEANuclearData.py new file mode 120000 index 0000000..33810fc --- /dev/null +++ b/Raphael/IAEANuclearData.py @@ -0,0 +1 @@ +../Cleopatra/IAEANuclearData.py \ No newline at end of file diff --git a/Raphael/README.md b/Raphael/README.md index ab23b9c..09ad857 100644 --- a/Raphael/README.md +++ b/Raphael/README.md @@ -23,6 +23,7 @@ The foundation of the code base are * assLegendreP.py - for associate Legendre polynomial for positive m * clebschGordan.py - for custom build CG, which is faster * opticalPotential.py - for optical potential, only have An & Cai for deuteron and Kronning for proton now +* reactionData.py - to digest and store the basic reaction data * ../Cleopatra/IAEANuclearData.py - for getting nuclear data like mass and spin-partiy * coulombWave.py - attemp to make fast CoulombWave.... diff --git a/Raphael/distortedWave.py b/Raphael/distortedWave.py index fdb57a8..a411dca 100755 --- a/Raphael/distortedWave.py +++ b/Raphael/distortedWave.py @@ -1,8 +1,6 @@ #!/usr/bin/env python3 -from boundState import BoundState - -from solveSE import WoodsSaxonPot, CoulombPotential, SpinOrbit_Pot, WS_SurfacePot, SolvingSE +from solveSE import SolvingSE from mpmath import coulombf, coulombg import numpy as np diff --git a/Raphael/dwba_zr.py b/Raphael/dwba_zr.py index 13c3bfe..8a72135 100755 --- a/Raphael/dwba_zr.py +++ b/Raphael/dwba_zr.py @@ -9,8 +9,8 @@ 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 +# 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 @@ -18,126 +18,8 @@ from boundState import BoundState from solveSE import WoodsSaxonPot, CoulombPotential, SpinOrbit_Pot, WS_SurfacePot from distortedWave import DistortedWave import opticalPotentials as op - +from reactionData import approximate_to_half_integer, ReactionData -def approximate_to_half_integer(value): - return round(value * 2) / 2 - -class ReactionData: - def __init__(self, nu_A:str, nu_a:str, nu_b:str, nu_B:str, JB:str, orbital:str, ExB:float, ELabPerU:float): - self.SpinBalanced = self.ReactionDigestion(nu_A, nu_a, nu_b, nu_B, JB, orbital, ExB, ELabPerU) - - - def ReactionDigestion(self, nu_A:str, nu_a:str, nu_b:str, nu_B:str, JB:str, orbital:str, ExB:float, ELabPerU:float): - iso = IsotopeClass() - - self.A_A, self.Z_A = iso.GetAZ(nu_A) - self.A_a, self.Z_a = iso.GetAZ(nu_a) - self.A_b, self.Z_b = iso.GetAZ(nu_b) - self.A_B, self.Z_B = iso.GetAZ(nu_B) - - self.ELab = self.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) - ExB = ExB - - # sym_A = iso.GetSymbol(A_A, Z_A) - # sym_B = iso.GetSymbol(A_B, Z_B) - - spin_A_str = iso.GetJpi(self.A_A, self.Z_A) - self.spin_A = float(eval(re.sub(r'[+-]', '', spin_A_str))) - self.spin_B = float(eval(re.sub(r'[+-]', '', JB))) - - if self.A_a == 2 and self.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 - self.A_x = abs(self.A_a - self.A_b) - self.Z_x = abs(self.Z_a - self.Z_b) - - mass_x = iso.GetMassFromAZ(self.A_x, self.Z_x) - - #======== core - if self.A_A < self.A_B : # (d,p) - self.A_c = self.A_A - self.Z_c = self.Z_A - self.BindingEnergy = mass_B - mass_A - mass_x + ExB - else: #(p,d) - self.A_c = self.A_B - self.Z_c = self.Z_B - self.BindingEnergy = mass_A - mass_B - 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 - - self.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) - - self.j = approximate_to_half_integer(self.j) - self.s = approximate_to_half_integer(self.s) - self.spin_a = approximate_to_half_integer(self.spin_a) - self.spin_b = approximate_to_half_integer(self.spin_b) - - 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 False - 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" - - self.Q_value = mass_A + mass_a - mass_b - mass_B - ExB - self.dwI = DistortedWave(nu_A, nu_a, self.ELab) - - Ecm_I = self.dwI.Ecm - Ecm_O = Ecm_I + self.Q_value - self.Eout = ((Ecm_O + mass_b + mass_B + ExB)**2 - (mass_b + mass_B + ExB)**2)/2/(mass_B + ExB) - self.dwO = DistortedWave(nu_B, nu_b, self.Eout) - - Eout2 = self.ELab + self.Q_value #this is incorrec, but used in ptolmey infileCreator - - print("==================================================") - print(self.reactionStr) - print(f"Transfer Orbtial : {orbital}") - print(f"Q-value : {self.Q_value:10.6f} MeV") - print(f"Binding : {self.BindingEnergy:10.6f} MeV") - print(f" Eout : {self.Eout} MeV | {Eout2}") - - return True - - 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): diff --git a/Raphael/reactionData.py b/Raphael/reactionData.py new file mode 100755 index 0000000..4ed4744 --- /dev/null +++ b/Raphael/reactionData.py @@ -0,0 +1,126 @@ +#!/usr/bin/env python3 +import re +# sys.path.append(os.path.join(os.path.dirname(__file__), '../Cleopatra')) +from IAEANuclearData import IsotopeClass + +from clebschGordan import obeys_triangle_rule +from distortedWave import DistortedWave +import opticalPotentials as op + + +def approximate_to_half_integer(value): + return round(value * 2) / 2 + +class ReactionData: + def __init__(self, nu_A:str, nu_a:str, nu_b:str, nu_B:str, JB:str, orbital:str, ExB:float, ELabPerU:float): + self.SpinBalanced = self.ReactionDigestion(nu_A, nu_a, nu_b, nu_B, JB, orbital, ExB, ELabPerU) + + + def ReactionDigestion(self, nu_A:str, nu_a:str, nu_b:str, nu_B:str, JB:str, orbital:str, ExB:float, ELabPerU:float): + iso = IsotopeClass() + + self.A_A, self.Z_A = iso.GetAZ(nu_A) + self.A_a, self.Z_a = iso.GetAZ(nu_a) + self.A_b, self.Z_b = iso.GetAZ(nu_b) + self.A_B, self.Z_B = iso.GetAZ(nu_B) + + self.ELab = self.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) + ExB = ExB + + # sym_A = iso.GetSymbol(A_A, Z_A) + # sym_B = iso.GetSymbol(A_B, Z_B) + + spin_A_str = iso.GetJpi(self.A_A, self.Z_A) + self.spin_A = float(eval(re.sub(r'[+-]', '', spin_A_str))) + self.spin_B = float(eval(re.sub(r'[+-]', '', JB))) + + if self.A_a == 2 and self.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 + self.A_x = abs(self.A_a - self.A_b) + self.Z_x = abs(self.Z_a - self.Z_b) + + mass_x = iso.GetMassFromAZ(self.A_x, self.Z_x) + + #======== core + if self.A_A < self.A_B : # (d,p) + self.A_c = self.A_A + self.Z_c = self.Z_A + self.BindingEnergy = mass_B - mass_A - mass_x + ExB + else: #(p,d) + self.A_c = self.A_B + self.Z_c = self.Z_B + self.BindingEnergy = mass_A - mass_B - 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 + + self.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) + + self.j = approximate_to_half_integer(self.j) + self.s = approximate_to_half_integer(self.s) + self.spin_a = approximate_to_half_integer(self.spin_a) + self.spin_b = approximate_to_half_integer(self.spin_b) + + 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 False + 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" + + self.Q_value = mass_A + mass_a - mass_b - mass_B - ExB + self.dwI = DistortedWave(nu_A, nu_a, self.ELab) + + Ecm_I = self.dwI.Ecm + Ecm_O = Ecm_I + self.Q_value + self.Eout = ((Ecm_O + mass_b + mass_B + ExB)**2 - (mass_b + mass_B + ExB)**2)/2/(mass_B + ExB) + self.dwO = DistortedWave(nu_B, nu_b, self.Eout) + + Eout2 = self.ELab + self.Q_value #this is incorrec, but used in ptolmey infileCreator + + print("==================================================") + print(self.reactionStr) + print(f"Transfer Orbtial : {orbital}") + print(f"Q-value : {self.Q_value:10.6f} MeV") + print(f"Binding : {self.BindingEnergy:10.6f} MeV") + print(f" Eout : {self.Eout} MeV | {Eout2}") + + return True diff --git a/Raphael/solveSE.py b/Raphael/solveSE.py index b8726f1..315e262 100755 --- a/Raphael/solveSE.py +++ b/Raphael/solveSE.py @@ -6,7 +6,7 @@ import re import sys, os import matplotlib.pyplot as plt -sys.path.append(os.path.join(os.path.dirname(__file__), '../Cleopatra')) +# sys.path.append(os.path.join(os.path.dirname(__file__), '../Cleopatra')) from IAEANuclearData import IsotopeClass class PotentialForm: