Program Listing for File DiffusionSDE.h

↰ Return to documentation for file (include/crpropa/module/DiffusionSDE.h)

#ifndef CRPROPA_DIFFUSIONSDE_H
#define CRPROPA_DIFFUSIONSDE_H

#include <iostream>
#include <vector>
#include <cmath>
#include <string>
#include <cstdlib>
#include <stdexcept>

#include "crpropa/Module.h"
#include "crpropa/magneticField/MagneticField.h"
#include "crpropa/advectionField/AdvectionField.h"
#include "crpropa/Units.h"
#include "crpropa/Random.h"

#include "kiss/logger.h"

namespace crpropa {
class DiffusionSDE : public Module{

private:
            ref_ptr<MagneticField> magneticField;
            ref_ptr<AdvectionField> advectionField;
            double minStep; // minStep/c_light is the minimum integration timestep
            double maxStep; // maxStep/c_light is the maximum integration timestep
            double tolerance; // tolerance is criterion for step adjustment. Step adjustment takes place when the tangential vector of the magnetic field line is calculated.
            double epsilon; // ratio of parallel and perpendicular diffusion coefficient D_par = epsilon*D_perp
            double alpha; // power law index of the energy dependent diffusion coefficient: D\propto E^alpha
            double scale; // scaling factor for the diffusion coefficient D = scale*D_0

public:
        DiffusionSDE(ref_ptr<crpropa::MagneticField> magneticField, double tolerance = 1e-4, double minStep = 10 * pc, double maxStep = 1 * kpc, double epsilon = 0.1);
        DiffusionSDE(ref_ptr<crpropa::MagneticField> magneticField, ref_ptr<crpropa::AdvectionField> advectionField, double tolerance = 1e-4, double minStep = 10 * pc, double maxStep = 1 * kpc, double epsilon = 0.1);

        void process(crpropa::Candidate *candidate) const;

        void tryStep(const Vector3d &Pos, Vector3d &POut, Vector3d &PosErr, double z, double propStep ) const;
        void driftStep(const Vector3d &Pos, Vector3d &LinProp, double h) const;
        void calculateBTensor(double rig, double BTen[], Vector3d pos, Vector3d dir, double z) const;

        void setMinimumStep(double minStep);
        void setMaximumStep(double maxStep);
        void setTolerance(double tolerance);
        void setEpsilon(double kappa);
        void setAlpha(double alpha);
        void setScale(double Scale);
        void setMagneticField(ref_ptr<crpropa::MagneticField> magneticField);
        void setAdvectionField(ref_ptr<crpropa::AdvectionField> advectionField);

        double getMinimumStep() const;
        double getMaximumStep() const;
        double getTolerance() const;
        double getEpsilon() const;
        double getAlpha() const;
        double getScale() const;
        std::string getDescription() const;

  ref_ptr<MagneticField> getMagneticField() const;
        Vector3d getMagneticFieldAtPosition(Vector3d pos, double z) const;
        ref_ptr<AdvectionField> getAdvectionField() const;
        Vector3d getAdvectionFieldAtPosition(Vector3d pos) const;

};
} //namespace crpropa

#endif // CRPROPA_DIFFUSIONSDE_H