"D-functions are according to hep-ph/9703315 and coincide with those of Julien's thesis (if we replace all the s_i with 1 in the latter)";

D1[x1_,x2_,x3_,x4_]=(4/Pi)Integrate[  A^(-2)  Sin[x1 A]  Sin[x2 A]  Sin [x3  A] Sin[x4 A] ,{A,0,Infinity},Assumptions-> x1∈Reals&&x2∈Reals&&x3∈Reals&&x4∈Reals]
(*Plot3D[D1[x1,x2,1,1],{x1,-2,2},{x2,-2,2}]*)
D2[x1_,x2_,x3_,x4_]= (4/Pi) (x3 x4) Integrate[  A^(-2)  Sin[x1 A] Sin[x2 A] (Cos[x3 A]- Sin[x3 A]/ (x3 A)) (Cos[x4 A]- Sin[x4 A]/ (x4 A)) ,{A,0,Infinity},Assumptions-> x1∈Reals&&x2∈Reals&&x3∈Reals&&x4∈Reals]
(*Plot3D[D2[x1,x2,1,1],{x1,-2,2},{x2,-2,2}]*)
"Dont know why this direct way of calculation is so time-consuming";
(*D3[x1_,x2_,x3_,x4_]= (4/Pi) (x1 x2 x3 x4) Integrate[  A^(-2)  (Cos[x1 A]- Sin[x1 A]/ (x1 A))(Cos[x2 A]- Sin[x2 A]/ (x2 A)) (Cos[x3 A]- Sin[x3 A]/ (x3 A)) (Cos[x4 A]- Sin[x4 A]/ (x4 A)) ,{A,0,Infinity}],Assumptions-> x1∈Reals&&x2∈Reals&&x3∈Reals&&x4∈Reals]*)
"So I use the answer from the Oleg's file that seem to plot the same result";
D3[x1_,x2_,x3_,x4_]=1/480 (-Abs[x1+x2-x3-x4]^5-Abs[x1-x2+x3-x4]^5+Abs[x1+x2+x3-x4]^5-Abs[x1-x2-x3+x4]^5+Abs[x1+x2-x3+x4]^5+Abs[x1-x2+x3+x4]^5+Abs[-x1+x2+x3+x4]^5-Abs[x1+x2+x3+x4]^5-120 x1 x2 x3 x4 (Abs[x1+x2-x3-x4]+Abs[x1-x2+x3-x4]+Abs[x1+x2+x3-x4]+Abs[x1-x2-x3+x4]+Abs[x1+x2-x3+x4]+Abs[x1-x2+x3+x4]+Abs[-x1+x2+x3+x4]+Abs[x1+x2+x3+x4])+20 x1 x2 (-Abs[x1+x2-x3-x4]^3+Abs[x1-x2+x3-x4]^3+Abs[x1+x2+x3-x4]^3+Abs[x1-x2-x3+x4]^3+Abs[x1+x2-x3+x4]^3-Abs[x1-x2+x3+x4]^3-Abs[-x1+x2+x3+x4]^3-Abs[x1+x2+x3+x4]^3)+20 x1 x3 (Abs[x1+x2-x3-x4]^3-Abs[x1-x2+x3-x4]^3+Abs[x1+x2+x3-x4]^3+Abs[x1-x2-x3+x4]^3-Abs[x1+x2-x3+x4]^3+Abs[x1-x2+x3+x4]^3-Abs[-x1+x2+x3+x4]^3-Abs[x1+x2+x3+x4]^3)+20 x1 x4 (Abs[x1+x2-x3-x4]^3+Abs[x1-x2+x3-x4]^3-Abs[x1+x2+x3-x4]^3-Abs[x1-x2-x3+x4]^3+Abs[x1+x2-x3+x4]^3+Abs[x1-x2+x3+x4]^3-Abs[-x1+x2+x3+x4]^3-Abs[x1+x2+x3+x4]^3)+20 x2 x3 (Abs[x1+x2-x3-x4]^3+Abs[x1-x2+x3-x4]^3+Abs[x1+x2+x3-x4]^3-Abs[x1-x2-x3+x4]^3-Abs[x1+x2-x3+x4]^3-Abs[x1-x2+x3+x4]^3+Abs[-x1+x2+x3+x4]^3-Abs[x1+x2+x3+x4]^3)+20 x2 x4 (Abs[x1+x2-x3-x4]^3-Abs[x1-x2+x3-x4]^3-Abs[x1+x2+x3-x4]^3+Abs[x1-x2-x3+x4]^3+Abs[x1+x2-x3+x4]^3-Abs[x1-x2+x3+x4]^3+Abs[-x1+x2+x3+x4]^3-Abs[x1+x2+x3+x4]^3)+20 x3 x4 (-Abs[x1+x2-x3-x4]^3+Abs[x1-x2+x3-x4]^3-Abs[x1+x2+x3-x4]^3+Abs[x1-x2-x3+x4]^3-Abs[x1+x2-x3+x4]^3+Abs[x1-x2+x3+x4]^3+Abs[-x1+x2+x3+x4]^3-Abs[x1+x2+x3+x4]^3)-60 x2 x3 x4 (-(-x1+x2+x3+x4)^2 Sign[x1-x2-x3-x4]-(x1+x2-x3-x4)^2 Sign[x1+x2-x3-x4]-(x1-x2+x3-x4)^2 Sign[x1-x2+x3-x4]-(x1+x2+x3-x4)^2 Sign[x1+x2+x3-x4]-(x1-x2-x3+x4)^2 Sign[x1-x2-x3+x4]-(x1+x2-x3+x4)^2 Sign[x1+x2-x3+x4]-(x1-x2+x3+x4)^2 Sign[x1-x2+x3+x4]-(x1+x2+x3+x4)^2 Sign[x1+x2+x3+x4])-60 x1 x2 x3 ((-x1+x2+x3+x4)^2 Sign[x1-x2-x3-x4]+(x1+x2-x3-x4)^2 Sign[x1+x2-x3-x4]+(x1-x2+x3-x4)^2 Sign[x1-x2+x3-x4]+(x1+x2+x3-x4)^2 Sign[x1+x2+x3-x4]-(x1-x2-x3+x4)^2 Sign[x1-x2-x3+x4]-(x1+x2-x3+x4)^2 Sign[x1+x2-x3+x4]-(x1-x2+x3+x4)^2 Sign[x1-x2+x3+x4]-(x1+x2+x3+x4)^2 Sign[x1+x2+x3+x4])-60 x1 x2 x4 ((-x1+x2+x3+x4)^2 Sign[x1-x2-x3-x4]+(x1+x2-x3-x4)^2 Sign[x1+x2-x3-x4]-(x1-x2+x3-x4)^2 Sign[x1-x2+x3-x4]-(x1+x2+x3-x4)^2 Sign[x1+x2+x3-x4]+(x1-x2-x3+x4)^2 Sign[x1-x2-x3+x4]+(x1+x2-x3+x4)^2 Sign[x1+x2-x3+x4]-(x1-x2+x3+x4)^2 Sign[x1-x2+x3+x4]-(x1+x2+x3+x4)^2 Sign[x1+x2+x3+x4])-60 x1 x3 x4 ((-x1+x2+x3+x4)^2 Sign[x1-x2-x3-x4]-(x1+x2-x3-x4)^2 Sign[x1+x2-x3-x4]+(x1-x2+x3-x4)^2 Sign[x1-x2+x3-x4]-(x1+x2+x3-x4)^2 Sign[x1+x2+x3-x4]+(x1-x2-x3+x4)^2 Sign[x1-x2-x3+x4]-(x1+x2-x3+x4)^2 Sign[x1+x2-x3+x4]+(x1-x2+x3+x4)^2 Sign[x1-x2+x3+x4]-(x1+x2+x3+x4)^2 Sign[x1+x2+x3+x4])-5 x2 (-(-x1+x2+x3+x4)^4 Sign[x1-x2-x3-x4]-(x1+x2-x3-x4)^4 Sign[x1+x2-x3-x4]+(x1-x2+x3-x4)^4 Sign[x1-x2+x3-x4]+(x1+x2+x3-x4)^4 Sign[x1+x2+x3-x4]+(x1-x2-x3+x4)^4 Sign[x1-x2-x3+x4]+(x1+x2-x3+x4)^4 Sign[x1+x2-x3+x4]-(x1-x2+x3+x4)^4 Sign[x1-x2+x3+x4]-(x1+x2+x3+x4)^4 Sign[x1+x2+x3+x4])-5 x3 (-(-x1+x2+x3+x4)^4 Sign[x1-x2-x3-x4]+(x1+x2-x3-x4)^4 Sign[x1+x2-x3-x4]-(x1-x2+x3-x4)^4 Sign[x1-x2+x3-x4]+(x1+x2+x3-x4)^4 Sign[x1+x2+x3-x4]+(x1-x2-x3+x4)^4 Sign[x1-x2-x3+x4]-(x1+x2-x3+x4)^4 Sign[x1+x2-x3+x4]+(x1-x2+x3+x4)^4 Sign[x1-x2+x3+x4]-(x1+x2+x3+x4)^4 Sign[x1+x2+x3+x4])-5 x4 (-(-x1+x2+x3+x4)^4 Sign[x1-x2-x3-x4]+(x1+x2-x3-x4)^4 Sign[x1+x2-x3-x4]+(x1-x2+x3-x4)^4 Sign[x1-x2+x3-x4]-(x1+x2+x3-x4)^4 Sign[x1+x2+x3-x4]-(x1-x2-x3+x4)^4 Sign[x1-x2-x3+x4]+(x1+x2-x3+x4)^4 Sign[x1+x2-x3+x4]+(x1-x2+x3+x4)^4 Sign[x1-x2+x3+x4]-(x1+x2+x3+x4)^4 Sign[x1+x2+x3+x4])-5 x1 ((-x1+x2+x3+x4)^4 Sign[x1-x2-x3-x4]-(x1+x2-x3-x4)^4 Sign[x1+x2-x3-x4]-(x1-x2+x3-x4)^4 Sign[x1-x2+x3-x4]+(x1+x2+x3-x4)^4 Sign[x1+x2+x3-x4]-(x1-x2-x3+x4)^4 Sign[x1-x2-x3+x4]+(x1+x2-x3+x4)^4 Sign[x1+x2-x3+x4]+(x1-x2+x3+x4)^4 Sign[x1-x2+x3+x4]-(x1+x2+x3+x4)^4 Sign[x1+x2+x3+x4]));
(*Plot3D[D3[x1,x2,1,1],{x1,-2,2},{x2,-2,2}]*)

