Add Continuous Flux Injection from Differential MCEq Profile as a new Body Class and Corresponding Updates

configure

@@ -709,6 +709,7 @@ OBJECTS = $(patsubst src/%.cpp,build/%.o,$(SOURCES))
 EXAMPLES := examples/Single_energy/single_energy \
             examples/Multiple_energy/multiple_energy \
             examples/Atm_default/atm_default \
+            examples/Emitting_Atmosphere/emitting_atmosphere \
             examples/Bodies/bodies \
             examples/Xsections/xsections \
             examples/NSI/nsi \
@@ -799,6 +800,10 @@ examples/Atm_default/atm_default : $(DYN_PRODUCT) examples/Atm_default/main.cpp
 	@echo Compiling atmospheric example
 	@$(CXX) $(EXAMPLES_FLAGS) examples/Atm_default/main.cpp -lnuSQuIDS $(LDFLAGS) -o $@
 
+examples/Emitting_Atmosphere/emitting_atmosphere : $(DYN_PRODUCT) examples/Emitting_Atmosphere/main.cpp
+	@echo Compiling emitting atmosphere example
+	@$(CXX) $(EXAMPLES_FLAGS) examples/Emitting_Atmosphere/main.cpp -lnuSQuIDS $(LDFLAGS) -o $@
+
 build/exBody.o : examples/Bodies/exBody.h examples/Bodies/exBody.cpp
 	@$(CXX) $(EXAMPLES_FLAGS) -c examples/Bodies/exBody.cpp -o $@
 

examples/Emitting_Atmosphere/main.cpp

