Line data Source code
1 : #include "crpropa/magneticField/UF23Field.h"
2 :
3 : #include <exception>
4 : #include <limits>
5 : #include <string>
6 : #include <cmath>
7 :
8 : // local helper functions and constants
9 : namespace uf23 {
10 :
11 : template<typename T>
12 : crpropa::Vector3d CylToCart(const T v, const double cosPhi, const double sinPhi)
13 : {
14 195 : return crpropa::Vector3d(v[0] * cosPhi - v[1] * sinPhi,
15 195 : v[0] * sinPhi + v[1] * cosPhi,
16 56 : v[2]);
17 : }
18 :
19 : template<typename T>
20 : crpropa::Vector3d CartToCyl(const T v, const double cosPhi, const double sinPhi)
21 : {
22 9 : return crpropa::Vector3d(v[0] * cosPhi + v[1] * sinPhi,
23 : -v[0] * sinPhi + v[1] * cosPhi,
24 : v[2]);
25 : }
26 :
27 : // logistic sigmoid function
28 : inline
29 : double
30 : Sigmoid(const double x, const double x0, const double w)
31 : {
32 185 : return 1 / (1 + exp(-(x-x0)/w));
33 : }
34 :
35 : // angle between v0 = (cos(phi0), sin(phi0)) and v1 = (cos(phi1), sin(phi1))
36 : inline
37 : double
38 6 : DeltaPhi(const double phi0, const double phi1)
39 : {
40 6 : return acos(cos(phi1)*cos(phi0) + sin(phi1)*sin(phi0));
41 : }
42 :
43 : // Interal units used in this code.
44 : // Convert to crpropa with e.g. uf23::kpc / crpropa::kpc.
45 : // The conversion is, however, only needed in the single non-private
46 : // getField() method, all other functions use uf23 units.
47 : const double kPi = 3.1415926535897932384626;
48 : const double kTwoPi = 2*kPi;
49 : const double degree = kPi/180.;
50 : const double kpc = 1;
51 : const double microgauss = 1;
52 : const double megayear = 1;
53 : const double Gpc = 1e6*kpc;
54 : const double pc = 1e-3*kpc;
55 : const double second = megayear / (1e6*60*60*24*365.25);
56 : const double kilometer = kpc / 3.0856775807e+16;
57 : }
58 :
59 : namespace crpropa {
60 8 : UF23Field::UF23Field(const ModelType mt) :
61 8 : fModelType(mt),
62 8 : fMaxRadiusSquared(pow(30*uf23::kpc, 2))
63 : {
64 :
65 : // all but expX model have a-->\infty, Eq.(38)
66 8 : fPoloidalA = 1 * Gpc;
67 :
68 8 : switch (fModelType) {
69 : // ---------------------------------------------
70 1 : case base: {
71 1 : fDiskB1 = 1.0878565e+00 * uf23::microgauss;
72 1 : fDiskB2 = 2.6605034e+00 * uf23::microgauss;
73 1 : fDiskB3 = 3.1166311e+00 * uf23::microgauss;
74 1 : fDiskH = 7.9408965e-01 * uf23::kpc;
75 1 : fDiskPhase1 = 2.6316589e+02 * uf23::degree;
76 1 : fDiskPhase2 = 9.7782269e+01 * uf23::degree;
77 1 : fDiskPhase3 = 3.5112281e+01 * uf23::degree;
78 1 : fDiskPitch = 1.0106900e+01 * uf23::degree;
79 1 : fDiskW = 1.0720909e-01 * uf23::kpc;
80 1 : fPoloidalB = 9.7775487e-01 * uf23::microgauss;
81 1 : fPoloidalP = 1.4266186e+00 * uf23::kpc;
82 1 : fPoloidalR = 7.2925417e+00 * uf23::kpc;
83 1 : fPoloidalW = 1.1188158e-01 * uf23::kpc;
84 1 : fPoloidalZ = 4.4597373e+00 * uf23::kpc;
85 1 : fStriation = 3.4557571e-01;
86 1 : fToroidalBN = 3.2556760e+00 * uf23::microgauss;
87 1 : fToroidalBS = -3.0914569e+00 * uf23::microgauss;
88 1 : fToroidalR = 1.0193815e+01 * uf23::kpc;
89 1 : fToroidalW = 1.6936993e+00 * uf23::kpc;
90 1 : fToroidalZ = 4.0242749e+00 * uf23::kpc;
91 1 : break;
92 : }
93 1 : case cre10: {
94 : // ---------------------------------------------
95 1 : fDiskB1 = 1.2035697e+00 * uf23::microgauss;
96 1 : fDiskB2 = 2.7478490e+00 * uf23::microgauss;
97 1 : fDiskB3 = 3.2104342e+00 * uf23::microgauss;
98 1 : fDiskH = 8.0844932e-01 * uf23::kpc;
99 1 : fDiskPhase1 = 2.6515882e+02 * uf23::degree;
100 1 : fDiskPhase2 = 9.8211313e+01 * uf23::degree;
101 1 : fDiskPhase3 = 3.5944588e+01 * uf23::degree;
102 1 : fDiskPitch = 1.0162759e+01 * uf23::degree;
103 1 : fDiskW = 1.0824003e-01 * uf23::kpc;
104 1 : fPoloidalB = 9.6938453e-01 * uf23::microgauss;
105 1 : fPoloidalP = 1.4150957e+00 * uf23::kpc;
106 1 : fPoloidalR = 7.2987296e+00 * uf23::kpc;
107 1 : fPoloidalW = 1.0923051e-01 * uf23::kpc;
108 1 : fPoloidalZ = 4.5748332e+00 * uf23::kpc;
109 1 : fStriation = 2.4950386e-01;
110 1 : fToroidalBN = 3.7308133e+00 * uf23::microgauss;
111 1 : fToroidalBS = -3.5039958e+00 * uf23::microgauss;
112 1 : fToroidalR = 1.0407507e+01 * uf23::kpc;
113 1 : fToroidalW = 1.7398375e+00 * uf23::kpc;
114 1 : fToroidalZ = 2.9272800e+00 * uf23::kpc;
115 1 : break;
116 : }
117 1 : case nebCor: {
118 : // ---------------------------------------------
119 1 : fDiskB1 = 1.4081935e+00 * uf23::microgauss;
120 1 : fDiskB2 = 3.5292400e+00 * uf23::microgauss;
121 1 : fDiskB3 = 4.1290147e+00 * uf23::microgauss;
122 1 : fDiskH = 8.1151971e-01 * uf23::kpc;
123 1 : fDiskPhase1 = 2.6447529e+02 * uf23::degree;
124 1 : fDiskPhase2 = 9.7572660e+01 * uf23::degree;
125 1 : fDiskPhase3 = 3.6403798e+01 * uf23::degree;
126 1 : fDiskPitch = 1.0151183e+01 * uf23::degree;
127 1 : fDiskW = 1.1863734e-01 * uf23::kpc;
128 1 : fPoloidalB = 1.3485916e+00 * uf23::microgauss;
129 1 : fPoloidalP = 1.3414395e+00 * uf23::kpc;
130 1 : fPoloidalR = 7.2473841e+00 * uf23::kpc;
131 1 : fPoloidalW = 1.4318227e-01 * uf23::kpc;
132 1 : fPoloidalZ = 4.8242603e+00 * uf23::kpc;
133 1 : fStriation = 3.8610837e-10;
134 1 : fToroidalBN = 4.6491142e+00 * uf23::microgauss;
135 1 : fToroidalBS = -4.5006610e+00 * uf23::microgauss;
136 1 : fToroidalR = 1.0205288e+01 * uf23::kpc;
137 1 : fToroidalW = 1.7004868e+00 * uf23::kpc;
138 1 : fToroidalZ = 3.5557767e+00 * uf23::kpc;
139 1 : break;
140 : }
141 1 : case neCL: {
142 : // ---------------------------------------------
143 1 : fDiskB1 = 1.4259645e+00 * uf23::microgauss;
144 1 : fDiskB2 = 1.3543223e+00 * uf23::microgauss;
145 1 : fDiskB3 = 3.4390669e+00 * uf23::microgauss;
146 1 : fDiskH = 6.7405199e-01 * uf23::kpc;
147 1 : fDiskPhase1 = 1.9961898e+02 * uf23::degree;
148 1 : fDiskPhase2 = 1.3541461e+02 * uf23::degree;
149 1 : fDiskPhase3 = 6.4909767e+01 * uf23::degree;
150 1 : fDiskPitch = 1.1867859e+01 * uf23::degree;
151 1 : fDiskW = 6.1162799e-02 * uf23::kpc;
152 1 : fPoloidalB = 9.8387831e-01 * uf23::microgauss;
153 1 : fPoloidalP = 1.6773615e+00 * uf23::kpc;
154 1 : fPoloidalR = 7.4084361e+00 * uf23::kpc;
155 1 : fPoloidalW = 1.4168192e-01 * uf23::kpc;
156 1 : fPoloidalZ = 3.6521188e+00 * uf23::kpc;
157 1 : fStriation = 3.3600213e-01;
158 1 : fToroidalBN = 2.6256593e+00 * uf23::microgauss;
159 1 : fToroidalBS = -2.5699466e+00 * uf23::microgauss;
160 1 : fToroidalR = 1.0134257e+01 * uf23::kpc;
161 1 : fToroidalW = 1.1547728e+00 * uf23::kpc;
162 1 : fToroidalZ = 4.5585463e+00 * uf23::kpc;
163 1 : break;
164 : }
165 1 : case spur: {
166 : // ---------------------------------------------
167 1 : fDiskB1 = -4.2993328e+00 * uf23::microgauss;
168 1 : fDiskH = 7.5019749e-01 * uf23::kpc;
169 1 : fDiskPhase1 = 1.5589875e+02 * uf23::degree;
170 1 : fDiskPitch = 1.2074432e+01 * uf23::degree;
171 1 : fDiskW = 1.2263120e-01 * uf23::kpc;
172 1 : fPoloidalB = 9.9302987e-01 * uf23::microgauss;
173 1 : fPoloidalP = 1.3982374e+00 * uf23::kpc;
174 1 : fPoloidalR = 7.1973387e+00 * uf23::kpc;
175 1 : fPoloidalW = 1.2262244e-01 * uf23::kpc;
176 1 : fPoloidalZ = 4.4853270e+00 * uf23::kpc;
177 1 : fSpurCenter = 1.5718686e+02 * uf23::degree;
178 1 : fSpurLength = 3.1839577e+01 * uf23::degree;
179 1 : fSpurWidth = 1.0318114e+01 * uf23::degree;
180 1 : fStriation = 3.3022369e-01;
181 1 : fToroidalBN = 2.9286724e+00 * uf23::microgauss;
182 1 : fToroidalBS = -2.5979895e+00 * uf23::microgauss;
183 1 : fToroidalR = 9.7536425e+00 * uf23::kpc;
184 1 : fToroidalW = 1.4210055e+00 * uf23::kpc;
185 1 : fToroidalZ = 6.0941229e+00 * uf23::kpc;
186 1 : break;
187 : }
188 1 : case synCG: {
189 : // ---------------------------------------------
190 1 : fDiskB1 = 8.1386878e-01 * uf23::microgauss;
191 1 : fDiskB2 = 2.0586930e+00 * uf23::microgauss;
192 1 : fDiskB3 = 2.9437335e+00 * uf23::microgauss;
193 1 : fDiskH = 6.2172353e-01 * uf23::kpc;
194 1 : fDiskPhase1 = 2.2988551e+02 * uf23::degree;
195 1 : fDiskPhase2 = 9.7388282e+01 * uf23::degree;
196 1 : fDiskPhase3 = 3.2927367e+01 * uf23::degree;
197 1 : fDiskPitch = 9.9034844e+00 * uf23::degree;
198 1 : fDiskW = 6.6517521e-02 * uf23::kpc;
199 1 : fPoloidalB = 8.0883734e-01 * uf23::microgauss;
200 1 : fPoloidalP = 1.5820957e+00 * uf23::kpc;
201 1 : fPoloidalR = 7.4625235e+00 * uf23::kpc;
202 1 : fPoloidalW = 1.5003765e-01 * uf23::kpc;
203 1 : fPoloidalZ = 3.5338550e+00 * uf23::kpc;
204 1 : fStriation = 6.3434763e-01;
205 1 : fToroidalBN = 2.3991193e+00 * uf23::microgauss;
206 1 : fToroidalBS = -2.0919944e+00 * uf23::microgauss;
207 1 : fToroidalR = 9.4227834e+00 * uf23::kpc;
208 1 : fToroidalW = 9.1608418e-01 * uf23::kpc;
209 1 : fToroidalZ = 5.5844594e+00 * uf23::kpc;
210 1 : break;
211 : }
212 1 : case twistX: {
213 : // ---------------------------------------------
214 1 : fDiskB1 = 1.3741995e+00 * uf23::microgauss;
215 1 : fDiskB2 = 2.0089881e+00 * uf23::microgauss;
216 1 : fDiskB3 = 1.5212463e+00 * uf23::microgauss;
217 1 : fDiskH = 9.3806180e-01 * uf23::kpc;
218 1 : fDiskPhase1 = 2.3560316e+02 * uf23::degree;
219 1 : fDiskPhase2 = 1.0189856e+02 * uf23::degree;
220 1 : fDiskPhase3 = 5.6187572e+01 * uf23::degree;
221 1 : fDiskPitch = 1.2100979e+01 * uf23::degree;
222 1 : fDiskW = 1.4933338e-01 * uf23::kpc;
223 1 : fPoloidalB = 6.2793114e-01 * uf23::microgauss;
224 1 : fPoloidalP = 2.3292519e+00 * uf23::kpc;
225 1 : fPoloidalR = 7.9212358e+00 * uf23::kpc;
226 1 : fPoloidalW = 2.9056201e-01 * uf23::kpc;
227 1 : fPoloidalZ = 2.6274437e+00 * uf23::kpc;
228 1 : fStriation = 7.7616317e-01;
229 1 : fTwistingTime = 5.4733549e+01 * uf23::megayear;
230 1 : break;
231 : }
232 1 : case expX: {
233 : // ---------------------------------------------
234 1 : fDiskB1 = 9.9258148e-01 * uf23::microgauss;
235 1 : fDiskB2 = 2.1821124e+00 * uf23::microgauss;
236 1 : fDiskB3 = 3.1197345e+00 * uf23::microgauss;
237 1 : fDiskH = 7.1508681e-01 * uf23::kpc;
238 1 : fDiskPhase1 = 2.4745741e+02 * uf23::degree;
239 1 : fDiskPhase2 = 9.8578879e+01 * uf23::degree;
240 1 : fDiskPhase3 = 3.4884485e+01 * uf23::degree;
241 1 : fDiskPitch = 1.0027070e+01 * uf23::degree;
242 1 : fDiskW = 9.8524736e-02 * uf23::kpc;
243 1 : fPoloidalA = 6.1938701e+00 * uf23::kpc;
244 1 : fPoloidalB = 5.8357990e+00 * uf23::microgauss;
245 1 : fPoloidalP = 1.9510779e+00 * uf23::kpc;
246 1 : fPoloidalR = 2.4994376e+00 * uf23::kpc;
247 : // internally, xi is fitted and z = tan(xi)*a
248 1 : fPoloidalXi = 2.0926122e+01 * uf23::degree;
249 1 : fPoloidalZ = fPoloidalA*tan(fPoloidalXi);
250 1 : fStriation = 5.1440500e-01;
251 1 : fToroidalBN = 2.7077434e+00 * uf23::microgauss;
252 1 : fToroidalBS = -2.5677104e+00 * uf23::microgauss;
253 1 : fToroidalR = 1.0134022e+01 * uf23::kpc;
254 1 : fToroidalW = 2.0956159e+00 * uf23::kpc;
255 1 : fToroidalZ = 5.4564991e+00 * uf23::kpc;
256 1 : break;
257 : }
258 0 : default: {
259 0 : throw std::runtime_error("unknown field model");
260 : break;
261 : }
262 : }
263 :
264 8 : fSinPitch = sin(fDiskPitch);
265 8 : fCosPitch = cos(fDiskPitch);
266 8 : fTanPitch = tan(fDiskPitch);
267 :
268 8 : }
269 :
270 : Vector3d
271 72 : UF23Field::getField(const Vector3d &position)
272 : const
273 : {
274 : Vector3d posInKpc = position / kpc;
275 72 : return (this->operator()(posInKpc)) / uf23::microgauss * microgauss;
276 : }
277 :
278 :
279 : Vector3d
280 72 : UF23Field::operator()(const Vector3d& posInKpc)
281 : const
282 : {
283 : const auto pos = posInKpc * uf23::kpc;
284 72 : if (pos.getR2() > fMaxRadiusSquared)
285 : return Vector3d(0, 0, 0);
286 : else {
287 72 : const auto diskField = getDiskField(pos);
288 72 : const auto haloField = getHaloField(pos);
289 : return (diskField + haloField) / uf23::microgauss;
290 : }
291 : }
292 :
293 : Vector3d
294 72 : UF23Field::getDiskField(const Vector3d& pos)
295 : const
296 : {
297 72 : if (fModelType == spur)
298 9 : return getSpurField(pos.x, pos.y, pos.z);
299 : else
300 63 : return getSpiralField(pos.x, pos.y, pos.z);
301 : }
302 :
303 :
304 : Vector3d
305 72 : UF23Field::getHaloField(const Vector3d& pos)
306 : const
307 : {
308 72 : if (fModelType == twistX)
309 9 : return getTwistedHaloField(pos.x, pos.y, pos.z);
310 : else
311 : return
312 63 : getToroidalHaloField(pos.x, pos.y, pos.z) +
313 126 : getPoloidalHaloField(pos.x, pos.y, pos.z);
314 : }
315 :
316 :
317 : Vector3d
318 9 : UF23Field::getTwistedHaloField(const double x, const double y, const double z)
319 : const
320 : {
321 9 : const double r = sqrt(x*x + y*y);
322 9 : const double cosPhi = r > std::numeric_limits<double>::min() ? x / r : 1;
323 8 : const double sinPhi = r > std::numeric_limits<double>::min() ? y / r : 0;
324 :
325 9 : const Vector3d bXCart = getPoloidalHaloField(x, y, z);
326 9 : const double bXCartTmp[3] = {bXCart.x, bXCart.y, bXCart.z};
327 : const Vector3d bXCyl = uf23::CartToCyl(bXCartTmp, cosPhi, sinPhi);
328 :
329 : const double bZ = bXCyl.z;
330 : const double bR = bXCyl.x;
331 :
332 : double bPhi = 0;
333 :
334 9 : if (fTwistingTime != 0 && r != 0) {
335 : // radial rotation curve parameters (fit to Reid et al 2014)
336 : const double v0 = -240 * uf23::kilometer/uf23::second;
337 : const double r0 = 1.6 * uf23::kpc;
338 : // vertical gradient (Levine+08)
339 : const double z0 = 10 * uf23::kpc;
340 :
341 : // Eq.(43)
342 8 : const double fr = 1 - exp(-r/r0);
343 : // Eq.(44)
344 8 : const double t0 = exp(2*std::abs(z)/z0);
345 8 : const double gz = 2 / (1 + t0);
346 :
347 : // Eq. (46)
348 8 : const double signZ = z < 0 ? -1 : 1;
349 8 : const double deltaZ = -signZ * v0 * fr / z0 * t0 * pow(gz, 2);
350 : // Eq. (47)
351 8 : const double deltaR = v0 * ((1-fr)/r0 - fr/r) * gz;
352 :
353 : // Eq.(45)
354 8 : bPhi = (bZ * deltaZ + bR * deltaR) * fTwistingTime;
355 :
356 : }
357 : const double bCylX[3] = {bR, bPhi , bZ};
358 9 : return uf23::CylToCart(bCylX, cosPhi, sinPhi);
359 : }
360 :
361 : Vector3d
362 63 : UF23Field::getToroidalHaloField(const double x, const double y, const double z)
363 : const
364 : {
365 63 : const double r2 = x*x + y*y;
366 63 : const double r = sqrt(r2);
367 : const double absZ = std::abs(z);
368 :
369 63 : const double b0 = z >= 0 ? fToroidalBN : fToroidalBS;
370 63 : const double rh = fToroidalR;
371 63 : const double z0 = fToroidalZ;
372 63 : const double fwh = fToroidalW;
373 : const double sigmoidR = uf23::Sigmoid(r, rh, fwh);
374 63 : const double sigmoidZ = uf23::Sigmoid(absZ, fDiskH, fDiskW);
375 :
376 : // Eq. (21)
377 63 : const double bPhi = b0 * (1. - sigmoidR) * sigmoidZ * exp(-absZ/z0);
378 :
379 : const double bCyl[3] = {0, bPhi, 0};
380 63 : const double cosPhi = r > std::numeric_limits<double>::min() ? x / r : 1;
381 56 : const double sinPhi = r > std::numeric_limits<double>::min() ? y / r : 0;
382 63 : return uf23::CylToCart(bCyl, cosPhi, sinPhi);
383 : }
384 :
385 : Vector3d
386 72 : UF23Field::getPoloidalHaloField(const double x, const double y, const double z)
387 : const
388 : {
389 72 : const double r2 = x*x + y*y;
390 72 : const double r = sqrt(r2);
391 :
392 72 : const double c = pow(fPoloidalA/fPoloidalZ, fPoloidalP);
393 72 : const double a0p = pow(fPoloidalA, fPoloidalP);
394 72 : const double rp = pow(r, fPoloidalP);
395 72 : const double abszp = pow(std::abs(z), fPoloidalP);
396 72 : const double cabszp = c*abszp;
397 :
398 : /*
399 : since $\sqrt{a^2 + b} - a$ is numerical unstable for $b\ll a$,
400 : we use $(\sqrt{a^2 + b} - a) \frac{\sqrt{a^2 + b} + a}{\sqrt{a^2
401 : + b} + a} = \frac{b}{\sqrt{a^2 + b} + a}$}
402 : */
403 :
404 72 : const double t0 = a0p + cabszp - rp;
405 72 : const double t1 = sqrt(pow(t0, 2) + 4*a0p*rp);
406 72 : const double ap = 2*a0p*rp / (t1 + t0);
407 :
408 : double a = 0;
409 72 : if (ap < 0) {
410 0 : if (r > std::numeric_limits<double>::min()) {
411 : // this should never happen
412 0 : throw std::runtime_error("ap = " + std::to_string(ap));
413 : }
414 : else
415 : a = 0;
416 : }
417 : else
418 72 : a = pow(ap, 1/fPoloidalP);
419 :
420 : // Eq.(29) and Eq.(32)
421 : const double radialDependence =
422 72 : fModelType == expX ?
423 9 : exp(-a/fPoloidalR) :
424 63 : 1 - uf23::Sigmoid(a, fPoloidalR, fPoloidalW);
425 :
426 : // Eq.(28)
427 72 : const double Bzz = fPoloidalB * radialDependence;
428 :
429 : // (r/a)
430 72 : const double rOverA = 1 / pow(2*a0p / (t1 + t0), 1/fPoloidalP);
431 :
432 : // Eq.(35) for p=n
433 72 : const double signZ = z < 0 ? -1 : 1;
434 : const double Br =
435 72 : Bzz * c * a / rOverA * signZ * pow(std::abs(z), fPoloidalP - 1) / t1;
436 :
437 : // Eq.(36) for p=n
438 72 : const double Bz = Bzz * pow(rOverA, fPoloidalP-2) * (ap + a0p) / t1;
439 :
440 72 : if (r < std::numeric_limits<double>::min())
441 : return Vector3d(0, 0, Bz);
442 : else {
443 : const double bCylX[3] = {Br, 0 , Bz};
444 64 : const double cosPhi = x / r;
445 64 : const double sinPhi = y / r;
446 : return uf23::CylToCart(bCylX, cosPhi, sinPhi);
447 : }
448 : }
449 :
450 : Vector3d
451 9 : UF23Field::getSpurField(const double x, const double y, const double z)
452 : const
453 : {
454 : // reference approximately at solar radius
455 : const double rRef = 8.2*uf23::kpc;
456 :
457 : // cylindrical coordinates
458 9 : const double r2 = x*x + y*y;
459 9 : const double r = sqrt(r2);
460 9 : if (r < std::numeric_limits<double>::min())
461 : return Vector3d(0, 0, 0);
462 :
463 8 : double phi = atan2(y, x);
464 8 : if (phi < 0)
465 4 : phi += uf23::kTwoPi;
466 :
467 8 : const double phiRef = fDiskPhase1;
468 : int iBest = -2;
469 : double bestDist = -1;
470 32 : for (int i = -1; i <= 1; ++i) {
471 24 : const double pphi = phi - phiRef + i*uf23::kTwoPi;
472 24 : const double rr = rRef*exp(pphi * fTanPitch);
473 24 : if (bestDist < 0 || std::abs(r-rr) < bestDist) {
474 11 : bestDist = std::abs(r-rr);
475 : iBest = i;
476 : }
477 : }
478 8 : if (iBest == 0) {
479 3 : const double phi0 = phi - log(r/rRef) / fTanPitch;
480 :
481 : // Eq. (16)
482 3 : const double deltaPhi0 = uf23::DeltaPhi(phiRef, phi0);
483 3 : const double delta = deltaPhi0 / fSpurWidth;
484 3 : const double B = fDiskB1 * exp(-0.5*pow(delta, 2));
485 :
486 : // Eq. (18)
487 : const double wS = 5*uf23::degree;
488 3 : const double phiC = fSpurCenter;
489 3 : const double deltaPhiC = uf23::DeltaPhi(phiC, phi);
490 3 : const double lC = fSpurLength;
491 3 : const double gS = 1 - uf23::Sigmoid(std::abs(deltaPhiC), lC, wS);
492 :
493 : // Eq. (13)
494 3 : const double hd = 1 - uf23::Sigmoid(std::abs(z), fDiskH, fDiskW);
495 :
496 : // Eq. (17)
497 3 : const double bS = rRef/r * B * hd * gS;
498 3 : const double bCyl[3] = {bS * fSinPitch, bS * fCosPitch, 0};
499 3 : const double cosPhi = x / r;
500 3 : const double sinPhi = y / r;
501 : return uf23::CylToCart(bCyl, cosPhi, sinPhi);
502 : }
503 : else
504 : return Vector3d(0, 0, 0);
505 :
506 : }
507 :
508 : Vector3d
509 63 : UF23Field::getSpiralField(const double x, const double y, const double z)
510 : const
511 : {
512 : // reference radius
513 : const double rRef = 5*uf23::kpc;
514 : // inner boundary of spiral field
515 : const double rInner = 5*uf23::kpc;
516 : const double wInner = 0.5*uf23::kpc;
517 : // outer boundary of spiral field
518 : const double rOuter = 20*uf23::kpc;
519 : const double wOuter = 0.5*uf23::kpc;
520 :
521 : // cylindrical coordinates
522 63 : const double r2 = x*x + y*y;
523 63 : if (r2 == 0)
524 : return Vector3d(0, 0, 0);
525 56 : const double r = sqrt(r2);
526 56 : const double phi = atan2(y, x);
527 :
528 : // Eq.(13)
529 56 : const double hdz = 1 - uf23::Sigmoid(std::abs(z), fDiskH, fDiskW);
530 :
531 : // Eq.(14) times rRef divided by r
532 : const double rFacI = uf23::Sigmoid(r, rInner, wInner);
533 56 : const double rFacO = 1 - uf23::Sigmoid(r, rOuter, wOuter);
534 : // (using lim r--> 0 (1-exp(-r^2))/r --> r - r^3/2 + ...)
535 56 : const double rFac = r > 1e-5*uf23::pc ? (1-exp(-r*r)) / r : r * (1 - r2/2);
536 56 : const double gdrTimesRrefByR = rRef * rFac * rFacO * rFacI;
537 :
538 : // Eq. (12)
539 56 : const double phi0 = phi - log(r/rRef) / fTanPitch;
540 :
541 : // Eq. (10)
542 56 : const double b =
543 56 : fDiskB1 * cos(1 * (phi0 - fDiskPhase1)) +
544 56 : fDiskB2 * cos(2 * (phi0 - fDiskPhase2)) +
545 56 : fDiskB3 * cos(3 * (phi0 - fDiskPhase3));
546 :
547 : // Eq. (11)
548 56 : const double fac = hdz * gdrTimesRrefByR;
549 : const double bCyl[3] =
550 56 : { b * fac * fSinPitch,
551 56 : b * fac * fCosPitch,
552 : 0};
553 :
554 56 : const double cosPhi = x / r;
555 56 : const double sinPhi = y / r;
556 : return uf23::CylToCart(bCyl, cosPhi, sinPhi);
557 : }
558 : }
|
