Line data Source code
1 : #include "crpropa/magneticField/KST24Field.h"
2 :
3 : /*
4 : The C++ implementation of the GMF model KST24 (A.Korochkin, D.Semikoz, P.Tinyakov 2024)
5 : The model was presented in arXiv:2407.02148 and published in A&A
6 : If you use the model, please cite A&A, 693, A284 (2025)
7 :
8 : In KST24 GMF model the position of the Solar System is at {-8.2 kpc, 0, 0}
9 : The Galactic Center is at {0, 0, 0}
10 : z-axis points to the North pole
11 : */
12 :
13 :
14 : namespace crpropa {
15 0 : KST24Field::KST24Field()
16 : {
17 0 : north_tor_B_gauss = -3.2e-6;
18 0 : north_tor_zmin_kpc = 1.185;
19 0 : north_tor_zmax_kpc = 2.1;
20 0 : north_tor_rmin_kpc = 1.0;
21 0 : north_tor_rmax_kpc = 9.1;
22 :
23 0 : south_tor_B_gauss = 3.2e-6;
24 0 : south_tor_zmin_kpc = -2.5;
25 0 : south_tor_zmax_kpc = -1.22;
26 0 : south_tor_rmin_kpc = 1.0;
27 0 : south_tor_rmax_kpc = 14.0;
28 :
29 0 : Xfield_B_gauss = 1.8e-6;
30 0 : Xfield_rmin_kpc = 1.0;
31 0 : Xfield_rmax_kpc = 6.2;
32 0 : Xfield_theta_deg = 28.0;
33 0 : Xfield_theta_rad = Xfield_theta_deg*M_PI/180.;
34 :
35 0 : LB_B_gauss = -3.5e-6;
36 0 : LB_lB_deg = 50.0;
37 0 : LB_bB_deg = 2.0;
38 0 : LB_rmin_kpc = 0.2;
39 0 : LB_dr_kpc = 0.03;
40 0 : LB_x0_kpc = -8.2;
41 0 : LB_y0_kpc = 0.095;
42 0 : LB_z0_kpc = -0.05;
43 0 : LB_Bdir_x = cos(LB_bB_deg*M_PI/180)*cos(LB_lB_deg*M_PI/180);
44 0 : LB_Bdir_y = cos(LB_bB_deg*M_PI/180)*sin(LB_lB_deg*M_PI/180);
45 0 : LB_Bdir_z = sin(LB_bB_deg*M_PI/180);
46 :
47 :
48 0 : ScutumArm_B_gauss = 4.9e-6;
49 0 : ScutumArm_pitch_deg = 20;
50 0 : ScutumArm_phi0_deg = -134.0;
51 0 : ScutumArm_x_shift_kpc = 1.0;
52 0 : ScutumArm_y_shift_kpc = 0.0;
53 0 : ScutumArm_arc_radius1_kpc = 0.8;
54 0 : ScutumArm_arc_radius2_kpc = 1.0;
55 0 : ScutumArm_arc_eps = 2.0;
56 0 : ScutumArm_arc_div_deg = 0.0;
57 0 : ScutumArm_rmin_kpc = 3.0;
58 0 : ScutumArm_rmax_kpc = 16.0;
59 0 : ScutumArm_zmin_kpc = -1.0;
60 0 : ScutumArm_zmax_kpc = 1.0;
61 :
62 0 : CarinaSagittariusArm_B_gauss = 1.3e-6;
63 0 : CarinaSagittariusArm_pitch_deg = 20;
64 0 : CarinaSagittariusArm_phi0_deg = -80.0;
65 0 : CarinaSagittariusArm_x_shift_kpc = 1.37;
66 0 : CarinaSagittariusArm_y_shift_kpc = 0.0;
67 0 : CarinaSagittariusArm_arc_radius1_kpc = 1.0;
68 0 : CarinaSagittariusArm_arc_radius2_kpc = 0.79;
69 0 : CarinaSagittariusArm_arc_eps = 2.3;
70 0 : CarinaSagittariusArm_arc_div_deg = 3.0;
71 0 : CarinaSagittariusArm_rmin_kpc = 3.0;
72 0 : CarinaSagittariusArm_rmax_kpc = 15.0;
73 0 : CarinaSagittariusArm_zmin_kpc = -1.2;
74 0 : CarinaSagittariusArm_zmax_kpc = 1.2;
75 :
76 0 : LocalArm_B_gauss = -3.5e-6;
77 0 : LocalArm_pitch_deg = 20;
78 0 : LocalArm_phi0_deg = -2.2;
79 0 : LocalArm_x_shift_kpc = -0.15;
80 0 : LocalArm_y_shift_kpc = 0.0;
81 0 : LocalArm_arc_radius1_kpc = 0.73;
82 0 : LocalArm_arc_radius2_kpc = 1.0;
83 0 : LocalArm_arc_eps = 1.45;
84 0 : LocalArm_arc_div_deg = 0.0;
85 0 : LocalArm_rmin_kpc = 3.3;
86 0 : LocalArm_rmax_kpc = 14.0;
87 0 : LocalArm_zmin_kpc = -1.0;
88 0 : LocalArm_zmax_kpc = 1.0;
89 :
90 0 : PerseusArm1_B_gauss = -2.0e-6;
91 0 : PerseusArm1_pitch_deg = 20;
92 0 : PerseusArm1_phi0_deg = 46.0;
93 0 : PerseusArm1_x_shift_kpc = -1.0;
94 0 : PerseusArm1_y_shift_kpc = 0.0;
95 0 : PerseusArm1_arc_radius1_kpc = 0.4;
96 0 : PerseusArm1_arc_radius2_kpc = 1.1;
97 0 : PerseusArm1_arc_eps = 2.0;
98 0 : PerseusArm1_arc_div_deg = 0.0;
99 0 : PerseusArm1_rmin_kpc = 3.0;
100 0 : PerseusArm1_rmax_kpc = 17.0;
101 0 : PerseusArm1_zmin_kpc = -1.0;
102 0 : PerseusArm1_zmax_kpc = 1.0;
103 :
104 0 : PerseusArm2_B_gauss = -3.5e-6;
105 0 : PerseusArm2_pitch_deg = 20;
106 0 : PerseusArm2_phi0_deg = 46.0;
107 0 : PerseusArm2_x_shift_kpc = -1.0;
108 0 : PerseusArm2_y_shift_kpc = 0.0;
109 0 : PerseusArm2_arc_radius1_kpc = 1.2;
110 0 : PerseusArm2_arc_radius2_kpc = 1.1;
111 0 : PerseusArm2_arc_eps = 2.0;
112 0 : PerseusArm2_arc_div_deg = 0.0;
113 0 : PerseusArm2_rmin_kpc = 3.0;
114 0 : PerseusArm2_rmax_kpc = 17.0;
115 0 : PerseusArm2_zmin_kpc = -2.0;
116 0 : PerseusArm2_zmax_kpc = 2.0;
117 0 : }
118 :
119 :
120 :
121 0 : Vector3d KST24Field::getField(const Vector3d& pos) const
122 : {
123 : Vector3d vals(0, 0, 0);
124 : Vector3d pos_kpc = pos/kpc;
125 :
126 0 : if (pos_kpc.z > north_tor_zmin_kpc) // northern halo height lower boundary
127 : {
128 0 : vals += get_toroidal(pos_kpc, north_tor_B_gauss, north_tor_zmin_kpc, north_tor_zmax_kpc,
129 0 : north_tor_rmin_kpc, north_tor_rmax_kpc);
130 0 : vals += get_Xfield(pos_kpc, Xfield_B_gauss, Xfield_rmin_kpc, Xfield_rmax_kpc, Xfield_theta_rad);
131 : return vals*gauss;
132 : }
133 :
134 0 : if (pos_kpc.z < south_tor_zmax_kpc) // southern halo height upper boundary
135 : {
136 0 : vals += get_toroidal(pos_kpc, south_tor_B_gauss, south_tor_zmin_kpc, south_tor_zmax_kpc,
137 0 : south_tor_rmin_kpc, south_tor_rmax_kpc);
138 0 : vals += get_Xfield(pos_kpc, Xfield_B_gauss, Xfield_rmin_kpc, Xfield_rmax_kpc, Xfield_theta_rad);
139 : return vals*gauss;
140 : }
141 :
142 0 : if (is_LB(pos_kpc, LB_rmin_kpc, LB_dr_kpc, LB_x0_kpc, LB_y0_kpc, LB_z0_kpc))
143 : {
144 0 : vals += get_LB(pos_kpc, LB_B_gauss, LB_lB_deg, LB_bB_deg,
145 0 : LB_rmin_kpc, LB_dr_kpc, LB_x0_kpc, LB_y0_kpc, LB_z0_kpc);
146 : return vals*gauss;
147 : }
148 :
149 0 : vals += get_logspiral(pos_kpc, ScutumArm_B_gauss, ScutumArm_pitch_deg, ScutumArm_phi0_deg,
150 0 : ScutumArm_x_shift_kpc, ScutumArm_y_shift_kpc, ScutumArm_arc_radius1_kpc,
151 0 : ScutumArm_arc_radius2_kpc, ScutumArm_arc_eps, ScutumArm_arc_div_deg,
152 0 : ScutumArm_rmin_kpc, ScutumArm_rmax_kpc, ScutumArm_zmin_kpc, ScutumArm_zmax_kpc);
153 :
154 0 : vals += get_logspiral(pos_kpc, CarinaSagittariusArm_B_gauss, CarinaSagittariusArm_pitch_deg, CarinaSagittariusArm_phi0_deg,
155 0 : CarinaSagittariusArm_x_shift_kpc, CarinaSagittariusArm_y_shift_kpc, CarinaSagittariusArm_arc_radius1_kpc,
156 0 : CarinaSagittariusArm_arc_radius2_kpc, CarinaSagittariusArm_arc_eps, CarinaSagittariusArm_arc_div_deg,
157 0 : CarinaSagittariusArm_rmin_kpc, CarinaSagittariusArm_rmax_kpc, CarinaSagittariusArm_zmin_kpc, CarinaSagittariusArm_zmax_kpc);
158 :
159 0 : vals += get_logspiral(pos_kpc, LocalArm_B_gauss, LocalArm_pitch_deg, LocalArm_phi0_deg,
160 0 : LocalArm_x_shift_kpc, LocalArm_y_shift_kpc, LocalArm_arc_radius1_kpc,
161 0 : LocalArm_arc_radius2_kpc, LocalArm_arc_eps, LocalArm_arc_div_deg,
162 0 : LocalArm_rmin_kpc, LocalArm_rmax_kpc, LocalArm_zmin_kpc, LocalArm_zmax_kpc);
163 :
164 0 : vals += get_logspiral(pos_kpc, PerseusArm1_B_gauss, PerseusArm1_pitch_deg, PerseusArm1_phi0_deg,
165 0 : PerseusArm1_x_shift_kpc, PerseusArm1_y_shift_kpc, PerseusArm1_arc_radius1_kpc,
166 0 : PerseusArm1_arc_radius2_kpc, PerseusArm1_arc_eps, PerseusArm1_arc_div_deg,
167 0 : PerseusArm1_rmin_kpc, PerseusArm1_rmax_kpc, PerseusArm1_zmin_kpc, PerseusArm1_zmax_kpc);
168 :
169 0 : vals += get_logspiral(pos_kpc, PerseusArm2_B_gauss, PerseusArm2_pitch_deg, PerseusArm2_phi0_deg,
170 0 : PerseusArm2_x_shift_kpc, PerseusArm2_y_shift_kpc, PerseusArm2_arc_radius1_kpc,
171 0 : PerseusArm2_arc_radius2_kpc, PerseusArm2_arc_eps, PerseusArm2_arc_div_deg,
172 0 : PerseusArm2_rmin_kpc, PerseusArm2_rmax_kpc, PerseusArm2_zmin_kpc, PerseusArm2_zmax_kpc);
173 :
174 0 : vals += get_Xfield(pos_kpc, Xfield_B_gauss, Xfield_rmin_kpc, Xfield_rmax_kpc, Xfield_theta_rad);
175 :
176 : return vals*gauss;
177 : }
178 :
179 0 : bool KST24Field::is_LB(const Vector3d pos_kpc, const double LB_rmin_kpc, const double LB_dr_kpc,
180 : const double LB_x0_kpc, const double LB_y0_kpc, const double LB_z0_kpc) const
181 : {
182 : double eRx, eRy, eRz;
183 : double cR;
184 :
185 0 : eRx = pos_kpc.x - LB_x0_kpc;
186 0 : eRy = pos_kpc.y - LB_y0_kpc;
187 0 : eRz = pos_kpc.z - LB_z0_kpc;
188 0 : cR = sqrt(eRx*eRx + eRy*eRy + eRz*eRz);
189 0 : if (cR < LB_rmin_kpc + LB_dr_kpc)
190 0 : return true;
191 : return false;
192 : }
193 :
194 0 : Vector3d KST24Field::get_toroidal(const Vector3d pos_kpc, const double tor_B_gauss,
195 : const double tor_zmin_kpc, const double tor_zmax_kpc,
196 : const double tor_rmin_kpc, const double tor_rmax_kpc) const
197 : {
198 : Vector3d vals_gauss(0, 0, 0);
199 0 : if ((pos_kpc.z <= tor_zmin_kpc) or (tor_zmax_kpc <= pos_kpc.z))
200 : return vals_gauss;
201 :
202 : double cR, theta;
203 0 : cR = sqrt(pow(pos_kpc.x, 2) + pow(pos_kpc.y, 2));
204 0 : if ((tor_rmin_kpc <= cR) and (cR <= tor_rmax_kpc))
205 : {
206 0 : theta = atan2(pos_kpc.y, pos_kpc.x);
207 0 : vals_gauss.x = tor_B_gauss*sin(theta);
208 0 : vals_gauss.y = -tor_B_gauss*cos(theta);
209 : vals_gauss.z = 0;
210 : }
211 : return vals_gauss;
212 : }
213 :
214 0 : Vector3d KST24Field::get_Xfield(const Vector3d pos_kpc, const double Xfield_B_gauss,
215 : const double Xfield_rmin_kpc, const double Xfield_rmax_kpc,
216 : const double Xfield_theta_rad) const
217 : {
218 : Vector3d vals_gauss(0, 0, 0);
219 :
220 : int sign = 1;
221 0 : if (pos_kpc.z < 0)
222 : sign = -1;
223 0 : if (fabs(pos_kpc.z) > 10)
224 : return vals_gauss;
225 :
226 : double cR, cR_0;
227 0 : cR = sqrt(pos_kpc.x*pos_kpc.x + pos_kpc.y*pos_kpc.y);
228 0 : cR_0 = cR - pos_kpc.z*sign*tan(Xfield_theta_rad);
229 :
230 0 : if ((cR_0 < Xfield_rmin_kpc) or (Xfield_rmax_kpc < cR_0))
231 : return vals_gauss;
232 :
233 0 : double mgn_frc = cR_0/cR;
234 0 : double xy_theta = atan2(pos_kpc.y, pos_kpc.x);
235 0 : vals_gauss.x = Xfield_B_gauss*mgn_frc*cos(xy_theta)*sign*sin(Xfield_theta_rad);
236 0 : vals_gauss.y = Xfield_B_gauss*mgn_frc*sin(xy_theta)*sign*sin(Xfield_theta_rad);
237 0 : vals_gauss.z = Xfield_B_gauss*mgn_frc*cos(Xfield_theta_rad);
238 :
239 0 : return vals_gauss;
240 : }
241 :
242 0 : Vector3d KST24Field::get_LB(const Vector3d pos_kpc, const double LB_B_gauss,
243 : const double LB_lB_deg, const double LB_bB_deg,
244 : const double LB_rmin_kpc, const double LB_dr_kpc,
245 : const double LB_x0_kpc, const double LB_y0_kpc, const double LB_z0_kpc) const
246 : {
247 : Vector3d vals_gauss(0, 0, 0);
248 :
249 : double eRx, eRy, eRz, cR;
250 0 : eRx = pos_kpc.x - LB_x0_kpc;
251 0 : eRy = pos_kpc.y - LB_y0_kpc;
252 0 : eRz = pos_kpc.z - LB_z0_kpc;
253 0 : cR = sqrt(eRx*eRx + eRy*eRy + eRz*eRz);
254 0 : if ((cR < LB_rmin_kpc) or (LB_rmin_kpc + LB_dr_kpc < cR))
255 : return vals_gauss;
256 0 : eRx /= cR;
257 0 : eRy /= cR;
258 0 : eRz /= cR;
259 :
260 : double fdir_x, fdir_y, fdir_z, cosTheta;
261 0 : cosTheta = eRx*LB_Bdir_x + eRy*LB_Bdir_y + eRz*LB_Bdir_z; // expanding cross product: b(ac) - c(ab)
262 0 : fdir_x = LB_Bdir_x - eRx*cosTheta;
263 0 : fdir_y = LB_Bdir_y - eRy*cosTheta;
264 0 : fdir_z = LB_Bdir_z - eRz*cosTheta;
265 :
266 : // magnetic field amplification factor
267 : double ampl, ampl_elec;
268 0 : ampl = (1 + pow(LB_rmin_kpc, 2)/(2*LB_rmin_kpc*LB_dr_kpc + LB_dr_kpc*LB_dr_kpc));
269 :
270 0 : vals_gauss.x = LB_B_gauss*ampl*fdir_x;
271 0 : vals_gauss.y = LB_B_gauss*ampl*fdir_y;
272 0 : vals_gauss.z = LB_B_gauss*ampl*fdir_z;
273 :
274 0 : return vals_gauss;
275 : }
276 :
277 0 : Vector3d KST24Field::get_logspiral(const Vector3d pos_kpc, const double B_gauss,
278 : const double pitch_deg, const double phi0_deg,
279 : const double x_shift_kpc, const double y_shift_kpc,
280 : const double arc_radius1_kpc, const double arc_radius2_kpc,
281 : const double arc_eps , const double arc_div_deg,
282 : const double rmin_kpc, const double rmax_kpc,
283 : const double zmin_kpc, const double zmax_kpc) const
284 : {
285 : Vector3d vals_gauss(0, 0, 0);
286 :
287 0 : if ((pos_kpc.z < zmin_kpc) or (zmax_kpc < pos_kpc.z))
288 : return vals_gauss;
289 :
290 : std::vector<double> pos_v;
291 0 : pos_v.push_back(pos_kpc.x + x_shift_kpc);
292 0 : pos_v.push_back(pos_kpc.y + y_shift_kpc);
293 0 : pos_v.push_back(pos_kpc.z);
294 :
295 0 : double r = sqrt(pos_v[0]*pos_v[0] + pos_v[1]*pos_v[1]);
296 0 : if ((r < rmin_kpc) or (rmax_kpc < r))
297 : return vals_gauss;
298 :
299 : double a_kpc = 3;
300 0 : double k = tan(pitch_deg*M_PI/180.);
301 0 : double cos_pitch = cos(pitch_deg*M_PI/180.);
302 0 : double sin_pitch = sin(pitch_deg*M_PI/180.);
303 0 : double phi0 = phi0_deg*M_PI/180.;
304 0 : double arc_div_rad = arc_div_deg*M_PI/180.;
305 :
306 : int nn;
307 0 : double phi = atan2(pos_v[1], pos_v[0]);
308 0 : nn = floor((log(r/a_kpc) - k*(phi + phi0))/(2*M_PI*k));
309 :
310 : double r1, r2;
311 0 : r1 = a_kpc*exp(k*(phi + phi0 + 2*M_PI*nn));
312 0 : r2 = r1*exp(k*2*M_PI);
313 :
314 0 : if (fabs(r - r1) > fabs(r - r2))
315 : {
316 : r1 = r2;
317 0 : nn++; // we need outer arm
318 : }
319 :
320 :
321 : // perpendicular to the spiral arm axis
322 : double xi, yi, dir1, dir2, dirtmp;
323 0 : xi = r1*cos(phi);
324 0 : yi = r1*sin(phi);
325 0 : dir1 = k*xi - yi;
326 0 : dir2 = k*yi + xi;
327 0 : dirtmp = sqrt(dir1*dir1 + dir2*dir2);
328 0 : dir1 /= dirtmp;
329 0 : dir2 /= dirtmp;
330 : // std::cout << ((xi - pos[0])*dir1 + (yi - pos[1])*dir2)/sqrt(pow((xi - pos[0]), 2) + pow((yi - pos[1]), 2)) << '\n';
331 :
332 :
333 : // first order correction
334 : double delta_phi;
335 0 : delta_phi = -((xi - pos_v[0])*dir1 + (yi - pos_v[1])*dir2)/(r1/cos_pitch);
336 0 : r1 = a_kpc*exp(k*(phi + phi0 + delta_phi + 2*M_PI*nn));
337 0 : xi = r1*cos(phi + delta_phi);
338 0 : yi = r1*sin(phi + delta_phi);
339 0 : dir1 = k*xi - yi;
340 0 : dir2 = k*yi + xi;
341 0 : dirtmp = sqrt(dir1*dir1 + dir2*dir2);
342 0 : dir1 /= dirtmp;
343 0 : dir2 /= dirtmp;
344 :
345 :
346 : double d1, L1, d_at_L1;
347 0 : d1 = sqrt(pow((xi - pos_v[0]), 2) + pow((yi - pos_v[1]), 2));
348 0 : L1 = r1/sin_pitch;
349 :
350 : double r_scale = 5;
351 : double L1_at_r_scale, arc_r_plane, arc_r_z;
352 0 : L1_at_r_scale = r_scale/sin_pitch;
353 0 : arc_r_plane = std::max(arc_radius2_kpc, arc_radius2_kpc*(1 + (L1 - L1_at_r_scale)*arc_div_rad));
354 : // arc_r_plane = arc_radius2_kpc;
355 0 : arc_r_z = std::max(arc_radius1_kpc, arc_radius1_kpc*(1 + (L1 - L1_at_r_scale)*arc_div_rad));
356 0 : if (arc_r_plane > 1.2)
357 : arc_r_plane = 1.2;
358 0 : if (arc_r_z > 1.2)
359 : arc_r_z = 1.2;
360 :
361 :
362 0 : d_at_L1 = sqrt(pow(fabs(d1/arc_r_plane), arc_eps) + pow(fabs(pos_v[2]/arc_r_z), arc_eps));
363 0 : if ((d_at_L1 > 1) or (d_at_L1 < 0.))
364 : return vals_gauss;
365 : else
366 : {
367 0 : double scaling_factor = (arc_radius1_kpc*arc_radius2_kpc)/(arc_r_plane*arc_r_z);
368 0 : vals_gauss.x = B_gauss*dir1*scaling_factor;
369 0 : vals_gauss.y = B_gauss*dir2*scaling_factor;
370 0 : vals_gauss.z = 0;
371 : }
372 :
373 0 : return vals_gauss;
374 0 : }
375 : }
|
