Merge branch 'master' into release

Update the relase branch.

.travis.yml

@@ -66,6 +66,7 @@ matrix:
         - CRPROPA_DATA_PATH=$TRAVIS_BUILD_DIR/build/data/
 
 before_install:
+  - echo "TRAVIS_OS_NAME=$TRAVIS_OS_NAME"
   - |
     if [[ "$TRAVIS_OS_NAME" == "osx" ]]; then
       HOMEBREW_NO_AUTO_UPDATE=1 brew install swig;
@@ -85,7 +86,7 @@ before_script:
   - cmake .. -DENABLE_PYTHON=True -DPYTHON_EXECUTABLE=$PYTHON_EXECUTABLE -DENABLE_TESTING=On
   - cp data.tar.gz $HOME/crpropa_cache
 script:
-  - make
+  - VERBOSE=1 make
   - make test
 after_failure:
   - cat Testing/Temporary/LastTest.log

CHANGELOG.md

@@ -3,6 +3,8 @@
 ### Bug fixes:
 
 ### New features:
+  * Vector3 now allows has index based access to its components. This improves
+    interoperability with third party libraries, in particular with numpy arrays.
 
 ### Features that are deprecated and will be removed after this release:
 

CMakeLists.txt

@@ -229,8 +229,8 @@ if (ENABLE_QUIMBY)
 		list(APPEND CRPROPA_EXTRA_LIBRARIES ${QUIMBY_LIBRARY})
 		add_definitions (-DCRPROPA_HAVE_QUIMBY)
 		list(APPEND CRPROPA_SWIG_DEFINES -DCRPROPA_HAVE_QUIMBY)
-		list(APPEND CRPROPA_SWIG_DEFINES -I${QUIMBY_INCLUDE_DIR}/../share/quimby)
-		list(APPEND CRPROPA_SWIG_DEFINES -I${QUIMBY_INCLUDE_DIR})
+		list(APPEND SWIG_INCLUDE_DIRECTORIES ${QUIMBY_INCLUDE_DIR}/../share/quimby)
+		list(APPEND SWIG_INCLUDE_DIRECTORIES ${QUIMBY_INCLUDE_DIR})
 	endif(QUIMBY_FOUND)
 endif(ENABLE_QUIMBY)
 
@@ -241,8 +241,8 @@ if(SAGA_FOUND)
 	list(APPEND CRPROPA_EXTRA_LIBRARIES ${SAGA_LIBRARY})
 	add_definitions (-DCRPROPA_HAVE_SAGA)
 	list(APPEND CRPROPA_SWIG_DEFINES -DCRPROPA_HAVE_SAGA)
-	list(APPEND CRPROPA_SWIG_DEFINES -I${SAGA_INCLUDE_DIR}/../python)
-	list(APPEND CRPROPA_SWIG_DEFINES -I${SAGA_INCLUDE_DIR})
+	list(APPEND SWIG_INCLUDE_DIRECTORIES ${SAGA_INCLUDE_DIR}/../python)
+	list(APPEND SWIG_INCLUDE_DIRECTORIES ${SAGA_INCLUDE_DIR})
 endif(SAGA_FOUND)
 
 # muparser (optional for generic source spectra)
@@ -252,7 +252,7 @@ if(MUPARSER_FOUND)
 	list(APPEND CRPROPA_EXTRA_LIBRARIES ${MUPARSER_LIBRARY})
 	add_definitions (-DCRPROPA_HAVE_MUPARSER)
 	list(APPEND CRPROPA_SWIG_DEFINES -DCRPROPA_HAVE_MUPARSER)
-	list(APPEND CRPROPA_SWIG_DEFINES -I${MUPARSER_INCLUDE_DIR})
+	list(APPEND SWIG_INCLUDE_DIRECTORIES ${MUPARSER_INCLUDE_DIR})
 endif(MUPARSER_FOUND)
 
 # zlib / zstream-cpp (optional for compressed output files)
@@ -263,7 +263,7 @@ if(ZLIB_FOUND)
 	list(APPEND CRPROPA_EXTRA_LIBRARIES ${ZLIB_LIBRARIES})
 	add_definitions (-DCRPROPA_HAVE_ZLIB)
 	list(APPEND CRPROPA_SWIG_DEFINES -DCRPROPA_HAVE_ZLIB)
-	list(APPEND CRPROPA_SWIG_DEFINES -I${ZLIB_INCLUDE_DIRS})
+	list(APPEND SWIG_INCLUDE_DIRECTORIES ${ZLIB_INCLUDE_DIRS})
 endif(ZLIB_FOUND)
 
 # HDF5 (optional for HDF5 output files)
@@ -277,7 +277,9 @@ if(ENABLE_HDF5)
 			list(APPEND CRPROPA_EXTRA_LIBRARIES ${HDF5_LIBRARIES})
 			add_definitions (-DCRPROPA_HAVE_HDF5)
 			list(APPEND CRPROPA_SWIG_DEFINES -DCRPROPA_HAVE_HDF5)
-			list(APPEND CRPROPA_SWIG_DEFINES -I${HDF5_INCLUDE_DIRS})
+			list(APPEND SWIG_INCLUDE_DIRECTORIES ${HDF5_INCLUDE_DIRS})
+			#string(REPLACE " " " -I" HDF5_INCLUDE_DIRS_SWIG ${HDF5_INCLUDE_DIRS})
+			#list(APPEND CRPROPA_SWIG_DEFINES -I${HDF5_INCLUDE_DIRS_SWIG})
 		endif(NOT HDF5_IS_PARALLEL)
 	endif(HDF5_FOUND)
 endif(ENABLE_HDF5)
