Line data Source code
1 : #include "crpropa/magneticField/TF17Field.h"
2 : #include "crpropa/Units.h"
3 :
4 : #include <algorithm>
5 : #include <string>
6 :
7 : namespace crpropa {
8 : using namespace std;
9 :
10 0 : TF17Field::TF17Field(TF17DiskModel disk_model_, TF17HaloModel halo_model_) {
11 0 : disk_model = disk_model_;
12 0 : halo_model = halo_model_;
13 :
14 0 : useHaloField = true;
15 0 : useDiskField = true;
16 :
17 0 : if ((halo_model == TF17HaloModel::C0) && (disk_model == TF17DiskModel::Ad1)) {
18 : // disk parameters
19 0 : set_r1_disk(3 * kpc);
20 0 : set_B1_disk(19.0 * muG);
21 0 : set_H_disk(0.055 * kpc);
22 0 : set_phi_star_disk(-54 * M_PI / 180);
23 0 : set_a_disk(0.9 / kpc / kpc);
24 : // halo parameters
25 0 : set_z1_halo(0);
26 0 : set_B1_halo(0.36 * muG);
27 0 : set_L_halo(3.0 * kpc);
28 0 : set_a_halo(1.17 / kpc / kpc);
29 : // shared parameters
30 0 : set_p0(-7.9 * M_PI / 180);
31 0 : set_Hp(5 * kpc);
32 0 : set_Lp(50 * kpc); // > 18 kpc
33 :
34 0 : } else if ((halo_model == TF17HaloModel::C0) && (disk_model == TF17DiskModel::Bd1)) {
35 : // disk parameters
36 0 : set_r1_disk(3 * kpc);
37 0 : set_B1_disk(2.0 * muG);
38 0 : set_H_disk(0.32 * kpc);
39 0 : set_phi_star_disk(153 * M_PI / 180);
40 : // halo parameters
41 0 : set_z1_halo(0);
42 0 : set_B1_halo(0.29 * muG);
43 0 : set_L_halo(3.4 * kpc);
44 0 : set_a_halo(0.88 / kpc / kpc);
45 : // shared parameters
46 0 : set_p0(-7.2 * M_PI / 180);
47 0 : set_Hp(9 * kpc);
48 0 : set_Lp(50 * kpc); // > 16 kpc
49 :
50 0 : } else if ((halo_model == TF17HaloModel::C0) && (disk_model == TF17DiskModel::Dd1)) {
51 : // disk parameters
52 0 : set_z1_disk(1.5 * kpc);
53 0 : set_B1_disk(0.065 * muG);
54 0 : set_L_disk(9.8 * kpc);
55 0 : set_phi_star_disk(14 * M_PI / 180);
56 : // halo parameters
57 0 : set_z1_halo(0);
58 0 : set_B1_halo(0.18 * muG);
59 0 : set_L_halo(4.8 * kpc);
60 0 : set_a_halo(0.61 / kpc / kpc);
61 : // shared parameters
62 0 : set_p0(-7.4 * M_PI / 180);
63 0 : set_Hp(4.2 * kpc);
64 0 : set_Lp(50 * kpc); // > 22 kpc
65 :
66 0 : } else if ((halo_model == TF17HaloModel::C1) && (disk_model == TF17DiskModel::Ad1)) {
67 : // disk parameters
68 0 : set_r1_disk(3 * kpc);
69 0 : set_B1_disk(32.0 * muG);
70 0 : set_H_disk(0.054 * kpc);
71 0 : set_phi_star_disk(-31 * M_PI / 180);
72 0 : set_a_disk(0.031 / kpc / kpc);
73 : // halo parameters
74 0 : set_z1_halo(0);
75 0 : set_B1_halo(9.0 * muG);
76 0 : set_L_halo(2.1 * kpc);
77 0 : set_phi_star_halo(198 * M_PI / 180);
78 0 : set_a_halo(0.33 / kpc / kpc);
79 : // shared parameters
80 0 : set_p0(-9.1 * M_PI / 180);
81 0 : set_Hp(1.2 * kpc);
82 0 : set_Lp(50 * kpc); // > 38 kpc
83 :
84 0 : } else if ((halo_model == TF17HaloModel::C1) && (disk_model == TF17DiskModel::Bd1)) {
85 : // disk parameters
86 0 : set_r1_disk(3 * kpc);
87 0 : set_B1_disk(24 * muG);
88 0 : set_H_disk(0.090 * kpc);
89 0 : set_phi_star_disk(-34 * M_PI / 180);
90 : // halo parameters
91 0 : set_z1_halo(0);
92 0 : set_B1_halo(8.2 * muG);
93 0 : set_L_halo(2.2 * kpc);
94 0 : set_phi_star_halo(197 * M_PI / 180);
95 0 : set_a_halo(0.38 / kpc / kpc);
96 : // shared parameters
97 0 : set_p0(-9.0 * M_PI / 180);
98 0 : set_Hp(1.2 * kpc);
99 0 : set_Lp(50 * kpc); // > 38 kpc
100 :
101 0 : } else if ((halo_model == TF17HaloModel::C1) && (disk_model == TF17DiskModel::Dd1)) {
102 : // disk parameters
103 0 : set_z1_disk(1.5 * kpc);
104 0 : set_B1_disk(0.40 * muG);
105 0 : set_L_disk(2.9 * kpc);
106 0 : set_phi_star_disk(120 * M_PI / 180);
107 : // halo parameters
108 0 : set_z1_halo(0);
109 0 : set_B1_halo(9.5 * muG);
110 0 : set_L_halo(2.1 * kpc);
111 0 : set_phi_star_halo(179 * M_PI / 180);
112 0 : set_a_halo(0.45 / kpc / kpc);
113 : // shared parameters
114 0 : set_p0(-8.4 * M_PI / 180);
115 0 : set_Hp(1.2 * kpc);
116 0 : set_Lp(50 * kpc); // > 30 kpc
117 : }
118 :
119 0 : epsilon = std::numeric_limits<double>::epsilon();
120 0 : }
121 :
122 0 : void TF17Field::setUseDiskField(bool use) { useDiskField = use; }
123 0 : void TF17Field::setUseHaloField(bool use) { useHaloField = use; }
124 :
125 0 : bool TF17Field::isUsingDiskField() { return useDiskField; }
126 0 : bool TF17Field::isUsingHaloField() { return useHaloField; }
127 :
128 0 : void TF17Field::set_B1_disk(const double B1){ B1_disk = B1; }
129 0 : void TF17Field::set_z1_disk(const double z1){ z1_disk = z1; }
130 0 : void TF17Field::set_r1_disk(const double r1){ r1_disk = r1; }
131 0 : void TF17Field::set_H_disk(const double H){ H_disk = H; }
132 0 : void TF17Field::set_L_disk(const double L){ L_disk = L; }
133 0 : void TF17Field::set_a_disk(const double a){ a_disk = a; }
134 0 : void TF17Field::set_phi_star_disk(const double phi){ phi_star_disk = phi; }
135 :
136 0 : void TF17Field::set_B1_halo(const double B1){ B1_halo = B1; }
137 0 : void TF17Field::set_z1_halo(const double z1){ z1_halo = z1; }
138 0 : void TF17Field::set_L_halo(const double L){ L_halo = L; }
139 0 : void TF17Field::set_a_halo(const double a){ a_halo = a; }
140 0 : void TF17Field::set_phi_star_halo(const double phi){ phi_star_halo = phi; }
141 :
142 0 : void TF17Field::set_Hp(const double H){ H_p = H; }
143 0 : void TF17Field::set_Lp(const double L){ L_p = L; }
144 0 : void TF17Field::set_p0(const double p0){
145 0 : p_0 = p0;
146 0 : cot_p0 = cos(p_0) / sin(p_0);
147 0 : }
148 :
149 0 : string TF17Field::getDiskModel() const {
150 : string model_name;
151 0 : switch (disk_model) {
152 : case TF17DiskModel::Ad1 : model_name = "Ad1"; break;
153 : case TF17DiskModel::Bd1 : model_name = "Bd1"; break;
154 : case TF17DiskModel::Dd1 : model_name = "Dd1"; break;
155 : }
156 0 : return model_name;
157 : }
158 0 : string TF17Field::getHaloModel() const {
159 : string model_name;
160 0 : switch (halo_model) {
161 : case TF17HaloModel::C0 : model_name = "C0"; break;
162 : case TF17HaloModel::C1 : model_name = "C1"; break;
163 : }
164 0 : return model_name;
165 : }
166 :
167 :
168 0 : Vector3d TF17Field::getField(const Vector3d& pos) const {
169 0 : double r = sqrt(pos.x * pos.x + pos.y * pos.y); // in-plane radius
170 0 : double phi = M_PI - pos.getPhi(); // azimuth in our convention
171 : // double cosPhi = pos.x / r;
172 0 : double cosPhi = cos(phi);
173 : // double sinPhi = pos.y / r;
174 0 : double sinPhi = sin(phi);
175 :
176 : Vector3d b(0.);
177 0 : if (useDiskField)
178 0 : b += getDiskField(r, pos.z, phi, sinPhi, cosPhi); // disk field
179 0 : if (useHaloField)
180 0 : b += getHaloField(r, pos.z, phi, sinPhi, cosPhi); // halo field
181 0 : return b;
182 : }
183 :
184 0 : Vector3d TF17Field::getDiskField(const double& r, const double& z, const double& phi, const double& sinPhi, const double& cosPhi) const {
185 : Vector3d b(0.);
186 0 : double B_r = 0;
187 : double B_phi = 0;
188 0 : double B_z = 0;
189 :
190 0 : if (disk_model == TF17DiskModel::Ad1) { // ==========================================================
191 0 : if (r > r1_disk) {
192 0 : double z1_disk_z = (1. + a_disk * r1_disk * r1_disk) / (1. + a_disk * r * r); // z1_disk / z
193 : // B components in (r, phi, z)
194 0 : double B_r0 = radialFieldScale(B1_disk, phi_star_disk, z1_disk_z*z, phi, r, z);
195 0 : B_r = (r1_disk / r) * z1_disk_z * B_r0;
196 0 : B_z = 2 * a_disk * r1_disk * z1_disk_z*z / (1+ a_disk * r * r) * B_r0;
197 0 : B_phi = azimuthalFieldComponent(r, z, B_r, B_z);
198 : } else {
199 : // within r = r1_disk, the field lines are straight in direction g_phi + phi_star_disk
200 : // and thus z = z1
201 0 : double phi1_disk = shiftedWindingFunction(r1_disk, z) + phi_star_disk;
202 0 : double B_amp = B1_disk * exp(-fabs(z) / H_disk);
203 0 : B_r = cos(phi1_disk - phi) * B_amp;
204 0 : B_phi = sin(phi1_disk - phi) * B_amp;
205 : }
206 :
207 0 : } else if (disk_model == TF17DiskModel::Bd1) { // ===================================================
208 : // for model Bd1, best fit for n = 2
209 0 : if ( r > epsilon ) {
210 0 : double r1_disk_r = r1_disk / r;
211 0 : double z1_disk_z = 5. / (r1_disk_r*r1_disk_r + 4./sqrt(r1_disk_r)); // z1_disk / z -> remove z dependancy
212 0 : double B_r0 = radialFieldScale(B1_disk, phi_star_disk, z1_disk_z*z, phi, r, z);
213 0 : B_r = r1_disk_r * z1_disk_z * B_r0;
214 0 : B_z = -0.4 * r1_disk_r / r * z1_disk_z* z1_disk_z * z * (r1_disk_r*r1_disk_r - 1./sqrt(r1_disk_r)) * B_r0;
215 : } else {
216 0 : double z1_disk_z = 5. * r*r / (r1_disk*r1_disk); // z1_disk / z -> remove z dependancy
217 0 : double B_r0 = radialFieldScale(B1_disk, phi_star_disk, z1_disk_z*z, phi, r, z);
218 0 : B_r = 5. * r / r1_disk * B_r0;
219 0 : B_z = -10. * z / r1_disk * B_r0;
220 : }
221 0 : B_phi = azimuthalFieldComponent(r, z, B_r, B_z);
222 :
223 0 : } else if (disk_model == TF17DiskModel::Dd1) { // ===================================================
224 : // for model Dd1, best fit for n = 0.5
225 0 : double z_sign = z >= 0 ? 1. : -1.;
226 0 : double z_abs = fabs(z);
227 0 : if ( z_abs > epsilon ) {
228 0 : double z1_disk_z = z1_disk / z_abs;
229 0 : double r1_disk_r = 1.5 / (sqrt(z1_disk_z) + 0.5/z1_disk_z); // r1_disk / r
230 0 : double F_r = r1_disk_r*r <= L_disk ? 1. : exp(1. - r1_disk_r*r/L_disk);
231 : // simplication of the equation in the cosinus
232 0 : double B_z0 = z_sign * B1_disk * F_r * cos(phi - shiftedWindingFunction(r, z) - phi_star_disk);
233 0 : B_r = -0.5/1.5 * r1_disk_r * r1_disk_r * r1_disk_r * r / z_abs * (sqrt(z1_disk_z) - 1/z1_disk_z) * B_z0;
234 0 : B_z = z_sign * r1_disk_r * r1_disk_r * B_z0;
235 : } else {
236 0 : double z_z1_disk = z_abs / z1_disk;
237 0 : double r1_disk_r = 1.5 * sqrt(z_abs / z1_disk); // r1_disk / r
238 0 : double F_r = r1_disk_r*r <= L_disk ? 1. : exp(1. - r1_disk_r*r/L_disk);
239 0 : double B_z0 = z_sign * B1_disk * F_r * cos(phi - shiftedWindingFunction(r, z) - phi_star_disk);
240 0 : B_r = -1.125 * r / z1_disk * (1 - 2.5 * z_z1_disk * sqrt(z_z1_disk)) * B_z0;
241 0 : B_z = z_sign * r1_disk_r * r1_disk_r * B_z0;
242 : }
243 0 : B_phi = azimuthalFieldComponent(r, z, B_r, B_z);
244 : }
245 :
246 : // Convert to (x, y, z) components
247 0 : b.x = - (B_r * cosPhi - B_phi * sinPhi); // flip x-component at the end
248 0 : b.y = B_r * sinPhi + B_phi * cosPhi;
249 0 : b.z = B_z;
250 0 : return b;
251 : }
252 :
253 0 : Vector3d TF17Field::getHaloField(const double& r, const double& z, const double& phi, const double& sinPhi, const double& cosPhi) const {
254 : int m;
255 :
256 : Vector3d b(0.);
257 0 : double r1_halo_r = (1. + a_halo * z1_halo * z1_halo) / (1. + a_halo * z * z);
258 : // B components in (r, phi, z)
259 : double B_z0;
260 :
261 0 : if (halo_model == TF17HaloModel::C0) { // m = 0
262 0 : B_z0 = B1_halo * exp(-r1_halo_r*r / L_halo);
263 0 : } else if (halo_model == TF17HaloModel::C1) { // m = 1
264 : // simplication of the equation in the cosinus
265 0 : double phi_prime = phi - shiftedWindingFunction(r, z) - phi_star_halo;
266 0 : B_z0 = B1_halo * exp(-r1_halo_r*r / L_halo) * cos(phi_prime);
267 : }
268 :
269 : // Contrary to article, Br has been rewriten to a little bit by replacing
270 : // (2 * a * r1**3 * z) / (r**2) by (2 * a * r1**2 * z) / (r * (1+a*z**2))
271 : // but that is strictly equivalent except we can reintroduce the z1 in the expression via r1
272 0 : double B_r = 2 * a_halo * r1_halo_r * r1_halo_r * r * z / (1. + a_halo * z * z) * B_z0;
273 0 : double B_z = r1_halo_r * r1_halo_r * B_z0;
274 0 : double B_phi = azimuthalFieldComponent(r, z, B_r, B_z);
275 :
276 : // Convert to (x, y, z) components
277 0 : b.x = - (B_r * cosPhi - B_phi * sinPhi); // flip x-component at the end
278 0 : b.y = B_r * sinPhi + B_phi * cosPhi;
279 0 : b.z = B_z;
280 :
281 0 : return b;
282 : }
283 :
284 0 : double TF17Field::azimuthalFieldComponent(const double& r, const double& z, const double& B_r, const double& B_z) const {
285 0 : double r_ = r / L_p;
286 0 : double rscale = r > epsilon ? r_ * exp(-r_) / (1 - exp(-r_)) : 1 - r_/2. - r_*r_/12. ;
287 0 : double B_phi = cot_p0 / zscale(z) * rscale * B_r;
288 0 : B_phi = B_phi - 2 * z * r / (H_p * H_p) / zscale(z) * shiftedWindingFunction(r, z) * B_z;
289 0 : return B_phi;
290 : }
291 :
292 0 : double TF17Field::radialFieldScale(const double& B1, const double& phi_star, const double& z1, const double& phi, const double& r, const double& z) const {
293 : // simplication of the equation in the cosinus
294 0 : double phi_prime = phi - shiftedWindingFunction(r, z) - phi_star;
295 : // This term occures is parameterizations of models A and B always bisymmetric (m = 1)
296 0 : return B1 * exp(-fabs(z1) / H_disk) * cos(phi_prime);
297 : }
298 :
299 0 : double TF17Field::shiftedWindingFunction(const double& r, const double& z) const {
300 0 : return cot_p0 * log(1 - exp(-r / L_p) + epsilon) / zscale(z);
301 : }
302 :
303 0 : double TF17Field::zscale(const double& z) const {
304 0 : return 1 + z * z / H_p / H_p;
305 : }
306 :
307 : } // namespace crpropa
|
