diff --git a/Raphael/DWBA.py b/Raphael/DWBA.py index a71d66a..348675f 100755 --- a/Raphael/DWBA.py +++ b/Raphael/DWBA.py @@ -30,38 +30,64 @@ from dwba_zr import DWBA_ZR ####################################### Simple distorted wave calculation -kaka = DistortedWave("11C", "p", 60) -kaka.SetRange(0, 0.1, 1000) -kaka.maxL = 14 +''' +kaka = DistortedWave("148Sm", "a", 50) kaka.PrintInput() kaka.ClearPotential() -kaka.AddPotential(WoodsSaxonPot(-34.714-6.749j, 1.122, 0.676), False) # False = only use 11C for radius calculation -kaka.AddPotential(WS_SurfacePot( -3.194j, 1.307, 0.524), False) -kaka.AddPotential(SpinOrbit_Pot( -4.532+0.477j, 0.894, 0.590), False) -kaka.AddPotential(CoulombPotential(1.578), False) -kaka.PrintPotentials() - -kaka.CalScatteringMatrix(True) -kaka.PrintScatteringMatrix() -kaka.PlotScatteringMatrix() - -# kaka.PlotDistortedWave(1, 1.5) - -# kaka.PlotDCSUnpolarized(180, 1, None, True) - -exit() - -kaka = DistortedWave("60Ni", "p", 30) -kaka.PrintInput() -kaka.ClearPotential() -kaka.AddPotential(WoodsSaxonPot(-47.937-2.853j, 1.120, 0.669), False) # False = only use 11C for radius calculation -kaka.AddPotential(WS_SurfacePot( -6.878j, 1.280, 0.550), False) -kaka.AddPotential(SpinOrbit_Pot( -5.250+0.162j, 1.020, 0.590), False) -kaka.AddPotential(CoulombPotential(1.258), False) +kaka.AddPotential(WoodsSaxonPot(-65.500 -29.800j, 1.427, 0.671), False) # False = only use 11C for radius calculation +kaka.AddPotential(CoulombPotential(1.4), False) kaka.PrintPotentials() kaka.CalScatteringMatrix() kaka.PrintScatteringMatrix() # kaka.PlotScatteringMatrix() -kaka.PlotDCSUnpolarized(180, 1, None, True) \ No newline at end of file +kaka.CalAngDistribution(180, 0.5, None, False) +kaka.PlotDCSUnpolarized() +''' + +''' +kaka = DistortedWave("60Ni", "p", 30) +# kaka.SetRange(0, 0.02, 1000) +kaka.PrintInput() +kaka.ClearPotential() +kaka.AddPotential(WoodsSaxonPot(-47.937 -2.853j, 1.200, 0.669), False) # False = only use 11C for radius calculation +kaka.AddPotential(WS_SurfacePot( -6.878j, 1.280, 0.550), False) +kaka.AddPotential(SpinOrbit_Pot( -5.250 +0.162j, 1.020, 0.590), False) +kaka.AddPotential(CoulombPotential(1.258), False) +kaka.PrintPotentials() + +# kaka.PlotPotential(0, 0.5, 10) + +kaka.CalScatteringMatrix(True) +kaka.PrintScatteringMatrix() +kaka.PlotScatteringMatrix() +# kaka.PlotDistortedWave(1, 1.5) + +kaka.CalAngDistribution(180, 0.5, None, False) +kaka.PlotDCSUnpolarized() +''' + +kaka = DistortedWave("11C", "p", 60) +# kaka.SetRange(0, 0.1, 1000) +kaka.maxL = 14 +kaka.PrintInput() +kaka.ClearPotential() +kaka.AddPotential(WoodsSaxonPot(-34.714 -6.749j, 1.122, 0.676), False) # False = only use 11C for radius calculation +kaka.AddPotential(WS_SurfacePot( -3.194j, 1.307, 0.524), False) +kaka.AddPotential(SpinOrbit_Pot( -4.532 +0.477j, 0.894, 0.590), False) +kaka.AddPotential(CoulombPotential(1.578), False) +kaka.PrintPotentials() + +# kaka.PlotPotential(10) + +kaka.CalScatteringMatrix(True) +kaka.PrintScatteringMatrix() +# kaka.PlotScatteringMatrix() + +# kaka.PlotDistortedWave(1, 1.5) +kaka.CalAngDistribution(180, 0.5, None, False) +kaka.PlotDCSUnpolarized() + + + diff --git a/Raphael/distortedWave.py b/Raphael/distortedWave.py index b833a73..849ed85 100755 --- a/Raphael/distortedWave.py +++ b/Raphael/distortedWave.py @@ -43,6 +43,13 @@ class DistortedWave(SolvingSE): if eta is None: eta = self.eta return np.angle(gamma(L+1+1j*eta)) + + def CoulombHankel(self, sign, rho, L = None, eta = None) -> complex: + if L is None: + L = self.L + if eta is None: + eta = self.eta + return float(coulombg(L, eta, rho)) + sign * 1j * float(coulombf(L, eta, rho)) def CalScatteringMatrix(self, normTo1 = False, maxL = None, verbose = False): start_time = time.time() # Start the timer @@ -56,7 +63,6 @@ class DistortedWave(SolvingSE): tempZeroList = np.zeros(len(self.rpos), dtype=np.complex128) for L in range(0, maxL+1): - # sigma = self.CoulombPhaseShift() temp_ScatMatrix = [] temp_distortedWaveU = [] @@ -70,22 +76,67 @@ class DistortedWave(SolvingSE): self.SetLJ(L, J) self.SolveByRK4() + #============ 2-point using Hankel + # r1 = self.rpos[-2] + # u1 = self.solU[-2] + # Hp1 = self.CoulombHankel( 1, self.k * r1) + # Hm1 = self.CoulombHankel(-1, self.k * r1) + + # r2 = self.rpos[-1] + # u2 = self.solU[-1] + # Hp2 = self.CoulombHankel( 1, self.k * r2) + # Hm2 = self.CoulombHankel(-1, self.k * r2) + + # ScatMatrix = (u1 * Hm2 - Hm1 * u2) / (u1 * Hp2 - Hp1 * u2) + # norm = 2j * (u1 * Hp2 - Hp1 * u2 ) / (Hp1 * Hm2 - Hm1 * Hp2) + + #=========== 2 point G, F, this is fastest r1 = self.rpos[-2] + u1 = self.solU[-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] + u2 = self.solU[-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] + print(f"{L:2d}, {J:4.1f} | {r1:.3f}, {f1:10.6f}, {g1:10.6f} | {r2:.3f}, {f2:10.6f}, {g2:10.6f}") + det = f2*g1 - f1*g2 - A = (f2*u1 - u2*f1) / det B = (u2*g1 - g2*u1) / det - ScatMatrix = (B + A * 1j)/(B - A * 1j) + norm = B - A * 1j + + #============ 5 points slopes + # r_list = self.rpos[-5:] + # u_list = self.solU[-5:] + # f_list = [float(coulombf(self.L, self.eta, self.k * r)) for r in r_list] + # g_list = [float(coulombg(self.L, self.eta, self.k * r)) for r in r_list] + # u0 = u_list[2] + # f0 = f_list[2] + # g0 = g_list[2] + # du = FivePointsSlope(u_list, 2) + # df = FivePointsSlope(f_list, 2) + # dg = FivePointsSlope(g_list, 2) + + # det = df*g0 - dg*f0 + # A = (df*u0 - du*f0) /det + # B = (du*g0 - dg*u0) /det + # ScatMatrix = (B + A * 1j)/(B - A * 1j) + # norm = B - A * 1j + + # #============ 100 points fitting + # r_list = self.rpos[-5:] + # u_list = self.solU[-5:] + # f_list = [float(coulombf(self.L, self.eta, self.k * r)) for r in r_list] + # g_list = [float(coulombg(self.L, self.eta, self.k * r)) for r in r_list] + + # # Fit u_list to A * g_list + B * f_list + # A, B = np.linalg.lstsq(np.vstack([g_list, f_list]).T, u_list, rcond=None)[0] + # ScatMatrix = (B + A * 1j) / (B - A * 1j) + # norm = B - A * 1j if verbose: print(f"{{{L},{J}, {np.real(ScatMatrix):10.6f} + {np.imag(ScatMatrix):10.6f}I}}") @@ -96,8 +147,7 @@ class DistortedWave(SolvingSE): if normTo1 : dwU /= self.maxSolU else: - #dwU *= np.exp(-1j*sigma)/(B-A*1j) - dwU *= 1./(B-A*1j) + dwU *= 1./norm temp_distortedWaveU.append(dwU) self.ScatMatrix.append(temp_ScatMatrix) @@ -239,25 +289,31 @@ class DistortedWave(SolvingSE): value = value / (2 * self.S + 1) return value - def PlotDCSUnpolarized(self, thetaRange = 180, thetaStepDeg = 0.2, maxL = None, verbose = False): - theta_values = np.linspace(0, thetaRange, int(thetaRange/thetaStepDeg)+1) + def CalAngDistribution(self, thetaRange = 180, thetaStepDeg = 0.2, maxL = None, verbose = False): + self.theta_list = np.linspace(0, thetaRange, int(thetaRange/thetaStepDeg)+1) + self.angDist = [] + for theta in self.theta_list: + if theta == 0: + self.angDist.append(1) + else: + self.angDist.append(self.DCSUnpolarized(theta, maxL)/ self.RutherFord(theta)) + if verbose : + print(f"{theta:6.2f}, {self.angDist[-1]:10.6f}") + def PrintAngDistribution(self): + print("======================= Angular Distribution") + for theta, dist in zip(self.theta_list, self.angDist): + print(f"{theta:5.1f}, {dist:10.3f}") + + def PlotDCSUnpolarized(self): thetaTick = 30 + thetaRange = max(self.theta_list) if thetaRange < 180: thetaTick = 10 - y_values = [] - for theta in theta_values: - if theta == 0: - y_values.append(1) - else: - y_values.append(self.DCSUnpolarized(theta, maxL)/ self.RutherFord(theta)) - if verbose : - print(f"{theta:6.2f}, {y_values[-1]:10.6f}") - plt.figure(figsize=(8, 6)) # plt.plot(theta_values, y_values, marker='o', linestyle='-', color='blue') - plt.plot(theta_values, y_values, linestyle='-', color='blue') + plt.plot(self.theta_list, self.angDist, linestyle='-', color='blue') plt.title("Differential Cross Section (Unpolarized)") plt.xlabel("Angle [deg]") plt.ylabel("D.C.S / Ruth") diff --git a/Raphael/solveSE.py b/Raphael/solveSE.py index efee108..b8726f1 100755 --- a/Raphael/solveSE.py +++ b/Raphael/solveSE.py @@ -4,6 +4,7 @@ import math import numpy as np import re import sys, os +import matplotlib.pyplot as plt sys.path.append(os.path.join(os.path.dirname(__file__), '../Cleopatra')) from IAEANuclearData import IsotopeClass @@ -22,7 +23,7 @@ class PotentialForm: # 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)) + self.R0= self.r0 * (A**(1./3.) +a**(1./3.)) def output(self, x): return 0 @@ -58,7 +59,7 @@ class WoodsSaxonPot(PotentialForm): self.id = 1 def output(self, x): - if self.V0 == 0.0: + if abs(self.V0) == 0.0: return 0 else: return self.V0/(1 + math.exp((x-self.R0)/self.a0)) @@ -75,13 +76,15 @@ class SpinOrbit_Pot(PotentialForm): self.id = 2 def output(self, x): - if self.V0 == 0.0 : + if abs(self.V0) == 0.0 : return 0 else: if x > 0 : - return 4*(self.V0 * math.exp((x-self.R0)/self.a0))/(self.a0*math.pow(1+math.exp((x-self.R0)/self.a0),2))/x + exponent = (x - self.R0) / self.a0 + return 4*(self.V0 * math.exp(exponent))/(self.a0*math.pow(1+math.exp(exponent),2))/x else : - return 4*1e+19 + # return 4*1e+19 + return 0 def printPot(self): return super().printPot("Spin-Orbit") @@ -94,7 +97,7 @@ class WS_SurfacePot(PotentialForm): self.id = 3 def output(self, x): - if self.V0 == 0 : + if abs(self.V0) == 0 : return 0 else: exponent = (x - self.R0) / self.a0 @@ -158,6 +161,8 @@ class SolvingSE: self.mu = (self.mass_A * self.mass_a)/(self.mass_A + self.mass_a) self.Ecm = 0.0 + self.Factor = 2*self.mu/math.pow(self.hbarc,2) + def ConstructUsingAZ(self, A, ZA, a, Za, ELabPerA): # print(f"ConstructUsingAZ : {A}, {ZA}, {a}, {Za}, {ELabPerA:.3f}") self.A_A = A @@ -183,7 +188,6 @@ class SolvingSE: self.Z = self.Z_A * self.Z_a self.Energy = ELabPerA - def CalCMConstants(self, useELabAsEcm = False): if useELabAsEcm: self.E_tot = self.Energy # total energy in CM @@ -240,10 +244,17 @@ class SolvingSE: value = 0 for pot in self.potential_List: if isinstance(pot, PotentialForm): - if pot.id == 2 and self.L > 0: - value = value + self.LS() * pot.output(x) + if pot.id == 2 : + if self.L > 0: + value = value + self.LS() * pot.output(x) else: value = value + pot.output(x) + + if x > 1e-6: + value = self.Factor*(value) + self.L*(1+self.L)/x/x + else: + value = 0 + return value def PrintPotentials(self): @@ -251,14 +262,38 @@ class SolvingSE: if isinstance(pot, PotentialForm): pot.printPot() + def PlotPotential(self, L, J, maxR = None): + self.SetLJ(L, J) + pot = [] + e_list = [] + for r in self.rpos: + pot.append(self.__PotentialValue(r)) + + pot_real = [np.real(self.__PotentialValue(r)) for r in self.rpos] + pot_imag = [np.imag(self.__PotentialValue(r)) for r in self.rpos] + + plt.figure() + plt.plot(self.rpos, pot_real, label='Real Part') + plt.plot(self.rpos, pot_imag, label='Imaginary Part') + plt.xlabel('r (fm)') + plt.ylabel('Potential (MeV)') + plt.title('Potential vs. r') + if maxR != None: + plt.xlim(-1, maxR) + plt.ylim([min(pot_real)*1.1, abs(min(pot_real)*1.1)]) + plt.legend() + plt.grid(True) + plt.show(block=False) + input("press anykey to continous....") + 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 + if x > 1e-6 : + return self.__PotentialValue(x)*y - self.Factor*self.Ecm*y else: return 0