@@ -0,0 +1,189 @@
+ /******************************************************************************
+ *    This program is free software: you can redistribute it and/or modify     *
+ *   it under the terms of the GNU General Public License as published by      *
+ *   the Free Software Foundation, either version 3 of the License, or         *
+ *   (at your option) any later version.                                       *
+ *                                                                             *
+ *   This program is distributed in the hope that it will be useful,           *
+ *   but WITHOUT ANY WARRANTY; without even the implied warranty of            *
+ *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the             *
+ *   GNU General Public License for more details.                              *
+ *                                                                             *
+ *   You should have received a copy of the GNU General Public License         *
+ *   along with this program.  If not, see <http://www.gnu.org/licenses/>.     *
+ *                                                                             *
+ *   Authors:                                                                  *
+ *      Carlos Arguelles (University of Wisconsin Madison)                     *
+ *         carguelles@icecube.wisc.edu                                         *
+ *      Jordi Salvado (University of Wisconsin Madison)                        *
+ *         jsalvado@icecube.wisc.edu                                           *
+ *      Christopher Weaver (University of Wisconsin Madison)                   *
+ *         chris.weaver@icecube.wisc.edu                                       *
+ ******************************************************************************/
+
+#include <vector>
+#include <iostream>
+#include <fstream>
+#include <nuSQuIDS/nuSQuIDS.h>
+
+/*
+ * The problem of solving the propagation of the atmospheric neutrinos is an
+ * energy and zenith dependent problem, for this we include the class nuSQUIDSAtm
+ * that allows to solve a set of nuSUIDS energy dependent objects to take in to account the 
+ * zenith dependence. The problem of flux insertion into the atmosphere due to cosmic ray sources
+ * during evolution requires a class derived from EarthAtm to be called into nuSQUIDSAtm
+ */
+
+
+using namespace nusquids;
+
+
+//Function that sets the initial flux if you want to add flux in at the top of the atmosphere
+//for this example we set it to 0
+double flux_function(double enu, double cz){
+  return 0;
+}
+
+
+int main()
+{
+  //Units and constants class
+  squids::Const units;
+  //Number of neutrinos (3) standard 4 1-sterile
+  //Number of neutrinos, 4 is the case with one sterile neutrino
+  std::cout << "Enter number of neutrino flavors to consider:\n";
+  std::cout << "(3) Three Active Neutrinos, " << "(4) 3+1 Three Active and One Sterile Neutrino" << std::endl;
+  unsigned int numneu;
+  std::cin >>numneu;
+  if( not(numneu==3 || numneu==4)){
+    throw std::runtime_error("Only (3) or (4) are valid options");
+  }
+  std::cout << "Consider Earth absorption and interactions? " << std::endl;
+  std::string use_int;
+  std::cin >> use_int;
+  bool interactions = false;
+  if(use_int=="yes" || use_int=="y")
+    interactions = true;
+
+  //Minimum and maximum values for the energy and cosine zenith, notice that the energy needs to have the 
+  //units, if they are omitted the input is in eV.
+  double Emin=1e-2*units.GeV;
+  double Emax=1.e6*units.GeV;
+  double czmin=-1;
+  double czmax=1;
+
+  //Declaration of the atmospheric object
+  //Must include the nuSQUIDS template or nuSQUIDSNSI user defined class derived from nuSQUIDS in order to 
+  //use the EmittingEarthAtm class or another user defined class derived from EarthAtm
+  std::cout << "Begin: constructing nuSQuIDS-Atm object" << std::endl;
+  nuSQUIDSAtm<nuSQUIDS,EmittingEarthAtm> nus_atm(linspace(czmin,czmax,40),logspace(Emin,Emax,100),numneu,both,interactions);
+  std::cout << "End: constructing nuSQuIDS-Atm object" << std::endl;
+
+  //Sets height of atmosphere in km (defaults to 22 km when not set)
+  nus_atm.Set_AtmHeight(60);
+
+  //Functions Necessary for injection
+  //Passes the energy range into the emitting bodies for interpolation
+  nus_atm.Set_AtmEmissionEnergies();
+  //Allows the nusquids objects to accept neutrino sources
+  nus_atm.Set_NeutrinoSources(true);
+
+  std::cout << "Begin: setting mixing angles." << std::endl;
+  // set mixing angles, mass differences and cp phases
+  nus_atm.Set_MixingAngle(0,1,0.563942);
+  nus_atm.Set_MixingAngle(0,2,0.154085);
+  nus_atm.Set_MixingAngle(1,2,0.785398);
+
+  nus_atm.Set_SquareMassDifference(1,7.65e-05);
+  nus_atm.Set_SquareMassDifference(2,0.00247);
+
+  nus_atm.Set_CPPhase(0,2,0);
+  if(numneu > 3){
+    nus_atm.Set_SquareMassDifference(3,-1.);
+    nus_atm.Set_MixingAngle(1,3,0.160875);
+  }
+  std::cout << "End: setting mixing angles." << std::endl;
+
+  //Setup integration precision
+  nus_atm.Set_rel_error(1.0e-6);
+  nus_atm.Set_abs_error(1.0e-6);
+  nus_atm.Set_GSL_step(gsl_odeiv2_step_rk4);  
+
+  //Array that contains the values of the energies and cosine of the zenith, is the same length for every zenith
+  auto e_range = nus_atm.GetERange();
+  auto cz_range = nus_atm.GetCosthRange();
+
+  std::cout << "Begin: setting initial state." << std::endl;
+
+  //Construct the initial state at the top of the atmosphere, we set it to 0 here
+  marray<double,4> inistate{nus_atm.GetNumCos(),nus_atm.GetNumE(),2,numneu};
+  std::fill(inistate.begin(),inistate.end(),0);
+  for ( int ci = 0 ; ci < nus_atm.GetNumCos(); ci++){
+    for ( int ei = 0 ; ei < nus_atm.GetNumE(); ei++){
+      for ( int rho = 0; rho < 2; rho ++ ){
+        for (int flv = 0; flv < numneu; flv++){
+          inistate[ci][ei][rho][flv] = (flv == 0) ? flux_function(e_range[ei], cz_range[ci]) : 0;//set the flux for the muon flavor: 1
+          inistate[ci][ei][rho][flv] = (flv == 1) ? flux_function(e_range[ei], cz_range[ci]) : 0;//set the flux for the muon flavor: 1
+        }
+      }
+    }
+  }
+
+  //Set the initial state in the atmSQuIDS object
+  nus_atm.Set_initial_state(inistate,flavor);
+  std::cout << "End: setting initial state." << std::endl;
+
+
+  //Set to true the monitoring prgress bar and the vacuum oscillations
+  nus_atm.Set_ProgressBar(true);
+  nus_atm.Set_IncludeOscillations(true);
+
+  //Here we do the evolution of all the states
+  std::cout << "Begin: Evolution" << std::endl;
+  nus_atm.EvolveState();
+  std::cout << "End: Evolution" << std::endl;
+
+  //We can save the current state in HDF5 format for future use.
+  //nus_atm.WriteStateHDF5("./atmospheric_emission_example_numneu_"+std::to_string(numneu)+".hdf5");
+
+  //This file will contain the final flux, since initially we set it to 1 for the muon, 
+  //this can be read as the muon ration F_final/F_initial in cos(zenith) and energy.
+  std::ofstream file("fluxes_flavor.txt");
+
+
+  //Set the resolution and the ranges for the ouput, remember that an interpolation in energy is used in 
+  //in the interaction picture, the vacuum oscillations are solve analytically with arbitrary Energy precision.
+  //For the zenith a linear interpolation is used.
+  int Nen=700;
+  int Ncz=100;
+  double lEmin=log10(Emin);
+  double lEmax=log10(Emax);;
+
+  //Writing to the file!  
+  file << "# log10(E) cos(zenith) E flux_i . . . ." << std::endl;
+  for(double cz=czmin;cz<czmax;cz+=(czmax-czmin)/(double)Ncz){
+    for(double lE=lEmin; lE<lEmax; lE+=(lEmax-lEmin)/(double)Nen){
+      double E=pow(10.0,lE);
+      file << lE - log10(units.GeV) << " " << cz << " " << E/units.GeV;
+      for(int fl=0; fl<numneu; fl++){
+        file << " " <<  nus_atm.EvalFlavor(fl,cz, E, 0);
+      }
+      for(int fl=0; fl<numneu; fl++){
+        file << " " <<  nus_atm.EvalFlavor(fl,cz, E, 1);
+      }
+      file << std::endl;
+    }
+    file << std::endl;
+  }
+
+  //This asks if you want to run the gnuplot plotting script.
+  std::string plt;
+  std::cout << std::endl <<  "Done! " << std::endl <<
+    "  Do you want to run the gnuplot script? yes/no" << std::endl;
+  std::cin >> plt;
+  if(plt=="yes" || plt=="y")
+    return system("./plot.plt");
+
+
+  return 0;
+}

