Athena: prob/cpaw3d.c Source File

Go to the documentation of this file.
00001 #include "copyright.h"
00002 /*============================================================================*/
00003 /*! \file cpaw3d.c
00004  *  \brief Problem generator for circularly polarized Alfven wave (CPAW) in 3D
00005  *   test.  
00006  *
00007  * PURPOSE: Problem generator for circularly polarized Alfven wave (CPAW) in 3D
00008  *   test.  The angles the wave propagates to the grid is automatically computed
00009  *   to be alpha12 = tan^{-1} (Y/X) and alpha23 = tan^{-1} (Z/Y). 
00010  *
00011  * The wave is defined with reference to a coordinate system (x,y,z) with
00012  *  transformation rules to the code coordinate system (x1,x2,x3)
00013  *
00014  *   First rotate about the y axis:
00015  *   - x' = x*cos(ang_2) - z*sin(ang_2)
00016  *   - y' = y
00017  *   - z' = x*sin(ang_2) + z*cos(ang_2) 
00018  *
00019  *   Next rotate about the z' axis:
00020  *   - x = x'*cos(ang_3) - y'*sin(ang_3)
00021  *   - y = x'*sin(ang_3) + y'*cos(ang_3)
00022  *   - z = z' 
00023  *
00024  *   Expanding this out we get:
00025  *   - x1 = x*cos(ang_2)*cos(ang_3) - y*sin(ang_3) - z*sin(ang_2)*cos(ang_3)
00026  *   - x2 = x*cos(ang_2)*sin(ang_3) - y*cos(ang_3) - z*sin(ang_2)*sin(ang_3)
00027  *   - x3 = x*sin(ang_2)                           + z*cos(ang_2)
00028  *
00029  *   This inverts to:
00030  *   - x =  x1*cos(ang_2)*cos(ang_3) + x2*cos(ang_2)*sin(ang_3) + x3*sin(ang_2)
00031  *   - y = -x1*sin(ang_3)            + x2*cos(ang_3)
00032  *   - z = -x1*sin(ang_2)*cos(ang_3) - x2*sin(ang_2)*sin(ang_3) + x3*cos(ang_2)
00033  *
00034  *   The magnetic field is given by:
00035  *   - B_x = b_par
00036  *   - B_y = b_perp*cos(k*x)
00037  *   - B_z = b_perp*sin(k*x)
00038  *   where k = 2.0*PI/lambda
00039  *
00040  * Note these transformations are the same used in linear_wave3d.c
00041  *
00042  * Can be used for either standing (problem/v_par=1.0) or travelling
00043  * (problem/v_par=0.0) waves.
00044  *
00045  * USERWORK_AFTER_LOOP function computes L1 error norm in solution by comparing
00046  *   to initial conditions.  Problem must be evolved for an integer number of
00047  *   wave periods for this to work.
00048  *
00049  * PRIVATE FUNCTION PROTOTYPES:
00050  * - A1() - 1-component of vector potential for initial conditions
00051  * - A2() - 2-component of vector potential for initial conditions
00052  * - A3() - 3-component of vector potential for initial conditions
00053  *
00054  * REFERENCE: G. Toth,  "The div(B)=0 constraint in shock capturing MHD codes",
00055  *   JCP, 161, 605 (2000)                                                     */
00056 /*============================================================================*/
00057 
00058 #include <math.h>
00059 #include <stdio.h>
00060 #include <stdlib.h>
00061 #include "defs.h"
00062 #include "athena.h"
00063 #include "globals.h"
00064 #include "prototypes.h"
00065 
00066 #ifndef MHD
00067 #error : The cpaw3d test only works for MHD.
00068 #endif
00069 
00070 /* Initial solution, shared with Userwork_after_loop to compute L1 error */
00071 static ConsS ***RootSoln=NULL;
00072 
00073 /* Parameters which define initial solution -- made global so that they can be
00074  * shared with functions A1,2,3 which compute vector potentials */
00075 static Real b_par, b_perp;
00076 static Real ang_2, ang_3; /* Rotation angles about the y and z' axis */
00077 static Real fac, sin_a2, cos_a2, sin_a3, cos_a3;
00078 static Real lambda, k_par; /* Wavelength, 2*PI/wavelength */
00079 
00080 /*==============================================================================
00081  * PRIVATE FUNCTION PROTOTYPES:
00082  * A1() - 1-component of vector potential for initial conditions
00083  * A2() - 2-component of vector potential for initial conditions
00084  * A3() - 3-component of vector potential for initial conditions
00085  *============================================================================*/
00086 
00087 static Real A1(const Real x1, const Real x2, const Real x3);
00088 static Real A2(const Real x1, const Real x2, const Real x3);
00089 static Real A3(const Real x1, const Real x2, const Real x3);
00090 
00091 /*=========================== PUBLIC FUNCTIONS ===============================*/
00092 /*----------------------------------------------------------------------------*/
00093 /* problem:  */
00094 
00095 void problem(DomainS *pDomain)
00096 {
00097   GridS *pGrid = pDomain->Grid;
00098   ConsS ***Soln;
00099   int i=0,j=0,k=0;
00100   int is,ie,js,je,ks,ke,nx1,nx2,nx3,dir;
00101   Real dx1 = pGrid->dx1;
00102   Real dx2 = pGrid->dx2;
00103   Real dx3 = pGrid->dx3;
00104   Real hdx1 = 0.5*pGrid->dx1;
00105   Real hdx2 = 0.5*pGrid->dx2;
00106   Real hdx3 = 0.5*pGrid->dx3;
00107   Real x1size, x2size, x3size;
00108   Real x,x1,x2,x3,cs,sn;
00109   Real v_par, v_perp, den, pres;
00110 #ifdef RESISTIVITY
00111   Real v_A, kva, omega_h, omega_l, omega_r;
00112 #endif
00113 
00114   is = pGrid->is; ie = pGrid->ie;
00115   js = pGrid->js; je = pGrid->je;
00116   ks = pGrid->ks; ke = pGrid->ke;
00117 
00118   nx1 = (ie-is)+1 + 2*nghost;
00119   nx2 = (je-js)+1 + 2*nghost;
00120   nx3 = (ke-ks)+1 + 2*nghost;
00121 
00122   if(pGrid->Nx[2] <= 1)
00123     ath_error("[problem]: cp_alfven3d assumes a 3D grid\n");
00124 
00125   if ((Soln = (ConsS***)calloc_3d_array(nx3,nx2,nx1,sizeof(ConsS)))==NULL)
00126     ath_error("[problem]: Error allocating memory for Soln\n");
00127 
00128   if (pDomain->Level == 0){
00129     if ((RootSoln = (ConsS***)calloc_3d_array(nx3,nx2,nx1,sizeof(ConsS)))
00130       == NULL) ath_error("[problem]: Error allocating memory for RootSoln\n");
00131   }
00132 
00133 /* Imposing periodicity and one wavelength along each grid direction */
00134 
00135   x1size = pDomain->RootMaxX[0] - pDomain->RootMinX[0];
00136   x2size = pDomain->RootMaxX[1] - pDomain->RootMinX[1];
00137   x3size = pDomain->RootMaxX[2] - pDomain->RootMinX[2];
00138 
00139   ang_3 = atan(x1size/x2size);
00140   sin_a3 = sin(ang_3);
00141   cos_a3 = cos(ang_3);
00142 
00143   ang_2 = atan(0.5*(x1size*cos_a3 + x2size*sin_a3)/x3size);
00144   sin_a2 = sin(ang_2);
00145   cos_a2 = cos(ang_2);
00146 
00147   x1 = x1size*cos_a2*cos_a3;
00148   x2 = x2size*cos_a2*sin_a3;
00149   x3 = x3size*sin_a2;
00150 
00151 /* For lambda choose the smaller of the 3 */
00152 
00153   lambda = x1;
00154   lambda = MIN(lambda,x2);
00155   lambda = MIN(lambda,x3);
00156 
00157 /* Initialize k_parallel */
00158 
00159   k_par = 2.0*PI/lambda;
00160   b_par = par_getd("problem","b_par");
00161   den = 1.0;
00162 
00163   ath_pout(0,"va_parallel = %g\n",b_par/sqrt(den));
00164 
00165   b_perp = par_getd("problem","b_perp");
00166   v_perp = b_perp/sqrt((double)den);
00167 
00168   dir    = par_geti_def("problem","dir",1); /* right(1)/left(2) polarization */
00169   if (dir == 1) /* right polarization */
00170     fac = 1.0;
00171   else          /* left polarization */
00172     fac = -1.0;
00173 
00174 #ifdef RESISTIVITY
00175   Q_Hall = par_getd("problem","Q_H");
00176   d_ind  = 0.0;
00177   v_A    = b_par/sqrt((double)den);
00178   if (Q_Hall > 0.0)
00179   {
00180     kva     = k_par*v_A;
00181     omega_h = 1.0/Q_Hall;
00182 
00183     omega_r = 0.5*SQR(kva)/omega_h*(sqrt(1.0+SQR(2.0*omega_h/kva)) + 1.0);
00184     omega_l = 0.5*SQR(kva)/omega_h*(sqrt(1.0+SQR(2.0*omega_h/kva)) - 1.0);
00185 
00186     if (dir == 1) /* right polarization (whistler wave) */
00187       v_perp = v_perp * kva / omega_r;
00188     else          /* left polarization */
00189       v_perp = v_perp * kva / omega_l;
00190   }
00191 #endif
00192 
00193 /* The gas pressure and parallel velocity are free parameters. */
00194 
00195   pres = par_getd("problem","pres");
00196   v_par = par_getd("problem","v_par");
00197 
00198 /* Use the vector potential to initialize the interface magnetic fields */
00199 
00200   for (k=ks; k<=ke+1; k++) {
00201     for (j=js; j<=je+1; j++) {
00202       for (i=is; i<=ie+1; i++) {
00203         cc_pos(pGrid,i,j,k,&x1,&x2,&x3);
00204         x1 -= 0.5*pGrid->dx1;
00205         x2 -= 0.5*pGrid->dx2;
00206         x3 -= 0.5*pGrid->dx3;
00207 
00208         pGrid->B1i[k][j][i] = (A3(x1,x2+dx2 ,x3+hdx3) - A3(x1,x2,x3+hdx3))/dx2 -
00209                               (A2(x1,x2+hdx2,x3+dx3 ) - A2(x1,x2+hdx2,x3))/dx3;
00210 
00211         pGrid->B2i[k][j][i] = (A1(x1+hdx1,x2,x3+dx3 ) - A1(x1+hdx1,x2,x3))/dx3 -
00212                               (A3(x1+dx1 ,x2,x3+hdx3) - A3(x1,x2,x3+hdx3))/dx1;
00213 
00214         pGrid->B3i[k][j][i] = (A2(x1+dx1,x2+hdx2,x3) - A2(x1,x2+hdx2,x3))/dx1 -
00215                               (A1(x1+hdx1,x2+dx2,x3) - A1(x1+hdx1,x2,x3))/dx2;
00216       }
00217     }
00218   }
00219 
00220 /* Now initialize the cell centered quantities */
00221 
00222   for (k=ks; k<=ke; k++) {
00223     for (j=js; j<=je; j++) {
00224       for (i=is; i<=ie; i++) {
00225         cc_pos(pGrid,i,j,k,&x1,&x2,&x3);
00226         x = cos_a2*(x1*cos_a3 + x2*sin_a3) + x3*sin_a2;
00227         sn = sin(k_par*x);
00228         cs = fac*cos(k_par*x);
00229 
00230         pGrid->U[k][j][i].d  = den;
00231         pGrid->U[k][j][i].M1 = den*(v_par*cos_a2*cos_a3 +
00232                                     v_perp*sn*sin_a3 +
00233                                     v_perp*cs*sin_a2*cos_a3);
00234 
00235         pGrid->U[k][j][i].M2 = den*(v_par*cos_a2*sin_a3 -
00236                                     v_perp*sn*cos_a3 +
00237                                     v_perp*cs*sin_a2*sin_a3);
00238 
00239         pGrid->U[k][j][i].M3 = den*(v_par*sin_a2 -
00240                                     v_perp*cs*cos_a2);
00241 
00242         pGrid->U[k][j][i].B1c = 0.5*(pGrid->B1i[k][j][i  ] +
00243                                      pGrid->B1i[k][j][i+1]);
00244         pGrid->U[k][j][i].B2c = 0.5*(pGrid->B2i[k][j  ][i] +
00245                                      pGrid->B2i[k][j+1][i]);
00246         pGrid->U[k][j][i].B3c = 0.5*(pGrid->B3i[k  ][j][i] +
00247                                      pGrid->B3i[k+1][j][i]);
00248 
00249 #ifndef ISOTHERMAL
00250         pGrid->U[k][j][i].E = pres/Gamma_1 
00251           + 0.5*(SQR(pGrid->U[k][j][i].B1c) +
00252                  SQR(pGrid->U[k][j][i].B2c) +
00253                  SQR(pGrid->U[k][j][i].B3c) )
00254           + 0.5*(SQR(pGrid->U[k][j][i].M1) +
00255                  SQR(pGrid->U[k][j][i].M2) +
00256                  SQR(pGrid->U[k][j][i].M3) )/den;
00257 #endif
00258 
00259 /* And store the initial state */
00260         Soln[k][j][i] = pGrid->U[k][j][i];
00261       }
00262     }
00263   }
00264 
00265 /* save solution on root grid */
00266 
00267   if (pDomain->Level == 0) {
00268     for (k=ks; k<=ke; k++) {
00269     for (j=js; j<=je; j++) {
00270     for (i=is; i<=ie; i++) {
00271       RootSoln[k][j][i].d  = Soln[k][j][i].d ;
00272       RootSoln[k][j][i].M1 = Soln[k][j][i].M1;
00273       RootSoln[k][j][i].M2 = Soln[k][j][i].M2;
00274       RootSoln[k][j][i].M3 = Soln[k][j][i].M3;
00275 #ifndef ISOTHERMAL
00276       RootSoln[k][j][i].E  = Soln[k][j][i].E ;
00277 #endif /* ISOTHERMAL */
00278       RootSoln[k][j][i].B1c = Soln[k][j][i].B1c;
00279       RootSoln[k][j][i].B2c = Soln[k][j][i].B2c;
00280       RootSoln[k][j][i].B3c = Soln[k][j][i].B3c;
00281 #if (NSCALARS > 0)
00282       for (n=0; n<NSCALARS; n++)
00283         RootSoln[k][j][i].s[n] = Soln[k][j][i].s[n];
00284 #endif
00285     }}}
00286   }
00287 
00288   return;
00289 }
00290 
00291 /*==============================================================================
00292  * PROBLEM USER FUNCTIONS:
00293  * problem_write_restart() - writes problem-specific user data to restart files
00294  * problem_read_restart()  - reads problem-specific user data from restart files
00295  * get_usr_expr()          - sets pointer to expression for special output data
00296  * get_usr_out_fun()       - returns a user defined output function pointer
00297  * get_usr_par_prop()      - returns a user defined particle selection function
00298  * Userwork_in_loop        - problem specific work IN     main loop
00299  * Userwork_after_loop     - problem specific work AFTER  main loop
00300  *----------------------------------------------------------------------------*/
00301         
00302 void problem_write_restart(MeshS *pM, FILE *fp)
00303 {
00304   return;
00305 }
00306 
00307 void problem_read_restart(MeshS *pM, FILE *fp)
00308 {
00309   return;
00310 }
00311 
00312 
00313 ConsFun_t get_usr_expr(const char *expr){
00314   return NULL;
00315 }
00316 
00317 VOutFun_t get_usr_out_fun(const char *name){
00318   return NULL;
00319 }
00320 
00321 #ifdef RESISTIVITY
00322 void get_eta_user(GridS *pG, int i, int j, int k,
00323                              Real *eta_O, Real *eta_H, Real *eta_A)
00324 {
00325   *eta_O = 0.0;
00326   *eta_H = 0.0;
00327   *eta_A = 0.0;
00328   
00329   return;
00330 }
00331 #endif
00332 
00333 void Userwork_in_loop(MeshS *pM)
00334 {
00335 }
00336 
00337 /*---------------------------------------------------------------------------
00338  * Userwork_after_loop: computes L1-error in CPAW,
00339  * ASSUMING WAVE HAS PROPAGATED AN INTEGER NUMBER OF PERIODS
00340  * Must set parameters in input file appropriately so that this is true
00341  */
00342 
00343 void Userwork_after_loop(MeshS *pM)
00344 {
00345   GridS *pGrid;
00346   int i=0,j=0,k=0;
00347   int is,ie,js,je,ks,ke;
00348   Real rms_error=0.0;
00349   ConsS error,total_error;
00350   FILE *fp;
00351   char *fname;
00352   int Nx1, Nx2, Nx3, count;
00353 #if defined MPI_PARALLEL
00354   double err[8], tot_err[8];
00355   int ierr,myID;
00356 #endif
00357 
00358   total_error.d = 0.0;
00359   total_error.M1 = 0.0;
00360   total_error.M2 = 0.0;
00361   total_error.M3 = 0.0;
00362   total_error.B1c = 0.0;
00363   total_error.B2c = 0.0;
00364   total_error.B3c = 0.0;
00365 #ifndef ISOTHERMAL
00366   total_error.E = 0.0;
00367 #endif /* ISOTHERMAL */
00368 
00369 /* Compute error only on root Grid, which is in Domain[0][0] */
00370 
00371   pGrid=pM->Domain[0][0].Grid;
00372   if (pGrid == NULL) return;
00373 
00374   is = pGrid->is; ie = pGrid->ie;
00375   js = pGrid->js; je = pGrid->je;
00376   ks = pGrid->ks; ke = pGrid->ke;
00377 
00378 /* compute L1 error in each variable, and rms total error */
00379 
00380   for (k=ks; k<=ke; k++) {
00381   for (j=js; j<=je; j++) {
00382     error.d = 0.0;
00383     error.M1 = 0.0;
00384     error.M2 = 0.0;
00385     error.M3 = 0.0;
00386     error.B1c = 0.0;
00387     error.B2c = 0.0;
00388     error.B3c = 0.0;
00389 #ifndef ISOTHERMAL
00390     error.E = 0.0;
00391 #endif /* ISOTHERMAL */
00392 
00393     for (i=is; i<=ie; i++) {
00394       error.d   += fabs(pGrid->U[k][j][i].d   - RootSoln[k][j][i].d);
00395       error.M1  += fabs(pGrid->U[k][j][i].M1  - RootSoln[k][j][i].M1);
00396       error.M2  += fabs(pGrid->U[k][j][i].M2  - RootSoln[k][j][i].M2);
00397       error.M3  += fabs(pGrid->U[k][j][i].M3  - RootSoln[k][j][i].M3); 
00398       error.B1c += fabs(pGrid->U[k][j][i].B1c - RootSoln[k][j][i].B1c);
00399       error.B2c += fabs(pGrid->U[k][j][i].B2c - RootSoln[k][j][i].B2c);
00400       error.B3c += fabs(pGrid->U[k][j][i].B3c - RootSoln[k][j][i].B3c);
00401 #ifndef ISOTHERMAL
00402       error.E   += fabs(pGrid->U[k][j][i].E   -  RootSoln[k][j][i].E);
00403 #endif /* ISOTHERMAL */
00404     }
00405 
00406     total_error.d += error.d;
00407     total_error.M1 += error.M1;
00408     total_error.M2 += error.M2;
00409     total_error.M3 += error.M3;
00410     total_error.B1c += error.B1c;
00411     total_error.B2c += error.B2c;
00412     total_error.B3c += error.B3c;
00413 #ifndef ISOTHERMAL
00414     total_error.E += error.E;
00415 #endif /* ISOTHERMAL */
00416   }}
00417 
00418 #if defined MPI_PARALLEL
00419   Nx1 = pM->Domain[0][0].Nx[0];
00420   Nx2 = pM->Domain[0][0].Nx[1];
00421   Nx3 = pM->Domain[0][0].Nx[2];
00422 #else
00423   Nx1 = ie - is + 1;
00424   Nx2 = je - js + 1;
00425   Nx3 = ke - ks + 1;
00426 #endif
00427   count = Nx1*Nx2*Nx3;
00428 
00429 #ifdef MPI_PARALLEL 
00430 /* Now we have to use an All_Reduce to get the total error over all the MPI
00431  * grids.  Begin by copying the error into the err[] array */
00432   
00433   err[0] = total_error.d;
00434   err[1] = total_error.M1;
00435   err[2] = total_error.M2;
00436   err[3] = total_error.M3;
00437   err[4] = total_error.B1c;
00438   err[5] = total_error.B2c;
00439   err[6] = total_error.B3c;
00440 #ifndef ISOTHERMAL
00441   err[7] = total_error.E;
00442 #endif /* ISOTHERMAL */
00443 
00444   ierr = MPI_Reduce(err,tot_err,8,MPI_DOUBLE,MPI_SUM,0,
00445     pM->Domain[0][0].Comm_Domain);
00446 
00447 /* If I'm the parent, copy the sum back to the total_error variable */
00448 
00449   ierr = MPI_Comm_rank(pM->Domain[0][0].Comm_Domain, &myID);
00450   if(myID == 0){ /* I'm the parent */
00451     total_error.d   = tot_err[0];
00452     total_error.M1  = tot_err[1];
00453     total_error.M2  = tot_err[2];
00454     total_error.M3  = tot_err[3];
00455     total_error.B1c = tot_err[4];
00456     total_error.B2c = tot_err[5];
00457     total_error.B3c = tot_err[6];
00458 #ifndef ISOTHERMAL
00459     total_error.E   = tot_err[7];
00460 #endif /* ISOTHERMAL */
00461   }
00462   else return; /* The child grids do not do any of the following code */
00463 
00464 #endif /* MPI_PARALLEL */
00465 
00466 /* Compute RMS error over all variables, and print out */
00467 
00468   rms_error = SQR(total_error.d) + SQR(total_error.M1) + SQR(total_error.M2)
00469                 + SQR(total_error.M3);
00470   rms_error += SQR(total_error.B1c) + SQR(total_error.B2c) 
00471                + SQR(total_error.B3c);
00472 #ifndef ISOTHERMAL
00473   rms_error += SQR(total_error.E);
00474 #endif /* ISOTHERMAL */
00475   rms_error = sqrt(rms_error)/(double)count;
00476 
00477 /* Print error to file "cpaw3d-errors.dat" */
00478   
00479 #ifdef MPI_PARALLEL
00480   fname = "../cpaw3d-errors.dat";
00481 #else
00482   fname = "cpaw3d-errors.dat";
00483 #endif
00484 /* The file exists -- reopen the file in append mode */
00485   if((fp=fopen(fname,"r")) != NULL){
00486     if((fp = freopen(fname,"a",fp)) == NULL){
00487       ath_error("[Userwork_after_loop]: Unable to reopen file.\n");
00488       return;
00489     }
00490   }
00491 /* The file does not exist -- open the file in write mode */
00492   else{
00493     if((fp = fopen(fname,"w")) == NULL){
00494       ath_perr(-1,"[Userwork_after_loop]: Unable to open file.\n");
00495       free(fname);
00496       return;
00497     }
00498 /* Now write out some header information */
00499     fprintf(fp,"# Nx1  Nx2  Nx3  RMS-Error  d  M1  M2  M3");
00500 #ifndef ISOTHERMAL
00501     fprintf(fp,"  E");
00502 #endif /* ISOTHERMAL */
00503     fprintf(fp,"  B1c  B2c  B3c");
00504     fprintf(fp,"\n#\n");
00505   }
00506 
00507   fprintf(fp,"%d  %d  %d  %e",Nx1,Nx2,Nx3,rms_error);
00508   fprintf(fp,"  %e  %e  %e  %e",
00509           (total_error.d/(double)count),
00510           (total_error.M1/(double)count),
00511           (total_error.M2/(double)count),
00512           (total_error.M3/(double)count) );
00513 #ifndef ISOTHERMAL
00514   fprintf(fp,"  %e",(total_error.E/(double)count) );
00515 #endif /* ISOTHERMAL */
00516   fprintf(fp,"  %e  %e  %e",
00517           (total_error.B1c/(double)count),
00518           (total_error.B2c/(double)count),
00519           (total_error.B3c/(double)count));
00520   fprintf(fp,"\n");
00521 
00522   fclose(fp);
00523 
00524   return;
00525 }
00526 
00527 /*=========================== PRIVATE FUNCTIONS ==============================*/
00528   
00529 /*----------------------------------------------------------------------------*/
00530 /*! \fn static Real A1(const Real x1, const Real x2, const Real x3)
00531  *  \brief A1: 1-component of vector potential, using a gauge such that Ax = 0,
00532  * and Ay, Az are functions of x and y alone.
00533  */
00534 
00535 static Real A1(const Real x1, const Real x2, const Real x3)
00536 {
00537   Real x, y;
00538   Real Ay, Az;
00539 
00540   x =  x1*cos_a2*cos_a3 + x2*cos_a2*sin_a3 + x3*sin_a2;
00541   y = -x1*sin_a3        + x2*cos_a3;
00542 
00543   Ay = fac*(b_perp/k_par)*sin(k_par*x);
00544   Az = (b_perp/k_par)*cos(k_par*x) + b_par*y;
00545 
00546   return -Ay*sin_a3 - Az*sin_a2*cos_a3;
00547 }
00548 
00549 /*----------------------------------------------------------------------------*/
00550 /*! \fn static Real A2(const Real x1, const Real x2, const Real x3)
00551  *  \brief A2: 2-component of vector potential
00552  */
00553 
00554 static Real A2(const Real x1, const Real x2, const Real x3)
00555 {
00556   Real x, y;
00557   Real Ay, Az;
00558 
00559   x =  x1*cos_a2*cos_a3 + x2*cos_a2*sin_a3 + x3*sin_a2;
00560   y = -x1*sin_a3        + x2*cos_a3;
00561 
00562   Ay = fac*(b_perp/k_par)*sin(k_par*x);
00563   Az = (b_perp/k_par)*cos(k_par*x) + b_par*y;
00564 
00565   return Ay*cos_a3 - Az*sin_a2*sin_a3;
00566 }
00567 
00568 /*----------------------------------------------------------------------------*/
00569 /*! \fn static Real A3(const Real x1, const Real x2, const Real x3)
00570  *  \brief A3: 3-component of vector potential
00571  */
00572 
00573 static Real A3(const Real x1, const Real x2, const Real x3)
00574 {
00575   Real x, y;
00576   Real Az;
00577 
00578   x =  x1*cos_a2*cos_a3 + x2*cos_a2*sin_a3 + x3*sin_a2;
00579   y = -x1*sin_a3        + x2*cos_a3;
00580 
00581   Az = (b_perp/k_par)*cos(k_par*x) + b_par*y;
00582 
00583   return Az*cos_a2;
00584 }