535 lines
16 KiB
Python
Executable File
535 lines
16 KiB
Python
Executable File
#!/usr/bin/env python3
|
|
import sys
|
|
|
|
filename = sys.argv[1]
|
|
|
|
if len(sys.argv) < 2:
|
|
print("Error: Not enough arguments provided.")
|
|
print("Usage: ./{sys.argv[0]} filename")
|
|
sys.exit(1)
|
|
|
|
#####################################################
|
|
import numpy as np
|
|
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):
|
|
r_data, y_data = [], []
|
|
|
|
start_line = None
|
|
|
|
with open(file_path, "r") as file:
|
|
for i, line in enumerate(file, start=1):
|
|
if "0 R RL,R" in line:
|
|
if start_line is None:
|
|
start_line = i+1
|
|
|
|
if len(r_data) > 1 and line[0] == "0" :
|
|
start_line = None
|
|
break
|
|
|
|
if start_line != None:
|
|
if start_line <= i :
|
|
columns = line.split()
|
|
if len(columns) >= 2:
|
|
r_data.append(float(columns[0])) # Convert to float
|
|
y_data.append(float(columns[1]))
|
|
|
|
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
|
|
|
|
with open(file_path, "r") as file:
|
|
for i, line in enumerate(file, start=1):
|
|
|
|
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 :
|
|
start_line = None
|
|
break
|
|
|
|
if start_line != None:
|
|
if start_line <= i :
|
|
columns = line.split()
|
|
temp = []
|
|
if len(columns) >= 2:
|
|
if line[0] == "+" : # these line are for outgoing
|
|
for i in range(int(2*(2*sb+1)+1)):
|
|
temp.append(float(columns[i+1]))
|
|
sAmpOut.append(temp)
|
|
else:
|
|
for i in range(int(2*(2*sa+1)+1)):
|
|
temp.append(float(columns[i]))
|
|
sAmpIn.append(temp)
|
|
|
|
return sAmpIn, sAmpOut
|
|
|
|
#-------------------------------------------------------
|
|
def extract_ElasticXsec(file_path=filename):
|
|
x_data, y_data = [], []
|
|
|
|
start_line = None
|
|
|
|
with open(file_path, "r") as file:
|
|
for i, line in enumerate(file, start=1):
|
|
if "0 Theta Sig(1)/Coul Sigma(1)" in line:
|
|
if start_line is None:
|
|
start_line = i+1
|
|
|
|
if len(x_data) > 1 and line[0] == "0" :
|
|
start_line = None
|
|
break
|
|
|
|
if start_line != None:
|
|
if start_line <= i :
|
|
columns = line.split()
|
|
if len(columns) >= 2:
|
|
x_data.append(float(columns[0])) # Convert to float
|
|
y_data.append(float(columns[1]))
|
|
|
|
return [x_data, y_data]
|
|
|
|
|
|
#-------------------------------------------------------
|
|
def extract_Xsec(file_path=filename):
|
|
x_data, y_data = [], []
|
|
|
|
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:
|
|
if start_line is None:
|
|
start_line = i+1
|
|
|
|
if len(x_data) > 1 and "0Tot-sig" in line :
|
|
start_line = None
|
|
break
|
|
|
|
if start_line != None:
|
|
if start_line <= i :
|
|
columns = line.split()
|
|
if len(columns) >= 2:
|
|
x_data.append(float(columns[0])) # Convert to float
|
|
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
|
|
|
|
with open(file_path, "r") as file:
|
|
for i, line in enumerate(file, start=1):
|
|
|
|
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 :
|
|
start_line = None
|
|
break
|
|
|
|
if start_line != None:
|
|
if start_line <= i :
|
|
l_data = int(line[:4])
|
|
real_data = float(line[4:14])
|
|
imag_data = float(line[14:])
|
|
data.append([l_data, real_data, imag_data])
|
|
|
|
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")
|
|
plt.xlabel("Radius [fm]")
|
|
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:
|
|
l_data.append(a[0])
|
|
real_data.append(a[1])
|
|
imag_data.append(a[2])
|
|
|
|
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 " + 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")
|
|
plt.xlabel("Angle [deg]")
|
|
if isRuth:
|
|
plt.ylabel("d.s.c / Ruth ")
|
|
else:
|
|
plt.ylabel("d.s.c [mb/sr]")
|
|
plt.yscale('log')
|
|
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, spin):
|
|
global plotSMID
|
|
|
|
l_data = []
|
|
nSpin = int(2*spin+1)
|
|
|
|
fig, axes = plt.subplots(1, nSpin, figsize=(6*nSpin, 4))
|
|
|
|
smIn = []
|
|
for a in data:
|
|
l_data.append(int(a[0]))
|
|
tempSM = []
|
|
for i in range(nSpin):
|
|
real = a[2*i+1]
|
|
imag = a[2*i+2]
|
|
tempSM.append((real + imag * 1j)*2j + 1 )
|
|
smIn.append(tempSM)
|
|
|
|
for i in range(nSpin):
|
|
real_parts = [entry[i].real for entry in smIn]
|
|
imag_parts = [entry[i].imag for entry in smIn]
|
|
|
|
start = 1
|
|
if i == nSpin -1 :
|
|
start = 0
|
|
|
|
# Plot real part vs L
|
|
axes[i].plot(l_data[start:], real_parts[start:], label=f'Re', marker='o')
|
|
|
|
# Plot imaginary part vs L
|
|
axes[i].plot(l_data[start:], imag_parts[start:], label=f'Im', marker='x')
|
|
|
|
# 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-spin/2.:+.1f}')
|
|
|
|
# Add grid lines
|
|
axes[i].set_xlim(-1, max(l_data)+1)
|
|
axes[i].set_ylim(-1.1, 1.1)
|
|
axes[i].set_xticks(np.arange(0, max(l_data)+3, 5)) # Set x-ticks from 0 to 20 with step 5
|
|
axes[i].grid(True)
|
|
|
|
# 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)
|
|
|
|
|
|
##############################################################################
|
|
##############################################################################
|
|
|
|
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)
|
|
|
|
# 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):
|
|
print(f"{{{r:7.3f}, {y_data[i]:10.7f}}},")
|
|
|
|
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"
|
|
|
|
radmat = extract_RadialMatrix(str_a, str_b)
|
|
|
|
if int(2*ma)%2 == 1 :
|
|
msg_a = f"Ja = La + {int(2*ma)}/2"
|
|
else:
|
|
msg_a = f"Ja = La + {int(ma)}"
|
|
|
|
if int(2*mb)%2 == 1 :
|
|
msg_b = f"Jb = Lb + {int(2*mb)}/2"
|
|
else:
|
|
msg_b = f"Jb = Jb + {int(mb)}"
|
|
|
|
if isPlot:
|
|
plot_RadialMatrix(radmat, f": {msg_a}, {msg_b}")
|
|
|
|
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 = 1 #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)
|
|
|
|
|
|
# radMat = plot_RadialMatrix2(-1, -0.5, True)
|
|
|
|
# for a in radMat:
|
|
# ll, real, imag = a
|
|
# print(f"{{{int(ll):2d}, {real:10.7f} + {imag:10.7f} I}},")
|
|
|
|
#================================================ 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.") |