From 257b80cae452beb583974649a956eaa144cb12ef Mon Sep 17 00:00:00 2001 From: "Ryan@Home" Date: Sun, 23 Feb 2025 20:53:50 -0500 Subject: [PATCH] scattering matrix very large for 11C(p,p).... --- Raphael/DWBA.py | 50 ++++++++++++++++++++++++++++++++++++ Raphael/DWBA_ZR.py | 39 ---------------------------- Raphael/README.md | 55 ++++++++++++++++++++++++++++++++++++++++ Raphael/distortedWave.py | 11 +++++--- 4 files changed, 112 insertions(+), 43 deletions(-) create mode 100755 Raphael/DWBA.py delete mode 100755 Raphael/DWBA_ZR.py create mode 100644 Raphael/README.md diff --git a/Raphael/DWBA.py b/Raphael/DWBA.py new file mode 100755 index 0000000..6bfe16a --- /dev/null +++ b/Raphael/DWBA.py @@ -0,0 +1,50 @@ +#!/usr/bin/env python3 + +import time +import matplotlib.pyplot as plt +from distortedWave import WoodsSaxonPot, CoulombPotential, SpinOrbit_Pot, WS_SurfacePot, DistortedWave +from dwba_zr import DWBA_ZR + +# haha = DWBA_ZR("16O", "d", "p", "17O", "1/2+", "1s1/2", 0.0, 10) +# haha = DWBA_ZR("16O", "d", "p", "17O", "5/2+", "0d5/2", 0.0, 10) + + +# haha.FindBoundState() +# haha.boundState.PlotBoundState() + +# haha.CalRadialIntegral() + +# haha.PlotScatteringMatrix(False) + +# haha.PlotIncomingDistortedWave(2, 1, 20) +# haha.PlotOutgoingDistortedWave(15, 15.5, 20) + +# haha.PrintRadialIntegral() +# haha.PlotRadialIntegral() +# haha.PlotRadialIntegralSigle(1, 1, -0.5) + + +# haha.CalAngDistribution(0, 180, 1) +# haha.PrintAngDist() +# haha.PlotAngDist() + + +####################################### Simple distorted wave calculation +kaka = DistortedWave("11C", "p", 60) +kaka.SetRange(0, 0.01, 2000) +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.500), 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) \ No newline at end of file diff --git a/Raphael/DWBA_ZR.py b/Raphael/DWBA_ZR.py deleted file mode 100755 index 194f369..0000000 --- a/Raphael/DWBA_ZR.py +++ /dev/null @@ -1,39 +0,0 @@ -#!/usr/bin/env python3 - -import time -import matplotlib.pyplot as plt -from dwba_zr import DWBA_ZR - -haha = DWBA_ZR("16O", "d", "p", "17O", "1/2+", "1s1/2", 0.0, 10) -# haha = DWBA_ZR("16O", "d", "p", "17O", "5/2+", "0d5/2", 0.0, 10) -haha.FindBoundState() - -# haha.boundState.PlotBoundState() - -haha.CalRadialIntegral() - -# haha.PlotScatteringMatrix(False) - -# haha.PlotIncomingDistortedWave(2, 1, 20) -# haha.PlotOutgoingDistortedWave(15, 15.5, 20) - -# haha.PrintRadialIntegral() -# haha.PlotRadialIntegral() -# haha.PlotRadialIntegralSigle(1, 1, -0.5) - - -haha.CalAngDistribution(0, 180, 1) -haha.PrintAngDist() -haha.PlotAngDist() - - -###################################### - -# print(haha.GetPreCalNineJ(1, 1, 3, 3.5)) - -# print(haha.Gamma(1, 1, 3, 3.5, 0, 1, 0.5)) - - - -# haha.PreCalLegendreP(5) -# haha.Beta(0, 1, 0.5) diff --git a/Raphael/README.md b/Raphael/README.md new file mode 100644 index 0000000..fe3d0e3 --- /dev/null +++ b/Raphael/README.md @@ -0,0 +1,55 @@ +# Background + +This is a python implementation of Distorted Wave Born Approximation calculation for direct nuclear reaction. The theory is stated on [here](https://nukephysik101.wordpress.com/2025/02/21/handbook-of-direct-nuclear-reaction-for-retarded-theorist-v1/) + +The motivation is two folds: 1) for my own curiousity, 2) All of DWBA code I know of is using Fortran and some of the code cannot be compiled anymore, and for those still can be compiled, the input card is a mess. So, i decided to have a mordern version of DWBA code. + +# Requirement +* numpy +* scipy - for gamma function, curve_fit, interp1d, simpson +* matplotlib +* mpmath - for coulombf, couombg and whittaker +* sympy - for winger_9j + +the version should not matter. + +# Start + +Open the DWBA.py, there are some examples + +# Code Components + +The foundation of the code base are +* assLegendreP.py - for associate Legendre polynomial for positive m +* clebschGordan.py - for custom build CG, which is faster +* opticalPotential.py - for optical potential, only have An & Cai for deuteron and Kronning for proton now +* ../Cleopatra/IAEANuclearData.py - for getting nuclear data like mass and spin-partiy + +Next, we have +* solveSE.py + +which have __PotentialFrom__ class, from this class, we have WoodsSaxon, WSSurface, Spin-Orbit, and Coulomb potential. Also, the __SolveSE__ class can get the mass from IAEA, calculated the Kinematics, and Solve the Radial wave function. + +based on the __SolveSE__ class, we have two classes +* boundState.py - for solving boudstate wave function +* distortedWave.py - for solving the distorted wave and extract the Scattering Matrix + +the __DistortedWave__ class can also calculated the elastics differential cross section in a ratio with RutherFord scattering. However, the code is not optimized and it may take minutes to have the angualr distribution. + +Finally, we have the +* dwba_zr.py + +which is a Zero-raneg DWBA class for calculating the radial integral and do the argular summation to obtain the differential cross section for transfer reaction. The code is kind of optimized. The optimization is done by precalculated the radial integral, the associated Legendre polaynomail, the Clebsch-Gordon coefficients, and the 9j-symbol. There may be some small tweets can be done, but it is python, no machines code level optimization. + +# Performanace + +for s-orbtial transfer, it takes 10 secs. for d-orbital, it takes like half minute. + +# Limitation & future work +* This is only for zero-range approximation, the angular distribution agree with experiment. +* Only for (d,p) or (p,d) (not and need test), should be work for all single nucleon transfer reaction. +* Need to work for inelastic scattering +* Need to implement the Finite-Range calcualtion +* Need to work on the 2-nucleons transfers +* Need to work on the Coupled-Channels +* Need to have a C++ version for speed sake (or use multiple cores, or both) diff --git a/Raphael/distortedWave.py b/Raphael/distortedWave.py index 4a9ff30..db74f20 100755 --- a/Raphael/distortedWave.py +++ b/Raphael/distortedWave.py @@ -6,7 +6,7 @@ from solveSE import WoodsSaxonPot, CoulombPotential, SpinOrbit_Pot, WS_SurfacePo from mpmath import coulombf, coulombg import numpy as np -from scipy.special import gamma, sph_harm, factorial +from scipy.special import gamma, factorial import matplotlib.pyplot as plt from assLegendreP import associated_legendre_array @@ -44,7 +44,7 @@ class DistortedWave(SolvingSE): eta = self.eta return np.angle(gamma(L+1+1j*eta)) - def CalScatteringMatrix(self, maxL = None, verbose = False): + def CalScatteringMatrix(self, normTo1 = False, maxL = None, verbose = False): start_time = time.time() # Start the timer if maxL is None: @@ -93,8 +93,11 @@ class DistortedWave(SolvingSE): temp_ScatMatrix.append(ScatMatrix) dwU = np.array(self.solU, dtype=np.complex128) - #dwU *= np.exp(-1j*sigma)/(B-A*1j) - dwU *= 1./(B-A*1j) + if normTo1 : + dwU /= self.maxSolU + else: + #dwU *= np.exp(-1j*sigma)/(B-A*1j) + dwU *= 1./(B-A*1j) temp_distortedWaveU.append(dwU) self.ScatMatrix.append(temp_ScatMatrix)