(* :Title:      kpsv_spfn
*)

(* :Summary:
Spectral function of He-3 by Salme et al
*)

(* :Author:	Sergey Kulagin
*)

(* :Package Version:    Dec 07, 2006
   :Updated:
*)

(* :Sources: *)

(* :Discussion:
This package reads the data file of
A. Kievsky et al Phys.Rev. C56 (1997) 64-75
He-3 spectral function and
returns the proton and neutron spectral functions.
*)

BeginPackage["kpsv`"]

protectkpsv = Unprotect[sfndat,dkpsv,ckpsv,pup,eup,nkmax,netot,ed,emn,pmn]
(*
sfndat::usage=
"sfndat[inuc,in] is {e,k,spfn} array for inuc=1,2 and in=1,2,3."
*)
dkpsv::usage=
"dkpsv[n,k] is the momentum distribution corresponding
to the deuteron pole contribution to n=0,1,2 spectral function."

ckpsv::usage=
"ckpsv[n,nuc,e,k] is the contribution from continous states to n=0,1,2
spectral function for proton (nuc=1) and neutron (nuc=2).
Energy in MeV, momentum in 1/fm and spectral function in fm^4."

pup::usage="Maximum momentum in the grid (in 1/fm)."
eup::usage="Maximum energy of spectator pair in the grid (in MeV)."
(*
nkmax::usage="Total number of momenta of the struck nucleon"
netot::usage="Total number of energies of spectator pair"
ed::usage="Deuteron binding energy (<0, in MeV)."
*)
emn::usage="Min energy in the grid"
pmn::usage="Min momentum in the grid"


(* Messages: *)
ckpsv::arg = "Illegal spectral function type `1` or nucleon isospin state `2` "
dkpsv::arg = "Illegal nucleon isospin state `1` "


Begin["`Private`"]
(* ENERGIES ARE IN MEV AND MOMENTA IN 1/FM *)

(* Options of numerical integration:
intopts = {PrecisionGoal -> 4, Method -> MultiDimensional, MaxRecursion -> 30};
*)
(* BEGIN reading the spectral function from data file *)
Print["Read KPSV data file and set up the spectral function..."];

(* eup   max energy of spectator pair in MeV.
   pup   max momentum of the struck nucleon in 1/fm
   netot total number of energies of spectator pair
   nkmax total number of momenta of the struck nucleon
   inuc=1 proton spectral function
   inuc=2 neutron spectral function
   pke=spl(inuc,it,3,ik,ie) unpolarized spectral function
       normalized to 4*Pi*\int dE/197.23 dk k^2 P(k,E) = 1
   f2ke=spl(inuc,it,2,ik,ie)
   f1ke=spl(inuc,it,1,ik,ie)
   pke, f1ke, f2ke in fm^4
*)

nema = {24,12,12,12};

Clear[inuc,it,ie,in,ik,e2f,e22,pki,eup,pup,netot,nkmax,sfndat,spl,dkpsv,ckpsv];
Clear[ddat,cdat];

sfndat[i_,j_] := {} /; (MemberQ[{1,2},i] && MemberQ[{1,2,3},j]);

rawdat = OpenRead["kpsv_spfn.dat"]

eup = Read[rawdat, Number];
pup = Read[rawdat, Number];
netot = Read[rawdat, Number];
nkmax = Read[rawdat, Number];


Do[(ne1=inuc;
  Do[(nep1=nema[[it]]+1;
      e2f = ReadList[rawdat, Number, nep1];
      pki = ReadList[rawdat, Number, nkmax];
        Do[(
           Do[(inn=Read[rawdat, Number];
               e22=Read[rawdat, Number];
                 Do[(spl=Read[rawdat, Number];
(* According to Salme we divide spl values in his data file by p *)
                     spl=spl/pki[[ik]];
                     If[( e22<0 && it==1 || e22>0 ), AppendTo[sfndat[inuc,in],{e22,pki[[ik]],spl}]]
                    ),{ik,1,nkmax} ]
              ),{in,1,3}]
           ),{ie,ne1,nep1}]
      )
  , {it, 1, 4}])
, {inuc, 1, 2}];

