rtang/PtolemyGUI
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity

completed distortedwave class, calculation should be optimized when calculate the elastic sum, but don't know how

Browse Source
This commit is contained in:
Ryan@Home 2025-02-20 00:02:12 -05:00
parent bd74ebdddd
commit db40a2bff0
3 changed files with 84 additions and 38 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

36
Raphael/DWBA_ZR.py
View File

@ -2,6 +2,7 @@
from boundState import BoundState from boundState import BoundState
from solveSE import WoodsSaxonPot, CoulombPotential, SpinOrbit_Pot, WS_SurfacePot from solveSE import WoodsSaxonPot, CoulombPotential, SpinOrbit_Pot, WS_SurfacePot
import matplotlib.pyplot as plt
# boundState = BoundState(16, 8, 1, 0, 1, 0, 0.5, -4.14) # boundState = BoundState(16, 8, 1, 0, 1, 0, 0.5, -4.14)
# boundState.SetPotential(1.10, 0.65, -6, 1.25, 0.65, 1.25) # boundState.SetPotential(1.10, 0.65, -6, 1.25, 0.65, 1.25)
@ -11,18 +12,37 @@ from solveSE import WoodsSaxonPot, CoulombPotential, SpinOrbit_Pot, WS_SurfacePo
from distortedWave import DistortedWave from distortedWave import DistortedWave
dw = DistortedWave("60Ni", "p", 30) # dw = DistortedWave("60Ni", "p", 30)
# dw.ClearPotential()
# dw.AddPotential(WoodsSaxonPot(-47.937-2.853j, 1.20, 0.669), False)
# dw.AddPotential(WS_SurfacePot(-6.878j, 1.28, 0.550), False)
# dw.AddPotential(SpinOrbit_Pot(-5.250 + 0.162j, 1.02, 0.590), False)
# dw.AddPotential(CoulombPotential(1.258), False)
dw = DistortedWave("60Ni", "d", 60)
dw.PrintInput()
dw.ClearPotential() dw.ClearPotential()
dw.AddPotential(WoodsSaxonPot(-47.937-2.853j, 1.20, 0.669), False) dw.AddPotential(WoodsSaxonPot(-81.919, 1.15, 0.768), False)
dw.AddPotential(WS_SurfacePot(-6.878j, 1.28, 0.550), False) dw.AddPotential(WoodsSaxonPot(-4.836j, 1.33, 0.464), False)
dw.AddPotential(SpinOrbit_Pot(-5.250 + 0.162j, 1.02, 0.590), False) dw.AddPotential(WS_SurfacePot(-8.994j, 1.373, 0.774), False)
dw.AddPotential(CoulombPotential(1.258), False) dw.AddPotential(SpinOrbit_Pot(-3.557, 0.972, 1.011), False)
dw.AddPotential(CoulombPotential(1.303), False)
dw.CalScatteringMatrix() dw.CalScatteringMatrix()
dw.PrintScatteringMatrix() dw.PrintScatteringMatrix()
# dw.PlotScatteringMatrix() dw.PlotDCSUnpolarized(180, 1)
dw.PlotDCSUnpolarized() # for i in range(1, 19):
# theta = 10*i
# # ruth = dw.RutherFord(theta)
# # coulAmp = dw.CoulombScatterintAmp(theta)
# dw.CalLegendre(theta)
# nuAmp1 = dw.NuclearScatteringAmp(-0.5, 0.5, 14)
# nuAmp2 = dw.NuclearScatteringAmp(0.5, -0.5, 14)
# # dsc = dw.DCSUnpolarized(theta, 14)
# # print(f"{theta:3.0f}, {nuAmp1:15.5f}, {nuAmp2:15.5f}, {dsc:10.6f}, {ruth:10.6f}")
# print(f"{theta:3.0f}, {nuAmp1:15.5f}, {nuAmp2:15.5f}")

18
Raphael/assLegendreP.py
View File

