(* ::Package:: *)

 BeginPackage["nuosc`"]
 
 
 nuosc::usage= "nuosc[enu, imh, cphase, len, inu, in, out] computes the oscillation probability if in is not same as out or probability of survival if in=out. enu is energy in GeV, Len length in km,  inu = 1, -1 for neutrinos antineutrino, in = start flavor (1 e, 2 mu , 3 tau) and out flavor.  imh = 1, -1 is normal, inverted mass hierarchy, cphase is CP phase. "
 
(* ::Title:: *)
(*Neutrino Oscillation Calculator*)
(*Oscillation calculation including matter effects.  Uses formulation by  Ohlsson and Snellman (J. Math Phys. 41, No 5 May 2000).*)
(*Calculation has several modes. neutrino/antinueutrino, 1-> 1, 2-> 1, 2-> 3, etc.    1: electron, 2: muon, 3: tau neutrino.    *)
(*Code was translated from the original Fortran to Mathematica , Milind Diwan,   Feb 16, 2023.  *)
(**)
(*This is for constant matter density (which can also be set to 0) *)
(*This formulation is expanded to 3 neutrinos and it then requires consideration of which mass is higher.  m_3 or m2.  m1 is always less than m_2 because of the results from the solar data. *)

 
 Begin["Private`"]

 
 nuosc[enu_, imh_ , cphase_ , len_, inu_,in_, out_]:= Module[ {},
 
(* ::Subsection:: *)
(*Initialize parameters 	*)


(* ::Input::Initialization:: *)

dmsq21 =7.5*10^-5; 
dmsq32=imh*2.4*10^-3;  (*eV^2*)
dmsq31=dmsq21+dmsq32;  (*eV^2*)


\[Rho]E = 3;  (* Earth's density  gm/cm3*)

                                                       
(* ::Subsection:: *)
(*mixing angles  *)


(* ::Input::Initialization:: *)
\[Theta]1 = 35; (* 35 deg *) 
\[Theta]2  = 45; 
\[Theta]3 = 10;  (* 10 deg *) 


(* ::Section:: *)
(*Now we enter energy and length dependent calculations.  We will put them all inside a module.  *)


(* ::Subsection:: *)
(*Delta m^2 phases of importance.  Matter potential,  Calculate the exact "c"'s using Umatrix, Al, and Eij  *)



(* ::Input::Initialization:: *)
s1 = Sin[\[Theta]1*Pi/180];
s2 = Sin[\[Theta]2*Pi/180];
s3 = Sin[\[Theta]3*Pi/180];
c1 = Sqrt[1-s1^2];
c2 = Sqrt[1-s2^2];
c3 = Sqrt[1-s3^2];
\[Delta]cp = cphase*Pi/180;
               
delta =inu*\[Delta]cp;


(* ::Input::Initialization:: *)
U2 = {{1,0,0},{0,c2,s2},{0,-s2,c2}};
(*U3alt = {{c3,0,s3},{0,1,0},{-s3*Exp[-\[ImaginaryI]*delta],0,c3*Exp[-\[ImaginaryI]*delta]}}; *)  (* This is an alternate form *)
U3 = {{c3,0,s3*Exp[-I*delta]},{0,1,0},{-s3*Exp[I*delta],0,c3}};
U1 = {{c1,s1,0},{-s1,c1,0},{0,0,1}};
Umatrix =U2.U3.U1;


(* ::Input::Initialization:: *)
Ue1 = Umatrix[[1]][[1]];
Ue2 = Umatrix[[1]][[2]];
Ue3 = Umatrix[[1]][[3]];
Ue1s = Conjugate[Ue1]*Ue1;
Ue2s = Conjugate[Ue2]*Ue2;
Ue3s = Conjugate[Ue3]*Ue3;

                                           
E21 = 2*(1.27*dmsq21*len/enu);
E31 = 2*(1.27*dmsq31*len/enu);
E32 = 2*(1.27*dmsq32*len/enu); 
E12 = -E21;
E13 = -E31;
E23 = -E32; 
(* matter potential *)  
amatter = inu*7.56*10^-5*\[Rho]E*enu; 
AL = amatter*E32/dmsq32; 
(* calculate exact C's *)  
Czero = -(2/27)*AL*AL*AL + 
            (AL*AL/9)*(Ue1s*(E21+E31) + Ue2s(E12+E32)+Ue3s*(E13+E23) ) - (AL/9)*(Ue1s*(E21+E23)*(E31+E32) + Ue2s*(E12+E13)*(E31+E32) + Ue3s*(E12+E13)*(E21+E23) +(1/3)*(E12^2+E13^2+E23^2+E12*E13+E21*E23+E31*E32)) - (1/27)*(E12+E13)*(E21+E23)*(E31+E32) ; 
C1one = -AL*AL/3 + (AL/3)*(Ue1s*(E21+E31) + Ue2s*(E12+E32) + Ue3s*(E13+E23)) - (1/9)*(E12^2+E13^2+E23^2+E12*E13+E21*E23+E31*E32); 
(* test if this is negative *)  
Ctest = (C1one/3)^3+ (Czero/2)^2; 
If[Ctest>= 0, Return[0]]; 
TL1 = { AL*Conjugate[Ue1]*Ue1 - AL/3 + (E12+E13)/3,  AL*Conjugate[Ue1]*Ue2,  AL*Conjugate[Ue1]*Ue3} ; 
TL2 = { AL*Conjugate[Ue2]*Ue1 ,  AL*Conjugate[Ue2]*Ue2 - AL/3 + (E21+E23)/3,  AL*Conjugate[Ue2]*Ue3} ; 
TL3 = { AL*Conjugate[Ue3]*Ue1 ,  AL*Conjugate[Ue3]*Ue2 ,  AL*Conjugate[Ue3]*Ue3- AL/3 + (E31+E32)/3} ; 
TL = {TL1, TL2, TL3}; 
TraceT = TL[[1]][[1]]+TL[[2]][[2]]+TL[[3]][[3]]; 
     y = -Czero*Czero - (4./27.)*C1one^3 ; 
z = Sqrt[y]/Czero;  
x = ArcTan[z]/3; 
If[Czero>=0, x = x+ Pi]; 
\[Lambda] = { -Sqrt[-C1one/3]*Cos[x] + Sqrt[-C1one]*Sin[x], -Sqrt[-C1one/3]*Cos[x] - Sqrt[-C1one]*Sin[x], 2*Sqrt[-C1one/3]*Cos[x]} ; 

TLs= Table[ 0,{i,1,3},{j,1,3}] ; 
TLs = Table[ Sum[TL[[i]][[k]]*TL[[k]][[j]],{k,1,3}], {i,1,3},{j,1,3}]; 
TauL = Table[0,{i,1,3},{j,1,3}]; 
TauL = Table[Sum[Umatrix[[i]][[m]]*Conjugate[Umatrix[[j]][[n]]]*TL[[m]][[n]],{n,1,3},{m,1,3}] ,{i,1,3},{j,1,3}] ;
TauLs = Table[0,{i,1,3},{j,1,3}];  
TauLs = Table[  Sum[Umatrix[[i]][[m]]*Conjugate[Umatrix[[j]][[n]]]*TLs[[m]][[n]],{n,1,3},{m,1,3}] ,{i,1,3},{j,1,3}] ;
kdelta = Table[If[m==n,1,0],{m,1,3},{n,1,3}];
ProbC =  Sum[Exp[-I*\[Lambda][[j]]] * ( ((\[Lambda][[j]]^2+C1one)*kdelta[[in]][[out]] +\[Lambda][[j]]*TauL[[in]][[out]]+TauLs[[in]][[out]])/(3.*\[Lambda][[j]]^2+C1one)),{j,1,3}]; 
Prob = Conjugate[ProbC]*ProbC;
Return[Prob] ;  ]   
End[]
 EndPackage[]
 
 
