Update python bindings install docs & fix typo in SetAtmosphereHeight binding

README.md

@@ -106,12 +106,7 @@ Compiling the Python Interface
 
 If you have activated the python interface by doing `--with-python-bindings` then proceed according to this instructions
 Even when configured the python interface is not built with the main library. 
-To compile it do the following:
+To compile and install it do the following:
 
-	cd resources/python/src/
-	make
-
-After successful compilation the bindings will be stored in `resources/python/bindings/`. 
-To make them available from within python, modify your PYTHONPATH:
-
-	export PYTHONPATH=$(PATH_TO_nuSQUIDS)/resources/python/bindings/:$PYTHONPATH
+	make python
+	make python-install

include/nuSQuIDS/nuSQuIDSLV.h

@@ -0,0 +1,77 @@
+#ifndef nusquidslv_H
+#define nusquidslv_H
+
+#include <vector>
+#include <nuSQuIDS/nuSQuIDS.h>
+
+namespace nusquids {
+
+struct LVParameters {
+  gsl_complex c_emu;
+  gsl_complex c_mutau;
+};
+
+class nuSQUIDSLV: public nuSQUIDS {
+  private:
+    bool lv_parameters_set = false;
+    int n_ = 1;
+    LVParameters c_params;
+    squids::SU_vector LVP;
+    std::vector<squids::SU_vector> LVP_evol;
+
+    void AddToPreDerive(double x);
+    void AddToWriteHDF5(hid_t hdf5_loc_id) const;
+    void AddToReadHDF5(hid_t hdf5_loc_id);
+    squids::SU_vector HI(unsigned int ie, unsigned int irho) const;
+    void initialize_LVP_evol_at_each_energy_node();
+
+  public:
+    nuSQUIDSLV() {}
+    nuSQUIDSLV(marray<double,1> E_vector,unsigned int numneu, NeutrinoType NT = both, bool iinteraction = false, std::shared_ptr<CrossSectionLibrary> ncs = nullptr);
+    nuSQUIDSLV(std::string hdf5_filename, std::string grp = "/", std::shared_ptr<InteractionStructure> int_struct = nullptr);
+    nuSQUIDSLV(unsigned int numneu, NeutrinoType NT = neutrino);
+
+    void Set_LV_OpMatrix(LVParameters &lv_params);
+    void Set_LV_OpMatrix(gsl_matrix_complex *cmatrix);
+    void Set_LV_Operator(squids::SU_vector op);
+    void Set_LV_EnergyPower(int n);
+    void Set_MixingAngle(unsigned int i, unsigned int j, double angle);
+    void Set_CPPhase(unsigned int i, unsigned int j, double angle);
+    void dump_probabilities() const;
+
+    void Set_initial_state(const marray<double, 1> &ini_state, Basis basis = flavor);
+    void Set_initial_state(const marray<double, 2> &ini_state, Basis basis = flavor);
+    void Set_initial_state(const marray<double, 3> &ini_state, Basis basis = flavor);
+};
+
+class nuSQUIDSLVAtm : public nuSQUIDSAtm<nuSQUIDSLV> {
+
+  public:
+
+    // Use the base class constructors
+    using nuSQUIDSAtm<nuSQUIDSLV>::nuSQUIDSAtm;
+
+    // nuSQUIDSLVAtm has an attribute nusq_array that contains one nuSQUIDSLV object
+    // for each zenith. This array can be obtained with this->GetnuSQuIDS(). The Setter
+    // methods have to be wrapped such that the Setter of each element of that array is called
+
+    void Set_LV_OpMatrix(LVParameters &lv_params){
+      for(nuSQUIDSLV& nsq : this->GetnuSQuIDS()) nsq.Set_LV_OpMatrix(lv_params);
+    }
+
+    void Set_LV_OpMatrix(gsl_matrix_complex *cmatrix){
+      for(nuSQUIDSLV& nsq : this->GetnuSQuIDS()) nsq.Set_LV_OpMatrix(cmatrix);
+    }
+
+    void Set_LV_Operator(squids::SU_vector op){
+      for(nuSQUIDSLV& nsq : this->GetnuSQuIDS()) nsq.Set_LV_Operator(op);
+    }
+
+    void Set_LV_EnergyPower(int n){
+      for(nuSQUIDSLV& nsq : this->GetnuSQuIDS()) nsq.Set_LV_EnergyPower(n);
+    }
+
+}; 
+} // namespace nusquids
+
+#endif // nusquidslv_H

resources/python/src/nuSQUIDSpy.cpp

@@ -22,6 +22,17 @@
  ******************************************************************************/
 
 #include "nuSQUIDSpy.h"
