From bd83983e40d68d78cb4ff1d4c56eb71ad73725d7 Mon Sep 17 00:00:00 2001 From: "Ryan@SOLARIS_testStation" Date: Wed, 12 Feb 2025 20:44:02 -0500 Subject: [PATCH] added comments on DWUCK4, update extractData for DWUCK4, add DWtest.DAT for 16O(d,p) --- dwuck4/BDWCK4.FOR | 4 +- dwuck4/DWtest.DAT | 16 ++ dwuck4/extractData.py | 343 ++++++++++++++++++++++++++++++++++-------- 3 files changed, 298 insertions(+), 65 deletions(-) create mode 100644 dwuck4/DWtest.DAT diff --git a/dwuck4/BDWCK4.FOR b/dwuck4/BDWCK4.FOR index 7622e2b..92978a2 100644 --- a/dwuck4/BDWCK4.FOR +++ b/dwuck4/BDWCK4.FOR @@ -534,8 +534,8 @@ C NORMALIZE RADIAL FUNCTIONS C CTEMP(1)=FR(IX1)*C(IX1)-FR(IX2)*C(IX2) CTEMP(2)=FR(IX1)*C(IX2)+FR(IX2)*C(IX1) - FR(IX1)=CTEMP(1) - FR(IX2)=CTEMP(2) + FR(IX1)=CTEMP(1) ! real + FR(IX2)=CTEMP(2) ! complex 30 CONTINUE IX=IX+LPL2 39 CONTINUE diff --git a/dwuck4/DWtest.DAT b/dwuck4/DWtest.DAT new file mode 100644 index 0000000..90b9dff --- /dev/null +++ b/dwuck4/DWtest.DAT @@ -0,0 +1,16 @@ +10001310500100000 16O(D,P)17O G.S. s1/2 orbital ++181. +00. +01.0 ++15+01+00+01 ++00.10 +12. ++20.00 +02. +01. +16. +08. +01.30 +02. ++01. -88.955 +01.149 +00.675 -02.348 +01.345 +00.603 ++02. +01.394 +00.687 +40.872 +01.394 +00.687 +-04. -14.228 +00.972 +01.011 +01.562 +00.477 ++01.92 +01. +01. +17. +08. +01.42 +01. ++01. -49.544 +01.146 +00.675 -02.061 +01.146 +00.675 ++02. +30.680 +01.302 +00.528 +-04. -21.184 +00.934 +00.590 +00.424 +00.934 +00.590 +-04.14 +01. +00. +16. +08. +01.30 +01. +-01. -01. +01.10 +00.65 +24. ++01. +00. +01. +01. +58. +9 END OF DATA DWUCK4 test cases \ No newline at end of file diff --git a/dwuck4/extractData.py b/dwuck4/extractData.py index 5d09090..e23b1b5 100755 --- a/dwuck4/extractData.py +++ b/dwuck4/extractData.py @@ -1,12 +1,11 @@ #!/usr/bin/env python3 import sys -s = float(sys.argv[1]) -filename = sys.argv[2] +filename = sys.argv[1] -if len(sys.argv) < 3: +if len(sys.argv) < 2: print("Error: Not enough arguments provided.") - print("Usage: ./{sys.argv[0]} spin filename") + print("Usage: ./{sys.argv[0]} filename") sys.exit(1) ##################################################### @@ -15,8 +14,34 @@ import re import numpy as np import matplotlib.pyplot as plt +Lmax = None +sa = None +sb = None + +def extract_LmaxSaSb(file_path=filename): + global Lmax, sa, sb + + with open(file_path, "r") as file: + for i, line in enumerate(file, start=1): + + if Lmax is None and "LMAX" in line : + columns = line.split() + Lmax = int(columns[6]) + print(f"LMAX = {Lmax}") + + if sa != None and sb is None and "2*STR" in line: + columns = line.split() + sb = int(float(columns[-1]))/2. + print(f"Spin out = {sb}") + + if sa is None and "2*STR" in line: + columns = line.split() + sa = int(float(columns[-1]))/2. + print(f"Spin In = {sa}") + + def extract_BoundState(file_path=filename): - x_data, y_data = [], [] + r_data, y_data = [], [] start_line = None @@ -26,7 +51,7 @@ def extract_BoundState(file_path=filename): if start_line is None: start_line = i+1 - if len(x_data) > 1 and line[0] == "0" : + if len(r_data) > 1 and line[0] == "0" : start_line = None break @@ -34,45 +59,84 @@ def extract_BoundState(file_path=filename): if start_line <= i : columns = line.split() if len(columns) >= 2: - x_data.append(float(columns[0])) # Convert to float + r_data.append(float(columns[0])) # Convert to float y_data.append(float(columns[1])) - return [x_data, y_data] + return [r_data, y_data] + +#------------------------------------------------------- +def extract_DistortedWave(file_path=filename): + + if Lmax is None or sa is None or sb is None : + print("Please run the extract_LmaxSaSb() first") + return + + r_list = [0.0] + distWaveIn = np.empty((Lmax+1, int(2*(2*sa+1))), dtype=object) # going to be Lmax * (2*sa+1) + distWaveOut = np.empty((Lmax+1, int(2*(2*sb+1))), dtype=object) # going to be Lmax * (2*sa+1) + + for i in range(Lmax+1): + for j in range(int(2*(2*sa+1))): + distWaveIn[i, j] = [0.0] + for j in range(int(2*(2*sb+1))): + distWaveOut[i, j] = [0.0] + + start_line = None + + l_list = [] + + with open(file_path, "r") as file: + for i, line in enumerate(file, start=1): + + if "0R1=" in line: + if start_line is None: + start_line = i+2 + columns = line.split() + r_list.append(float(columns[1])) + l_list = [] + + if len(l_list) > Lmax: + start_line = None + + if start_line != None: + if start_line <= i : + # print(line) + if line[0] == "+" : + l = int(line[70:72]) + for j in range(int(2*sb+1)): + distWaveOut[l,2*j].append(float(line[20*j+72: 20*j+82])) + distWaveOut[l,2*j+1].append(float(line[20*j+82:20*j+92])) + else: + l = int(line[:4]) + l_list.append(l) + for j in range(int(2*sa+1)): + # print(float(line[20*j+4: 20*j+14])) + distWaveIn[l,2*j].append(float(line[20*j+4: 20*j+14])) + distWaveIn[l,2*j+1].append(float(line[20*j+14:20*j+24])) + + return r_list, distWaveIn, distWaveOut + #------------------------------------------------------- def extract_ScatAmp(file_path=filename): + + if Lmax is None or sa is None or sb is None : + print("Please run the extract_LmaxSaSb() first") + return + sAmpIn = [] sAmpOut = [] start_line = None - LMAX = None - Spin2In = None - Spin2Out = None with open(file_path, "r") as file: for i, line in enumerate(file, start=1): - - if LMAX is None and "LMAX" in line : - columns = line.split() - LMAX = int(columns[6]) - print(f"LMAX = {LMAX}") - - if Spin2In != None and Spin2Out is None and "2*STR" in line: - columns = line.split() - Spin2Out = int(float(columns[-1])) - print(f"2*Spin out = {Spin2Out}") - - if Spin2In is None and "2*STR" in line: - columns = line.split() - Spin2In = int(float(columns[-1])) - print(f"2*Spin In = {Spin2In}") - if "L REAL D1 IMAG D1 REAL D2" in line: if start_line is None: start_line = i+1 - if LMAX != None and len(sAmpIn) > LMAX : + if Lmax != None and len(sAmpIn) > Lmax : start_line = None break @@ -82,15 +146,15 @@ def extract_ScatAmp(file_path=filename): temp = [] if len(columns) >= 2: if line[0] == "+" : # these line are for outgoing - for i in range(2*(Spin2Out+1)+1): + for i in range(int(2*(2*sb+1)+1)): temp.append(float(columns[i+1])) sAmpOut.append(temp) else: - for i in range(2*(Spin2In+1)+1): + for i in range(int(2*(2*sa+1)+1)): temp.append(float(columns[i])) sAmpIn.append(temp) - return Spin2In, sAmpIn, Spin2Out, sAmpOut + return sAmpIn, sAmpOut #------------------------------------------------------- def extract_ElasticXsec(file_path=filename): @@ -124,7 +188,6 @@ def extract_Xsec(file_path=filename): start_line = None - with open(file_path, "r") as file: for i, line in enumerate(file, start=1): if "0 Theta Inelsig,fm**2" in line: @@ -140,29 +203,29 @@ def extract_Xsec(file_path=filename): columns = line.split() if len(columns) >= 2: x_data.append(float(columns[0])) # Convert to float - y_data.append(float(columns[1])) + y_data.append(float(columns[1])*10) ## factor convert fm^2 to mb return [x_data, y_data] #------------------------------------------------------- def extract_RadialMatrix(ma:str, mb:str, file_path=filename): + + if Lmax is None or sa is None or sb is None : + print("Please run the extract_LmaxSaSb() first") + return + data = [] start_line = None - LMAX = None with open(file_path, "r") as file: for i, line in enumerate(file, start=1): - if LMAX is None and "LMAX" in line : - columns = line.split() - LMAX = int(columns[6]) - if "0 RADIAL MATRIX ELEMENTS," in line and ma in line and mb in line: if start_line is None: start_line = i+2 - if LMAX != None and len(data) > LMAX : + if Lmax != None and len(data) > Lmax : start_line = None break @@ -176,7 +239,12 @@ def extract_RadialMatrix(ma:str, mb:str, file_path=filename): return data #------------------------------------------------------- +plotID = 0 +plotSMID = 0 + def plot_BoundState(data): + global plotID + x_data, y_data = data plt.figure(figsize=(8, 5)) plt.plot(x_data, y_data, marker="o", linestyle="-", color="b", label="Extracted Data") @@ -184,9 +252,16 @@ def plot_BoundState(data): plt.ylabel("Value") plt.title("Bound state radial function") plt.grid(True) + + plotID = plotID + 1 + manager = plt.get_current_fig_manager() + manager.window.wm_geometry(f"+{850*plotID}+50") + plt.show(block=False) def plot_RadialMatrix(data, msg:str = ""): + global plotID + l_data, real_data, imag_data = [], [], [] for a in data: @@ -201,10 +276,18 @@ def plot_RadialMatrix(data, msg:str = ""): plt.ylabel("Value") plt.title("Radial Matrix " + msg) plt.grid(True) + plt.legend() + + plotID = plotID + 1 + manager = plt.get_current_fig_manager() + manager.window.wm_geometry(f"+{850*plotID}+50") + plt.show(block=False) def plot_Xsec(data, isRuth = False): + global plotID + x_data, y_data = data plt.figure(figsize=(8, 5)) plt.plot(x_data, y_data, linestyle="-", color="b", label="Extracted Data") @@ -217,13 +300,19 @@ def plot_Xsec(data, isRuth = False): plt.xlim(-5, 185) plt.xticks(np.arange(0, 181, 20)) plt.grid(True) + + plotID = plotID + 1 + manager = plt.get_current_fig_manager() + manager.window.wm_geometry(f"+{850*plotID}+50") + plt.show(block=False) -def plot_SMatrix(data, spin2): +def plot_SMatrix(data, spin): + global plotSMID l_data = [] - nSpin = spin2+1 + nSpin = int(2*spin+1) fig, axes = plt.subplots(1, nSpin, figsize=(6*nSpin, 4)) @@ -254,7 +343,7 @@ def plot_SMatrix(data, spin2): # Adding labels and title axes[i].set_xlabel('L') axes[i].set_ylabel('Value') - axes[i].set_title(f'Real and Imaginary Parts vs L for Spin {i-spin2/2.:+.1f}') + axes[i].set_title(f'Real and Imaginary Parts vs L for Spin {i-spin/2.:+.1f}') # Add grid lines axes[i].set_xlim(-1, max(l_data)+1) @@ -265,44 +354,172 @@ def plot_SMatrix(data, spin2): # Show the legend axes[i].legend() + plotSMID = plotSMID + 1 + manager = plt.get_current_fig_manager() + manager.window.wm_geometry(f"+850+{650*plotSMID}") + # Adjust layout to prevent overlapping subplots plt.tight_layout() plt.show(block=False) + +def plot_DistortWave(r_list, dw_real, dw_imag, title:str = ""): + global plotID + + plt.figure(figsize=(8, 5)) + plt.plot(r_list, dw_real, marker="o", linestyle="-", color="b", label="Re") + plt.plot(r_list, dw_imag, marker="x", linestyle="-", color="r", label="Im") + plt.xlabel("Radius [fm]") + plt.ylabel("Value") + plt.title("Distorted Wave radial function: " + title) + plt.legend() + plt.grid(True) + + plotID = plotID + 1 + manager = plt.get_current_fig_manager() + manager.window.wm_geometry(f"+{850*plotID}+50") + + plt.show(block=False) + + ############################################################################## ############################################################################## -# bs_data = extract_BoundState() -# plot_BoundState(bs_data) +extract_LmaxSaSb() ## must be run first -# spin2In, sAmpIn, spin2Out, sAmpOut = extract_ScatAmp() +bs_data = extract_BoundState() +plot_BoundState(bs_data) -# plot_SMatrix(sAmpIn, spin2In) -# plot_SMatrix(sAmpOut, spin2Out) +sAmpIn, sAmpOut = extract_ScatAmp() +plot_SMatrix(sAmpIn, sa) +plot_SMatrix(sAmpOut, sb) -# elXsec_data = extract_ElasticXsec() -# plot_Xsec(elXsec_data) +elXsec_data = extract_ElasticXsec() +plot_Xsec(elXsec_data) -# xsec_data = extract_Xsec() -# plot_Xsec(xsec_data) +xsec_data = extract_Xsec() +plot_Xsec(xsec_data) -radmat0 = extract_RadialMatrix("+ 2/2", "+ 1/2") -plot_RadialMatrix(radmat0, "+1,+1/2") +def plot_RadialMatrix2(ma:float, mb:float, isPlot:bool=True): + str_a = f"+{int(2*ma):2.0f}/2" + str_b = f"+{int(2*mb):2.0f}/2" -radmat1 = extract_RadialMatrix("+ 2/2", "+-1/2") -plot_RadialMatrix(radmat1, "+1,-1/2") + radmat = extract_RadialMatrix(str_a, str_b) -radmat2 = extract_RadialMatrix("+ 0/2", "+ 1/2") -plot_RadialMatrix(radmat2, "+0,+1/2") + if int(2*ma)%2 == 1 : + msg_a = f"Ja = La + {int(2*ma)}/2" + else: + msg_a = f"Ja = La + {int(ma)}" -radmat3 = extract_RadialMatrix("+ 0/2", "+-1/2") -plot_RadialMatrix(radmat3, "+0,-1/2") + if int(2*mb)%2 == 1 : + msg_b = f"Jb = Lb + {int(2*mb)}/2" + else: + msg_b = f"Jb = Jb + {int(mb)}" -radmat4 = extract_RadialMatrix("+-2/2", "+ 1/2") -plot_RadialMatrix(radmat4, "-1,+1/2") + if isPlot: + plot_RadialMatrix(radmat, f": {msg_a}, {msg_b}") -radmat5 = extract_RadialMatrix("+-2/2", "+-1/2") -plot_RadialMatrix(radmat5, "-1,-1/2") + return radmat + + + +rList, dwIn, dwOut = extract_DistortedWave() +def plot_DW(isIncoming:bool, L:int, m:float): + if isIncoming : + k = int(m + sa) + plot_DistortWave(rList, dwIn[L][2*k], dwIn[L][2*k+1], f"Incoming L={L}, J=L+{m}") + else: + k = int(m + sb) + plot_DistortWave(rList, dwOut[L][2*k], dwOut[L][2*k+1], f"Outgoing L={L}, J=L+{m}") + +import scipy.special as sp +import scipy.interpolate as interp +def CoulombPS(L, eta): + return np.angle(sp.gamma(1 + L + 1j * 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 : + plot_DW(True, L, ma) + plot_DW(False, L, mb) + + etaI = 0.3997 + etaO = 0.276894 + + radmat = plot_RadialMatrix2(ma, mb, isPlot) + + prod_re, prod_im = [], [] + + total = 0 + for i, r in enumerate(rList): + ia = int(2*(ma + sa)) + ib = int(2*(mb + sb)) + dw_a = dwIn[L][ia][i] + dwIn[L][ia+1][i] * 1j + dw_b = dwOut[L][ib][i] + dwOut[L][ib+1][i] * 1j + bound = 0 + if rList[i] <= max(r_list) and rList[i] >= min(r_list) : + bound = interp_radial(r) + total += dw_a * dw_b * bound + prod_re.append(np.real(dw_a * dw_b * bound)) + prod_im.append(np.imag(dw_a * dw_b * bound)) + if verbose >= 1: + print(f"{i:3d} {bound:8.5f}, {rList[i]:4.1f}, {np.real(dw_a):8.5f}, {np.imag(dw_a):8.5f}, ({np.abs(dw_a):8.5f}), {np.real(dw_b):8.5f}, {np.imag(dw_b):8.5f} | {np.real(dw_a * dw_b * bound):9.6f}, {np.imag(dw_a * dw_b * bound):9.6f}") + + total = total * 0.1 * 17./16. + phase = np.exp( 1j * (CoulombPS(L, etaI)- CoulombPS(L, etaO)) ) + + if verbose >= 2: + print("-------------------------") + print(f" Radial integral before CoulombPS : {total}") + print(f" CoulombPS : {phase}") + + total = phase * total + print(f" total : [{L}, {np.real(total):8.5f}, {np.imag(total):8.6f}]") + print(f" DWUCK4 : {radmat[L]}") + + # plt.figure(figsize=(8, 5)) + # plt.plot(rList, prod_re, marker="o", linestyle="-", color="b", label="Real") + # plt.plot(rList, prod_im, marker="x", linestyle="-", color="r", label="Imag") + # plt.xlabel("L") + # plt.ylabel("Value") + # plt.title("Product of Radial") + # plt.grid(True) + # plt.legend() + # plt.show(block=False) + + return total + +CalRadialIntgeral(3, 1, 0.5) + + +#================================================ cal Radial matrix and Plot +# radialIn = [] + +# for ll in range(15): +# radialIn.append([ll, CalRadialIntgeral(ll, 1, 0.5, False)]) + + +# l_data, real_data, imag_data = [], [], [] + +# for a in radialIn: +# l_data.append(a[0]) +# real_data.append(np.real(a[1])) +# imag_data.append(np.imag(a[1])) + +# plt.figure(figsize=(8, 5)) +# plt.plot(l_data, real_data, marker="o", linestyle="-", color="b", label="Real") +# plt.plot(l_data, imag_data, marker="x", linestyle="-", color="r", label="Imag") +# plt.xlabel("L") +# plt.ylabel("Value") +# plt.title("Radial Matrix (cal)") +# plt.grid(True) +# plt.show(block=False) + +# plot_RadialMatrix2(1, 1/2) input("Press Enter to exit.") \ No newline at end of file