Program Listing for File JF12Field.h

↰ Return to documentation for file (include/crpropa/magneticField/JF12Field.h)

#ifndef CRPROPA_JF12FIELD_H
#define CRPROPA_JF12FIELD_H

#include "crpropa/magneticField/MagneticField.h"
#include "crpropa/Grid.h"
#include "kiss/logger.h"

namespace crpropa {
class JF12Field: public MagneticField {
protected:
        bool useRegularField;
        bool useStriatedField;
        bool useTurbulentField;
        bool useDiskField;
        bool useToroidalHaloField;
        bool useXField;

        // disk spiral arms
        double rArms[8];       // radii where each arm crosses the negative x-axis
        double pitch;          // pitch angle
        double sinPitch, cosPitch, tanPitch, cotPitch, tan90MinusPitch;

        // Regular field ----------------------------------------------------------
        // disk
        double bDisk[11];      // field strengths of the 8 arms at r=5 kpc; additional entries added for periodic closure in JF12FieldSolenoidal
        double bRing;          // ring field strength 3<r<5 kpc
        double hDisk, wDisk;   // disk/halo transistion and width
        // toroidal halo
        double bNorth, bSouth; // northern, southern halo field strength
        double rNorth, rSouth; // northern, southern transistion radius
        double wHalo, z0;      // transistion width and vertical scale height
        // poloidal halo
        double bX;             // field strength at origin
        double thetaX0;        // constant elevation angle at r > rXc, z = 0
        double sinThetaX0, cosThetaX0, tanThetaX0, cotThetaX0;
        double rXc;            // radius of varying elevation angle region
        double rX;             // exponential scale height

        // Striated field ---------------------------------------------------------
        double sqrtbeta;       // relative strength of striated field
        ref_ptr<Grid1f> striatedGrid;

        // Turbulent field --------------------------------------------------------
        ref_ptr<Grid3f> turbulentGrid;
        // disk
        double bDiskTurb[8]; // field strengths in arms at r=5 kpc
        double bDiskTurb5;   // field strength at r<5kpc
        double zDiskTurb;        // Gaussian scale height of disk
        // halo
        double bHaloTurb; // halo field strength
        double rHaloTurb; // exponential scale length
        double zHaloTurb; // Gaussian scale height

public:
        JF12Field();

        // Create and set a random realization for the striated field
        void randomStriated(int seed = 0);

#ifdef CRPROPA_HAVE_FFTW3F
        // Create a random realization for the turbulent field
        void randomTurbulent(int seed = 0);
#endif

        void setStriatedGrid(ref_ptr<Grid1f> grid);

        void setTurbulentGrid(ref_ptr<Grid3f> grid);

        ref_ptr<Grid1f> getStriatedGrid();
        ref_ptr<Grid3f> getTurbulentGrid();

        void setUseRegularField(bool use);
        virtual void setUseStriatedField(bool use);
        virtual void setUseTurbulentField(bool use);
        void setUseDiskField(bool use);
        void setUseToroidalHaloField(bool use);
        void setUseXField(bool use);

        bool isUsingRegularField();
        bool isUsingStriatedField();
        bool isUsingTurbulentField();
        bool isUsingDiskField();
        bool isUsingToroidalHaloField();
        bool isUsingXField();

        double logisticFunction(const double& x, const double& x0, const double& w) const;

        // Regular field components
        Vector3d getRegularField(const Vector3d& pos) const;
        virtual Vector3d getDiskField(const double& r, const double& z, const double& phi, const double& sinPhi, const double& cosPhi) const;
        Vector3d getToroidalHaloField(const double& r, const double& z, const double& sinPhi, const double& cosPhi) const;
        virtual Vector3d getXField(const double& r, const double& z, const double& sinPhi, const double& cosPhi) const;

        // Regular and striated field component
        Vector3d getStriatedField(const Vector3d& pos) const;

        // Brms of the turbulent field
        double getTurbulentStrength(const Vector3d& pos) const;

        // Turbulent field component
        Vector3d getTurbulentField(const Vector3d& pos) const;

        // All set field components
        Vector3d getField(const Vector3d& pos) const;
};
class PlanckJF12bField: public JF12Field {
        public:
                PlanckJF12bField();
};
} // namespace crpropa

#endif // CRPROPA_JF12FIELD_H