@@ -393,6 +395,7 @@ add_library(crpropa SHARED
 	src/magneticField/MagneticField.cpp
 	src/magneticField/MagneticFieldGrid.cpp
 	src/magneticField/PT11Field.cpp
+    src/magneticField/TF17Field.cpp
 	src/magneticField/ArchimedeanSpiralField.cpp
 	src/advectionField/AdvectionField.cpp
 	src/massDistribution/Nakanishi.cpp

include/CRPropa.h

@@ -61,6 +61,7 @@
 #include "crpropa/magneticField/MagneticField.h"
 #include "crpropa/magneticField/MagneticFieldGrid.h"
 #include "crpropa/magneticField/PT11Field.h"
+#include "crpropa/magneticField/TF17Field.h"
 #include "crpropa/magneticField/QuimbyMagneticField.h"
 #include "crpropa/magneticField/ArchimedeanSpiralField.h"
 

include/crpropa/Common.h

@@ -3,7 +3,6 @@
 
 #include <string>
 #include <vector>
-
 /**
  @file
  @brief Common helper functions
@@ -67,7 +66,21 @@ inline double pow_integer<0>(double base)
   return 1;
 }
 
-
+// - input:  function over which to integrate, integration limits A and B
+// - output: 8-points Gauß-Legendre integral
+static const double X[8] = {.0950125098, .2816035507, .4580167776, .6178762444, .7554044083, .8656312023, .9445750230, .9894009349};
+static const double W[8] = {.1894506104, .1826034150, .1691565193, .1495959888, .1246289712, .0951585116, .0622535239, .0271524594};
+template<typename Integrand>
+double gaussInt(Integrand&& integrand, double A, double B) {
+	const double XM = 0.5 * (B + A);
+	const double XR = 0.5 * (B - A);
+	double SS = 0.;
+	for (int i = 0; i < 8; ++i) {
+		double DX = XR * X[i];
+		SS += W[i] * (integrand(XM + DX) + integrand(XM - DX));
+	}
+	return XR * SS;
+}
 
 } // namespace crpropa
 

include/crpropa/PhotonBackground.h

@@ -1,6 +1,8 @@
 #ifndef CRPROPA_PHOTONBACKGROUND_H
 #define CRPROPA_PHOTONBACKGROUND_H
 
+#include "crpropa/Common.h"
+
 #include <string>
 
 namespace crpropa {
@@ -31,7 +33,68 @@ double photonFieldScaling(PhotonField photonField, double z);
 // Returns a string representation of the field
 std::string photonFieldName(PhotonField photonField);
 
+
+/**
+ @class photonFieldSampling
+ @brief Reimplementation of SOPHIA photon sampling. Naming and unit conventions are taken from SOPHIA to ease comparisions.
+ */
+class PhotonFieldSampling {
+public:
+	PhotonFieldSampling();
+
+	/**
+	 Constructor to mimic SOPHIA structure.
+	  @param bgFlag		1: CMB | 2: IRB_Kneiske04
+	 */
+	explicit PhotonFieldSampling(int bgFlag);
+
+	/**
+	 SOPHIA's photon sampling method. Returns energy [J] of a photon of the photon field.
+	 @param onProton	particle type: proton or neutron
+	 @param E_in		energy of incoming nucleon
+	 @param z_in		redshift of incoming nucleon
+	 */
+	double sample_eps(bool onProton, double E_in, double z_in) const;
+protected:
+	int bgFlag;
+
+	// called by: sample_eps
+	// - input: photon energy [eV], redshift
+	// - output: photon density per unit energy [#/(eVcm^3)]
+	double getPhotonDensity(double eps, double z_in) const;
+
+	// called by: sample_eps  
+	// - input: s [GeV^2] 
+	// - output: (s-p^2) * sigma_(nucleon/gamma) [GeV^2 * mubarn]
+	double functs(double s, bool onProton) const;
+
+	// called by: sample_eps, gaussInt
+	// - input: photon energy eps [eV], E_in [GeV]
+	// - output: probability to encounter photon of energy eps
+	double prob_eps(double eps, bool onProton, double E_in, double z_in) const;
+
+	// called by: functs
+	// - input: photon energy [eV]
+	// - output: crossection of nucleon-photon-interaction [mubarn]
+	double crossection(double eps, bool onProton) const;
+
+	// called by: crossection
+	// - input: photon energy [eV], threshold [eV], max [eV], unknown [no unit]
+	// - output: unknown [no unit]
+	double Pl(double x, double xth, double xmax, double alpha) const;
+
+	// called by: crossection
+	// - input: photon energy [eV], threshold [eV], unknown [eV]
+	// - output: unknown [no unit]
+	double Ef(double x, double th, double w) const;
+
+	// called by: crossection
+	// - input: cross section [µbarn], width [GeV], mass [GeV/c^2], rest frame photon energy [GeV]
+	// - output: Breit-Wigner crossection of a resonance of width Gamma
+	double breitwigner(double sigma_0, double Gamma, double DMM, double epsPrime, bool onProton) const;
+};
 /** @}*/
+
 } // namespace crpropa
 
 #endif // CRPROPA_PHOTONBACKGROUND_H