examples/Emitting_Atmosphere/plot.plt

@@ -0,0 +1,23 @@
+#!/usr/bin/env gnuplot
+if ( GPVAL_VERSION >= 4.4 && strstrt(GPVAL_TERMINALS, 'wxt') > 0 ) set terminal wxt persist
+if ( GPVAL_VERSION >= 4.4 && strstrt(GPVAL_TERMINALS, 'qt') > 0 ) set terminal qt persist
+if ( GPVAL_VERSION >= 4.4 && strstrt(GPVAL_TERMINALS, 'wxt') == 0 && strstrt(GPVAL_TERMINALS, 'qt') == 0 ) print "wxt and qt terminals not available, proceeding with default"
+if ( GPVAL_VERSION < 4.4 ) print "gnuplot is too old to check for available terminals" ; print "attempting to use wxt terminal and hoping for the best" ; set terminal wxt persist
+
+#set cbrange [0:1]
+set logscale cb
+set xlabel "log_{10}(E/GeV)"
+set ylabel "Cos(zenith)"
+set pm3d map
+splot "fluxes_flavor.txt" u 1:2:5
+
+set xrange [2.:6.]
+set terminal png size 800,600 enhanced font 'Verdana,10'
+set output "plot.png"
+#set terminal postscript enhance eps color
+#set output "plot.eps"
+#set terminal svg size 350,262 fname 'Verdana' fsize 10
+#set output "plot.svg"
+replot
+
+

include/nuSQuIDS/body.h