@ -2,36 +2,36 @@
import numpy as np import numpy as np
def associated_legendre_polynomials_single_angle(L, M, theta_deg): def associated_legendre_array(maxL, maxM, theta_deg):
# Convert theta from degrees to radians # Convert theta from degrees to radians
theta = np.radians(theta_deg) theta = np.radians(theta_deg)
x = np.cos(theta) x = np.cos(theta)
P = np.zeros((L + 1, M + 1)) P = np.zeros((maxL + 1, maxM + 1))
# P^m_l for m = 0, l = 0 # P^m_l for m = 0, l = 0
P[0, 0] = 1.0 P[0, 0] = 1.0
# P^m_l for m = 0, l > 0 # P^m_l for m = 0, l > 0
for l in range(1, L + 1): for l in range(1, maxL + 1):
P[l, 0] = ((2*l - 1) * x * P[l-1, 0] - (l - 1) * P[l-2, 0]) / l P[l, 0] = ((2*l - 1) * x * P[l-1, 0] - (l - 1) * P[l-2, 0]) / l
# P^m_l for m > 0 (using recursion) # P^m_l for m > 0 (using recursion)
for m in range(1, M + 1): for m in range(1, maxM + 1):
# P^m_m # P^m_m
P[m, m] = (1 - 2*m) * P[m-1, m-1] P[m, m] = (1 - 2*m) * P[m-1, m-1]
P[m, m] *= np.sqrt(1 - x**2) P[m, m] *= np.sqrt(1 - x**2)
for l in range(m + 1, L + 1): for l in range(m + 1, maxL + 1):
# P^m_l for m < l # P^m_l for m < l
P[l, m] = ((2*l - 1) * x * P[l-1, m] - (l + m - 1) * P[l-2, m]) / (l - m) P[l, m] = ((2*l - 1) * x * P[l-1, m] - (l + m - 1) * P[l-2, m]) / (l - m)
return P return P
# Example usage # Example usage
#L = 15 # Maximum l degree # L = 15 # Maximum l degree
#M = 5 # Maximum m order # M = 3 # Maximum m order
#legendre_polynomials = associated_legendre_polynomials_single_angle(L, M, 45) # legendre_polynomials = associated_legendre_array(L, M, 10)
#print(legendre_polynomials) # print(legendre_polynomials)

68
Raphael/distortedWave.py
View File

