diff --git a/Raphael/DWBA_ZR.py b/Raphael/DWBA_ZR.py index 2d86212..3584c16 100755 --- a/Raphael/DWBA_ZR.py +++ b/Raphael/DWBA_ZR.py @@ -1,12 +1,7 @@ #!/usr/bin/env python3 from boundState import BoundState - -from solveSE import WoodsSaxonPot, CoulombPotential, SpinOrbit_Pot, WS_SurfacePot, SolvingSE -from mpmath import coulombf, coulombg - -import numpy as np -from scipy.special import gamma +from solveSE import WoodsSaxonPot, CoulombPotential, SpinOrbit_Pot, WS_SurfacePot # 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) @@ -14,19 +9,9 @@ from scipy.special import gamma # # boundState.PrintWF() # boundState.PlotBoundState() -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 +from distortedWave import DistortedWave -def FivePointsSlope(data, n): - return ( data[n + 2] - 8 * data[n + 1] + 8 * data[n - 1] - data[n - 2] ) / 12 - - -dw = SolvingSE("60Ni", "p", 30) -dw.SetRange(0, 0.1, 300) -dw.SetLJ(0, 0.5) -dw.dsolu0 = 1 -dw.CalCMConstants() -dw.PrintInput() +dw = DistortedWave("60Ni", "p", 30) dw.ClearPotential() dw.AddPotential(WoodsSaxonPot(-47.937-2.853j, 1.20, 0.669), False) @@ -34,47 +19,10 @@ 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) -rpos = dw.rpos -solU = dw.SolveByRK4() +dw.CalScatteringMatrix() -solU /= dw.maxSolU +dw.PrintScatteringMatrix() -# for r, u in zip(rpos, solU): -# print(f"{r:.3f} {np.real(u):.6f} {np.imag(u):.6f}") +# dw.PlotScatteringMatrix() -def CoulombPhaseShift(L, eta): - return np.angle(gamma(L+1+1j*eta)) - -sigma = CoulombPhaseShift(dw.L, dw.eta) -# find pahse shift by using the asymptotic behavior of the wave function -r1 = rpos[-2] -f1 = float(coulombf(dw.L, dw.eta, dw.k*r1)) -g1 = float(coulombg(dw.L, dw.eta, dw.k*r1)) -u1 = solU[-2] - -r2 = rpos[-1] -f2 = float(coulombf(dw.L, dw.eta, dw.k*r2)) -g2 = float(coulombg(dw.L, dw.eta, dw.k*r2)) -u2 = solU[-1] - -det = f2*g1 - f1*g2 -A = (f2*u1 - u2*f1) / det -B = (u2*g1 - g2*u1) / det - -print(f"A = {np.real(A):.6f} + {np.imag(A):.6f} I") -print(f"B = {np.real(B):.6f} + {np.imag(B):.6f} I") - -ScatMatrix = (B + A * 1j)/(B - A * 1j) -print(f"Scat Matrix = {np.real(ScatMatrix):.6f} + {np.imag(ScatMatrix):.6f} I") - -solU = np.array(solU, dtype=np.complex128) -solU *= np.exp(1j * sigma)/(B-A*1j) - -from matplotlib import pyplot as plt - -plt.plot(rpos, np.real(solU), label="Real") -plt.plot(rpos, np.imag(solU), label="Imaginary") -plt.legend() -plt.show(block=False) - -input("Press Enter to continue...") \ No newline at end of file +dw.PlotDCSUnpolarized() \ No newline at end of file diff --git a/Raphael/coulombWave.py b/Raphael/coulombWave.py index bdc9090..11507bd 100755 --- a/Raphael/coulombWave.py +++ b/Raphael/coulombWave.py @@ -1,14 +1,56 @@ #!/usr/bin/env python3 import numpy as np -from mpmath import coulombf, coulombg +from scipy.special import gamma -L = 20 -eta = 2.0 -rho = 30 +def pochhammer(x, n): + """Compute Pochhammer symbol (x)_n""" + if n == 0: + return 1.0 + result = 1.0 + for i in range(n): + result *= (x + i) + return result -f_values = coulombf(L, eta, rho) -g_values = coulombg(L, eta, rho) +def hyp1f1_series(a, b, z, max_terms=1000, tol=1e-10): + """ + Compute _1F_1(a, b, z) using series expansion. + + :param a, b: Parameters of the hypergeometric function + :param z: Complex argument + :param max_terms: Maximum number of terms to sum + :param tol: Tolerance for convergence + :return: Approximation of _1F_1(a, b, z) + """ + sum = 0.0 + 0.0j + last_term = 0.0 + 0.0j + for n in range(max_terms): + term = pochhammer(a, n) / pochhammer(b, n) * z**n / np.math.factorial(n) + sum += term + if abs(term - last_term) < tol * abs(sum): # Check for convergence + return sum + last_term = term + return sum # If we reach here, we've used all terms without converging to desired tolerance -print(f_values) -print(g_values) \ No newline at end of file +def coulomb_wave_function(L, eta, rho): + """ + Compute the regular Coulomb wave function F_L(eta, rho). + + :param L: Angular momentum quantum number + :param eta: Sommerfeld parameter + :param rho: Radial coordinate scaled by k (wavenumber) + :return: F_L(eta, rho) + """ + # Compute normalization constant C_L(eta) + C_L = (2**L * np.exp(-np.pi * eta / 2) * np.abs(gamma(L + 1 + 1j * eta))) / gamma(2*L + 2) + + # Compute _1F_1 using our series expansion + hyp1f1_value = hyp1f1_series(L + 1 - 1j * eta, 2*L + 2, 2j * rho) + + # Return the Coulomb wave function + return C_L * (rho ** (L+1)) * np.exp(-1j * rho) * hyp1f1_value + +# Example usage +#L, eta, rho = 1, 0.5, 1.0 +#result = coulomb_wave_function(L, eta, rho) +#print(f"F_{L}({eta}, {rho}) = {result}") \ No newline at end of file diff --git a/Raphael/distortedWave.py b/Raphael/distortedWave.py new file mode 100755 index 0000000..c90bd4b --- /dev/null +++ b/Raphael/distortedWave.py @@ -0,0 +1,227 @@ +#!/usr/bin/env python3 + +from boundState import BoundState + +from solveSE import WoodsSaxonPot, CoulombPotential, SpinOrbit_Pot, WS_SurfacePot, SolvingSE +from mpmath import coulombf, coulombg + +import numpy as np +from scipy.special import gamma, sph_harm, factorial +import matplotlib.pyplot as plt + + +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 + +def FivePointsSlope(data, n): + return ( data[n + 2] - 8 * data[n + 1] + 8 * data[n - 1] - data[n - 2] ) / 12 + +from sympy.physics.quantum.cg import CG +from sympy import S +def clebsch_gordan(j1, m1, j2, m2, j, m): + cg = CG(S(j1), S(m1), S(j2), S(m2), S(j), S(m)) + result = cg.doit() + return np.complex128(result) + +def KroneckerDelta(i, j): + if i == j: + return 1 + else: + return 0 + +class DistortedWave(SolvingSE): + def __init__(self, target, projectile, ELab): + super().__init__(target, projectile, ELab) + self.SetRange(0, 0.1, 300) + self.CalCMConstants() + + self.ScatMatrix = [] + self.distortedWaveU = [] + + def SetLJ(self, L, J): + self.L = L + self.J = J + self.dsolu0 = pow(0.1, 2*L+1) + + def CoulombPhaseShift(self, L = None, eta = None): + if L is None: + L = self.L + if eta is None: + eta = self.eta + return np.angle(gamma(L+1+1j*eta)) + + def CalScatteringMatrix(self, maxL = None, verbose = False): + if maxL is None: + maxL = self.maxL + + self.ScatMatrix = [] + self.distortedWaveU = [] + + for L in range(0, maxL+1): + sigma = self.CoulombPhaseShift() + + temp_ScatMatrix = [] + temp_distortedWaveU = [] + + for J in np.arange(L-self.S, L + self.S+1, 1): + if J < 0: + temp_ScatMatrix.append(0) + temp_distortedWaveU.append(0) + continue + + self.SetLJ(L, J) + self.SolveByRK4() + + r1 = self.rpos[-2] + f1 = float(coulombf(self.L, self.eta, self.k*r1)) + g1 = float(coulombg(self.L, self.eta, self.k*r1)) + u1 = self.solU[-2] + + r2 = self.rpos[-1] + f2 = float(coulombf(self.L, self.eta, self.k*r2)) + g2 = float(coulombg(self.L, self.eta, self.k*r2)) + u2 = self.solU[-1] + + det = f2*g1 - f1*g2 + + A = (f2*u1 - u2*f1) / det + B = (u2*g1 - g2*u1) / det + + ScatMatrix = (B + A * 1j)/(B - A * 1j) + + if verbose: + print(f"{{{L},{J}, {np.real(ScatMatrix):10.6f} + {np.imag(ScatMatrix):10.6f}I}}") + + temp_ScatMatrix.append(ScatMatrix) + + dwU = np.array(self.solU, dtype=np.complex128) + dwU *= np.exp(-1j*sigma)/(B-A*1j) + temp_distortedWaveU.append(dwU) + + self.ScatMatrix.append(temp_ScatMatrix) + self.distortedWaveU.append(temp_distortedWaveU) + + return [self.ScatMatrix, self.distortedWaveU] + + def PrintScatteringMatrix(self): + for L in range(0, len(self.ScatMatrix)): + 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() + + def GetScatteringMatrix(self, L, J): + return self.ScatMatrix[L][J-L+self.S] + + def GetDistortedWave(self, L, J): + return self.distortedWaveU[L][J-L+self.S] + + def PlotDistortedWave(self, L, J): + plt.plot(self.rpos, np.real(self.GetDistortedWave(L, J)), label="Real") + plt.plot(self.rpos, np.imag(self.GetDistortedWave(L, J)), label="Imaginary") + plt.legend() + plt.grid() + plt.show(block=False) + input("Press Enter to continue...") + + def PlotScatteringMatrix(self): + nSpin = int(self.S*2+1) + fig, axes = plt.subplots(1, nSpin, figsize=(6*nSpin, 4)) + + for i in range(0, nSpin): + sm = [] + l_list = [] + for L in range(0, len(self.ScatMatrix)): + if i == 0 and L == 0 : + continue + l_list.append(L) + sm.append(self.ScatMatrix[L][i]) + + axes[i].plot(l_list, np.real(sm), label="Real", marker='o') + axes[i].plot(l_list, np.imag(sm), label="Imaginary", marker='x') + axes[i].legend() + axes[i].set_xlabel('L') + axes[i].set_ylabel('Value') + if self.S*2 % 2 == 0 : + str = f'{int(i-self.S):+d}' + else: + str = f'{int(2*(i-self.S)):+d}/2' + axes[i].set_title(f'Real and Imaginary Parts vs L for Spin J = L{str}') + axes[i].grid() + + plt.tight_layout() + plt.show(block=False) + input("Press Enter to continue...") + + def RutherFord(self, theta): + sin_half_theta = np.sin(theta / 2) + result = self.eta**2 / (4 * (self.k**2) * (sin_half_theta**4)) + return result + + def CoulombScatterintAmp(self, theta_deg): + sin_sq = pow(np.sin(np.radians(theta_deg)/2), 2) + coulPS = self.CoulombPhaseShift(0) + return - self.eta/ (2*self.k*sin_sq) * np.exp(2j*(coulPS - self.eta*np.log(sin_sq))) + + def GMatrix1Spin(self, v, v0, l ) -> complex: + if self.S == 0 : + return self.ScatMatrix[l][0] - KroneckerDelta(v, v0) + else: + Jmin = l - self.S + Jmax = l + self.S + value = 0 + for J in np.arange(Jmin, Jmax + 1, 1): + index = int(J - Jmin) + value += clebsch_gordan(l, 0, self.S, v0, J, v0) * clebsch_gordan(l, v0-v, self.S, v, J, v0) * self.ScatMatrix[l][index] + + return value - KroneckerDelta(v, v0) + + def NuclearScatteringAmp(self, v, v0, theta, phi, maxL = None ) -> complex: + value = 0 + if maxL is None: + maxL = self.maxL + for l in range(0, maxL+1): + if abs(v0-v) > l : + value += 0 + else: + 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) + + return value * np.sqrt(4*np.pi)/ 2 / 1j / self.k + + def DCSUnpolarized(self, theta, phi, maxL = None): + value = 0 + for v 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 = value / (2 * self.S + 1) + return value + + def PlotDCSUnpolarized(self, thetaRange = 180, thetaStepDeg = 0.2, maxL = None): + theta_values = np.linspace(0, thetaRange, int(thetaRange/thetaStepDeg)) + + thetaTick = 30 + if thetaRange < 180: + thetaTick = 10 + + y_values = [] + for theta in theta_values: + theta = np.deg2rad(theta) + if theta == 0: + y_values.append(1) + else: + y_values.append(self.DCSUnpolarized(theta , 0, maxL)/ self.RutherFord(theta)) + print(f"{np.rad2deg(theta):6.2f}, {y_values[-1]:10.6f}") + + plt.figure(figsize=(8, 6)) + plt.plot(theta_values, y_values, marker='o', linestyle='-', color='blue', label='Real Part') + plt.title("Differential Cross Section (Unpolarized)") + plt.xlabel("Theta (degrees)") + plt.ylabel("Value") + plt.yscale("log") + plt.xticks(np.arange(0, thetaRange + 1, thetaTick)) + plt.xlim(0, thetaRange) + plt.legend() + plt.grid() + plt.show(block=False) + input("Press Enter to continue...") \ No newline at end of file diff --git a/Raphael/solveSE.py b/Raphael/solveSE.py index cd943ab..f9c86c1 100755 --- a/Raphael/solveSE.py +++ b/Raphael/solveSE.py @@ -89,7 +89,7 @@ class SolvingSE: dr = 0.05 nStep = 600*5 rpos = np.arange(rStart, rStart+nStep*dr, dr) - SolU = [] # raidal wave function + solU = [] # raidal wave function maxSolU = 0.0 #constant @@ -199,7 +199,7 @@ class SolvingSE: self.dr = dr self.nStep = nStep self.rpos = np.arange(self.rStart, self.rStart+self.nStep*dr, self.dr) - self.SolU = [] + self.solU = [] self.maxSolU = 0.0 def ClearPotential(self): @@ -237,7 +237,7 @@ class SolvingSE: # Using Rungu-Kutta 4th method to solve u''[r] = G[r, u[r], u'[r]] def SolveByRK4(self): #initial condition - self.SolU = [self.solu0] + self.solU = [self.solu0] dSolU = [self.dsolu0] dyy = np.array([1., 0., 0., 0., 0.], dtype= complex) @@ -247,7 +247,7 @@ class SolvingSE: for i in range(self.nStep-1): r = self.rStart + self.dr * i - y = self.SolU[i] + y = self.solU[i] z = dSolU[i] for j in range(4): @@ -257,13 +257,13 @@ class SolvingSE: dy = sum(self.parD[j] * dyy[j + 1] for j in range(4)) dz = sum(self.parD[j] * dzz[j + 1] for j in range(4)) - self.SolU.append(y + dy) + self.solU.append(y + dy) dSolU.append(z + dz) - if np.real(self.SolU[-1]) > self.maxSolU: - self.maxSolU = abs(self.SolU[-1]) + if np.real(self.solU[-1]) > self.maxSolU: + self.maxSolU = abs(self.solU[-1]) - return self.SolU + return self.solU def NearestPosIndex(self, r): return min(len(self.rpos)-1, int((r - self.rStart) / self.dr))