From b1433f3d983aa4d43645235cf4bd2fa67ab83377 Mon Sep 17 00:00:00 2001 From: "Ryan@SOLARIS_testStation" Date: Tue, 25 Feb 2025 18:19:09 -0500 Subject: [PATCH] snapshot. xsec scaling factor problem... --- Raphael/IAEANuclearData.py | 1 - Raphael/dwba_zr.py | 51 ++++-- Raphael/reactionData.py | 36 ++-- Raphael/solveSE.py | 2 +- dwuck4/extractData.py | 15 +- dwuck4/inFileCreatorDW.py | 331 ++++++++----------------------------- 6 files changed, 140 insertions(+), 296 deletions(-) delete mode 120000 Raphael/IAEANuclearData.py diff --git a/Raphael/IAEANuclearData.py b/Raphael/IAEANuclearData.py deleted file mode 120000 index 33810fc..0000000 --- a/Raphael/IAEANuclearData.py +++ /dev/null @@ -1 +0,0 @@ -../Cleopatra/IAEANuclearData.py \ No newline at end of file diff --git a/Raphael/dwba_zr.py b/Raphael/dwba_zr.py index 8a72135..1c21146 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 @@ -21,9 +21,9 @@ import opticalPotentials as op from reactionData import approximate_to_half_integer, ReactionData 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): + def __init__(self, nu_A:str, nu_a:str, nu_b:str, JB:str, orbital:str, ExB:float, ELabPerU:float): - self.reactDigest = ReactionData(nu_A, nu_a, nu_b, nu_B, JB, orbital, ExB, ELabPerU) + self.reactDigest = ReactionData(nu_A, nu_a, nu_b, JB, orbital, ExB, ELabPerU) if self.reactDigest.SpinBalanced == False : return @@ -106,20 +106,37 @@ class DWBA_ZR: self.dwO.PrintPotentials() #---------------------------------------- other constants + print("========================================") D0 = 1.55e+4 # for (d,p) mass_I = self.dwI.mu - mass_O = self.dwO.mu k_I = self.dwI.k + mass_O = self.dwO.mu k_O = self.dwO.k # wave number of outgoing channel + # print(f" mu(I) : {mass_I}") + # print(f" k(I) : {k_I}") + # print(f" mu(O) : {mass_O}") + # print(f" k(O) : {k_O}") + self.massBoverMassA = A_B/A_A self.ffactor = np.sqrt(4*np.pi)/k_I /k_O - self.xsecScalingfactor = 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) + # print(f"spin A : {self.spin_A}") + # print(f"spin a : {self.spin_a}") + # print(f"spin B : {self.spin_B}") + + # self.spinFactor = (2*self.spin_B + 1) / (2*self.spin_A+1) / (2*self.s +1) + self.spinFactor = (2*self.spin_B + 1) / (2*self.spin_A+1) / (2*self.spin_a +1) + + # print(f" spin factor : {self.spinFactor}") + + self.xsecScalingfactor = D0 * mass_I * mass_O / np.pi / self.dwI.hbarc**4 / k_I**3 / k_O * self.spinFactor self.radialInt = None + print(f"Xsec Scaling factor : {self.xsecScalingfactor:.6f}") + self.PreCalNineJ() self.PreCalClebschGordan() @@ -132,7 +149,7 @@ class DWBA_ZR: return f"{int(2*spin):+d}/2" def FindBoundState(self): - self.boundState.FindPotentialDepth(-70, -45, 0.5) + self.boundState.FindPotentialDepth(-80, -45, 0.5) def ConvertLJ2RadialIndex(self, L1:int, J1, L2:int, J2): index1 = int(J1 - L1 + self.spin_a) @@ -200,6 +217,7 @@ class DWBA_ZR: pf2 = np.exp(1j*self.dwO.CoulombPhaseShift(L2)) integral = simpson (bs*wf1*wf2, dx=self.boundState.dr) indexL2 = int(L2 - L1 + self.l) + # product = integral * pf1 * pf2 product = integral * pf1 * pf2 * self.massBoverMassA self.radialInt[L1][index1][indexL2][index2] = product # if J1 == L1 + self.spin_a and L2 == L1 + 1 and J2 == L2 - self.spin_b: @@ -279,10 +297,10 @@ class DWBA_ZR: plt.show(block=False) input("Press Enter to continue...") - def PlotScatteringMatrix(self, isIncoming): - if isIncoming : + def PlotIncomingScatteringMatrix(self): self.dwI.PlotScatteringMatrix() - else: + + def PlotOutgoingScatteringMatrix(self): self.dwO.PlotScatteringMatrix() def PlotIncomingDistortedWave(self, L, J, maxR = None): @@ -344,8 +362,8 @@ class DWBA_ZR: stop_time = time.time() print(f"Total time for pre-cal all CG {(stop_time - start_time) * 1000:.2f} msec") - print(f"self.maxL1, maxJ1, maxJ2, maxJ3") - print(self.CG.shape) + print(f"max(L1 J1, L2, J2) = {self.maxL1}, {maxJ1}, {maxJ2}, {maxJ3}") + print("CG shape : ",self.CG.shape) def GetPreCalCG(self, j1, m1, j2, m2, j3, m3): return self.CG[int(2*j1), int(2*m1 + 2*self.maxJ1+1), @@ -462,9 +480,14 @@ class DWBA_ZR: stop_time = time.time() print(f"\nTotal time {(stop_time - start_time) :.2f} sec") - def PrintAngDist(self): + def PrintAngDist(self, step:int = 1): + count = 0 for th, xs in zip(self.angList, self.angDist): - print(f"{th:6.1f}, {xs:13.10f}") + if step > 1 and count % step != 0: + count += 1 + continue + print(f"{{{th:6.1f}, {xs:13.10f}}},") + count += 1 def PlotAngDist(self, angMin = None, angMax = None): plt.plot(self.angList, self.angDist) diff --git a/Raphael/reactionData.py b/Raphael/reactionData.py index 4ed4744..17e5f02 100755 --- a/Raphael/reactionData.py +++ b/Raphael/reactionData.py @@ -1,6 +1,7 @@ #!/usr/bin/env python3 import re -# sys.path.append(os.path.join(os.path.dirname(__file__), '../Cleopatra')) +import sys, os +sys.path.append(os.path.join(os.path.dirname(__file__), '../Cleopatra')) from IAEANuclearData import IsotopeClass from clebschGordan import obeys_triangle_rule @@ -12,33 +13,39 @@ 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 __init__(self, nu_A:str, nu_a:str, nu_b:str, JB:str, orbital:str, ExB:float, ELabPerU:float): + self.SpinBalanced = self.ReactionDigestion(nu_A, nu_a, 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): + def ReactionDigestion(self, nu_A:str, nu_a: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.A_B = self.A_A + self.A_a - self.A_b + self.Z_B = self.Z_A + self.Z_a - self.Z_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) + mass_B = iso.GetMassFromAZ(self.A_B, self.Z_B) ExB = ExB - # sym_A = iso.GetSymbol(A_A, Z_A) - # sym_B = iso.GetSymbol(A_B, Z_B) + self.sym_A = iso.GetSymbol(self.A_A, self.Z_A) + self.sym_B = iso.GetSymbol(self.A_B, self.Z_B) + + nu_B = f"{self.A_B}{self.sym_B}" + # print(nu_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))) + print("-------- spin_B",self.spin_B) + if self.A_a == 2 and self.Z_a == 1: self.spin_a = 1.0 self.spin_b = 0.5 @@ -69,10 +76,10 @@ class ReactionData: index = match.start() # Get position of the first letter self.node = int(orbital[:index]) - l_sym = orbital[index:index+1] + self.l_sym = orbital[index:index+1] j_sym = orbital[index+1:] self.j = eval(j_sym) - self.l = op.ConvertLSym(l_sym) + self.l = op.ConvertLSym(self.l_sym) self.j = approximate_to_half_integer(self.j) self.s = approximate_to_half_integer(self.s) @@ -109,11 +116,16 @@ class ReactionData: self.Q_value = mass_A + mass_a - mass_b - mass_B - ExB self.dwI = DistortedWave(nu_A, nu_a, self.ELab) + self.mass_I = self.dwI.mu + self.k_I = self.dwI.k + Ecm_I = self.dwI.Ecm - Ecm_O = Ecm_I + self.Q_value + 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) + self.mass_O = self.dwO.mu + Eout2 = self.ELab + self.Q_value #this is incorrec, but used in ptolmey infileCreator print("==================================================") diff --git a/Raphael/solveSE.py b/Raphael/solveSE.py index 315e262..b8726f1 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: diff --git a/dwuck4/extractData.py b/dwuck4/extractData.py index c285d8b..4fa56c1 100755 --- a/dwuck4/extractData.py +++ b/dwuck4/extractData.py @@ -289,7 +289,8 @@ def plot_Xsec(data, isRuth = False): global plotID x_data, y_data = data - plt.figure(figsize=(8, 5)) + # plt.figure(figsize=(8, 5)) + plt.figure() plt.plot(x_data, y_data, linestyle="-", color="b", label="Extracted Data") plt.xlabel("Angle [deg]") if isRuth: @@ -390,18 +391,19 @@ extract_LmaxSaSb() ## must be run first bs_data = extract_BoundState() # plot_BoundState(bs_data) -# sAmpIn, sAmpOut = extract_ScatAmp() -# plot_SMatrix(sAmpIn, sa) -# plot_SMatrix(sAmpOut, sb) +sAmpIn, sAmpOut = extract_ScatAmp() +plot_SMatrix(sAmpIn, sa) +plot_SMatrix(sAmpOut, sb) # elXsec_data = extract_ElasticXsec() # plot_Xsec(elXsec_data) xsec_data = extract_Xsec() plot_Xsec(xsec_data) - x_data, y_data = xsec_data for i, r in enumerate(x_data): + if i % 5 != 0: + continue print(f"{{{r:7.3f}, {y_data[i]:10.7f}}},") def plot_RadialMatrix2(ma:float, mb:float, isPlot:bool=True): @@ -425,8 +427,6 @@ def plot_RadialMatrix2(ma:float, mb:float, isPlot:bool=True): return radmat - - rList, dwIn, dwOut = extract_DistortedWave() def plot_DW(isIncoming:bool, L:int, m:float): if isIncoming : @@ -444,7 +444,6 @@ def CoulombPS(L, eta): r_list, bsW = bs_data interp_radial = interp.interp1d(r_list, bsW, kind='cubic') - def CalRadialIntgeral(L, ma, mb, isPlot:bool = True, verbose:int = 1): if isPlot : diff --git a/dwuck4/inFileCreatorDW.py b/dwuck4/inFileCreatorDW.py index c839eeb..c8eda80 100755 --- a/dwuck4/inFileCreatorDW.py +++ b/dwuck4/inFileCreatorDW.py @@ -2,6 +2,11 @@ import sys import os +if len(sys.argv) < 7: + print("Error: Not enough arguments provided.") + print("Usage: ./{sys.argv[0]} reaction target_gs-spin orbital spin-pi Ex ELab[Mev/u]") + sys.exit(1) + reaction = sys.argv[1] JA_pi = sys.argv[2] orbital = sys.argv[3] @@ -9,176 +14,18 @@ JB_pi = sys.argv[4] Ex = float(sys.argv[5]) ELab = float(sys.argv[6]) -if len(sys.argv) < 7: - print("Error: Not enough arguments provided.") - print("Usage: ./{sys.argv[0]} reaction target_gs-spin orbital spin-pi Ex ELab[Mev/u]") - sys.exit(1) - sys.path.append(os.path.join(os.path.dirname(__file__), '../Cleopatra')) +sys.path.append(os.path.join(os.path.dirname(__file__), '../Raphael')) from IAEANuclearData import IsotopeClass +import opticalPotentials as op +from reactionData import ReactionData ##################################################### # only for (d,p) or (p,d) using An & Cai, Kronning ##################################################### -import numpy as np import re -import matplotlib.pyplot as plt - -# 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 #================== digest reaction @@ -189,82 +36,46 @@ nu_a = nuclei[1] nu_b = nuclei[2] nu_B = nuclei[3] -iso = IsotopeClass() +reactionData = ReactionData(nu_A, nu_a, nu_b, JB_pi, orbital, Ex, ELab) -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) +sym_A = reactionData.sym_A +A_A = reactionData.A_A +Z_A = reactionData.Z_A +A_a = reactionData.A_a +Z_a = reactionData.Z_a +A_B = reactionData.A_B +Z_B = reactionData.Z_B +A_b = reactionData.A_b +Z_b = reactionData.Z_b +A_x = reactionData.A_x +Z_x = reactionData.Z_x +A_c = reactionData.A_c +Z_c = reactionData.Z_c +node = reactionData.node +l_sym = reactionData.l_sym -A_x = abs(A_a - A_b) -Z_x = abs(Z_a - Z_b) +spin_a = reactionData.spin_a +spin_b = reactionData.spin_b -#---- check mass number and charge number is balnaced -if A_A + A_a - A_b - A_B != 0 or Z_A + Z_a - Z_b - Z_B != 0 : - print("reaction is incorrect, mass or charge not balanced.") - exit() +l = reactionData.l +j = reactionData.j -#---- check is (d,p) or (p, d) -if (Z_a !=1 or Z_b != 1) or (A_a + A_b != 3) : - print("not (d,p) or (p,d) reaction. stop.") - exit() - -mass_A = iso.GetMassFromSym(nu_A) -mass_a = iso.GetMassFromSym(nu_a) -mass_b = iso.GetMassFromSym(nu_b) -mass_B = iso.GetMassFromSym(nu_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 + 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) -sym_B = iso.GetSymbol(A_B, Z_B) - -if A_a == 2 and Z_a == 1: - spin_a = 1.0 - spin_b = 0.5 -else: - spin_a = 0.5 - spin_b = 1.0 - -Q_value = mass_A + mass_a - mass_b - mass_B - Ex - -print(f"Q-value : {Q_value:10.6f} MeV") -print(f"Binding : {BindingEnergy:10.6f} MeV") - -#=================== 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:] -j = eval(j_sym) -l = ConvertLSym(l_sym) +Q_value = reactionData.Q_value +BindingEnergy = reactionData.BindingEnergy #=================== outfile name fileOutName = str(sym_A) + str(A_A) + "_" + str(nu_a) + str(nu_b) + "_" \ + str(node) + str(l_sym) + str(int(2*j)) + "_" + str(Ex) + "_" + str(ELab) + ".in" -print(fileOutName) #=================== 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 * ELab)/hbarc # wave number of incoming channel +mass_I = reactionData.mass_I # reduced mass of incoming channel +k_I = reactionData.k_I # 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}") + +print(f"file out : {fileOutName}") +print(f" max L : {maxL}") #=================== create outfile with open(fileOutName, "w") as file: @@ -274,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 : - AnCai(A_A, Z_A, A_a*ELab) + op.AnCai(A_A, Z_A, A_a*ELab) else: - Koning(A_A, Z_A, A_a*ELab, Z_a) + 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"{rc0:+08.4f}") + file.write(f"{op.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"{-v:+08.4f}") # real - file.write(f"{r0:+08.4f}") # - file.write(f"{a:+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"{"":8s}") # spin-orbit skipped - file.write(f"{-vi:+08.4f}") # imag - file.write(f"{ri0:+08.4f}") # - file.write(f"{ai:+08.4f}\n") # + file.write(f"{-op.vi:+08.4f}") # imag + file.write(f"{op.ri0:+08.4f}") # + file.write(f"{op.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*vsi:+08.4f}") # imag - file.write(f"{rsi0:+08.4f}") # - file.write(f"{asi:+08.4f}\n") # + file.write(f"{4*op.vsi:+08.4f}") # imag + file.write(f"{op.rsi0:+08.4f}") # + file.write(f"{op.asi:+08.4f}\n") # # Spin-Orbit file.write(f"{-4:+08.4f}") - file.write(f"{-4*vso:+08.4f}") # real - file.write(f"{rso0:+08.4f}") # - file.write(f"{aso:+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"{"":8s}") # spin-orbit skipped - file.write(f"{-4*vsoi:+08.4f}") # imag - file.write(f"{rsoi0:+08.4f}") # - file.write(f"{asoi:+08.4f}\n") # + file.write(f"{-4*op.vsoi:+08.4f}") # imag + file.write(f"{op.rsoi0:+08.4f}") # + file.write(f"{op.asoi:+08.4f}\n") # #===== Block 6 if A_a == 2 : - Koning(A_B, Z_B, A_a*ELab + Q_value - Ex, Z_b) + op.Koning(A_B, Z_B, A_a*ELab + Q_value - Ex, Z_b) else: - AnCai(A_B, Z_B, A_a*ELab + Q_value - Ex) + 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"{rc0:+08.4f}") + file.write(f"{op.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"{-v:+08.4f}") # real - file.write(f"{r0:+08.4f}") # - file.write(f"{a:+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"{"":8s}") # spin-orbit skipped - file.write(f"{-vi:+08.4f}") # imag - file.write(f"{ri0:+08.4f}") # - file.write(f"{ai:+08.4f}\n") # + file.write(f"{-op.vi:+08.4f}") # imag + file.write(f"{op.ri0:+08.4f}") # + file.write(f"{op.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*vsi:+08.4f}") # imag - file.write(f"{rsi0:+08.4f}") # - file.write(f"{asi:+08.4f}\n") # + file.write(f"{4*op.vsi:+08.4f}") # imag + file.write(f"{op.rsi0:+08.4f}") # + file.write(f"{op.asi:+08.4f}\n") # # Spin-Orbit file.write(f"{-4:+08.4f}") - file.write(f"{-4*vso:+08.4f}") # real - file.write(f"{rso0:+08.4f}") # - file.write(f"{aso:+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"{"":8s}") # spin-orbit skipped - file.write(f"{-4*vsoi:+08.4f}") # imag - file.write(f"{rsoi0:+08.4f}") # - file.write(f"{asoi:+08.4f}\n") # + file.write(f"{-4*op.vsoi:+08.4f}") # imag + file.write(f"{op.rsoi0:+08.4f}") # + file.write(f"{op.asoi:+08.4f}\n") # #====== bound state file.write(f"{BindingEnergy:+08.4f}") file.write(f"{A_x:+08.4f}")