@ -9,6 +9,7 @@ import numpy as np
from scipy.special import gamma, sph_harm, factorial from scipy.special import gamma, sph_harm, factorial
import matplotlib.pyplot as plt import matplotlib.pyplot as plt
from assLegendreP import associated_legendre_array
def SevenPointsSlope(data, n): def SevenPointsSlope(data, n):
return (-data[n + 3] + 9 * data[n + 2] - 45 * data[n + 1] + 45 * data[n - 1] - 9 * data[n - 2] + data[n - 3]) / 60 return (-data[n + 3] + 9 * data[n + 2] - 45 * data[n + 1] + 45 * data[n - 1] - 9 * data[n - 2] + data[n - 3]) / 60
@ -29,6 +30,7 @@ def KroneckerDelta(i, j):
else: else:
return 0 return 0
############################################################
class DistortedWave(SolvingSE): class DistortedWave(SolvingSE):
def __init__(self, target, projectile, ELab): def __init__(self, target, projectile, ELab):
super().__init__(target, projectile, ELab) super().__init__(target, projectile, ELab)
@ -38,6 +40,8 @@ class DistortedWave(SolvingSE):
self.ScatMatrix = [] self.ScatMatrix = []
self.distortedWaveU = [] self.distortedWaveU = []
self.legendreArray = []
def SetLJ(self, L, J): def SetLJ(self, L, J):
self.L = L self.L = L
self.J = J self.J = J
@ -105,9 +109,15 @@ class DistortedWave(SolvingSE):
def PrintScatteringMatrix(self): def PrintScatteringMatrix(self):
for L in range(0, len(self.ScatMatrix)): for L in range(0, len(self.ScatMatrix)):
print("{", end="")
for i in range(0, len(self.ScatMatrix[L])): for i in range(0, len(self.ScatMatrix[L])):
print(f"{{{L:2d},{i-self.S:4.1f}, {np.real(self.ScatMatrix[L][i]):10.6f} + {np.imag(self.ScatMatrix[L][i]):10.6f}I}}", end=" ") print("{", end="")
print() print(f"{L:2d},{L+i-self.S:4.1f}, {np.real(self.ScatMatrix[L][i]):10.6f} + {np.imag(self.ScatMatrix[L][i]):10.6f}I", end="")
if i < len(self.ScatMatrix[L])-1 :
print("}, ", end="")
else:
print("}", end="")
print("},")
def GetScatteringMatrix(self, L, J): def GetScatteringMatrix(self, L, J):
return self.ScatMatrix[L][J-L+self.S] return self.ScatMatrix[L][J-L+self.S]
@ -152,15 +162,15 @@ class DistortedWave(SolvingSE):
plt.show(block=False) plt.show(block=False)
input("Press Enter to continue...") input("Press Enter to continue...")
def RutherFord(self, theta): def RutherFord(self, theta_deg):
sin_half_theta = np.sin(theta / 2) sin_half_theta = np.sin(np.radians(theta_deg + 1e-20) / 2)
result = self.eta**2 / (4 * (self.k**2) * (sin_half_theta**4)) result = self.eta**2 / (4 * (self.k**2) * (sin_half_theta**4))
return result return result
def CoulombScatterintAmp(self, theta_deg): def CoulombScatterintAmp(self, theta_deg):
sin_sq = pow(np.sin(np.radians(theta_deg)/2), 2) sin_sq = pow(np.sin(np.radians(theta_deg + 1e-20)/2), 2)
coulPS = self.CoulombPhaseShift(0) coulPS = self.CoulombPhaseShift(0)
return - self.eta/ (2*self.k*sin_sq) * np.exp(2j*(coulPS - self.eta*np.log(sin_sq))) return - self.eta / (2 * self.k * sin_sq) * np.exp(1j * (2*coulPS - self.eta * np.log(sin_sq)))
def GMatrix1Spin(self, v, v0, l ) -> complex: def GMatrix1Spin(self, v, v0, l ) -> complex:
if self.S == 0 : if self.S == 0 :
@ -175,7 +185,21 @@ class DistortedWave(SolvingSE):
return value - KroneckerDelta(v, v0) return value - KroneckerDelta(v, v0)
def NuclearScatteringAmp(self, v, v0, theta, phi, maxL = None ) -> complex: def CalLegendre(self, theta_deg, maxL = None, maxM = None):
if maxL is None:
maxL = self.maxL
if maxM is None:
maxM = int(2*self.S)
self.legendreArray = associated_legendre_array(maxL, maxM, theta_deg)
return self.legendreArray
def GetPreCalLegendre(self, L, M):
if abs (M) <= int(2*self.S):
return self.legendreArray[L][int(abs(M))]
else :
return 0
def NuclearScatteringAmp(self, v, v0, maxL = None ) -> complex:
value = 0 value = 0
if maxL is None: if maxL is None:
maxL = self.maxL maxL = self.maxL
@ -184,21 +208,24 @@ class DistortedWave(SolvingSE):
value += 0 value += 0
else: else:
coulPS = self.CoulombPhaseShift(l) coulPS = self.CoulombPhaseShift(l)
value += np.sqrt(2*l+1) * sph_harm(v0 - v, l, phi, theta) * np.exp(2j * coulPS)* self.GMatrix1Spin(v, v0, l) fact = pow(-1, v0-v) * np.sqrt(factorial(l - abs(v0-v))/factorial(l + abs(v0-v)))
value += (2*l+1) * fact * self.GetPreCalLegendre(l, v0-v) * np.exp(2j * coulPS)* self.GMatrix1Spin(v, v0, l)
return value * np.sqrt(4*np.pi)/ 2 / 1j / self.k return value / 2j / self.k
def DCSUnpolarized(self, theta, phi, maxL = None): def DCSUnpolarized(self, theta_deg, maxL = None):
value = 0 value = 0
self.CalLegendre(theta_deg)
jaja = self.CoulombScatterintAmp(theta_deg)
for v in np.arange(-self.S, self.S + 1, 1): for v in np.arange(-self.S, self.S + 1, 1):
for v0 in np.arange(-self.S, self.S + 1, 1): for v0 in np.arange(-self.S, self.S + 1, 1):
value += abs(self.CoulombScatterintAmp(theta) * KroneckerDelta(v, v0) + self.NuclearScatteringAmp(v, v0, theta, phi, maxL))**2 value += abs( jaja * KroneckerDelta(v, v0) + self.NuclearScatteringAmp(v, v0, maxL))**2
value = value / (2 * self.S + 1) value = value / (2 * self.S + 1)
return value return value
def PlotDCSUnpolarized(self, thetaRange = 180, thetaStepDeg = 0.2, maxL = None): def PlotDCSUnpolarized(self, thetaRange = 180, thetaStepDeg = 0.2, maxL = None):
theta_values = np.linspace(0, thetaRange, int(thetaRange/thetaStepDeg)) theta_values = np.linspace(0, thetaRange, int(thetaRange/thetaStepDeg)+1)
thetaTick = 30 thetaTick = 30
if thetaRange < 180: if thetaRange < 180:
@ -206,22 +233,21 @@ class DistortedWave(SolvingSE):
y_values = [] y_values = []
for theta in theta_values: for theta in theta_values:
theta = np.deg2rad(theta)
if theta == 0: if theta == 0:
y_values.append(1) y_values.append(1)
else: else:
y_values.append(self.DCSUnpolarized(theta , 0, maxL)/ self.RutherFord(theta)) y_values.append(self.DCSUnpolarized(theta, maxL)/ self.RutherFord(theta))
print(f"{np.rad2deg(theta):6.2f}, {y_values[-1]:10.6f}") print(f"{theta:6.2f}, {y_values[-1]:10.6f}")
plt.figure(figsize=(8, 6)) plt.figure(figsize=(8, 6))
plt.plot(theta_values, y_values, marker='o', linestyle='-', color='blue', label='Real Part') # plt.plot(theta_values, y_values, marker='o', linestyle='-', color='blue')
plt.plot(theta_values, y_values, linestyle='-', color='blue')
plt.title("Differential Cross Section (Unpolarized)") plt.title("Differential Cross Section (Unpolarized)")
plt.xlabel("Theta (degrees)") plt.xlabel("Angle [deg]")
plt.ylabel("Value") plt.ylabel("D.C.S / Ruth")
plt.yscale("log") plt.yscale("log")
plt.xticks(np.arange(0, thetaRange + 1, thetaTick)) plt.xticks(np.arange(0, thetaRange + 1, thetaTick))
plt.xlim(0, thetaRange) plt.xlim(0, thetaRange)
plt.legend()
plt.grid() plt.grid()
plt.show(block=False) plt.show(block=False)
input("Press Enter to continue...") input("Press Enter to continue...")