scattering matrix very large for 11C(p,p)....

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)

Raphael/DWBA_ZR.py

@@ -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)

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)

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,6 +93,9 @@ class DistortedWave(SolvingSE):
         temp_ScatMatrix.append(ScatMatrix)
 
         dwU = np.array(self.solU, dtype=np.complex128)
+        if normTo1 :
+          dwU /= self.maxSolU
+        else:
           #dwU *= np.exp(-1j*sigma)/(B-A*1j)
           dwU *= 1./(B-A*1j)
         temp_distortedWaveU.append(dwU)