Line data Source code
1 : #include "crpropa/magneticField/CMZField.h"
2 : #include "crpropa/Units.h"
3 :
4 : namespace crpropa {
5 :
6 5 : CMZField::CMZField() {
7 5 : useMCField = false;
8 5 : useICField = true;
9 5 : useNTFField = false;
10 5 : useRadioArc = false;
11 5 : }
12 :
13 8 : double CMZField::getA(double a1) const {
14 8 : return 4*log(2)/a1/a1;
15 : }
16 :
17 8 : double CMZField::getL(double a2) const {
18 8 : return a2/2*log(2);
19 : }
20 :
21 9 : Vector3d CMZField::BPol(const Vector3d& position,const Vector3d& mid, double B1, double a, double L) const{
22 : // cylindircal coordinates
23 : Vector3d pos = position - mid;
24 9 : double r = sqrt(pos.x*pos.x + pos.y*pos.y);
25 9 : double phi = std::atan2(pos.y, pos.x);
26 :
27 9 : double r1 = 1/(1+a*pos.z*pos.z);
28 9 : double Bs = B1*exp(-r1*r/L);
29 9 : double Br = 2*a*pow_integer<3>(r1)*r*pos.z*Bs;
30 :
31 : Vector3d b = Vector3d(0.);
32 9 : b.z = r1*r1*Bs;
33 : // transform to cartesian coordinates
34 9 : b.x = Br*cos(phi);
35 9 : b.y = Br*sin(phi);
36 9 : return b;
37 : }
38 :
39 14 : Vector3d CMZField::BAz(const Vector3d& position, const Vector3d& mid, double B1, double eta, double R) const {
40 : // cylindrical coordinates
41 : Vector3d pos = position - mid;
42 14 : double r = sqrt(pos.x*pos.x + pos.y*pos.y);
43 14 : double phi = std::atan2(pos.y,pos.x);
44 :
45 : Vector3d bVec(0.);
46 14 : double Hc = R/sqrt(log(2));
47 : double b = 1.;
48 : double m = 1;
49 14 : double r1 = R/10;
50 14 : double v = m/eta*log((r+b)/(R+b));
51 14 : double cosV = cos(v + m*phi);
52 :
53 : double Br=0;
54 : double Bphi=0;
55 :
56 14 : if(r>r1){
57 13 : double Pre = B1*cosV*exp(-pos.z*pos.z/Hc/Hc);
58 13 : Br = Pre*R/r;
59 13 : Bphi=-Pre/eta*R/(r+b);
60 : }
61 : else{
62 1 : double Pre = B1*exp(-pos.z*pos.z/Hc/Hc)*R/r1*(3*r/r1 - 2*r*r/r1/r1)*cosV;
63 : Br = Pre;
64 1 : Bphi = 1 + 6*(r-r1)/(2*r-3*r1)*(sin(v+m*phi)-sin(v))/cosV;
65 1 : Bphi *= -Pre*r/eta/(r+b);
66 : }
67 :
68 14 : bVec.x = Br*cos(phi) - Bphi*sin(phi);
69 14 : bVec.y = Br*sin(phi) + Bphi*cos(phi);
70 :
71 14 : return bVec;
72 : }
73 :
74 2 : bool CMZField::getUseMCField() const {
75 2 : return useMCField;
76 : }
77 2 : bool CMZField::getUseICField() const {
78 2 : return useICField;
79 : }
80 2 : bool CMZField::getUseNTFField() const {
81 2 : return useNTFField;
82 : }
83 2 : bool CMZField::getUseRadioArc() const {
84 2 : return useRadioArc;
85 : }
86 :
87 1 : void CMZField::setUseMCField(bool use) {
88 1 : useMCField = use;
89 1 : }
90 1 : void CMZField::setUseICField(bool use) {
91 1 : useICField = use;
92 1 : }
93 1 : void CMZField::setUseNTFField(bool use) {
94 1 : useNTFField = use;
95 1 : }
96 1 : void CMZField::setUseRadioArc(bool use) {
97 1 : useRadioArc = use;
98 1 : }
99 :
100 1 : Vector3d CMZField::getMCField(const Vector3d& pos) const {//Field in molecular clouds
101 : Vector3d b(0.);
102 : double eta=0.01;
103 : double N=59; // normalization factor, depends on eta
104 :
105 : // azimuthal component in dense clouds
106 : //A=SgrC
107 : Vector3d mid(0,-81.59*pc, -16.32*pc);
108 : double R = 1.7*pc;
109 : double B1=2.1e-3/N;
110 1 : b += BAz(pos, mid, B1, eta, R);
111 :
112 : // A=G0.253+0.016 Dust Ridge A
113 : mid = Vector3d(0, 37.53*pc, -2.37*pc);
114 : R=2.4*pc;
115 : B1=2.5e-3/N;
116 1 : b += BAz(pos, mid, B1, eta, R);
117 :
118 : //A=Dust Ridge B
119 : mid = Vector3d(0, 50.44, 8.16)*pc;
120 : R=1.9*pc;
121 : B1=0.9e-3/N;
122 1 : b += BAz(pos, mid, B1, eta, R);
123 :
124 : //A=Dust Ridge C
125 : mid = Vector3d(0,56.37,7.71)*pc;
126 : R=1.9*pc;
127 : B1=1.2e-3/N;
128 1 : b += BAz(pos, mid, B1, eta, R);
129 :
130 : //A=Dust Ridge D
131 : mid = Vector3d(0, 60.82, 7.42)*pc;
132 : R=3.3*pc;
133 : B1=1.7e-3/N;
134 1 : b += BAz(pos, mid, B1, eta, R);
135 :
136 : //A=Dust Ridge E
137 : mid = Vector3d(0, 70.91, 0.74)*pc;
138 : R=3.5*pc;
139 : B1=4.1e-3/N;
140 1 : b += BAz(pos, mid, B1, eta, R);
141 :
142 : //A=Dust Ridge F
143 : mid = Vector3d(0, 73.58, 2.97)*pc;
144 : R=2.4*pc;
145 : B1=3.9e-3/N;
146 1 : b += BAz(pos, mid, B1, eta, R);
147 :
148 : //Sgr D
149 : mid = Vector3d(0, 166.14, -10.38)*pc;
150 : R=1.8*pc;
151 : B1=0.8e-3/N;
152 1 : b += BAz(pos, mid, B1, eta, R);
153 :
154 : //Sgr B2
155 : mid = Vector3d(0, 97.01, -5.93)*pc;
156 : R=14*pc;
157 : B1=1.0e-3/N;
158 1 : b += BAz(pos, mid, B1, eta, R);
159 :
160 : //A=Inner R=5pc
161 : mid = Vector3d(0, -8.3, -6.9)*pc;
162 : R=5*pc;
163 : B1=3.0e-3/0.91;
164 1 : b += BAz(pos, mid, B1, 0.77, R); // different eta value!
165 :
166 : //20 km s^-1
167 : mid = Vector3d(0, -19.29,-11.87)*pc;
168 : R=9.4*pc;
169 : B1=2.7e-3/N;
170 1 : b += BAz(pos, mid, B1, eta, R);
171 :
172 : //50 km s^-1
173 : mid = Vector3d(0, -2.97, -10.38)*pc;
174 : R=9.4*pc;
175 : B1=3.7e-3/N;
176 1 : b += BAz(pos, mid, B1, eta, R);
177 :
178 : //SgrA* is different orrientated!
179 : //only phi component
180 1 : double x = pos.x;
181 1 : double y = pos.y + 8.3*pc;
182 1 : double z = pos.z + 6.9*pc;
183 : R=1.2e12;
184 : B1=65./3.07;
185 : double Hc = R/sqrt(log(2));
186 1 : double r = sqrt(x*x + y*y);
187 : double r1 = R/10;
188 1 : double phi= std::atan2(y,x);
189 : double Bphi;
190 :
191 1 : if(r>r1){
192 1 : Bphi = - B1*exp(-z*z/Hc/Hc)*R/r;
193 : }
194 : else{
195 0 : - B1*exp(-z*z/Hc/Hc)*R/r1*(3*r/r1- 2*r*r/r1*r1);
196 : }
197 :
198 1 : b.x -= Bphi*sin(phi);
199 1 : b.y += Bphi*cos(phi);
200 :
201 1 : return b*gauss;
202 : }
203 :
204 1 : Vector3d CMZField::getICField(const Vector3d& pos) const {//Field in intercloud medium--> poloidal field
205 : Vector3d mid(0.,-8.3*pc,-6.9*pc);
206 :
207 : double eta = 0.85;
208 : double B1 = 1e-5*gauss;
209 : double B2 = B1/eta;
210 : double a = 4*log(2)/pow(70*pc, 2);
211 : double L = 158*pc/log(2);
212 :
213 1 : return BPol(pos, mid, B2, a, L);
214 : }
215 :
216 1 : Vector3d CMZField::getNTFField(const Vector3d& pos) const {//Field in the non-thermal filaments--> predominantly poloidal field (except "pelical"-> azimuthal)
217 : Vector3d b(0.);
218 : Vector3d mid(0.);
219 :
220 : //A=SgrC
221 : mid = Vector3d(0., -81.59,-1.48)*pc;
222 : double a1=27.44*pc;
223 : double a2=1.73*pc;
224 : double eta=0.48;
225 : double B1=1.e-4;
226 1 : b += BPol(pos, mid, B1/eta, getA(a1), getL(a2));
227 :
228 : //A=G359.15-0.2 The Snake
229 : mid = Vector3d(0, -126.1,-25.22)*pc;
230 : a1=12.86*pc;
231 : a2=2.22*pc;
232 : B1=88.e-6;
233 1 : b += BPol(pos, mid, B1/eta, getA(a1), getL(a2));
234 :
235 : //A=G359.54+0.18 Nonthermal Filament
236 : mid = Vector3d(0, -68.24,25.22)*pc;
237 : a1=15.08*pc;
238 : a2=2.72;
239 : B1=1.e-3;
240 1 : b += BPol(pos, mid, B1/eta, getA(a1), getL(a2));
241 :
242 : //A=G359.79 +17 Nonthermal Filament
243 : mid = Vector3d(0,-31.15,23.74)*pc;
244 : a1=16.07*pc;
245 : a2=3.46*pc;
246 : B1=1.e-3;
247 1 : b += BPol(pos, mid, B1/eta, getA(a1), getL(a2));
248 :
249 : //A=G359.96 +0.09 Nonthermal Filament Southern Thread
250 : mid = Vector3d(0, 5.93, 16.32)*pc;
251 : a1=28.68*pc;
252 : a2=1.73*pc;
253 : B1=1.e-4;
254 1 : b += BPol(pos, mid, B1/eta, getA(a1), getL(a2));
255 :
256 : //A=G0.09 +0.17 Nonthermal Filament Northern thread
257 : mid = Vector3d(0, 13.35, 25.22)*pc;
258 : a1=29.42*pc;
259 : a2=2.23*pc;
260 : B1=140.e-6;
261 1 : b += BPol(pos, mid, B1/eta, getA(a1), getL(a2));
262 :
263 : //A=G359.85+0.47 Nonthermal Filament The Pelican is not poloidal but azimuthal
264 : mid = Vector3d(0, -22.25, 69.73)*pc;
265 : a1=11.37*pc;
266 : a2=2.23*pc;
267 : B1=70.e-6 /eta;
268 : // by and bz switched because pelican is differently oriented
269 1 : Vector3d bPelican = BPol(pos, mid, B1/eta, getA(a1), getL(a2));
270 1 : b.x += bPelican.x;
271 1 : b.y += bPelican.z;
272 1 : b.z += bPelican.y;
273 :
274 1 : return b*gauss;
275 : }
276 :
277 1 : Vector3d CMZField::getRadioArcField(const Vector3d& pos) const {//Field in the non-thermal filaments--> predominantly poloidal field
278 : //poloidal field in the non-thermal filament region A=RadioArc
279 : double eta=0.48;
280 : Vector3d mid(0,26.7*pc,10.38*pc);
281 : double a1=70.47*pc;// arcmin-> deg->cm
282 : double a2=9.89*pc;// arcmin-> deg-> cm
283 : double B1=1.e-3;
284 1 : return BPol(pos, mid, B1/eta, getA(a1), getL(a2))*gauss;
285 : }
286 :
287 1 : Vector3d CMZField::getField(const Vector3d& pos) const{
288 : Vector3d b(0.);
289 :
290 1 : if(useMCField){
291 0 : b += getMCField(pos);
292 : }
293 1 : if(useICField){
294 1 : b += getICField(pos);
295 : }
296 1 : if(useNTFField){
297 0 : b += getNTFField(pos);
298 : }
299 1 : if(useRadioArc){
300 0 : b += getRadioArcField(pos);
301 : }
302 :
303 1 : return b;
304 : }
305 :
306 : } // namespace crpropa
|
