added comments on DWUCK4, update extractData for DWUCK4, add DWtest.DAT for 16O(d,p)
This commit is contained in:
parent
04505f81c1
commit
bd83983e40
|
|
@ -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
|
||||
|
|
|
|||
16
dwuck4/DWtest.DAT
Normal file
16
dwuck4/DWtest.DAT
Normal file
|
|
@ -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
|
||||
|
|
@ -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.")
|
||||
Loading…
Reference in New Issue
Block a user