+#include <nuSQuIDS/nuSQuIDSLV.h>
+
+// Helper function to convert `gsl_complex` to Python `complex`
+std::complex<double> GSLComplexToPython(gsl_complex c) {
+    return std::complex<double>(GSL_REAL(c), GSL_IMAG(c));
+}
+
+// Helper function to convert Python `complex` to `gsl_complex`
+gsl_complex PythonToGSLComplex(const std::complex<double>& c) {
+    return gsl_complex_rect(c.real(), c.imag());
+}
 
 BOOST_PYTHON_MODULE(nuSQuIDS)
 {
@@ -376,7 +387,7 @@ BOOST_PYTHON_MODULE(nuSQuIDS)
     .def(init<std::string>())
     .def("GetRadius",&EarthAtm::GetRadius)
     .def("GetAtmosphereHeight",&EarthAtm::GetAtmosphereHeight)
-    .def("SetAtmosphereHeight",&EarthAtm::GetAtmosphereHeight)
+    .def("SetAtmosphereHeight",&EarthAtm::SetAtmosphereHeight)
     .def("MakeTrack",&EarthAtm::MakeTrack)
     .def("MakeTrackWithCosine",&EarthAtm::MakeTrackWithCosine)
     ;
@@ -409,4 +420,148 @@ BOOST_PYTHON_MODULE(nuSQuIDS)
   marray_from_python<2>();
   marray_from_python<3>();
   marray_from_python<4>();
+
+  // Bind LVParameters
+    boost::python::class_<nusquids::LVParameters>("LVParameters", "Class representing Lorentz Violation parameters")
+      .add_property(
+          "c_emu",
+          +[](const nusquids::LVParameters& self) { return GSLComplexToPython(self.c_emu); },
+          +[](nusquids::LVParameters& self, const std::complex<double>& value) { self.c_emu = PythonToGSLComplex(value); }
+      )
+      .add_property(
+          "c_mutau",
+          +[](const nusquids::LVParameters& self) { return GSLComplexToPython(self.c_mutau); },
+          +[](nusquids::LVParameters& self, const std::complex<double>& value) { self.c_mutau = PythonToGSLComplex(value); }
+      );
+
+  // Bind nuSQUIDSLV
+  auto nusquids_lv_register = RegisterBasicNuSQuIDSPythonBindings<nuSQUIDSLV>("nuSQUIDSLV");
+  auto nusquids_lv_class_object = nusquids_lv_register.GetClassObject();
+  nusquids_lv_class_object->def(
+      "Set_LV_OpMatrix",
+      (void (nuSQUIDSLV::*)(nusquids::LVParameters&)) &nuSQUIDSLV::Set_LV_OpMatrix
+  );
+  nusquids_lv_class_object->def(
+      "Set_LV_OpMatrix",
+      (void (nuSQUIDSLV::*)(gsl_matrix_complex*)) &nuSQUIDSLV::Set_LV_OpMatrix
+  );
+  nusquids_lv_class_object->def(
+      "Set_LV_Operator",
+      (void (nuSQUIDSLV::*)(squids::SU_vector)) &nuSQUIDSLV::Set_LV_Operator
+  );
+  nusquids_lv_class_object->def(
+      "Set_LV_EnergyPower",
+      (void (nuSQUIDSLV::*)(int)) &nuSQUIDSLV::Set_LV_EnergyPower
+  );
+  nusquids_lv_class_object->def(
+      "Set_MixingAngle",
+      (void (nuSQUIDSLV::*)(unsigned int, unsigned int, double)) &nuSQUIDSLV::Set_MixingAngle
+  );
+  nusquids_lv_class_object->def(
+      "Set_CPPhase",
+      (void (nuSQUIDSLV::*)(unsigned int, unsigned int, double)) &nuSQUIDSLV::Set_CPPhase
+  );
+  nusquids_lv_class_object->def(
+      "dump_probabilities",
+      (void (nuSQUIDSLV::*)() const) &nuSQUIDSLV::dump_probabilities
+  );
+  nusquids_lv_class_object->def(
+      "Set_initial_state",
+      (void (nuSQUIDSLV::*)(const marray<double, 1>&, Basis)) &nuSQUIDSLV::Set_initial_state,
+      (bp::arg("ini_state"), bp::arg("basis") = flavor)
+  );
+  nusquids_lv_class_object->def(
+      "Set_initial_state",
+      (void (nuSQUIDSLV::*)(const marray<double, 2>&, Basis)) &nuSQUIDSLV::Set_initial_state,
+      (bp::arg("ini_state"), bp::arg("basis") = flavor)
+  );
+  nusquids_lv_class_object->def(
+      "Set_initial_state",
+      (void (nuSQUIDSLV::*)(const marray<double, 3>&, Basis)) &nuSQUIDSLV::Set_initial_state,
+      (bp::arg("ini_state"), bp::arg("basis") = flavor)
+  );
+
+  // Bind nuSQUIDSLVAtm
+  // Ideally this would be done with RegisterBasicAtmNuSQuIDSPythonBindings. I tried this and it compiles.
+  // However, if I afterwards define LV-specific methods not defined in RegisterBasicAtmNuSQuIDSPythonBindings 
+  // I get an error at runtime in python. The error complains that the first argument to the method, which is
+  // in python the 'self', is not recognized as a valid C++ nuSQUIDSLVAtm class. Below is a workaround which
+  // is just more verbose than using RegisterBasicAtmNuSQuIDSPythonBindings.
+
+  class_<nuSQUIDSLVAtm, boost::noncopyable, std::shared_ptr<nuSQUIDSLVAtm> >("nuSQUIDSLVAtm", no_init)
+    .def(init<marray<double,1>,marray<double,1>,unsigned int,NeutrinoType>(args("CosZenith_vector","E_vector","numneu","NT")))
+    .def(init<marray<double,1>,marray<double,1>,unsigned int,NeutrinoType,bool>(args("CosZenith_vector","E_vector","numneu","NT","iinteraction")))
+    .def(init<marray<double,1>,marray<double,1>,unsigned int,NeutrinoType,bool,std::shared_ptr<CrossSectionLibrary>>(args("CosZenith_vector","E_vector","numneu","NT","iinteraction","ncs")))
+    .def(init<std::string>(args("filename")))
+    .def("EvolveState",&nuSQUIDSLVAtm::EvolveState)
+    .def("Set_TauRegeneration",&nuSQUIDSLVAtm::Set_TauRegeneration)
+    .def("EvalFlavor",(double(nuSQUIDSLVAtm::*)(unsigned int,double,double,unsigned int,bool) const)&nuSQUIDSLVAtm::EvalFlavor,
+        nuSQUIDSAtm_EvalFlavor_overload<nuSQUIDSLVAtm>(args("Flavor","cos(theta)","Neutrino Energy","NeuType","BoolToRandomzeProdutionHeight"),
+          "nuSQuIDSAtm evaluate flux.."))
+    .def("Set_EvalThreads",&nuSQUIDSLVAtm::Set_EvalThreads)
+    .def("Get_EvalThreads",&nuSQUIDSLVAtm::Get_EvalThreads)
+    .def("Set_EarthModel",&nuSQUIDSLVAtm::Set_EarthModel)
+    .def("WriteStateHDF5",&nuSQUIDSLVAtm::WriteStateHDF5)
+    .def("ReadStateHDF5",&nuSQUIDSLVAtm::ReadStateHDF5)
+    .def("Set_MixingAngle",&nuSQUIDSLVAtm::Set_MixingAngle)
+    .def("Get_MixingAngle",&nuSQUIDSLVAtm::Get_MixingAngle)
+    .def("Set_CPPhase",&nuSQUIDSLVAtm::Set_CPPhase)
+    .def("Get_CPPhase",&nuSQUIDSLVAtm::Get_CPPhase)
+    .def("Set_SquareMassDifference",&nuSQUIDSLVAtm::Set_SquareMassDifference)
+    .def("Get_SquareMassDifference",&nuSQUIDSLVAtm::Get_SquareMassDifference)
+    .def("Set_h",(void(nuSQUIDSLVAtm::*)(double))&nuSQUIDSLVAtm::Set_h)
+    .def("Set_h",(void(nuSQUIDSLVAtm::*)(double,unsigned int))&nuSQUIDSLVAtm::Set_h)
+    .def("Set_h_max",(void(nuSQUIDSLVAtm::*)(double))&nuSQUIDSLVAtm::Set_h_max)
+    .def("Set_h_max",(void(nuSQUIDSLVAtm::*)(double,unsigned int))&nuSQUIDSLVAtm::Set_h_max)
+    .def("Set_h_min",(void(nuSQUIDSLVAtm::*)(double))&nuSQUIDSLVAtm::Set_h_min)
+    .def("Set_h_min",(void(nuSQUIDSLVAtm::*)(double,unsigned int))&nuSQUIDSLVAtm::Set_h_min)
+    .def("Set_ProgressBar",&nuSQUIDSLVAtm::Set_ProgressBar)
+    .def("Set_MixingParametersToDefault",&nuSQUIDSLVAtm::Set_MixingParametersToDefault)
+    .def("Set_GSL_step",wrap_nusqatm_Set_GSL_STEP<nuSQUIDSLV>)
+    .def("Set_rel_error",(void(nuSQUIDSLVAtm::*)(double))&nuSQUIDSLVAtm::Set_rel_error)
+    .def("Set_rel_error",(void(nuSQUIDSLVAtm::*)(double, unsigned int))&nuSQUIDSLVAtm::Set_rel_error)
+    .def("Set_abs_error",(void(nuSQUIDSLVAtm::*)(double))&nuSQUIDSLVAtm::Set_abs_error)
+    .def("Set_abs_error",(void(nuSQUIDSLVAtm::*)(double, unsigned int))&nuSQUIDSLVAtm::Set_abs_error)
+    .def("Set_EvolLowPassCutoff",&nuSQUIDSLVAtm::Set_EvolLowPassCutoff)
+    .def("Set_EvolLowPassScale",&nuSQUIDSLVAtm::Set_EvolLowPassScale)
+    .def("GetNumE",&nuSQUIDSLVAtm::GetNumE)
+    .def("GetNumCos",&nuSQUIDSLVAtm::GetNumCos)
+    .def("GetNumNeu",&nuSQUIDSLVAtm::GetNumNeu)
+    .def("GetNumRho",&nuSQUIDSLVAtm::GetNumRho)
+    .def("GetnuSQuIDS",(std::vector<nuSQUIDSLV>&(nuSQUIDSLVAtm::*)())&nuSQUIDSLVAtm::GetnuSQuIDS,boost::python::return_internal_reference<>())
+    .def("GetnuSQuIDS",(nuSQUIDSLV&(nuSQUIDSLVAtm::*)(unsigned int))&nuSQUIDSLVAtm::GetnuSQuIDS,boost::python::return_internal_reference<>())
+    .def("Set_initial_state",(void(nuSQUIDSLVAtm::*)(const marray<double,3>&, Basis))&nuSQUIDSLVAtm::Set_initial_state,nuSQUIDSAtm_Set_initial_state<nuSQUIDSLVAtm>())
+    .def("Set_initial_state",(void(nuSQUIDSLVAtm::*)(const marray<double,4>&, Basis))&nuSQUIDSLVAtm::Set_initial_state,nuSQUIDSAtm_Set_initial_state<nuSQUIDSLVAtm>())
+    .def("GetStates", (marray<double,2>(nuSQUIDSLVAtm::*)(unsigned int))&nuSQUIDSLVAtm::GetStates,
+      nuSQUIDSAtm_GetStates_overload<nuSQUIDSLVAtm>(args("rho"), "Get evolved states of all nodes"))
+    .def("GetERange",&nuSQUIDSLVAtm::GetERange)
+    .def("GetCosthRange",&nuSQUIDSLVAtm::GetCosthRange)
+    .def("Set_IncludeOscillations",&nuSQUIDSLVAtm::Set_IncludeOscillations)
+    .def("Set_GlashowResonance",&nuSQUIDSLVAtm::Set_GlashowResonance)
+    .def("Set_TauRegeneration",&nuSQUIDSLVAtm::Set_TauRegeneration)
+    .def("Set_AllowConstantDensityOscillationOnlyEvolution",&nuSQUIDSLVAtm::Set_AllowConstantDensityOscillationOnlyEvolution)
+    .def("Set_PositivyConstrain",&nuSQUIDSLVAtm::Set_PositivityConstrain)
+    .def("Set_PositivyConstrainStep",&nuSQUIDSLVAtm::Set_PositivityConstrainStep)
+    .def("Get_EvalThreads",&nuSQUIDSLVAtm::Get_EvalThreads)
+    .def("Set_EvalThreads",&nuSQUIDSLVAtm::Set_EvalThreads)
+    .def("Set_EarthModel",&nuSQUIDSLVAtm::Set_EarthModel)
+    .def("SetNeutrinoCrossSections",&nuSQUIDSLVAtm::SetNeutrinoCrossSections)
+    .def("GetNeutrinoCrossSections",&nuSQUIDSLVAtm::GetNeutrinoCrossSections)
+    .def(
+        "Set_LV_OpMatrix",
+        (void (nuSQUIDSLVAtm::*)(nusquids::LVParameters&)) &nuSQUIDSLVAtm::Set_LV_OpMatrix
+    )
+    .def(
+        "Set_LV_OpMatrix",
+        (void (nuSQUIDSLVAtm::*)(gsl_matrix_complex*)) &nuSQUIDSLVAtm::Set_LV_OpMatrix
+    )
+    .def(
+        "Set_LV_Operator",
+        (void (nuSQUIDSLVAtm::*)(squids::SU_vector)) &nuSQUIDSLVAtm::Set_LV_Operator
+    )
+    .def(
+        "Set_LV_EnergyPower",
+        (void (nuSQUIDSLVAtm::*)(int)) &nuSQUIDSLVAtm::Set_LV_EnergyPower
+    );
+
 }