ASCOT5
Loading...
Searching...
No Matches
gctransform.c
Go to the documentation of this file.
1
32#include "offload_acc_omp.h"
33#include <stdlib.h>
34#include <stdio.h>
35#include <math.h>
36#include "ascot5.h"
37#include "math.h"
38#include "consts.h"
39#include "physlib.h"
40#include "gctransform.h"
41
42
44//static int GCTRANSFORM_ORDER = 1;
45#ifdef _OPENACC
46static int GCTRANSFORM_ORDER = 1;
47#pragma acc declare copyin(GCTRANSFORM_ORDER)
48#elif defined(_OPENMP)
49DECLARE_TARGET
50static int GCTRANSFORM_ORDER = 1;
51DECLARE_TARGET_END
52#else
53static int GCTRANSFORM_ORDER = 1;
54#endif
55
56
62void gctransform_setorder(int order) {
63 GCTRANSFORM_ORDER = order;
64}
65
90void gctransform_particle2guidingcenter(
91 real mass, real charge, real* B_dB,
92 real r, real phi, real z, real pr, real pphi, real pz,
93 real* R, real* Phi, real* Z, real* ppar, real* mu, real* zeta) {
94
95 /* |B| */
96 real Bnorm = sqrt(B_dB[0]*B_dB[0] + B_dB[4]*B_dB[4] + B_dB[8]*B_dB[8]);
97
98 /* Guiding center transformation is more easily done in cartesian
99 * coordinates so we switch to using those */
100 real rpz[3] = {r, phi, z};
101 real prpz[3] = {pr, pphi, pz};
102 real xyz[3];
103 real pxyz[3];
104 real B_dBxyz[12];
105 math_rpz2xyz(rpz, xyz);
106 math_vec_rpz2xyz(prpz, pxyz, phi);
107 math_jac_rpz2xyz(B_dB, B_dBxyz, r, phi);
108
109 /* bhat = Unit vector of B */
110 real bhat[3] = { B_dBxyz[0]/Bnorm, B_dBxyz[4]/Bnorm, B_dBxyz[8]/Bnorm };
111
112 /* Magnetic field norm gradient */
113 real gradB[3];
114 gradB[0] = (bhat[0]*B_dBxyz[1] + bhat[1]*B_dBxyz[5] + bhat[2]*B_dBxyz[9]);
115 gradB[1] = (bhat[0]*B_dBxyz[2] + bhat[1]*B_dBxyz[6] + bhat[2]*B_dBxyz[10]);
116 gradB[2] = (bhat[0]*B_dBxyz[3] + bhat[1]*B_dBxyz[7] + bhat[2]*B_dBxyz[11]);
117
118 /* nabla x |B| */
119 real curlB[3] = {B_dBxyz[10]-B_dBxyz[7],
120 B_dBxyz[3]-B_dBxyz[9],
121 B_dBxyz[5]-B_dBxyz[2]};
122
123 /* Magnetic field torsion = bhat dot ( nabla X bhat ) which is equivalent to
124 bhat dot ( nabla x |B| ) / |B|
125 because nabla x |B| = |B| nabla x bhat + (nabla |B|) x bhat */
126 real tau = math_dot(bhat, curlB)/Bnorm;
127
128 /* Gradient of magnetic field unit vector */
129 real nablabhat[9];
130 nablabhat[0] = (B_dBxyz[1] - gradB[0] * bhat[0]) / Bnorm;
131 nablabhat[1] = (B_dBxyz[2] - gradB[0] * bhat[1]) / Bnorm;
132 nablabhat[2] = (B_dBxyz[3] - gradB[0] * bhat[2]) / Bnorm;
133 nablabhat[3] = (B_dBxyz[5] - gradB[1] * bhat[0]) / Bnorm;
134 nablabhat[4] = (B_dBxyz[6] - gradB[1] * bhat[1]) / Bnorm;
135 nablabhat[5] = (B_dBxyz[7] - gradB[1] * bhat[2]) / Bnorm;
136 nablabhat[6] = (B_dBxyz[9] - gradB[2] * bhat[0]) / Bnorm;
137 nablabhat[7] = (B_dBxyz[10] - gradB[2] * bhat[1]) / Bnorm;
138 nablabhat[8] = (B_dBxyz[11] - gradB[2] * bhat[2]) / Bnorm;
139
140 /* Magnetic field curvature vector = bhat dot nablabhat.
141 Note that nabla x bhat = tau bhat + bhat x kappa */
142 real kappa[3];
143 kappa[0] =
144 ( bhat[0] * B_dBxyz[1] +
145 bhat[1] * B_dBxyz[5] +
146 bhat[2] * B_dBxyz[9] - math_dot(bhat, gradB) * bhat[0] ) / Bnorm;
147 kappa[1] =
148 ( bhat[0] * B_dBxyz[2] +
149 bhat[1] * B_dBxyz[6] +
150 bhat[2] * B_dBxyz[10] - math_dot(bhat, gradB) * bhat[1] ) / Bnorm;
151 kappa[2] =
152 ( bhat[0] * B_dBxyz[3] +
153 bhat[1] * B_dBxyz[7] +
154 bhat[2] * B_dBxyz[11] - math_dot(bhat, gradB) * bhat[2] ) / Bnorm;
155
156 /* Zeroth order mu and ppar */
157 real ppar0 =
158 pxyz[0] * bhat[0] + pxyz[1] * bhat[1] + pxyz[2] * bhat[2];
159
160 /* Gyrovector rhohat = bhat X perphat */
161 real pperp[3] = {pxyz[0] - ppar0 * bhat[0],
162 pxyz[1] - ppar0 * bhat[1],
163 pxyz[2] - ppar0 * bhat[2]};
164 real mu0 = math_dot(pperp, pperp) / ( 2 * mass * Bnorm );
165 real perphat[3];
166 math_unit(pperp, perphat);
167 if( charge < 0 ) {
168 perphat[0] = - perphat[0];
169 perphat[1] = - perphat[1];
170 perphat[2] = - perphat[2];
171 }
172 real rhohat[3];
173 math_cross(bhat, perphat, rhohat);
174
175 /* Double product of dyadic a1 = -(1/2) * (rhohat perphat + perphat rhohat)
176 and gradb (gradient of magnetic field unit vector) */
177 real a1ddotgradb =
178 -0.5 * ( 2 * ( rhohat[0] * perphat[0] * nablabhat[0] +
179 rhohat[1] * perphat[1] * nablabhat[4] +
180 rhohat[2] * perphat[2] * nablabhat[8] ) +
181 ( rhohat[0] * perphat[1] + rhohat[1] * perphat[0] ) *
182 ( nablabhat[1] + nablabhat[3] ) +
183 ( rhohat[0] * perphat[2] + rhohat[2] * perphat[0] ) *
184 ( nablabhat[2] + nablabhat[6] ) +
185 ( rhohat[1] * perphat[2] + rhohat[2] * perphat[1] ) *
186 ( nablabhat[5] + nablabhat[7] ) );
187
188 /* Double product of dyadic a2 = (1/4) * (perphat perphat - rhohat rhohat)
189 and gradb (gradient of magnetic field unit vector) */
190 real a2ddotgradb =
191 0.25 * (
192 perphat[0] * perphat[0] * nablabhat[0] +
193 perphat[1] * perphat[0] * nablabhat[3] +
194 perphat[2] * perphat[0] * nablabhat[6] +
195 perphat[0] * perphat[1] * nablabhat[1] +
196 perphat[1] * perphat[1] * nablabhat[4] +
197 perphat[2] * perphat[1] * nablabhat[7] +
198 perphat[0] * perphat[2] * nablabhat[2] +
199 perphat[1] * perphat[2] * nablabhat[5] +
200 perphat[2] * perphat[2] * nablabhat[8]
201 )
202 -0.25 * (
203 rhohat[0] * rhohat[0] * nablabhat[0] +
204 rhohat[0] * rhohat[1] * nablabhat[3] +
205 rhohat[0] * rhohat[2] * nablabhat[6] +
206 rhohat[1] * rhohat[0] * nablabhat[1] +
207 rhohat[1] * rhohat[1] * nablabhat[4] +
208 rhohat[1] * rhohat[2] * nablabhat[7] +
209 rhohat[2] * rhohat[0] * nablabhat[2] +
210 rhohat[2] * rhohat[1] * nablabhat[5] +
211 rhohat[2] * rhohat[2] * nablabhat[8]
212 );
213
214 /* Choose a fixed basis so that e2 = bhat x e1. We choose e1 = bhat x z */
215 real e1[3] = {0, 0, 1};
216 real e2[3];
217 math_cross(bhat, e1, e2);
218 math_unit(e2, e1);
219 math_cross(bhat, e1, e2);
220
221 /* Gyrolength */
222 real rho0 = 0.0;
223 if(charge != 0) {
224 rho0 = sqrt( ( 2* mass * mu0 ) / Bnorm ) / fabs(charge);
225 }
226
227 /* First order position */
228 real XYZ[3];
229 XYZ[0] = xyz[0] - rho0 * rhohat[0];
230 XYZ[1] = xyz[1] - rho0 * rhohat[1];
231 XYZ[2] = xyz[2] - rho0 * rhohat[2];
232
233 real RPZ[3];
234 math_xyz2rpz(XYZ, RPZ);
235 *R = RPZ[0];
236 *Phi = RPZ[1];
237 *Z = RPZ[2];
238
239 /* Zeroth order gyroangle is directly defined from basis {e1, e2} and
240 gyrovector as tan(zeta) = -rhohat dot e2 / rhohat dot e1 */
241 real zeta0 =
242 atan2( -math_dot(rhohat, e2), math_dot(rhohat, e1) );
243
244 /* First order momentum terms ppar, mu1, and zeta1 */
245 real ppar1 = 0.0;
246 if(charge != 0) {
247 ppar1 = -ppar0 * rho0 * math_dot(rhohat, kappa) +
248 ( mass * mu0 / charge ) * ( tau + a1ddotgradb );
249 }
250
251 real ppar2 = ppar0 * ppar0; // Power, not second order component ;)
252 real temp[3] = { mu0 * gradB[0] + ppar2 * kappa[0] / mass,
253 mu0 * gradB[1] + ppar2 * kappa[1] / mass,
254 mu0 * gradB[2] + ppar2 * kappa[2] / mass };
255 real mu1 = 0.0;
256 if(charge != 0) {
257 mu1 = ( rho0 / Bnorm ) * math_dot(rhohat, temp) -
258 ( ppar0 * mu0 / ( charge * Bnorm ) ) *
259 ( tau + a1ddotgradb );
260 }
261
262 /* This monster is Littlejohn's gyro gauge vector (nabla e1) dot e2 */
263 real bx2by2 = bhat[0] * bhat[0] + bhat[1] * bhat[1];
264 real b1x = -bhat[0]*bhat[2]/bx2by2;
265 real b2x = -bhat[1]/(bx2by2*bx2by2);
266 real b1y = -bhat[1]*bhat[2]/bx2by2;
267 real b2y = -bhat[0]/(bx2by2*bx2by2);
268 real Rvec[3];
269 Rvec[0] =
270 b1x * ( nablabhat[1] + b2x * ( nablabhat[0] + nablabhat[1] ) ) +
271 b1y * ( nablabhat[0] + b2y * ( nablabhat[0] + nablabhat[1] ) );
272 Rvec[1] =
273 b1x * ( nablabhat[4] + b2x * ( nablabhat[3] + nablabhat[4] ) ) +
274 b1y * ( nablabhat[3] + b2y * ( nablabhat[3] + nablabhat[4] ) );
275 Rvec[2] =
276 b1x * ( nablabhat[7] + b2x * ( nablabhat[6] + nablabhat[7] ) ) +
277 b1y * ( nablabhat[6] + b2y * ( nablabhat[6] + nablabhat[7] ) );
278
279 temp[0] = gradB[0] + kappa[0] * ( ppar2 / mass ) / (2 * mu0);
280 temp[1] = gradB[1] + kappa[1] * ( ppar2 / mass ) / (2 * mu0);
281 temp[2] = gradB[2] + kappa[2] * ( ppar2 / mass ) / (2 * mu0);
282 real zeta1 = 0.0;
283 if(charge != 0) {
284 zeta1 = -rho0 * math_dot(rhohat, Rvec) +
285 ( ppar0 / ( charge * Bnorm ) ) * a2ddotgradb +
286 ( rho0 / Bnorm ) * math_dot(perphat, temp);
287 }
288
289 /* Choose whether to use first order transformation in momentum space */
290 if(GCTRANSFORM_ORDER) {
291 *ppar = ppar0 + ppar1;
292 *mu = fabs(mu0 + mu1);
293 *zeta = zeta0 + zeta1;
294 }
295 else {
296 *ppar = ppar0;
297 *mu = mu0;
298 *zeta = zeta0;
299 }
300
301 /* zeta is defined to be in interval [0, 2pi] */
302 *zeta = fmod(CONST_2PI + (*zeta), CONST_2PI);
303}
304
329void gctransform_guidingcenter2particle(
330 real mass, real charge, real* B_dB,
331 real R, real Phi, real Z, real ppar, real mu, real zeta,
332 real* r, real* phi, real* z, real* pparprt, real* muprt, real* zetaprt) {
333
334 /* |B| */
335 real Bnorm = sqrt(B_dB[0]*B_dB[0] + B_dB[4]*B_dB[4] + B_dB[8]*B_dB[8]);
336
337 /* Guiding center transformation is more easily done in cartesian
338 * coordinates so we switch to using those */
339 real RPZ[3] = {R, Phi, Z};
340 real XYZ[3];
341 real B_dBxyz[12];
342 math_rpz2xyz(RPZ, XYZ);
343 math_jac_rpz2xyz(B_dB, B_dBxyz, R, Phi);
344
345 /* bhat = Unit vector of B */
346 real bhat[3] = { B_dBxyz[0]/Bnorm, B_dBxyz[4]/Bnorm, B_dBxyz[8]/Bnorm };
347
348 /* Magnetic field norm gradient */
349 real gradB[3];
350 gradB[0] = (bhat[0]*B_dBxyz[1] + bhat[1]*B_dBxyz[5] + bhat[2]*B_dBxyz[9]);
351 gradB[1] = (bhat[0]*B_dBxyz[2] + bhat[1]*B_dBxyz[6] + bhat[2]*B_dBxyz[10]);
352 gradB[2] = (bhat[0]*B_dBxyz[3] + bhat[1]*B_dBxyz[7] + bhat[2]*B_dBxyz[11]);
353
354 /* nabla x |B| */
355 real curlB[3] = {B_dBxyz[10]-B_dBxyz[7],
356 B_dBxyz[3]-B_dBxyz[9],
357 B_dBxyz[5]-B_dBxyz[2]};
358
359 /* Magnetic field torsion = bhat dot ( nabla X bhat ) which is equivalent to
360 bhat dot ( nabla x |B| ) / |B|
361 because nabla x |B| = |B| nabla x bhat + (nabla |B|) x bhat */
362 real tau = math_dot(bhat, curlB)/Bnorm;
363
364 /* Gradient of magnetic field unit vector */
365 real nablabhat[9];
366 nablabhat[0] = (B_dBxyz[1] - gradB[0] * bhat[0]) / Bnorm;
367 nablabhat[1] = (B_dBxyz[2] - gradB[0] * bhat[1]) / Bnorm;
368 nablabhat[2] = (B_dBxyz[3] - gradB[0] * bhat[2]) / Bnorm;
369 nablabhat[3] = (B_dBxyz[5] - gradB[1] * bhat[0]) / Bnorm;
370 nablabhat[4] = (B_dBxyz[6] - gradB[1] * bhat[1]) / Bnorm;
371 nablabhat[5] = (B_dBxyz[7] - gradB[1] * bhat[2]) / Bnorm;
372 nablabhat[6] = (B_dBxyz[9] - gradB[2] * bhat[0]) / Bnorm;
373 nablabhat[7] = (B_dBxyz[10] - gradB[2] * bhat[1]) / Bnorm;
374 nablabhat[8] = (B_dBxyz[11] - gradB[2] * bhat[2]) / Bnorm;
375
376 /* Magnetic field curvature vector = bhat dot nablabhat.
377 Note that nabla x bhat = tau bhat + bhat x kappa */
378 real kappa[3];
379 kappa[0] =
380 ( bhat[0] * B_dBxyz[1] +
381 bhat[1] * B_dBxyz[5] +
382 bhat[2] * B_dBxyz[9] - math_dot(bhat, gradB) * bhat[0] ) / Bnorm;
383 kappa[1] =
384 ( bhat[0] * B_dBxyz[2] +
385 bhat[1] * B_dBxyz[6] +
386 bhat[2] * B_dBxyz[10] - math_dot(bhat, gradB) * bhat[1] ) / Bnorm;
387 kappa[2] =
388 ( bhat[0] * B_dBxyz[3] +
389 bhat[1] * B_dBxyz[7] +
390 bhat[2] * B_dBxyz[11] - math_dot(bhat, gradB) * bhat[2] ) / Bnorm;
391
392 /* Choose a fixed basis so that e2 = bhat x e1. We choose e1 = bhat x z */
393 real e1[3] = {0, 0, 1};
394 real e2[3];
395 math_cross(bhat, e1, e2);
396 math_unit(e2, e1);
397 math_cross(bhat, e1, e2);
398
399 /* Gyrolength */
400 real rho0 = 0.0;
401 if(charge != 0) {
402 rho0 = sqrt( ( 2* mass * mu ) / Bnorm ) / fabs(charge);
403 }
404
405 /* Gyrovector rhohat and pperphat*/
406 real rhohat[3];
407 real perphat[3];
408 real c = cos(zeta);
409 real s = sin(zeta);
410
411 rhohat[0] = c * e1[0] - s * e2[0];
412 rhohat[1] = c * e1[1] - s * e2[1];
413 rhohat[2] = c * e1[2] - s * e2[2];
414
415 perphat[0] = -s * e1[0] - c * e2[0];
416 perphat[1] = -s * e1[1] - c * e2[1];
417 perphat[2] = -s * e1[2] - c * e2[2];
418
419 /* Double product of dyadic a1 = -(1/2) * (rhohat perphat + perphat rhohat)
420 and gradb (gradient of magnetic field unit vector) */
421 real a1ddotgradb =
422 -0.5 * ( 2 * ( rhohat[0] * perphat[0] * nablabhat[0] +
423 rhohat[1] * perphat[1] * nablabhat[4] +
424 rhohat[2] * perphat[2] * nablabhat[8] ) +
425 ( rhohat[0] * perphat[1] + rhohat[1] * perphat[0] ) *
426 ( nablabhat[1] + nablabhat[3] ) +
427 ( rhohat[0] * perphat[2] + rhohat[2] * perphat[0] ) *
428 ( nablabhat[2] + nablabhat[6] ) +
429 ( rhohat[1] * perphat[2] + rhohat[2] * perphat[1] ) *
430 ( nablabhat[5] + nablabhat[7] ) );
431
432 /* Double product of dyadic a2 = (1/4) * (perphat perphat - rhohat rhohat)
433 and gradb (gradient of magnetic field unit vector) */
434 real a2ddotgradb =
435 0.25 * (
436 perphat[0] * perphat[0] * nablabhat[0] +
437 perphat[1] * perphat[0] * nablabhat[3] +
438 perphat[2] * perphat[0] * nablabhat[6] +
439 perphat[0] * perphat[1] * nablabhat[1] +
440 perphat[1] * perphat[1] * nablabhat[4] +
441 perphat[2] * perphat[1] * nablabhat[7] +
442 perphat[0] * perphat[2] * nablabhat[2] +
443 perphat[1] * perphat[2] * nablabhat[5] +
444 perphat[2] * perphat[2] * nablabhat[8]
445 )
446 -0.25 * (
447 rhohat[0] * rhohat[0] * nablabhat[0] +
448 rhohat[0] * rhohat[1] * nablabhat[3] +
449 rhohat[0] * rhohat[2] * nablabhat[6] +
450 rhohat[1] * rhohat[0] * nablabhat[1] +
451 rhohat[1] * rhohat[1] * nablabhat[4] +
452 rhohat[1] * rhohat[2] * nablabhat[7] +
453 rhohat[2] * rhohat[0] * nablabhat[2] +
454 rhohat[2] * rhohat[1] * nablabhat[5] +
455 rhohat[2] * rhohat[2] * nablabhat[8]
456 );
457
458 /* First order terms */
459 real ppar1 = 0.0;
460 if(charge != 0) {
461 ppar1 = -ppar * rho0 * math_dot(rhohat, kappa) +
462 ( mass * mu / charge ) * ( tau + a1ddotgradb );
463 }
464
465 real ppar2 = ppar * ppar; // Second power, not second order term
466 real temp[3] =
467 { mu * gradB[0] + ppar2 * kappa[0] / mass,
468 mu * gradB[1] + ppar2 * kappa[1] / mass,
469 mu * gradB[2] + ppar2 * kappa[2] / mass };
470 real mu1 = 0.0;
471 if(charge != 0) {
472 mu1 = ( rho0 / Bnorm ) * math_dot(rhohat, temp) -
473 ( ppar * mu / ( charge * Bnorm ) ) *
474 ( tau + a1ddotgradb );
475 }
476
477 temp[0] = gradB[0] + kappa[0] * ppar2 / (2 * mass * mu);
478 temp[1] = gradB[1] + kappa[1] * ppar2 / (2 * mass * mu);
479 temp[2] = gradB[2] + kappa[2] * ppar2 / (2 * mass * mu);
480
481 /* This monster is Littlejohn's gyro gauge vector (nabla e1) dot e2 */
482 real bx2by2 = bhat[0] * bhat[0] + bhat[1] * bhat[1];
483 real b1x = -bhat[0]*bhat[2]/bx2by2;
484 real b2x = -bhat[1]/(bx2by2*bx2by2);
485 real b1y = -bhat[1]*bhat[2]/bx2by2;
486 real b2y = -bhat[0]/(bx2by2*bx2by2);
487 real Rvec[3];
488 Rvec[0] =
489 b1x * ( nablabhat[1] + b2x * ( nablabhat[0] + nablabhat[1] ) ) +
490 b1y * ( nablabhat[0] + b2y * ( nablabhat[0] + nablabhat[1] ) );
491 Rvec[1] =
492 b1x * ( nablabhat[4] + b2x * ( nablabhat[3] + nablabhat[4] ) ) +
493 b1y * ( nablabhat[3] + b2y * ( nablabhat[3] + nablabhat[4] ) );
494 Rvec[2] =
495 b1x * ( nablabhat[7] + b2x * ( nablabhat[6] + nablabhat[7] ) ) +
496 b1y * ( nablabhat[6] + b2y * ( nablabhat[6] + nablabhat[7] ) );
497 real zeta1 = 0.0;
498 if(charge != 0) {
499 zeta1 = -rho0 * math_dot(rhohat, Rvec) +
500 ( ppar / ( charge * Bnorm ) ) * a2ddotgradb +
501 ( rho0 / Bnorm ) * math_dot(perphat, temp);
502 }
503
504 /* Choose whether to use first or zeroth order momentum transform */
505 if(GCTRANSFORM_ORDER) {
506 mu = fabs(mu - mu1);
507 ppar -= ppar1;
508 zeta -= zeta1;
509 mu = fabs(mu);
510
511 /* Calculate new unit vector for position */
512 c = cos(zeta);
513 s = sin(zeta);
514
515 rhohat[0] = c * e1[0] - s * e2[0];
516 rhohat[1] = c * e1[1] - s * e2[1];
517 rhohat[2] = c * e1[2] - s * e2[2];
518 }
519
520 /* First order position */
521 real xyz[3];
522 xyz[0] = XYZ[0] + rho0 * rhohat[0];
523 xyz[1] = XYZ[1] + rho0 * rhohat[1];
524 xyz[2] = XYZ[2] + rho0 * rhohat[2];
525
526 real rpz[3];
527 math_xyz2rpz(xyz, rpz);
528
529 *r = rpz[0];
530 *phi = rpz[1];
531 *z = rpz[2];
532 *pparprt = ppar;
533 *muprt = mu;
534 *zetaprt = zeta;
535}
536
554void gctransform_pparmuzeta2prpphipz(real mass, real charge, real* B_dB,
555 real phi, real ppar, real mu, real zeta,
556 real* pr, real* pphi, real* pz) {
557 /* Find magnetic field norm and unit vector */
558 real Brpz[3] = {B_dB[0], B_dB[4], B_dB[8]};
559 real Bxyz[3];
560 math_vec_rpz2xyz(Brpz, Bxyz, phi);
561
562 real bhat[3];
563 math_unit(Bxyz, bhat);
564 real Bnorm = math_norm(Bxyz);
565
566 /* Find the basis vectors e1 and e2 */
567 real e1[3] = {0, 0, 1};
568 real e2[3];
569 math_cross(bhat, e1, e2);
570 math_unit(e2, e1);
571 math_cross(bhat, e1, e2);
572
573 /* Perpendicular basis vector */
574 real c = cos(zeta);
575 real s = sin(zeta);
576 real perphat[3];
577 perphat[0] = -s * e1[0] - c * e2[0];
578 perphat[1] = -s * e1[1] - c * e2[1];
579 perphat[2] = -s * e1[2] - c * e2[2];
580
581 /* Perpendicular momentum, negative particles travel opposite to perphat */
582 real pperp = sqrt(2.0 * mass * Bnorm * mu );
583 if( charge < 0 ) {
584 pperp = -pperp;
585 }
586
587 /* Evaluate the momentum vector from ppar and pperp */
588 real pxyz[3];
589 pxyz[0] = ppar * bhat[0] + pperp * perphat[0];
590 pxyz[1] = ppar * bhat[1] + pperp * perphat[1];
591 pxyz[2] = ppar * bhat[2] + pperp * perphat[2];
592
593 /* Back to cylindrical coordinates */
594 real prpz[3];
595 math_vec_xyz2rpz(pxyz, prpz, phi);
596
597 *pr = prpz[0];
598 *pphi = prpz[1];
599 *pz = prpz[2];
600}
ascot5.h
Main header file for ASCOT5.
real
double real
Definition ascot5.h:85
consts.h
Header file containing physical and mathematical constants.
CONST_2PI
#define CONST_2PI
2*pi
Definition consts.h:14
gctransform_guidingcenter2particle
void gctransform_guidingcenter2particle(real mass, real charge, real *B_dB, real R, real Phi, real Z, real ppar, real mu, real zeta, real *r, real *phi, real *z, real *pparprt, real *muprt, real *zetaprt)
Transform guiding center to particle phase space.
Definition gctransform.c:329
GCTRANSFORM_ORDER
static int GCTRANSFORM_ORDER
Definition gctransform.c:53
gctransform_pparmuzeta2prpphipz
void gctransform_pparmuzeta2prpphipz(real mass, real charge, real *B_dB, real phi, real ppar, real mu, real zeta, real *pr, real *pphi, real *pz)
Transform particle ppar, mu, and zeta to momentum vector.
Definition gctransform.c:554
gctransform_particle2guidingcenter
void gctransform_particle2guidingcenter(real mass, real charge, real *B_dB, real r, real phi, real z, real pr, real pphi, real pz, real *R, real *Phi, real *Z, real *ppar, real *mu, real *zeta)
Transform particle to guiding center phase space.
Definition gctransform.c:90
gctransform_setorder
void gctransform_setorder(int order)
Set the order of the transformation.
Definition gctransform.c:62
gctransform.h
Header file for gctransform.c.
fmod
real fmod(real x, real y)
Compute the modulus of two real numbers.
Definition math.c:22
math_jac_rpz2xyz
void math_jac_rpz2xyz(real *rpz, real *xyz, real r, real phi)
Convert a Jacobian from cylindrical to cartesian coordinates.
Definition math.c:44
math.h
Header file for math.c.
math_dot
#define math_dot(a, b)
Calculate dot product a[3] dot b[3].
Definition math.h:28
math_vec_xyz2rpz
#define math_vec_xyz2rpz(vxyz, vrpz, phi)
Transform vector from cartesian to cylindrical basis: vxyz -> vrpz, phi is the toroidal angle in radi...
Definition math.h:90
math_unit
#define math_unit(a, b)
Calculate unit vector b from a 3D vector a.
Definition math.h:70
math_xyz2rpz
#define math_xyz2rpz(xyz, rpz)
Convert cartesian coordinates xyz to cylindrical coordinates rpz.
Definition math.h:74
math_vec_rpz2xyz
#define math_vec_rpz2xyz(vrpz, vxyz, phi)
Transform vector from cylindrical to cartesian basis: vrpz -> vxyz, phi is the toroidal angle in radi...
Definition math.h:83
math_rpz2xyz
#define math_rpz2xyz(rpz, xyz)
Convert cylindrical coordinates rpz to cartesian coordinates xyz.
Definition math.h:78
math_cross
#define math_cross(a, b, c)
Calculate cross product for 3D vectors c = a x b.
Definition math.h:31
math_norm
#define math_norm(a)
Calculate norm of 3D vector a.
Definition math.h:64
physlib.h
Methods to evaluate elementary physical quantities.
Generated on Wed Sep 25 2024 07:10:18 for ASCOT5 by doxygen 1.10.0