From db7943aa9a58c2ccfc9dbec390c2408c6f9375e6 Mon Sep 17 00:00:00 2001 From: "Ryan@SOLARIS_testStation" Date: Tue, 18 Feb 2025 18:46:09 -0500 Subject: [PATCH] add Rapheal, a python code for DWBA --- Raphael/assLegendreP.py | 37 ++++++ Raphael/coulombWave.py | 14 +++ Raphael/solveSE.py | 246 ++++++++++++++++++++++++++++++++++++++++ dwuck4/LGNDR.py | 66 ----------- 4 files changed, 297 insertions(+), 66 deletions(-) create mode 100755 Raphael/assLegendreP.py create mode 100755 Raphael/coulombWave.py create mode 100755 Raphael/solveSE.py delete mode 100755 dwuck4/LGNDR.py diff --git a/Raphael/assLegendreP.py b/Raphael/assLegendreP.py new file mode 100755 index 0000000..6da6a7e --- /dev/null +++ b/Raphael/assLegendreP.py @@ -0,0 +1,37 @@ +#!/usr/bin/env python3 + +import numpy as np + +def associated_legendre_polynomials_single_angle(L, M, theta_deg): + # Convert theta from degrees to radians + theta = np.radians(theta_deg) + x = np.cos(theta) + + P = np.zeros((L + 1, M + 1)) + + # P^m_l for m = 0, l = 0 + P[0, 0] = 1.0 + + # P^m_l for m = 0, l > 0 + for l in range(1, L + 1): + 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) + for m in range(1, M + 1): + # P^m_m + P[m, m] = (1 - 2*m) * P[m-1, m-1] + P[m, m] *= np.sqrt(1 - x**2) + + for l in range(m + 1, L + 1): + # 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) + + return P + +# Example usage +#L = 15 # Maximum l degree +#M = 5 # Maximum m order + +#legendre_polynomials = associated_legendre_polynomials_single_angle(L, M, 45) + +#print(legendre_polynomials) \ No newline at end of file diff --git a/Raphael/coulombWave.py b/Raphael/coulombWave.py new file mode 100755 index 0000000..bdc9090 --- /dev/null +++ b/Raphael/coulombWave.py @@ -0,0 +1,14 @@ +#!/usr/bin/env python3 + +import numpy as np +from mpmath import coulombf, coulombg + +L = 20 +eta = 2.0 +rho = 30 + +f_values = coulombf(L, eta, rho) +g_values = coulombg(L, eta, rho) + +print(f_values) +print(g_values) \ No newline at end of file diff --git a/Raphael/solveSE.py b/Raphael/solveSE.py new file mode 100755 index 0000000..2efd270 --- /dev/null +++ b/Raphael/solveSE.py @@ -0,0 +1,246 @@ +#!/usr/bin/env python3 + +import math +import numpy as np + +import sys, os + +sys.path.append(os.path.join(os.path.dirname(__file__), '../Cleopatra')) +from IAEANuclearData import IsotopeClass + +class Coulomb: + def __init__(self, rc): + self.rc = rc + self.id = 0 + self.ee = 1.43996 # MeV.fm + + def setA(self, A): + self.Rc = self.rc * math.pow(A, 1/3) + + def setAa(self, A, a): + self.Rc = self.rc * (math.pow(A, 1/3) + math.pow(a, 1/3)) + + def setCharge(self, Z): + self.Charge = Z + + def output(self, x): + if x >self.Rc: + return (self.Charge * self.ee) / (x + 1e-20) # Add a small value to avoid division by zero + else: + return (self.Charge * self.ee) / (2 * self.Rc) * (3 - (x / self.Rc)**2) + +class WS: + def __init__(self, V0, r0, a0) : + self.V0 = V0 + self.r0 = r0 + self.a0 = a0 + self.id = 1 + + def setA(self, A): + self.R0 = self.r0 * math.pow(A, 1/3) + + def setAa(self, A, a): + self.R0 = self.r0 * (math.pow(A, 1/3) + math.pow(a, 1/3)) + + def output(self, x): + return self.V0/(1 + math.exp((x-self.R0)/self.a0)) + +class SO: + def __init__(self, VSO, rSO, aSO) : + # the LS factor is put in the SolvingSE Class + self.VSO = VSO + self.rSO = rSO + self.aSO = aSO + self.id = 2 + + def setA(self, A): + self.RSO = self.rSO * math.pow(A, 1/3) + + def setAa(self, A, a): + self.RSO = self.rSO * (math.pow(A, 1/3) + math.pow(a, 1/3)) + + def output(self, x): + if x > 0 : + return 4*(self.VSO * math.exp((x-self.RSO)/self.aSO))/(self.aSO*math.pow(1+math.exp((x-self.RSO)/self.aSO),2))/x + else : + return 4*1e+19 + +class WSSurface: + def __init__(self, V0, r0, a0): + self.V0 = V0 + self.r0 = r0 + self.a0 = a0 + self.id = 3 + + def setA(self, A): + self.R0 = self.r0 * math.pow(A, 1/3) + + def setAa(self, A, a): + self.R0 = self.r0 * (math.pow(A, 1/3) + math.pow(a, 1/3)) + + def output(self, x): + exponent = (x - self.R0) / self.a0 + return self.V0 * math.exp(exponent) / (1 + math.exp(exponent))**2 + +#======================================== +class SolvingSE: + #grid setting + rStart = 0.0 + dr = 0.05 + nStep = 600*5 + rpos = np.arange(rStart, rStart+nStep*dr, dr) + SolU = [] + maxSolU = 0.0 + WF = np.empty([nStep], dtype=float) # radial wave function + maxWF = 0.0 + + #constant + mn = 939.56539 #MeV/c2 + amu = 931.494 #MeV/c2 + hbarc = 197.326979 #MeV.fm + ee = 1.43996 # MeV.fm + + #RK4 constants + parC = [0, 1./2, 1./2, 1.] + parD = [1./6, 2./6, 2./6, 1./6] + + #inital condition + solu0 = 0.0 + dsolu0 = 0.0001 + + potential_List = [] + + def PrintInput(self): + print(f" A : ({self.A:3d}, {self.ZA:3d})") + print(f" a : ({self.a:3d}, {self.Za:3d})") + print(f" Elab : {self.Energy : 10.3f} MeV") + print(f" mu : {self.mu: 10.3f} MeV/c2") +# print(f" Ecm : {self.Ecm: 10.3f} MeV") +# print(f" k : {self.k: 10.3f} MeV/c") +# print(f" eta : {self.eta: 10.3f}") +# print(f" L : {self.L}, maxL : {self.maxL}") + print(f" dr : {self.dr} fm, nStep : {self.nStep}") + print(f"rStart : {self.rStart} fm, rMax : {self.nStep * self.dr} fm") + print(f"spin-A : {self.sA}, spin-a : {self.sa} ") + + + def __init__(self, A, ZA, a, Za, Energy): + self.A = A + self.a = a + self.ZA = ZA + self.Za = Za + self.Z = ZA * Za + + self.sA = 0 + self.sa = 0 + + self.L = 0 + self.S = 0 + self.J = 0 + + haha = IsotopeClass() + self.mass_A = haha.GetMassFromAZ(self.A, self.ZA) + self.mass_a = haha.GetMassFromAZ(self.a, self.Za) + + #self.mu = (A * a)/(A + a) * self.amu + self.mu = (self.mass_A * self.mass_a)/(self.mass_A + self.mass_a) + + self.Energy = Energy + self.Ecm = self.Energy +# self.E_tot = math.sqrt(math.pow((a+A)*self.amu,2) + 2 * A * self.amu * eng_Lab) +# self.Ecm = self.E_tot - (a + A) * self.amu +# self.k = math.sqrt(self.mu * 2 * abs(self.Ecm)) / self.hbarc +# self.eta = self.Z * self.ee * math.sqrt( self.mu/2/self.Ecm ) / self.hbarc + +# self.maxL = int(self.k * (1.4 * (self.A**(1/3) + self.a**(1/3)) + 3)) + + def SetSpin(self, sA, sa): + self.sA = sA + self.sa = sa + self.S = self.sa + + def SetLJ(self, L, J): + self.L = L + self.J = J + + def LS(self, L = None, J = None) : + if L is None: + L = self.L + if J is None: + J = self.J + return (J*(J+1)-L*(L+1)-self.S*(self.S))/2. + + # set the range in fm + def SetRange(self, rStart, dr, nStep): + self.rStart = rStart + self.dr = dr + self.nStep = nStep + self.WF=np.empty([nStep], dtype=float) + self.rpos = np.arange(self.rStart, self.rStart+self.nStep*dr, self.dr) + self.WF = np.empty([self.nStep], dtype=float) + self.maxWF = 0.0 + self.SolU = [] + self.maxSolU = 0.0 + + def ClearPotential(self): + self.potential_List = [] + + def AddPotential(self, pot, useBothMass : bool = False): + if pot.id == 0: + pot.setCharge(self.Z) + if useBothMass: + pot.setAa(self.A, self.a) + else: + pot.setA(self.A) + self.potential_List.append(pot) + + def __PotentialValue(self, x): + value = 0 + for pot in self.potential_List: + if pot.id == 2: + value = value + self.LS() * pot.output(x) + else: + value = value + pot.output(x) + return value + + def GetPotentialValue(self, x): + return self.__PotentialValue(x) + + # The G-function, u''[r] = G[r, u[r], u'[r]] + def __G(self, x, y, dy): + #return -2*x*dy -2*y # solution of gaussian + if x > 0 : + return 2*self.mu/math.pow(self.hbarc,2)*(self.__PotentialValue(x) - self.Ecm)*y + self.L*(1+self.L)/x/x*y + else: + return 0 + + # Using Rungu-Kutta 4th method to solve u''[r] = G[r, u[r], u'[r]] + def SolveByRK4(self): + #initial condition + self.SolU = [self.solu0] + dSolU = [self.dsolu0] + + dyy = np.array([1., 0., 0., 0., 0.], dtype= complex) + dzz = np.array([1., 0., 0., 0., 0.], dtype= complex) + + self.maxSolU = 0.0 + + for i in range(self.nStep-1): + r = self.rStart + self.dr * i + y = self.SolU[i] + z = dSolU[i] + + for j in range(4): + dyy[j + 1] = self.dr * (z + self.parC[j] * dzz[j]) + dzz[j + 1] = self.dr * self.__G(r + self.parC[j] * self.dr, y + self.parC[j] * dyy[j], z + self.parC[j] * dzz[j]) + + 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) + dSolU.append(z + dz) + + return self.SolU + + def normalize_boundState(self): + pass diff --git a/dwuck4/LGNDR.py b/dwuck4/LGNDR.py deleted file mode 100755 index 6bb09f7..0000000 --- a/dwuck4/LGNDR.py +++ /dev/null @@ -1,66 +0,0 @@ -#!/usr/bin/env python3 - -import numpy as np -from scipy.special import lpmv - -def lgndr(mplus, lplus, thet): - """ - Calculates Legendre polynomials Plm - - Parameters: - mplus : int - Number of m's > 0 - lplus : int - Number of l's > 0 - thet : float - Angle in degrees - - Returns: - plm : list - List containing Legendre polynomials - """ - - theta = np.radians(thet) - y = np.cos(theta) - z = np.sin(theta) - plm = np.zeros(459, dtype=np.float64) - - ix = 0 - for m in range(1, mplus + 1): # For MPLUS = 1, LPLUS = 16 - lx = m - 1 # LX = 0 - l2 = 0 # L2 = 0 - p3 = 1.0 # P3 = 1.0 - fl1 = float(lx) # FL1 = 0 - - if lx != 0: - for lt in range(1, lx + 1): - fl1 += 1.0 - p3 *= fl1 * z / 2.0 - - p2 = 0.0 # P2 = 0.0 - fl2 = fl1 + 1.0 # FL2 = 1.0 - fl3 = 1.0 # FL3 = 1.0 - - for lt in range(1, lplus + 1): # Loop Lb - ix1 = ix + lt - - if l2 < lx: - plm[ix1] = 0.0 - else: - if l2 > lx: - p3 = (fl2 * y * p2 - fl1 * p1) / fl3 - fl1 += 1.0 - fl2 += 2.0 - fl3 += 1.0 - plm[ix1] = p3 - print(f'PLM, {lx:3d}, {l2:3d}, {ix1:3d}, {plm[ix1]:15.10f}') - p1, p2 = p2, p3 - - l2 += 1 - - ix += lplus - - return plm - - -plm = lgndr(3, 16, 1) \ No newline at end of file