Quit[]

Cosmological parameters, time-redshift relation

(*Here we define the values of various cosmological parameters, and the relationship between the age of the universe and redshift, neglecting the cosmological constant. Cosmological parameters are taken from http://planck.caltech.edu/pub/2015results/Planck_2015_Results_XIII_Cosmological_Parameters.pdf, TT, TE, EE+low P 68% limits. *)
melec = 510.998903 10^3; (* in eV *)
OmegaM = 0.3156; 
OmegaLambda = 0.6844; 
OmegaRad = 8 10^(-5);
kmperMpc = 3.0857 10^19;
hval = 0.6727; 
H0 = 100 hval / kmperMpc;
OmegaB = 0.02225/hval^2;
OmegaC = 0.1198/hval^2;
Time[redshift_,OmegaM_,OmegaRad_,H0_] := 2 (OmegaM Sqrt[OmegaRad + OmegaM/redshift] + 2 OmegaRad^(3/2) redshift - 2 OmegaRad Sqrt[redshift (OmegaM + OmegaRad redshift)])/(3 OmegaM^2 redshift H0)
(*Approximate relationship neglecting dark energy, will fail at low redshifts - difference is <0.1% at 1+z=10. Note 'redshift' here means 1+z.*)
StandardTime[redshift_] := Time[redshift,OmegaM, OmegaRad, H0]
rhoDM = OmegaC*hval^2* 1.05375 10^(-5) ; (*present-day DM density in (GeV/c^2)/cm^3, taken from the PDG*)
HubbleParam[redshift_] := H0 Sqrt[OmegaM redshift^3 + OmegaRad redshift^4 + OmegaLambda];
ConversionFactor[redshift_] := HubbleParam[redshift]/redshift^3; (*"Conversion factor", defined as H(z)/(1+z)^3 where H(z) is the Hubble parameter*)

Import data, set up plotting defaults

(*Here we import the data files we will use in the remainder of the calculation*)

(**NEW RESULTS ON GENERAL ENERGY DEPOSITION FRACTIONS**)
DataFolder = "data/"; (*Insert the path to the folder where the .fits files are stored*)
Species = "elec"; (*adjust species to "phot" or "elec" here*)
mass = If[Species=="elec",melec,0];
DataFile = DataFolder <>Species <>"_processed_results.fits";
TempGrid=Import[DataFile,"Data"];
{outputredshift, energy, inputredshift,channels,origfractions,origf,newfractions,contcorrfractions,fion,fcorr,floss,fhigh,photenggrid,elecenggrid,fphotspec,felecspec} = { TempGrid[[1,1,1]],TempGrid[[1,1,2]], TempGrid[[1,1,3]], TempGrid[[1,1,4]],Partition[Partition[TempGrid[[1,1,5]],63],40],
Partition[TempGrid[[1,1,6]],63],Partition[Partition[Partition[TempGrid[[1,1,7]],63],40],63],Partition[Partition[TempGrid[[1,1,8]],63],40],
Partition[TempGrid[[1,1,9]],63],Partition[TempGrid[[1,1,10]],63],Partition[TempGrid[[1,1,11]],63],Partition[Partition[TempGrid[[1,1,12]],63],40], Partition[TempGrid[[1,1,13]],40], Partition[TempGrid[[1,1,14]],40], Partition[Partition[TempGrid[[1,1,15]],63],40], Partition[Partition[TempGrid[[1,1,16]],63],40]};
(**SUPPLEMENTARY DATA ON LOW-ENERGY PARTITION**)
SupplementFile = DataFolder <> "deposition_fractions_supplement.fits";
TempGrid2=Import[SupplementFile,"Data"];
{redshift,xh, basefracs, calibrationeng, fractable} =  {TempGrid2[[1,1,1]], TempGrid2[[1,1,2]], Partition[TempGrid2[[1,1,4]],63],TempGrid2[[1,1,5]],Partition[Partition[TempGrid2[[1,1,6]],6],63]};

(*PLOTTING DEFAULTS*)
PlotFolder = "plots/"; (*Insert the folder where you want the plots to go*)
PrettyExp[x_,MAX_:0]:=If[Abs[x]>MAX,Superscript[10,x],If[x≥0,(10^x),(10.^x)]]
fticks[min_,max_,AMAX_:0]:=Flatten[Table[Prepend[Flatten[Table[{x+Log[10,i]," ",{.01,0}},{i,2,9}],{1}],{x,PrettyExp[x,AMAX],{.03,0}}],{x,Ceiling[min],Floor[max],1}],1]
ticksb={fticks[0,12],fticks[1,4],None,None}; (*ticks on left and bottom only*)
ColorScheme = "Rainbow";

Plots for energy deposition fractions

(*Plots showing various integrals over the energy absorption function. Not necessary for later sections*)
(*Here we plot the aggregate fraction of the particle's energy that is deposited into various channels by redshift 10. Our best estimate is obtained by setting Prescription = ""; two versions of an alternate simplified method are implemented by setting Prescription to "SSCK" or "3keV".*)
SuffixList = {"Hion","Heion","LyA","Heat","Cont"};
Prescription = "";
fractionstable = Switch[Prescription,"",newfractions, "3keV",Table[fractable[[i,l,6]]*(origfractions[[j,k,l]] - contcorrfractions[[j,k,l]]) + If[i == 5, contcorrfractions[[j,k,l]],0],{i,1,Length[SuffixList]},{j,1,Length[inputredshift]},{k,1,Length[energy]},{l,1,Length[outputredshift]}],"SSCK",Table[basefracs[[i,l]]*(origfractions[[j,k,l]] - contcorrfractions[[j,k,l]]) + If[i == 5, contcorrfractions[[j,k,l]],0],{i,1,Length[SuffixList]},{j,1,Length[inputredshift]},{k,1,Length[energy]},{l,1,Length[outputredshift]}]];
ChannelPlotTotal[index_]:= Module[{Aggregate,thisPlot},Aggregate = Table[{energy[[j]], Log[10,inputredshift[[k]]],Sum[fractionstable[[index,k,j,i]],{i,1,Length[outputredshift]}]},{j,1,Length[energy]},{k,1,Length[inputredshift]}];
thisPlot = ListDensityPlot[Flatten[Aggregate,1], PlotRange ->{All,Log[10,{10,3000}],All},FrameLabel ->{"Energy (eV)","Redshift (1+z)"},LabelStyle->16, FrameTicks->ticksb, ColorFunction -> ColorScheme,PlotLegends -> BarLegend[Automatic,LegendMarkerSize -> 290,LabelStyle->16]];
Export[PlotFolder<>"totaldep_"<>Species<>"_"<>SuffixList[[index]]<>Prescription <>".pdf",thisPlot,"AllowRasterization"->True,ImageSize->360,ImageResolution->100];thisPlot]

ChannelPlotTotal[5]

(*Here we plot the amount of energy that is deposited into various channels as a function of redshift, normalized to the total power injected at that redshift, for the energy injection history corresponding to DM annihilation (that is, an equivalent "f-curve" separated by channel)*)
SuffixList = {"Hion","Heion","LyA","Heat","Cont"};
ChannelPlotF[index_]:= Module[{ChannelF,thisPlot},ChannelF = Table[{energy[[j]], Log[10,outputredshift[[i]]],ConversionFactor[outputredshift[[i]]] Sum[fractionstable[[index,k,j,i]] / ConversionFactor[inputredshift[[k]]],{k,1,Length[inputredshift]}]},{j,1,Length[energy]},{i,1,Length[outputredshift]}];
thisPlot = ListDensityPlot[Flatten[ChannelF,1], PlotRange ->{All,Log[10,{10,3000}],All},FrameLabel ->{"Energy (eV)","Redshift (1+z)"},LabelStyle->16, FrameTicks->ticksb, ColorFunction -> ColorScheme,PlotLegends -> BarLegend[Automatic,LegendMarkerSize -> 290,LabelStyle->16]];
Export[PlotFolder<>"fbychannel_"<>Species<>"_"<>SuffixList[[index]]<>Prescription<>".pdf",thisPlot,"AllowRasterization"->True,ImageSize->360,ImageResolution->100];thisPlot]

ChannelPlotF[1]

(*Here we plot the equivalent f-curve for total ionization, equivalent to summing channels 1 (H ionization) and 2 (He ionization) above*)
FIonPlot=
ListDensityPlot[Flatten[Table[{energy[[j]], Log[10,outputredshift[[i]]],fion[[j,i]]},{j,1,Length[energy]},{i,1,Length[outputredshift]}],1],  PlotRange ->{All,Log[10,{10,3000}],All},FrameLabel ->{"Energy (eV)","Redshift (1+z)"},LabelStyle->16, FrameTicks->ticksb, ColorFunction -> ColorScheme,PlotLegends -> BarLegend[Automatic,LegendMarkerSize -> 290,LabelStyle->16]]

(*Here we plot the equivalent overall f-curve, obtained by rescaling the deposited ionization to previous prescriptions for the partition of deposited energy into ionization*)
IonizationPrescription = "SSCK";
baseline = If[IonizationPrescription == "SSCK", basefracs[[1,;;]] +basefracs[[2,;;]],Module[{calibrationposition,newtable},calibrationposition = Flatten[Position[calibrationeng, 3000.]][[1]]; newtable=fractable[[;;,;;,calibrationposition]];newtable[[1,;;]] + newtable[[2,;;]]]];
(*Under the SSCK prescription, the power into ionization is (1 - xH)/3; it is stored in the basefracs array for each redshift (based on a background ionization history for xH, stored in the variable xh). Under the updated prescription used by the Planck Collaboration in 1502.01589, the ionization fractions are taken to be equal to those for 3 keV electrons (stored in the row of the array fractable corresponding to calibrationeng=3000). The appropriate equivalent f-curve will depend on the analysis you are using to set constraints -- analyses using the latter prescription will set slightly stronger limits on effective f, but effective f will be smaller for any individual model.*)
AggregateTotalF = Table[{energy[[j]], Log[10,outputredshift[[i]]],If[baseline[[i]] > 10^(-3),fion[[j,i]]/baseline[[i]],1000]},{j,1,Length[energy]},{i,1,Length[outputredshift]}];
thisPlot = ListDensityPlot[Flatten[AggregateTotalF,1], PlotRange ->{All,Log[10,{10,1200}],{0,2.5}},FrameLabel ->{"Energy (eV)","Redshift (1+z)"},LabelStyle->16, FrameTicks->ticksb, ColorFunction -> ColorScheme,ColorFunctionScaling -> False,PlotLegends -> BarLegend[ColorScheme,LegendMarkerSize -> 290,LabelStyle->16]]
Export[PlotFolder<>"fcurve_"<>Species<>"_"<>IonizationPrescription<>".pdf",thisPlot,"AllowRasterization"->True,ImageSize->360,ImageResolution->100];

(*Here we plot an alternate prescription for the equivalent overall f-curve, obtained by subtracting a correction from the total deposited energy, to account for losses to continuum photons by electrons above 3 keV. This gives slightly weaker constraints and is a conservative choice; it is more appropriate than the curves above when adapting constraints (computed using older prescriptions for the partition of deposited energy) from observables that depend significantly on contributions other than ionization*)
AggregateTotalFCorr = Table[{energy[[j]], Log[10,outputredshift[[i]]],fcorr[[j,i]]},{j,1,Length[energy]},{i,1,Length[outputredshift]}];
thisPlot = ListDensityPlot[Flatten[AggregateTotalFCorr,1], PlotRange ->{All,Log[10,{10,1200}],{0,2}},FrameLabel ->{"Energy (eV)","Redshift (1+z)"},LabelStyle->16, FrameTicks->ticksb, ColorFunction -> ColorScheme,ColorFunctionScaling -> False,PlotLegends -> BarLegend[ColorScheme,LegendMarkerSize -> 290,LabelStyle->16]]
Export[PlotFolder<>"fcurve_"<>Species<>"_continuumcorrected.pdf",thisPlot,"AllowRasterization"->True,ImageSize->360,ImageResolution->100];

Cross-checks, computing continuum losses

(*As a cross-check, fcorr can be computed from the contcorrfractions table, as follows:*)
fractionlost = Table[{energy[[j]], Log[10,outputredshift[[i]]],ConversionFactor[outputredshift[[i]]] Sum[contcorrfractions[[k,j,i]] / ConversionFactor[inputredshift[[k]]],{k,1,Length[inputredshift]}]},{j,1,Length[energy]},{i,1,Length[outputredshift]}];
fmodified = Table[{energy[[j]], Log[10,outputredshift[[i]]], origf[[j,i]] - fractionlost[[j,i,3]]},{j,1,Length[energy]},{i,1,Length[outputredshift]}];
(*We compare against the result from the panel above*)
GraphicsRow[{ListDensityPlot[Flatten[fmodified,1],  PlotRange ->{All,Log[10,{10,1200}],{0,2}},FrameLabel ->{"Energy (eV)","Redshift (1+z)"},FrameTicks->ticksb, ColorFunction -> ColorScheme],ListDensityPlot[Flatten[AggregateTotalFCorr,1], PlotRange ->{All,Log[10,{10,1200}],{0,2}},FrameLabel ->{"Energy (eV)","Redshift (1+z)"}, FrameTicks->ticksb,ColorFunction -> ColorScheme]}]
(*Likewise we can plot the fractional correction to the original f-curve, measuring total deposited energy, due to this effect; this is the fraction of "deposited" energy by more-than-3-keV electrons that ends up in low-energy continuum photons.*)
fractionalcontcorrection = Table[{energy[[j]], Log[10,outputredshift[[i]]],  Max[{fractionlost[[j,i,3]]/origf[[j,i]],0}]},{j,1,Length[energy]},{i,1,Length[outputredshift]}];
thisPlot =ListDensityPlot[Flatten[fractionalcontcorrection,1], PlotRange ->{All,Log[10,{10,1200}],{0,1}},ColorFunction -> ColorScheme,PlotLegends -> BarLegend[Automatic,LegendMarkerSize -> 290,LabelStyle->16],FrameLabel ->{"Energy (eV)","Redshift (1+z)"},LabelStyle->16, FrameTicks->ticksb]
Export[PlotFolder<>"f_contlossfraction_"<>Species <>".pdf",thisPlot/.{EdgeForm[],r_?(MemberQ[{RGBColor,Hue,CMYKColor,GrayLevel},Head[#]]&),i___}:>{EdgeForm[r],r,i}]
Export[PlotFolder<>"f_contlossfraction_"<>Species<>".pdf",thisPlot,"AllowRasterization"->True,ImageSize->360,ImageResolution->100];