add broken powerlaw

include/crpropa/diffusionTensor/DiffusionTensor.h

@@ -85,6 +85,33 @@ class DiffusionBrokenPowerlaw: public DiffusionTensor {
 };
 
 
+class DiffusionMultipleBreaks: public DiffusionTensor {
+private: 
+    double initialIndex;            //< power law index before the first break
+    double norm;                    //< normalisation at a given energy
+    double Enorm;                   //< normalisation energy
+    double eps;                     //< anisotropy of the diffusion tensor
+
+    std::vector<double> index_s;    //< power law slope after break s
+    std::vector<double> break_s;    //< energy of the break
+    std::vector<double> norm_s;  //< normalisation at a given break 
+
+public:
+    DiffusionMultipleBreaks(double initialIndex, double eps, double norm = 6.1e24 , double Enorm = 4 * GeV);
+
+    double getDiffusionCoefficient(double E) const;
+    Vector3d getTensorDiagonal(double E, int id = 0, double B = 0) const; 
+
+    /**
+        add a break to the diffusion tensor
+        @param index    power law slope after the break
+        @param E        Energy of the break
+    */
+    void addBreak(double index, double E);
+
+
+};
+
 class DiffusionCRINGE : public DiffusionTensor {
 private: 
     const double norm = 5.18e24; // norm of the diffusion tensor in m2/s
@@ -101,6 +128,40 @@ class DiffusionCRINGE : public DiffusionTensor {
 
 };
 
+class DiffusionBPLsoft : public DiffusionTensor {
+  private: 
+    double norm = 7.03e25; // norm of the diffusion tensor in m2/s
+    double breakEnergy = 63 * GeV; // energy of the break
+    double index1 = 0.334; // spectral index below the break
+    double index2 = -0.324; // spectral index above the break
+    double epsilon = 0.1; // anisotropy of the tensor
+
+  public:
+
+    DiffusionBPLsoft(double norm = 7.03e25, double breakEnergy = 63 * GeV, double index1 = 0.334, double index2 = -0.324, double epsilon = 0.1);
+
+    Vector3d getTensorDiagonal(double Energy, int id = 0, double B = 0) const;
+
+    std::string getDescription() {
+        return "Diffusion Tensor with a broken power law with a soft break.\n";
+    }
+
+    void setEpsilon(double eps);
+    double getEpsilon() const;
+
+    void setIndex1(double a);
+    double getIndex1() const;
+
+    void setIndex2(double a);
+    double getIndex2() const;
+
+    void setBreakEnergy(double E);
+    double getBreakEnergy() const;
+
+    void setNorm(double n);
+    double getNorm() const;
+};
+
 } // namespace crpropa
 
 #endif // CRPROPA_DIFFUSIONTENSOR_H

src/diffusionTensor/DiffusionTensor.cpp

@@ -91,6 +91,47 @@ double DiffusionBrokenPowerlaw::getEpsilon() const {
     return epsilon;
 }
 
+// ------------------------------------------------------------------------
+
+DiffusionMultipleBreaks::DiffusionMultipleBreaks(double index, double eps, double norm, double Enorm) : 
+    initialIndex(index), eps(eps), norm(norm), Enorm(Enorm) {
+        index_s.push_back(index);
+        break_s.push_back(Enorm);
+        norm_s.push_back(norm);
+    }
+
+double DiffusionMultipleBreaks::getDiffusionCoefficient(double E) const {
+    // find break iB = 0 is the first reference
+    int iB = 0;
+
+    while((iB + 1 < index_s.size())) { // check if next break exists
+        if (E < break_s[iB + 1]) 
+            break;
+        iB++;
+    }
+
+    double r = E / break_s[iB];
+    double D = norm_s[iB]; 
+    double i = index_s[iB];
+
+
+    return D * pow(r, i);
+}
+
+void DiffusionMultipleBreaks::addBreak(double index, double E) {
+    if (E < break_s.back())
+        throw std::runtime_error("next break must have a larger energy than all breaks before");
+
+    norm_s.push_back(getDiffusionCoefficient(E));
+    index_s.push_back(index);
+    break_s.push_back(E);
+}
+
+Vector3d DiffusionMultipleBreaks::getTensorDiagonal(double E, int id, double B) const {
+    return Vector3d(1, eps, eps) * getDiffusionCoefficient(E);
+}
+
+
 // ------------------------------------------------------------------------
 
 Vector3d DiffusionCRINGE::getTensorDiagonal(double E, int id, double B) const {
@@ -104,4 +145,51 @@ Vector3d DiffusionCRINGE::getTensorDiagonal(double E, int id, double B) const {
     }
 
     return Vector3d(D);
-}
+}
+
+// ------------------------------------------------------------------------
+
+DiffusionBPLsoft::DiffusionBPLsoft(double norm, double breakEnergy, double index1, double index2, double epsilon) : 
+    norm(norm), breakEnergy(breakEnergy), index1(index1), index2(index2), epsilon(epsilon) {}
+
+Vector3d DiffusionBPLsoft::getTensorDiagonal(double E, int id, double B) const {
+    double x = E / breakEnergy;
+    return norm * pow(x, index1) / (1 + pow(x, index2 - index1)) * Vector3d(1, epsilon, epsilon);
+}
+
+void DiffusionBPLsoft::setEpsilon(double eps) {
+    epsilon = eps;
+}
+double DiffusionBPLsoft::getEpsilon() const {
+    return epsilon;
+}
+
+void DiffusionBPLsoft::setNorm(double n) {
+    norm = n;
+}
+double DiffusionBPLsoft::getNorm() const {
+    return norm;
+}
+
+void DiffusionBPLsoft::setBreakEnergy(double E) {
+    breakEnergy = E;
+}
+double DiffusionBPLsoft::getBreakEnergy() const {
+    return breakEnergy;
+}
+
+void DiffusionBPLsoft::setIndex1(double a) {
+    index1 = a;
+}
+double DiffusionBPLsoft::getIndex1() const {
+    return index1;
+}
+
+void DiffusionBPLsoft::setIndex2(double a) {
+    index2 = a;
+}
+double DiffusionBPLsoft::getIndex2() const {
+    return index2;
+}
+
+// ------------------------------------------------------------------------