(* :Title:	DeuteronWF *)

(* :Summary:
Parameterization of the deuteron wave function
calculated with the Bonn and Paris potentials
*)

(* :Author:	Sergey Kulagin *)

(* :Package Version: 	March 14, 2006 
*)

(* :Sources:
[1] Paris WF: Lacombe, Loiseau et al. Phys.Lett. 101B (1981) 139.
[2] Bonn WF:  R. Machleidt et al., Phys. Rept. 149, 1 (1987).
*)

(* :Discussion:
This package defines the s-wave and d-wave deutron wave functions in both
configuration, u(r) and w(r), and momentum, psi0(q) and psi2(q), spaces using
the parameterization discussed in [1,2]. For definition of the s- and d-wave
functions and normalization constants see any textbook on nuclear physics
(summary is provided in DeuteronWF.nb).
Normalization adopted here:
Integrate[u^2+w^2,{r,0,Infinity}]=1
Integrate[q^2 (psi0^2+psi2^2),{q,0,Infinity}]=1
*)

BeginPackage["Deuteron`"]

zaschita=Unprotect[ed,setdeuteron,psi]

setdeuteron::usage = "setdeuteron[type] defines the definitions for the Paris (type=1) and the Bonn (type=2) deuteron wave function"

psi::usage = "psi[space,l,r] is the s-wave (l=0) or d-wave (l=2) deuteron wave function in configuration space (space=1, in this case r is relative nucleon separation in Fm) or momentum space (space=0, in this case r is momentum in 1/Fm)."

ed::usage = "Deuteron binding energy in 1/Fm."

Begin["`Private`"]

hc=197.327053; (* MeV/Fm conversion constant *)
ed = -2.2246/hc; (* Deuteron binding energy in 1/Fm *)

setdeuteron[type_?IntegerQ] := (
(Which[				type==1,
Print["Setting up the Paris wave function..."]; 
(* Paris WF definitions	*)
(* Coefficients C and D of the deuteron WF parametrization *)
   c = {0.88688076, -0.34717093, -0.30502380*10.^1,
        0.56207766*10.^2, -0.74957334*10.^3,  0.53365279*10.^4,
       -0.22706863*10.^5,  0.60434469*10.^5, -0.10292058*10.^6,
        0.11223357*10.^6, -0.75925226*10.^5,  0.29059715*10.^5 };
               d = { 0.23135193*10.^-1, -0.85604572*10.^0, 0.56068193*10.^1,
		    -0.69462922*10.^2,  0.41631118*10.^3, -0.12546621*10.^4,
		     0.12387830*10.^4,  0.33739172*10.^4, -0.13041151*10.^5,
		     0.19512524*10.^5 };
n = 1+Length[c]; (* number of terms in parameterization *)
(* Mass parameters in 1/fm:
m0 = 1.0 ; alpha = 0.23162461;
*)
m = Table[0.23162461 + (j - 1), {j, 1, n}],
				type==2, 
Print["Setting up the Bonn wave function..."];
(* Bonn WF definitions *)
   c = { 0.88628672, -0.27591814, -0.11610727, -12.975243, 77.490155, 
      -272.98039, 534.02693, -563.28069, 302.14616, -64.920925 };
	       d = { 0.023237078, -0.52115578, -0.57197401, 
                     2.7570246,	-26.157324, 84.419883, -98.308997, 38.498490 };
n = 1+Length[c]; (* number of terms in parameterization *)
(* Mass parameters in 1/fm:
m0 = 0.9 ; alpha = 0.231607;
*)
m = Table[0.231607 + 0.9*(j - 1), {j, 1, n}]
 ]);
    m2 = m^2;
    c = Join[ c, {-Sum[c[[i]], {i, n - 1}]} ];
    sum1 = Sum[ d[[i]], {i, n - 3}];
    sum2 = d.Drop[1/m2, -3];
    sum3 = d.Drop[m2, -3];
    mlast3 = Take[m2, -3];
d = Join[ d,
{#1/(#3 - #1)/(#2 - #1)(-#2 #3 sum2 + (#2 + #3) sum1 - sum3) & @@  mlast3,
  #1/(#3 - #1)/(#2 - #1)(-#2 #3 sum2 + (#2 + #3) sum1 - sum3) & @@  RotateLeft[mlast3],
  #1/(#3 - #1)/(#2 - #1)(-#2 #3 sum2 + (#2 + #3) sum1 - sum3) & @@  RotateLeft[mlast3, 2]}
	];
Print["Done."]
);

psi[space_?IntegerQ, l_?IntegerQ,q_]:=(
   Which[space==0, (* Momentum space wave functions: *)
	           Which[l==0,Sqrt[2/Pi]*(c.(1/(q^2 + m2))),
	                 l==2,Sqrt[2/Pi]*(d.(1/(q^2 + m2)))],
	 space==1, (* Coordinate space wave functions: *)
	           r=q;
	           Which[l==0, c.Exp[-m r],
	                 l==2, d.(Exp[-m r] (1+3/(m r)+3/(m2 r^2))) ] 
        ]
);

End[]
Protect[Evaluate[zaschita]]
EndPackage[]

Print["If you see no error messages the package DeuteronWF.m was loaded successfully.\nThis package defines the deuteron wave function calculated with Paris or Bonn OBEP (energy-independent version) in both configuration and momentum spaces. For more information type ?setdeuteron, psi, ed."];

(* End *)