LCOV - code coverage report
Current view: top level - src/magneticField - JF12Field.cpp (source / functions) Hit Total Coverage
Test: coverage.info.cleaned Lines: 0 231 0.0 %
Date: 2024-04-29 14:43:01 Functions: 0 29 0.0 %

          Line data    Source code
       1             : #include "crpropa/magneticField/JF12Field.h"
       2             : #include "crpropa/Units.h"
       3             : #include "crpropa/magneticField/turbulentField/SimpleGridTurbulence.h"
       4             : #include "crpropa/Random.h"
       5             : 
       6             : namespace crpropa {
       7             : 
       8           0 : JF12Field::JF12Field() {
       9           0 :         useRegularField = true;
      10           0 :         useStriatedField = false;
      11           0 :         useTurbulentField = false;
      12           0 :         useDiskField = true;
      13           0 :         useToroidalHaloField = true;
      14           0 :         useXField = true;
      15             : 
      16             :         // spiral arm parameters
      17           0 :         pitch = 11.5 * M_PI / 180;
      18           0 :         sinPitch = sin(pitch);
      19           0 :         cosPitch = cos(pitch);
      20           0 :         tanPitch = tan(pitch);
      21           0 :         cotPitch =  1. / tanPitch;
      22           0 :         tan90MinusPitch = tan(M_PI / 2 - pitch);
      23             : 
      24           0 :         rArms[0] = 5.1 * kpc;
      25           0 :         rArms[1] = 6.3 * kpc;
      26           0 :         rArms[2] = 7.1 * kpc;
      27           0 :         rArms[3] = 8.3 * kpc;
      28           0 :         rArms[4] = 9.8 * kpc;
      29           0 :         rArms[5] = 11.4 * kpc;
      30           0 :         rArms[6] = 12.7 * kpc;
      31           0 :         rArms[7] = 15.5 * kpc;
      32             : 
      33             :         // regular field parameters
      34           0 :         bRing = 0.1 * muG;
      35           0 :         hDisk = 0.40 * kpc;
      36           0 :         wDisk = 0.27 * kpc;
      37             : 
      38           0 :         bDisk[0] = 0.1 * muG;
      39           0 :         bDisk[1] = 3.0 * muG;
      40           0 :         bDisk[2] = -0.9 * muG;
      41           0 :         bDisk[3] = -0.8 * muG;
      42           0 :         bDisk[4] = -2.0 * muG;
      43           0 :         bDisk[5] = -4.2 * muG;
      44           0 :         bDisk[6] = 0.0 * muG;
      45           0 :         bDisk[7] = 2.7 * muG;
      46             : 
      47           0 :         bNorth = 1.4 * muG;
      48           0 :         bSouth = -1.1 * muG;
      49           0 :         rNorth = 9.22 * kpc;
      50           0 :         rSouth = 17 * kpc;
      51           0 :         wHalo = 0.20 * kpc;
      52           0 :         z0 = 5.3 * kpc;
      53             : 
      54           0 :         bX = 4.6 * muG;
      55           0 :         thetaX0 = 49.0 * M_PI / 180;
      56           0 :         sinThetaX0 = sin(thetaX0);
      57           0 :         cosThetaX0 = cos(thetaX0);
      58           0 :         tanThetaX0 = tan(thetaX0);
      59           0 :         cotThetaX0 = 1. / tanThetaX0;
      60           0 :         rXc = 4.8 * kpc;
      61           0 :         rX = 2.9 * kpc;
      62             : 
      63             :         // striated field parameter
      64           0 :         sqrtbeta = sqrt(1.36);
      65             : 
      66             :         // turbulent field parameters
      67           0 :         bDiskTurb[0] = 10.81 * muG;
      68           0 :         bDiskTurb[1] = 6.96 * muG;
      69           0 :         bDiskTurb[2] = 9.59 * muG;
      70           0 :         bDiskTurb[3] = 6.96 * muG;
      71           0 :         bDiskTurb[4] = 1.96 * muG;
      72           0 :         bDiskTurb[5] = 16.34 * muG;
      73           0 :         bDiskTurb[6] = 37.29 * muG;
      74           0 :         bDiskTurb[7] = 10.35 * muG;
      75             : 
      76           0 :         bDiskTurb5 = 7.63 * muG;
      77           0 :         zDiskTurb = 0.61 * kpc;
      78             : 
      79           0 :         bHaloTurb = 4.68 * muG;
      80           0 :         rHaloTurb = 10.97 * kpc;
      81           0 :         zHaloTurb = 2.84 * kpc;
      82           0 : }
      83             : 
      84           0 : void JF12Field::randomStriated(int seed) {
      85           0 :         useStriatedField = true;
      86             :         int N = 100;
      87           0 :         striatedGrid = new Grid1f(Vector3d(0.), N, 0.1 * kpc);
      88             : 
      89           0 :         Random random;
      90           0 :         if (seed != 0)
      91           0 :                 random.seed(seed);
      92             : 
      93           0 :         for (int ix = 0; ix < N; ix++)
      94           0 :                 for (int iy = 0; iy < N; iy++)
      95           0 :                         for (int iz = 0; iz < N; iz++) {
      96           0 :                                 float &f = striatedGrid->get(ix, iy, iz);
      97           0 :                                 f = round(random.rand()) * 2 - 1;
      98             :                         }
      99           0 : }
     100             : 
     101             : #ifdef CRPROPA_HAVE_FFTW3F
     102           0 : void JF12Field::randomTurbulent(int seed) {
     103           0 :         useTurbulentField = true;
     104             :         // turbulent field with Kolmogorov spectrum, B_rms = 1 (will be scaled) and Lc = 60 parsec, and 256 grid points.
     105             :         // Note that the inertial range of the turbulence is less than 2 orders of magnitude.
     106             :     const double lMin = 8 * parsec;
     107             :     const double lMax = 272 * parsec;
     108             :     const double Brms = 1;
     109             :     const double spacing = 4 * parsec;
     110             :     const double grid_n = 256;
     111             : 
     112             :     auto spectrum = SimpleTurbulenceSpectrum(Brms, lMin, lMax);
     113             :     auto gp = GridProperties(Vector3d(0.), grid_n, spacing);
     114           0 :     auto tf = SimpleGridTurbulence(spectrum, gp, seed);
     115           0 :     turbulentGrid = tf.getGrid();
     116             : 
     117           0 : }
     118             : #endif
     119             : 
     120           0 : void JF12Field::setStriatedGrid(ref_ptr<Grid1f> grid) {
     121           0 :         useStriatedField = true;
     122           0 :         striatedGrid = grid;
     123           0 : }
     124             : 
     125           0 : void JF12Field::setTurbulentGrid(ref_ptr<Grid3f> grid) {
     126           0 :         useTurbulentField = true;
     127           0 :         turbulentGrid = grid;
     128           0 : }
     129             : 
     130           0 : ref_ptr<Grid1f> JF12Field::getStriatedGrid() {
     131           0 :         return striatedGrid;
     132             : }
     133             : 
     134           0 : ref_ptr<Grid3f> JF12Field::getTurbulentGrid() {
     135           0 :         return turbulentGrid;
     136             : }
     137             : 
     138           0 : void JF12Field::setUseRegularField(bool use) {
     139           0 :         useRegularField = use;
     140           0 : }
     141             : 
     142           0 : void JF12Field::setUseDiskField(bool use) {
     143           0 :         useDiskField = use;
     144           0 : }
     145             : 
     146           0 : void JF12Field::setUseToroidalHaloField(bool use) {
     147           0 :         useToroidalHaloField = use;
     148           0 : }
     149             : 
     150           0 : void JF12Field::setUseXField(bool use) {
     151           0 :         useXField = use;
     152           0 : }
     153             : 
     154           0 : void JF12Field::setUseStriatedField(bool use) {
     155           0 :         if ((use) and !(striatedGrid)) {
     156           0 :                 KISS_LOG_WARNING << "JF12Field: No striated field set: ignored. Run e.g. randomStriated().";
     157           0 :                 return;
     158             :         }
     159           0 :         useStriatedField = use;
     160             : }
     161             : 
     162           0 : void JF12Field::setUseTurbulentField(bool use) {
     163           0 :         if ((use) and !(turbulentGrid)) {
     164           0 :                 KISS_LOG_WARNING << "JF12Field: No turbulent field set: ignored. Run e.g. randomTurbulent().";
     165           0 :                 return;
     166             :         }
     167           0 :         useTurbulentField = use;
     168             : }
     169             : 
     170           0 : bool JF12Field::isUsingRegularField() {
     171           0 :         return useRegularField;
     172             : }
     173             : 
     174           0 : bool JF12Field::isUsingDiskField() {
     175           0 :         return useDiskField;
     176             : }
     177             : 
     178           0 : bool JF12Field::isUsingToroidalHaloField() {
     179           0 :         return useToroidalHaloField;
     180             : }
     181             : 
     182           0 : bool JF12Field::isUsingXField() {
     183           0 :         return useXField;
     184             : }
     185             : 
     186           0 : bool JF12Field::isUsingStriatedField() {
     187           0 :         return useStriatedField;
     188             : }
     189             : 
     190           0 : bool JF12Field::isUsingTurbulentField() {
     191           0 :         return useTurbulentField;
     192             : }
     193             : 
     194           0 : double JF12Field::logisticFunction(const double& x, const double& x0, const double& w) const {
     195           0 :         return 1. / (1. + exp(-2. * (fabs(x) - x0) / w));
     196             : }
     197             : 
     198           0 : Vector3d JF12Field::getRegularField(const Vector3d& pos) const {
     199             :         Vector3d b(0.);
     200             : 
     201             :         double d = pos.getR(); // distance to galactic center
     202             : 
     203           0 :         if (d < 20 * kpc) {
     204           0 :                 double r = sqrt(pos.x * pos.x + pos.y * pos.y); // in-plane radius
     205           0 :                 double phi = pos.getPhi(); // azimuth
     206           0 :                 double sinPhi = sin(phi);
     207           0 :                 double cosPhi = cos(phi);
     208             : 
     209           0 :                 b += getDiskField(r, pos.z, phi, sinPhi, cosPhi);
     210           0 :                 b += getToroidalHaloField(r, pos.z, sinPhi, cosPhi);
     211           0 :                 b += getXField(r, pos.z, sinPhi, cosPhi);
     212             :         }
     213             : 
     214           0 :         return b;
     215             : }
     216             : 
     217           0 : Vector3d JF12Field::getDiskField(const double& r, const double& z, const double& phi, const double& sinPhi, const double& cosPhi) const {
     218             :         Vector3d b(0.);
     219           0 :         if (useDiskField) {
     220           0 :                 double lfDisk = logisticFunction(z, hDisk, wDisk);
     221           0 :                 if (r > 3 * kpc) {
     222             :                         double bMag;
     223           0 :                         if (r < 5 * kpc) {
     224             :                                 // molecular ring
     225           0 :                                 bMag = bRing * (5 * kpc / r) * (1 - lfDisk);
     226           0 :                                 b.x += -bMag * sinPhi;
     227           0 :                                 b.y += bMag * cosPhi;
     228             :                         } else {
     229             :                                 // spiral region
     230           0 :                                 double r_negx = r * exp(-(phi - M_PI) / tan90MinusPitch);
     231           0 :                                 if (r_negx > rArms[7])
     232           0 :                                         r_negx = r * exp(-(phi + M_PI) / tan90MinusPitch);
     233           0 :                                 if (r_negx > rArms[7])
     234           0 :                                         r_negx = r * exp(-(phi + 3 * M_PI) / tan90MinusPitch);
     235             : 
     236           0 :                                 for (int i = 7; i >= 0; i--)
     237           0 :                                         if (r_negx < rArms[i])
     238           0 :                                                 bMag = bDisk[i];
     239             : 
     240           0 :                                 bMag *= (5 * kpc / r) * (1 - lfDisk);
     241           0 :                                 b.x += bMag * (sinPitch * cosPhi - cosPitch * sinPhi);
     242           0 :                                 b.y += bMag * (sinPitch * sinPhi + cosPitch * cosPhi);
     243             :                         }
     244             :                 }
     245             :         }
     246           0 :         return b;
     247             : }
     248             : 
     249           0 : Vector3d JF12Field::getToroidalHaloField(const double& r, const double& z, const double& sinPhi, const double& cosPhi) const {
     250             :         Vector3d b(0.);
     251             : 
     252           0 :         if (useToroidalHaloField && (r * r + z * z > 1 * kpc * kpc)){
     253             : 
     254           0 :                 double lfDisk = logisticFunction(z, hDisk, wDisk);
     255           0 :                 double bMagH = exp(-fabs(z) / z0) * lfDisk;
     256             : 
     257           0 :                 if (z >= 0)
     258           0 :                         bMagH *= bNorth * (1 - logisticFunction(r, rNorth, wHalo));
     259             :                 else
     260           0 :                         bMagH *= bSouth * (1 - logisticFunction(r, rSouth, wHalo));
     261           0 :                 b.x += -bMagH * sinPhi;
     262           0 :                 b.y += bMagH * cosPhi;
     263             :         }
     264           0 :         return b;
     265             : }
     266             : 
     267           0 : Vector3d JF12Field::getXField(const double& r, const double& z, const double& sinPhi, const double& cosPhi) const {
     268             :         Vector3d b(0.);
     269             : 
     270           0 :         if (useXField && (r * r + z * z > 1 * kpc * kpc)){
     271             :                 double bMagX;
     272             :                 double sinThetaX, cosThetaX;
     273             :                 double rp;
     274           0 :                 double rc = rXc + fabs(z) / tanThetaX0;
     275           0 :                 if (r < rc) {
     276             :                         // varying elevation region
     277           0 :                         rp = r * rXc / rc;
     278           0 :                         bMagX = bX * exp(-1 * rp / rX) * pow(rXc / rc, 2.);
     279           0 :                         double thetaX = atan2(fabs(z), (r - rp));
     280           0 :                         if (z == 0)
     281             :                                 thetaX = M_PI / 2.;
     282           0 :                         sinThetaX = sin(thetaX);
     283           0 :                         cosThetaX = cos(thetaX);
     284             :                 } else {
     285             :                         // constant elevation region
     286           0 :                         rp = r - fabs(z) / tanThetaX0;
     287           0 :                         bMagX = bX * exp(-rp / rX) * (rp / r);
     288           0 :                         sinThetaX = sinThetaX0;
     289           0 :                         cosThetaX = cosThetaX0;
     290             :                 }
     291           0 :                 double zsign = z < 0 ? -1 : 1;
     292           0 :                 b.x += zsign * bMagX * cosThetaX * cosPhi;
     293           0 :                 b.y += zsign * bMagX * cosThetaX * sinPhi;
     294           0 :                 b.z += bMagX * sinThetaX;
     295             :         }
     296           0 :         return b;
     297             : }
     298             : 
     299           0 : Vector3d JF12Field::getStriatedField(const Vector3d& pos) const {
     300           0 :         return (getRegularField(pos)
     301           0 :                         * (1. + sqrtbeta * striatedGrid->closestValue(pos)));
     302             : }
     303             : 
     304           0 : double JF12Field::getTurbulentStrength(const Vector3d& pos) const {
     305           0 :         if (pos.getR() > 20 * kpc)
     306             :                 return 0;
     307             : 
     308           0 :         double r = sqrt(pos.x * pos.x + pos.y * pos.y); // in-plane radius
     309           0 :         double phi = pos.getPhi(); // azimuth
     310             : 
     311             :         // disk
     312             :         double bDisk = 0;
     313           0 :         if (r < 5 * kpc) {
     314           0 :                 bDisk = bDiskTurb5;
     315             :         } else {
     316             :                 // spiral region
     317           0 :                 double r_negx = r * exp(-(phi - M_PI) / tan90MinusPitch);
     318           0 :                 if (r_negx > rArms[7])
     319           0 :                         r_negx = r * exp(-(phi + M_PI) / tan90MinusPitch);
     320           0 :                 if (r_negx > rArms[7])
     321           0 :                         r_negx = r * exp(-(phi + 3 * M_PI) / tan90MinusPitch);
     322             : 
     323           0 :                 for (int i = 7; i >= 0; i--)
     324           0 :                         if (r_negx < rArms[i])
     325           0 :                                 bDisk = bDiskTurb[i];
     326             : 
     327           0 :                 bDisk *= (5 * kpc) / r;
     328             :         }
     329           0 :         bDisk *= exp(-0.5 * pow(pos.z / zDiskTurb, 2));
     330             : 
     331             :         // halo
     332           0 :         double bHalo = bHaloTurb * exp(-r / rHaloTurb)
     333           0 :                         * exp(-0.5 * pow(pos.z / zHaloTurb, 2));
     334             : 
     335             :         // modulate turbulent field
     336           0 :         return sqrt(pow(bDisk, 2) + pow(bHalo, 2));
     337             : }
     338             : 
     339           0 : Vector3d JF12Field::getTurbulentField(const Vector3d& pos) const {
     340           0 :         return (turbulentGrid->interpolate(pos) * getTurbulentStrength(pos));
     341             : }
     342             : 
     343           0 : Vector3d JF12Field::getField(const Vector3d& pos) const {
     344             :         Vector3d b(0.);
     345           0 :         if (useTurbulentField)
     346           0 :                 b += getTurbulentField(pos);
     347           0 :         if (useStriatedField)
     348           0 :                 b += getStriatedField(pos);
     349           0 :         else if (useRegularField)
     350           0 :                 b += getRegularField(pos);
     351           0 :         return b;
     352             : }
     353             : 
     354             : 
     355             : 
     356           0 : PlanckJF12bField::PlanckJF12bField() : JF12Field::JF12Field(){
     357             :         // regular field parameters
     358           0 :         bDisk[5] = -3.5 * muG;
     359           0 :         bX = 1.8 * muG;
     360             : 
     361             :         // turbulent field parameters;
     362           0 :         bDiskTurb[0] = 3.12 * muG;
     363           0 :         bDiskTurb[1] = 6.24 * muG;
     364           0 :         bDiskTurb[2] = 3.12 * muG;
     365           0 :         bDiskTurb[3] = 6.24 * muG;
     366           0 :         bDiskTurb[4] = 3.12 * muG;
     367           0 :         bDiskTurb[5] = 6.24 * muG;
     368           0 :         bDiskTurb[6] = 3.12 * muG;
     369           0 :         bDiskTurb[7] = 6.24 * muG;
     370             : 
     371           0 :         bDiskTurb5 = 3.90 * muG;
     372             : 
     373           0 :         bHaloTurb = 7.332 * muG;
     374           0 : }
     375             : 
     376             : } // namespace crpropa

Generated by: LCOV version 1.14