(* :Title:      ss_spfn
*)

(* :Summary:
Spectral function of He-3
*)

(* :Author:	Sergey Kulagin
*)

(* :Package Version:  May 29, 2003
   :Updated:          Dec 08, 2006
*)

(* :Sources: *)

(* :Discussion:
This package reads the data file of Sch\"ultze-Sauer spectral function
[R. Sch\"ultze and P. Sauer, Phys. Rev. C48 (1993) 38]
and returns the functions characterizing different channels. These functions
are constructed as interpolation between the grid points given in the SS
data file.
C               The spectral function is based on the three-body    *
C               wave function calculated on a finite grid with      *
C               40*30 mesh points for the (p,q=p_N)-momenta. It is  *
C               again calculated on a finite grid with (30*30) mesh *
C               points for the momentum/energy mesh (p_N,E).        * 
*)

BeginPackage["ss`"]

prtcted = Unprotect[d0,d1,d2,c0t0,c1t0,c2t0,c0t1,c1t1,c2t1]

d0::usage=
"d0[p] the deuteron pole contribution to the spectral function f0."
d1::usage=
"d1[p] the deuteron pole contribution to the spectral function f1."
d2::usage=
"d2[p] the deuteron pole contribution to the spectral function f2."

c0t0::usage=
"c0t0[e,p] the contribution from the two-body continuum with isospin 0
to the spectral function f0."
c1t0::usage=
"c1t0[e,p] the contribution from the two-body continuum with isospin 0
to the spectral function f1."
c2t0::usage=
"c2t0[e,p] the contribution from the two-body continuum with isospin 0
to the spectral function f2."

c0t1::usage=
"c0t1[e,p] the contribution from the two-body continuum with isospin 1
to the spectral function f0."
c1t1::usage=
"c1t1[e,p] the contribution from the two-body continuum with isospin 1
to the spectral function f1."
c2t1::usage=
"c2t1[e,p] the contribution from the two-body continuum with isospin 1
to the spectral function f2."

pmin::usage=
"Minimal momentum in the momentum mesh."
pmax::usage=
"Maximum momentum in the momentum mesh."
emin0::usage=
"Minimum energy of the isospin 0 state in the energy mesh
(deuteron binding energy)."
emin1::usage=
"Minimum energy of the isospin 1 state in the energy mesh."
emax::usage=
"Maximum energy in the energy mesh."

(* Messages: *)

Begin["`Private`"]

Print["Setting up the SS spectral function of He-3..."];

(*
Read spectral function from data file. Note the following rules:

  The spectral function is given for the pair-isospin T=0 and T=1.
  The two-body breakup is given in the first columne, i.e. R0T0(1,J) etc.
  R0T0(2..NERG,J) contains the three-body breakup.
  The normalization is chosen such that the proton and  neutron
  contribution is given by:

      R0P = (F0T0 + F0T1)*4.*PI/3.
      R0N = 2.*F0T1*4.*PI/3.

  with R0P = 2/3
  and  R0N = 1/3

  Since the spectral function is given on a finite mesh the sumrule is
  not exactly fullfilled !
     The mesh was chosen for a particular kinematic situation.
     The momentum mesh is given for Gaussian points.
*)
(* ENERGIES ARE IN MEV AND MOMENTA IN FM^-1 *)

(* Options of numerical integration:
intopts = {PrecisionGoal -> 4, Method -> MultiDimensional, MaxRecursion -> 30};
*)
nmesh = 30;
r0t0 = Table[0, {i, 1, nmesh}, {j, 1, nmesh}];
r1t0 = r0t0; r2t0 = r0t0;
r0t1 = r0t0; r1t1 = r0t0; r2t1 = r0t0;

(* BEGIN reading the spectral function from data file *)
(* ener   is the energy mesh points in MeV.
   nerg   is the number of energy mesh points.
   qmesh  are the momentum mesh points in Fm^-1.
   nq     is the number of momentum mesh points.
   qweigh are the Gassian weights (not used in this code).
   r0t0, r1t0, r2t0 are the points of the corresponding spectral functions
                    with isospin 0.
   r0t1, r1t1, r2t1 are the points of the corresponding spectral functions
                    with isospin 1.
*)
rawdat = OpenRead["ss_spfn.dat"]
nerg = Read[rawdat, Number];
nq = Read[rawdat, Number];
ener = ReadList[rawdat, Number, nerg];
qmesh = ReadList[rawdat, Number, nq];
qweigh = ReadList[rawdat, Number, nq];
Do[(
      Do[ (
          r0t0[[i, j]] = Read[rawdat, Number];
          r1t0[[i, j]] = Read[rawdat, Number];
          r2t0[[i, j]] = Read[rawdat, Number]
          )
        , {j, 1, nq}])
    , {i, 1, nerg}];
Do[(
      Do[ (
          r0t1[[i, j]] = Read[rawdat, Number];
          r1t1[[i, j]] = Read[rawdat, Number];
          r2t1[[i, j]] = Read[rawdat, Number]
          )
        , {j, 1, nq}])
    , {i, 1, nerg}];
Close[rawdat]
(* END reading the spectral function from data file *)


(*
        Interpolation of spectral functions
Since the spectral functions change for a few orders of magnitude,
technically it is more convenient to interpolate the Log of the spectral
functions.
*)
(* Interpolation order *)
intord = 1;

(* Isospin 0 *)
(*
The deuteron in the final state (2-body breakup) corresponds to the 1st
row in r0t0, r1t0, and r2t0.
We treat this state separately.
*)

deut0 = Table[{qmesh[[j]], r0t0[[1, j]] // Log}, {j, 1, nq}];
deut1 = Table[{qmesh[[j]], r1t0[[1, j]] // Log}, {j, 1, nq}];
deut2 = Table[{qmesh[[j]], r2t0[[1, j]] // Log}, {j, 1, nq}];

logd0 = Interpolation[deut0, InterpolationOrder -> intord]
logd1 = Interpolation[deut1, InterpolationOrder -> intord]
logd2 = Interpolation[deut2, InterpolationOrder -> intord]

d0[q_]:=(logd0[q] //Exp) //Re
d1[q_]:=(logd1[q] //Exp) //Re
d2[q_]:=(logd2[q] //Exp) //Re

(* Two-nucleon continuum in the final state (3-body breakup): *)

con0t0 = Table[{ener[[i]], qmesh[[j]], r0t0[[i, j]] // Log},
	{i, 2, nerg}, {j,1, nq}];
con1t0 = Table[{ener[[i]], qmesh[[j]], r1t0[[i, j]] // Log},
	{i, 2, nerg}, {j,1, nq}];
con2t0 = Table[{ener[[i]], qmesh[[j]], r2t0[[i, j]] // Log},
	{i, 2, nerg}, {j,1, nq}];
con0t0 = Flatten[con0t0, 1];
con1t0 = Flatten[con1t0, 1];
con2t0 = Flatten[con2t0, 1];

logc0t0 = Interpolation[con0t0, InterpolationOrder -> intord]
logc1t0 = Interpolation[con1t0, InterpolationOrder -> intord]
logc2t0 = Interpolation[con2t0, InterpolationOrder -> intord]

c0t0[e_,q_]:=(logc0t0[e,q] //Exp) //Re
c1t0[e_,q_]:=(logc1t0[e,q] //Exp) //Re
c2t0[e_,q_]:=(logc2t0[e,q] //Exp) //Re

(* Isospin 1 *)
(*
The deuteron (two-body breakup) does not contribute in this channel.
Two-nucleon continuum in the final state (3-body breakup):
*)

con0t1 = Table[{ener[[i]], qmesh[[j]], r0t1[[i, j]] // Log},
	{i, 2, nerg}, {j, 1, nq}];
con1t1 = Table[{ener[[i]], qmesh[[j]], r1t1[[i, j]] // Log},
	{i, 2, nerg}, {j, 1, nq}];
con2t1 = Table[{ener[[i]], qmesh[[j]], r2t1[[i, j]] // Log},
	{i, 2, nerg}, {j, 1, nq}];
con0t1 = Flatten[con0t1, 1];
con1t1 = Flatten[con1t1, 1];
con2t1 = Flatten[con2t1, 1];

logc0t1 = Interpolation[con0t1, InterpolationOrder -> intord]
logc1t1 = Interpolation[con1t1, InterpolationOrder -> intord]
logc2t1 = Interpolation[con2t1, InterpolationOrder -> intord]

c0t1[e_,q_]:=(logc0t1[e,q] //Exp) //Re
c1t1[e_,q_]:=(logc1t1[e,q] //Exp) //Re
c2t1[e_,q_]:=(logc2t1[e,q] //Exp) //Re


(*
The proton and neutron momentum distributions mdp[p] and mdn[p].
These distributions are given by the energy integrals of the spectral
functions with definite isospin:

mdp[p] = 4*Pi*(md0[p]+md1[p])
mdn[p] = 4*Pi*(2*md1[p])

The distributions are normalized as

\int_pmin^pmax dp p^2 mdp(p) = 2
                      mdn(p) = 1
*)

(*
        The edge points of the momentum and energy mesh.
*)
pmin=First[qmesh];
pmax=Last[qmesh]; (* pmax=6.5; *)
emin0=ener[[1]];
emin1=ener[[2]];
emax=Last[ener];


End[ ]
Protect[Evaluate[prtcted]]
EndPackage[ ]

Print["If no error messages, the SS spectral function was successfully loaded."];

(* End of ss_spfn *)