Close[rawdat]

(* END reading from data file *)

Clear[inuc,in,e,k];

ede = sfndat[1,3][[1,1]]; (* Deuteron binding energy *)
emn = sfndat[2,3][[1,1]]; (* Min energy in the grid *)
pmn = sfndat[2,3][[1,2]]; (* Min momentum in the grid *)

(* Separation of the deuteron pole from data file *)

 ddat[0] = (sfndat[1,3]// Take[#, nkmax] &)// Take[#, All, -2] &;
 ddat[1] = (sfndat[1,1]// Take[#, nkmax] &)// Take[#, All, -2] &;
 ddat[2] = (sfndat[1,2]// Take[#, nkmax] &)// Take[#, All, -2] &;


(*
ddat[in_] := (Which[ in==0, sfndat[1, in+3],
                    in==1 || in==2, sfndat[1, in]
                  ]// Take[#, nkmax] &)// Take[#, All, -2] &;
*)
(* Separation of continuum from data file *)
(*
cdat[1,in_] := Which[in==0, sfndat[1,in+3],
                     in==1 || in==2, sfndat[1,in]
                    ]// Drop[#,nkmax] &;

cdat[2,in_] := Which[in==0, sfndat[2,in+3],
                     in==1 || in==2, sfndat[2,in]
                    ];
*)
(* proton *)

cdat[1,0] = sfndat[1,3]// Drop[#,nkmax] &;
cdat[1,1] = sfndat[1,1]// Drop[#,nkmax] &;
cdat[1,2] = sfndat[1,2]// Drop[#,nkmax] &;

(* neutron *)

cdat[2,0] = sfndat[2,3];
cdat[2,1] = sfndat[2,1];
cdat[2,2] = sfndat[2,2];

(* Make an array of spin f_1 spectral function instead of B_1 *)

(* Deuteron pole term, only relevant for proton *)

ddat[1] = Table[{ddat[1][[k,1]], ddat[1][[k,2]] + ddat[2][[k,2]]/3}, {k,1,nkmax}];

(* Continuum *)

cdat[1,1] = Table[{cdat[1,1][[k,1]], cdat[1,1][[k,2]], cdat[1,1][[k,3]] + cdat[1,2][[k,3]]/3}, {k, 1, Length[cdat[1,1]]}];
cdat[2,1] = Table[{cdat[2,1][[k,1]], cdat[2,1][[k,2]], cdat[2,1][[k,3]] + cdat[2,2][[k,3]]/3}, {k, 1, Length[cdat[1,1]]}];

(*    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 *)

intopt2={InterpolationOrder -> {3,3}};
intopt1={InterpolationOrder -> 3};

(* Deuteron contribution to spectral function *)

dkpsv[in_,k_] := If[ MemberQ[{0,1,2},in]
                  , Interpolation[ddat[in], intopt1// Evaluate][k]
                  , (Message[dkpsv::arg, in]; Abort[])
                  ];

(* Construct the interpolation explicitly *)
(*
dkpsv[0,k_] = Interpolation[ddat[0], intopt1// Evaluate][k];
dkpsv[1,k_] = Interpolation[ddat[1], intopt1// Evaluate][k];
dkpsv[2,k_] = Interpolation[ddat[2], intopt1// Evaluate][k];
dkpsv[in_,k_] := If[ FreeQ[{0,1,2},in], (Message[dkpsv::arg, in]; Abort[]) ];
*)

(* Continuous spectrum contribution to spectral function *)

ckpsv[in_,inuc_,e_,k_] := Which[
     MemberQ[{0,1,2},in] && MemberQ[{1,2},inuc],Interpolation[cdat[inuc,in], intopt2// Evaluate][e,k]
    ,True, (Message[ckpsv::arg, in,inuc]; Abort[])
                              ];


(*
Clear[ddat];
Clear[cdat];
*)

End[ ]
Protect[Evaluate[protectkpsv]]
EndPackage[ ]

Print["If no error messages spectral function was successfully loaded.
For more information on defined functions type ?dkpsv ckpsv pup eup emn pmn."];

(* End of package *)