Line data Source code
1 : #include "crpropa/magneticField/JF12FieldSolenoidal.h"
2 : #include "crpropa/Units.h"
3 : #include "crpropa/GridTools.h"
4 : #include "crpropa/Random.h"
5 :
6 : namespace crpropa {
7 :
8 0 : JF12FieldSolenoidal::JF12FieldSolenoidal(double delta, double zs) {
9 0 : zS = zs; // set scale heigth for the parabolic X field lines
10 0 : r1 = 5 * kpc; // inner boundary of the disk field
11 0 : r2 = 20 * kpc; // outer boudary of the disk field
12 0 : r1s = r1 + delta; // the magnetic flux of the spirals is redirected for r in [r1,r1s]
13 0 : r2s = r2 - delta; // same here at outer boundary between r2s and r2
14 0 : phi0 = 0.; // somewhat arbitrary choice, has to be chosen in [-pi,pi]
15 :
16 0 : for (int i = 1;i < 9; i++){
17 : // fill the array with angles in [-pi,pi] where the 8 spiral arms intersect the r1 - ring
18 : // indexing starts at 1 to match the indexing in the papers on the JF12 field!
19 0 : phi0Arms[i] = M_PI - cotPitch * log(rArms[i-1] / r1);
20 : }
21 :
22 : // cyclic closure of the array, with next values periodically continued
23 : // outside [-pi,pi] to simplify looping and searching for correct spiral arms
24 0 : phi0Arms[0] = phi0Arms[8] + 2 * M_PI;
25 0 : phi0Arms[9] = phi0Arms[1] - 2 * M_PI;
26 0 : phi0Arms[10] = phi0Arms[2] - 2 *M_PI;
27 :
28 : // determine the position of phi0 in the array, i.e. find the correct spiral arm.
29 : int idx0 = 1; // corresponding index in phi0Arms such that phi0Arms[idx0] < phi0 < phi0Arms[idx0-1]
30 0 : while (phi0 < phi0Arms[idx0]){
31 0 : idx0 += 1; // search clockwise, starting with the check if phi0Arms[1] < phi0 < phi0Arms[0]
32 : }
33 :
34 : // fill the bDisk array with spiral field strengths at r = r1.
35 : // note the indexing starting with 1 here to match the indexing in the JF12 papers!
36 : // for a position (r1,phi), phi in [-pi,pi], the correct field strength is given by
37 : // bDisk[i] if phi0Arms[i] < phi0 < phi0Arms[i-1].
38 0 : bDisk[1] = 0.1 * muG;
39 0 : bDisk[2] = 3.0 * muG;
40 0 : bDisk[3] = -0.9 * muG;
41 0 : bDisk[4] = -0.8 * muG;
42 0 : bDisk[5] = -2.0 * muG;
43 0 : bDisk[6] = -4.2 * muG;
44 0 : bDisk[7] = 0.0 * muG;
45 :
46 : // re-compute b_8 for actual (net flux = 0)-correction of the spiral field with minimal round-off errors
47 : double flux1to7 = 0.;
48 0 : for (int i = 1; i < 8; i++){
49 0 : flux1to7 += (phi0Arms[i-1] - phi0Arms[i]) * bDisk[i];
50 : }
51 0 : bDisk[8] = -flux1to7 / (phi0Arms[7] - phi0Arms[8]);
52 :
53 0 : bDisk[0] = bDisk[8]; // again close the array periodically
54 0 : bDisk[9] = bDisk[1];
55 0 : bDisk[10] = bDisk[2];
56 :
57 : // set coefficients for the evaluation of the phi-integral over the piecewise constant field strengths at r=r1
58 : // such that it may be evaluated as H(phi) = phiCoeff[j] + bDisk[j] * phi later on
59 : // start integration at phi0Arms[0] first, shift to lower integration boundary phi0 later
60 0 : phiCoeff[0] = 0;
61 0 : for (int i = 1; i < 10; i++){
62 0 : phiCoeff[i] = phiCoeff[i-1] + (bDisk[i-1] - bDisk[i]) * phi0Arms[i-1];
63 : }
64 :
65 : // correct for H(phi0) = 0
66 0 : corr = phiCoeff[idx0] + bDisk[idx0] * phi0;
67 0 : for (int i = 1; i < 10; i++){
68 0 : phiCoeff[i] = phiCoeff[i] - corr;
69 : }
70 0 : }
71 :
72 0 : void JF12FieldSolenoidal::setDiskTransitionWidth(double delta) {
73 0 : r1s = r1 + delta;
74 0 : r2s = r2 - delta;
75 0 : }
76 :
77 0 : void JF12FieldSolenoidal::setXScaleHeight(double zs) {
78 0 : zS = zs;
79 0 : }
80 :
81 0 : double JF12FieldSolenoidal::getDiskTransitionWidth() const {
82 0 : return (r1s - r1);
83 : }
84 :
85 0 : double JF12FieldSolenoidal::getXScaleHeight() const {
86 0 : return zS;
87 : }
88 :
89 0 : void JF12FieldSolenoidal::deactivateOuterTransition() {
90 0 : r2s = r2;
91 0 : }
92 :
93 0 : void JF12FieldSolenoidal::setUseStriatedField(bool use) {
94 0 : if ((use) and (striatedGrid)) {
95 0 : KISS_LOG_WARNING << "JF12FieldSolenoidal: No striated field set: ignored.";
96 0 : return;
97 : }
98 0 : useStriatedField = use;
99 : }
100 :
101 0 : void JF12FieldSolenoidal::setUseTurbulentField(bool use) {
102 0 : if ((use) and (turbulentGrid)) {
103 0 : KISS_LOG_WARNING << "JF12FieldSolenoidal: No turbulent field set: ignored.";
104 0 : return;
105 : }
106 0 : useTurbulentField = use;
107 : }
108 :
109 0 : Vector3d JF12FieldSolenoidal::getDiskField(const double& r, const double& z, const double& phi, const double& sinPhi, const double& cosPhi) const {
110 : Vector3d b(0.);
111 :
112 0 : if (useDiskField){
113 0 : double lfDisk = logisticFunction(z, hDisk, wDisk); // for vertical scaling as in initial JF12
114 :
115 0 : double hint = getHPhiIntegral(r, phi); // phi integral to restore solenoidality in transition region, only enters if r is in [r1,r1s] or [r2s,r2]
116 0 : double mag1 = getSpiralFieldStrengthConstant(r, phi); // returns bDisk[j] for the current spiral arm
117 :
118 0 : if ((r1 < r) && (r < r2)) {
119 0 : double pdelta = getDiskTransitionPolynomial(r);
120 0 : double qdelta = getDiskTransitionPolynomialDerivative(r);
121 0 : double br = pdelta * mag1 * sinPitch;
122 0 : double bphi = pdelta * mag1 * cosPitch - qdelta * hint * sinPitch;
123 :
124 0 : b.x += br * cosPhi - bphi * sinPhi;
125 0 : b.y += br * sinPhi + bphi * cosPhi;
126 :
127 0 : b *= (1 - lfDisk);
128 : }
129 : }
130 0 : return b;
131 : }
132 :
133 0 : Vector3d JF12FieldSolenoidal::getXField(const double& r, const double& z, const double& sinPhi, const double& cosPhi) const {
134 : Vector3d b(0.);
135 :
136 0 : if (useXField){
137 : double bMagX;
138 : double sinThetaX, cosThetaX;
139 : double rp; // radius where current intial field line passes z = 0
140 0 : double rc = rXc + fabs(z) / tanThetaX0;
141 0 : double r0c = rXc + zS / tanThetaX0; // radius where field line through rXc passes z = zS
142 : double f, r0, br0, bz0;
143 : bool inner = true; // distinguish between inner and outer region
144 :
145 : // return intial field if z>=zS
146 0 : if (fabs(z) > zS){
147 0 : if ((r == 0.)){
148 0 : b.z = bX / ((1. + fabs(z) * cotThetaX0 / rXc) * (1. + fabs(z) * cotThetaX0 / rXc));
149 0 : return b;
150 : }
151 :
152 0 : if (r < rc) {
153 : // inner varying elevation region
154 0 : rp = r * rXc / rc;
155 0 : bMagX = bX * exp(-1 * rp / rX) * (rXc / rc) * (rXc / rc);
156 :
157 0 : double thetaX = atan(fabs(z) / (r - rp));
158 :
159 0 : if (z == 0)
160 : thetaX = M_PI / 2.;
161 :
162 0 : sinThetaX = sin(thetaX);
163 0 : cosThetaX = cos(thetaX);
164 : }
165 : else {
166 : // outer constant elevation region
167 0 : rp = r - fabs(z) / tanThetaX0;
168 0 : bMagX = bX * exp(-rp / rX) * (rp / r);
169 :
170 0 : sinThetaX = sinThetaX0;
171 0 : cosThetaX = cosThetaX0;
172 : }
173 0 : double zsign = z < 0 ? -1 : 1;
174 0 : b.x += zsign * bMagX * cosThetaX * cosPhi;
175 0 : b.y += zsign * bMagX * cosThetaX * sinPhi;
176 0 : b.z += bMagX * sinThetaX;
177 : }
178 : // parabolic field lines for z<zS
179 : else {
180 : // determine r at which parabolic field line through (r,z) passes z = zS
181 0 : r0 = r * 1. / (1.- 1./ (2. * (zS + rXc * tanThetaX0)) * (zS - z * z / zS));
182 :
183 : // determine correct region (inner/outer)
184 : // and compute factor F for solenoidality
185 0 : if (r0 >= r0c){
186 0 : r0 = r + 1. / (2. * tanThetaX0) * (zS - z * z / zS);
187 0 : f = 1. + 1/ (2 * r * tanThetaX0/ zS) * (1. - (z / zS) * (z / zS));
188 : }
189 : else
190 : {
191 0 : f = 1. / ((1. - 1./( 2. + 2. * (rXc * tanThetaX0/ zS)) * (1. - (z / zS) * (z / zS))) * (1. - 1./( 2. + 2. * (rXc * tanThetaX0/ zS)) * (1. - (z / zS) * (z / zS))));
192 : }
193 :
194 : // field strength at that position
195 0 : if (r0 < r0c){
196 0 : rp = r0 * rXc / r0c;
197 0 : double thetaX = atan(zS / (r0 - rp));
198 :
199 : // field strength at (r0,zS) for inner region
200 0 : br0 = bX * exp(- rp / rX) * (rXc/ r0c) * (rXc/ r0c) * cos(thetaX);
201 0 : bz0 = bX * exp(- rp / rX) * (rXc/ r0c) * (rXc/ r0c) * sin(thetaX);
202 : }
203 : else {
204 : // field strength at (r0,zS) for outer region
205 0 : rp = r0 - zS / tanThetaX0;
206 0 : br0 = bX * exp(- rp / rX) * (rp/r0) * cosThetaX0;
207 0 : bz0 = bX * exp(- rp / rX) * (rp/r0) * sinThetaX0;
208 : }
209 :
210 0 : double br = z / zS * f * br0;
211 0 : double bz = bz0 * f;
212 :
213 0 : b.x += br * cosPhi;
214 0 : b.y += br * sinPhi;
215 0 : b.z += bz;
216 : }
217 : }
218 : return b;
219 : }
220 :
221 0 : double JF12FieldSolenoidal::getDiskTransitionPolynomial(const double& r) const {
222 : // 0 disk field outside
223 0 : if ((r < r1) || (r > r2)) {
224 : return 0.;
225 : }
226 : // unchanged field
227 0 : if ((r > r1s) && (r < r2s)) {
228 0 : return r1/r;
229 : }
230 : // transitions region parameters
231 : double r_a = r1;
232 : double r_b = r1s;
233 :
234 0 : if (r >= r2s) {
235 : r_a = r2;
236 : r_b = r2s;
237 : }
238 : // differentiable transition at r_s, continous at r_a
239 0 : double fakt = (r_a / r_b - 2.) / ((r_a - r_b) * (r_a - r_b));
240 0 : return (r1/r_b) * (2. - r / r_b + fakt * (r-r_b) * (r-r_b));
241 : }
242 :
243 0 : double JF12FieldSolenoidal::getDiskTransitionPolynomialDerivative(const double& r) const {
244 : // 0 disk field outside
245 0 : if ((r < r1) || (r > r2)) {
246 : return 0.;
247 : }
248 : // unchanged field
249 0 : if ((r > r1s) && (r < r2s)) {
250 : return 0.;
251 : }
252 : // transitions region parameters
253 : double r_a = r1;
254 : double r_b = r1s;
255 :
256 0 : if (r >= r2s) {
257 : r_a = r2;
258 : r_b = r2s;
259 : }
260 : // differentiable transition polynomial at r_s, continous at r_a
261 0 : double fakt = (r_a / r_b - 2.) / ((r_a - r_b) * (r_a - r_b));
262 0 : return (r1/r_b) * (2. - 2. * r/r_b + fakt * (3. * r * r - 4. * r * r_b + r_b * r_b));
263 : }
264 :
265 0 : double JF12FieldSolenoidal::getHPhiIntegral(const double& r, const double& phi) const {
266 : // Evaluates the H(phi1) integral for solenoidality for the position (r,phi) which is mapped back to (r1=5kpc,phi1)
267 : // along the spiral field line.
268 : double H_ret = 0.;
269 :
270 0 : if ((r1 < r) && (r < r2)){
271 : // find index of the correct spiral arm for (r1,phi1) just like in getSpiralFieldStrengthConstant
272 : int idx = 1;
273 0 : double phi1 = phi - log(r/r1) * cotPitch;
274 0 : phi1 = atan2(sin(phi1), cos(phi1));
275 0 : while (phi1 < phi0Arms[idx]){
276 0 : idx += 1;
277 : }
278 0 : H_ret = phi1 * bDisk[idx] + phiCoeff[idx];
279 : }
280 0 : return H_ret;
281 : }
282 :
283 0 : double JF12FieldSolenoidal::getSpiralFieldStrengthConstant(const double& r, const double& phi) const {
284 : // For a given position (r, phi) in polar coordinates, this method returns the field strength
285 : // of the spiral field at r1 = 5 kpc for the magnetic spiral arm where (r, phi) is located.
286 : // The method first computes the angle phi1 at which the spiral field line passing through (r, phi) intersects
287 : // the circle with radius r1 = 5 kpc. Afterwards, the correct spiral arm is found by searching the index idx
288 : // such that phi0Arms[idx] < phi1 < phi0Arms[idx-1]. The correct field strength of the respective spiral arm
289 : // where (r, phi) is located is then given as bDisk[idx].
290 : double b_ret = 0.;
291 : int idx = 1;
292 0 : if ((r1 < r) && (r < r2)){
293 0 : double phi1 = phi - log(r/r1) * cotPitch; // map the position (r, phi) to (5 kpc, phi1) along the logarithmic spiral field line
294 0 : phi1 = atan2(sin(phi1), cos(phi1)); // map this angle to [-pi,+pi]
295 0 : while (phi1 < phi0Arms[idx]){
296 0 : idx += 1; // run clockwise through the spiral arms; the cyclic closure of phi0Arms[9] = phi0Arms[1] - 2 pi is needed if -pi <= phi1 <= phi0Arms[8].
297 : }
298 0 : b_ret = bDisk[idx];
299 : }
300 0 : return b_ret;
301 : }
302 : } // namespace crpropa
|
