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
|