grid = Import["~/Repositories/pyBBN/src/tests/grid.json"]
outer = Outer[{#1,#2,#3,#4}&,grid,grid,grid,grid];

Length@Gather@Flatten@(Apply[D1,outer,{4}]- Import["~/Repositories/pyBBN/src/tests/D1.json"])/(((Length@grid))^4)
Union@Flatten@(Apply[D1,outer,{4}]/Import["~/Repositories/pyBBN/src/tests/D1.json"])

Length@Gather@Flatten@(Apply[D2,outer,{4}]- Import["~/Repositories/pyBBN/src/tests/D2.json"])/(((Length@grid))^4)
Apply[D2,outer,{4}]/ Import["~/Repositories/pyBBN/src/tests/D2.json"]

Length@Gather@Flatten@(Apply[D3,outer,{4}]- Import["~/Repositories/pyBBN/src/tests/D3.json"])/(((Length@grid))^4)
Apply[D3,outer,{4}]/ Import["~/Repositories/pyBBN/src/tests/D3.json"]

SetDirectory["/media/data/anderl/Documents/Work/programming/bbn_code"];
Print[Directory[]];
Print[FileExistsQ["./Input/ew_nofluctuations.dat"]];

"Standard quantities";
"All energy quantiiies are measured in MeV";
GF=1.17 10^-11;
geff=11.;
MPl=1.22 10^22;
sectoMeV=1.52 10^21;
(Subscript[θ, W])=ArcSin[Sqrt[0.23]];gL=.5+Sin[(Subscript[θ, W])]^2 ;gR=Sin[(Subscript[θ, W])]^2;gtL=-.5+Sin[(Subscript[θ, W])]^2;
nF[EoT_]=1/(Exp[EoT]+1);
(*SetDirectory[""];*)
EWdata=Import["./Input/ew_nofluctuations.dat"];
ge=EWdata[[All,{1,8}]];
he=EWdata[[All,{1,9}]];
geffEW=Interpolation[ge,InterpolationOrder->1];
heffEW=Interpolation[he,InterpolationOrder->1];
"Hubble rate";
H[TTemp_]:=Sqrt[4Pi^3/45 geffEW[TTemp]]TTemp^2/MPl;
"a is the scale factor, (Subscript[a, norm])(Subscript[T, norm])=1 is the normalization";
Tnorm=70;
a[TTemp_]:=(heffEW[Tnorm]/heffEW[TTemp])^(1/3)/TTemp;
(*LogPlot[96Pi^3/GF^2/Mass^5/sectoMeV,{Mass,20,1},PlotRange->{{1,20},{0.005,10^4}}, Frame->True,FrameLabel->{"(Subscript[M, s]) [MeV]","(Subscript[τ, N]) [sec]"},Filling->{1->{Bottom,RGBColor[0.8,0.8,0.95]}}]
Export["Lifetime-theta=1.eps",%];*)

(*"We assume that all particles but sterile neutrino are relativistic so that their masses are negligible";
KfuncIntegrand[ps_,p2_,p3_,p4_,Ms_]:=128*96 Pi^3*(3*Es p2 p3 p4 D1[ps,p2,p3,p4]+3D3[ps,p2,p3,p4]+Es p4 D2[ps,p4,p2,p3]-p2 p3 D2[p2,p3,ps,p4]+Es p2 D2[ps,p2,p3,p4]+p3 p4 D2[p3,p4,ps,p2])/.{Es->Sqrt[ps^2+Ms^2]}
"Momentum ranges are taken finite to reduce computation time";
Kfunc[EoT_,MoT_]:=If[EoT-MoT>0,(1/64/Pi^3/EoT/Sqrt[EoT^2-MoT^2])NIntegrate[HeavisideTheta[EoT+p2oT-p3oT](1-nF[p2oT])nF[p3oT]nF[EoT+p2oT-p3oT]/nF[EoT]KfuncIntegrand[Sqrt[EoT^2-MoT^2],p2oT,p3oT,EoT+p2oT-p3oT,MoT],{p2oT,0,15},{p3oT,0,15},Method->"AdaptiveMonteCarlo"],0]; 
"Kfunc is a function of only two arguments and hence can be interpolated only ONCE and used in advance for ALL computations";
KfuncExport={};
Do[AppendTo[KfuncExport,{poT,MoT,Kfunc[Sqrt[poT^2+MoT^2],MoT]}],{poT,0,10,1},{MoT,0,20,2}];
Export["Input/Kfunc2.dat",KfuncExport];*)

KfInt=Import["Input/Kfunc2.dat"];
KfInter=KfInt[[All,{1,2,3}]];
KfuncImport=Interpolation[KfInter,InterpolationOrder->{2,2}];
Plot3D[KfuncImport[poT,MoT],{poT,0,5},{MoT,0,20}];
"Rate of sterile neutrino scatterings and decays";
Γscat[ps_,Ms_,T_,τ_]:=If[(ps/T-5)(ps/T)<0&&(Ms/T-20)(Ms/T)<0,T^5KfuncImport[ps/T,Ms/T]/τ/Ms^5 ,0]
Γdec[ps_,Ms_,τ_]:= Ms/Sqrt[ps^2+Ms^2]/τ
"τ here is the sterile neutrino lifetime measured in MeV";

Msinglet=140;τsinglet=0.07*sectoMeV;
(*Kinetic equation integration range*)Tin=100;Tfin=3.;
"Solution of Boltzmann equation";FSinglet[pconfTemp_,TTemp_]:=fS[T]/.NDSolve[{fS'[T]==(Γscat[pconfTemp/a[T],Msinglet,T,τsinglet]+Γdec[pconfTemp/a[T],Msinglet,τsinglet])/T/H[T]*( fS[T]-nF[Sqrt[pconfTemp^2/a[T]^2+Msinglet^2]/T] ),fS[Tin]==nF[Sqrt[pconfTemp^2/a[Tin]^2+Msinglet^2]/Tin]},fS[T],{T,Tin,Tfin}][[1]]/.{T->TTemp};
(*Export at temperature*)TExp=5.;FSingletExport={};Do[AppendTo[FSingletExport,{pPhys,FSinglet[pPhys*a[TExp],TExp]}],{pPhys,0,10.TExp,TExp/10.}];
Export["Output/Distribution_Ms-140_tau-0.07s_TExp-5MeV.dat",FSingletExport];