@@ -669,6 +669,7 @@ class EarthAtm: public Body{
     double x_ye_min;
     /// \brief Electron fraction at maximum radius.
     double x_ye_max;
+
   public:
     /// \brief Default constructor using supplied PREM.
     EarthAtm();
@@ -699,6 +700,7 @@ class EarthAtm: public Body{
     /// \brief EarthAtm trajectory
     class Track: public Body::Track{
       friend class EarthAtm;
+      friend class EmittingEarthAtm;
       private:
         /// \brief Cosine of the zenith angle.
         double cosphi;
@@ -762,6 +764,90 @@ class EarthAtm: public Body{
     Track MakeTrackWithCosine(double cosphi);
 };
 
+
+
+/// \class EmittingEarthAtm
+/// \brief A model of the production of neutrino cosmic ray flux production within the atmosphere.
+class EmittingEarthAtm: public EarthAtm{
+  protected:
+
+    /// \brief Min/Max coszen.
+    double czen_min;
+    double czen_max;
+    /// \brief Min/Max height.
+    double height_min;
+    double height_max;
+    /// \brief Min/Max energy.
+    double energy_min;
+    double energy_max;
+
+    /// \brief vector containing interpolators for the differential flux
+    std::vector<std::vector<TriCubicInterpolator>> Interpolators;
+
+    /// \brief indices for the electron and muon flavors
+    /// defaulted to the 0 and 1 states (first and second0
+    /// WARNING!!!: careful when changing these as other parts of nuSQuIDS are hard coded to
+    /// default to 0, 1, & 2 for nu_e, nu_mu, and nu_tau respectively
+    /// Do this only if you have defined your own HI matrix by overriding squids::SU_vector HI
+    /// and if you have explicitly initialized all non-used mixing angles to 0
+    int nu_e_index = 0;
+    int nu_mu_index = 1;
+
+  public:
+
+    /// \brief Default constructor using supplied MCEQ neutrino flux.
+    EmittingEarthAtm();
+
+    /// \brief Constructor from a user supplied production model.
+    /// @param prodmodel Path to the production model data file.    
+    /// \details The input file should have 3+2n columns for n=the number of produced flavors.
+    /// The first one is cosine(zenith) in descending order from or between 1 and -1
+    /// representing the range of coszen values corresponding to tracks relevant to your propagation
+    /// The second column must contain the height within the atmosphere in cm
+    /// in descending order from or between 6000000 (60 km) and 0
+    /// The third column must correspond to the energy in GeV 
+    /// in descending order over any relevant range
+    /// The remaining columns must correspond to the differential produced neutrino fluxes for each flavor
+    /// The order is assumed nu_e, nu_e_bar, nu_mu, nu_mu_bar with additional flavors added 
+    /// with the corresponding antineutrino differential produced flux in the following column
+    /// the produced flux corresponds to the coszen, height, and energy in the row
+    EmittingEarthAtm(std::string prodmodel);
+
+    /// \brief Deconstructor
+    ~EmittingEarthAtm();
+
+    /// \brief Function that sets the height of the atmosphere
+    /// \details Used by nuSQUIDSATM to set the height for every
+    /// instance of the body
+    void SetAtmosphereHeight(double height){
+      atm_height = height;
+      earth_with_atm_radius = radius + atm_height;
+    }
+    /// \brief set index of produced flavors
+    /// \details Set the flavor index of the electron and muon neutrinos
+    /// or to different flavor indices as user input data requires
+    /// WARNING!!!: careful when using this as other parts of nuSQuIDS are hard coded to
+    /// default to 0, 1, & 2 for nu_e, nu_mu, and nu_tau respectively
+    /// Do this only if you have defined your own HI matrix by overriding squids::SU_vector HI
+    /// and if you have explicitly initialized all non-used mixing angles to 0
+    void Set_ProducedFlavors(int nu_e_index_, int nu_mu_index_){
+      nu_e_index = nu_e_index_;
+      nu_mu_index = nu_mu_index_;
+    }
+
+    //Important function that overrides the default flux injection with ATM data
+    void injected_neutrino_flux(marray<double, 3>& flux, const GenericTrack& track, const nuSQUIDS& nusquids) override;
+
+    /// \brief Returns the body identifier.
+    static unsigned int GetId() {return 8;}
+    /// \brief Returns the name of the body.
+    static std::string GetName() {return "EmittingEarthAtm";}
+};
+
+
+
+
+
   // type defining
   typedef Body::Track Track;