Athena: utils.c Source File

Go to the documentation of this file.
00001 #include "copyright.h"
00002 /*============================================================================*/
00003 /*! \file utils.c
00004  *  \brief A variety of useful utility functions.
00005  *
00006  * PURPOSE: A variety of useful utility functions.
00007  *
00008  * CONTAINS PUBLIC FUNCTIONS: 
00009  * - ath_strdup()     - not supplied by fancy ANSI C, but ok in C89 
00010  * - ath_gcd()        - computes greatest common divisor by Euler's method
00011  * - ath_big_endian() - run-time detection of endianism of the host cpu
00012  * - ath_bswap()      - fast byte swapping routine
00013  * - ath_error()      - fatal error routine
00014  * - minmax1()        - fast Min/Max for a 1d array using registers
00015  * - minmax2()        - fast Min/Max for a 2d array using registers
00016  * - minmax3()        - fast Min/Max for a 3d array using registers
00017  *============================================================================*/
00018 
00019 #include <stdio.h>
00020 #include <stdlib.h>
00021 #include <string.h>
00022 #include <stdarg.h>
00023 #include <math.h>
00024 #include "defs.h"
00025 #include "athena.h"
00026 #include "prototypes.h"
00027 
00028 /*----------------------------------------------------------------------------*/
00029 /*! \fn char *ath_strdup(const char *in)
00030  *  \brief This is really strdup(), but strdup is not available in 
00031  *   ANSI  (-pendantic or -ansi will leave it undefined in gcc)
00032  *   much like allocate.
00033  */
00034 
00035 char *ath_strdup(const char *in)
00036 {
00037   char *out = (char *)malloc((1+strlen(in))*sizeof(char));
00038   if(out == NULL) {
00039     ath_perr(-1,"ath_strdup: failed to alloc %d\n",(int)(1+strlen(in)));
00040     return NULL; /* malloc failed */
00041   }
00042   return strcpy(out,in);
00043 }
00044 
00045 /*----------------------------------------------------------------------------*/
00046 /*! \fn int ath_gcd(int a, int b)
00047  *  \brief Calculate the Greatest Common Divisor by Euler's method
00048  */
00049 
00050 int ath_gcd(int a, int b)
00051 {
00052   int c;
00053   if(b>a) {c=a; a=b; b=c;} 
00054   while((c=a%b)) {a=b; b=c;}
00055   return b;
00056 }
00057 
00058 /*----------------------------------------------------------------------------*/
00059 /*! \fn int ath_big_endian(void)
00060  *  \brief Return 1 if the machine is big endian (e.g. Sun, PowerPC)
00061  * return 0 if not (e.g. Intel)
00062  */
00063 
00064 int ath_big_endian(void)
00065 {
00066   short int n = 1;
00067   char *ep = (char *)&n;
00068 
00069   return (*ep == 0); /* Returns 1 on a big endian machine */
00070 }
00071 
00072 /*----------------------------------------------------------------------------*/
00073 /*! \fn void ath_bswap(void *vdat, int len, int cnt)
00074  *  \brief Swap bytes, code stolen from NEMO  
00075  */
00076  
00077 void ath_bswap(void *vdat, int len, int cnt)
00078 {
00079   char tmp, *dat = (char *) vdat;
00080   int k;
00081  
00082   if (len==1)
00083     return;
00084   else if (len==2)
00085     while (cnt--) {
00086       tmp = dat[0];  dat[0] = dat[1];  dat[1] = tmp;
00087       dat += 2;
00088     }
00089   else if (len==4)
00090     while (cnt--) {
00091       tmp = dat[0];  dat[0] = dat[3];  dat[3] = tmp;
00092       tmp = dat[1];  dat[1] = dat[2];  dat[2] = tmp;
00093       dat += 4;
00094     }
00095   else if (len==8)
00096     while (cnt--) {
00097       tmp = dat[0];  dat[0] = dat[7];  dat[7] = tmp;
00098       tmp = dat[1];  dat[1] = dat[6];  dat[6] = tmp;
00099       tmp = dat[2];  dat[2] = dat[5];  dat[5] = tmp;
00100       tmp = dat[3];  dat[3] = dat[4];  dat[4] = tmp;
00101       dat += 8;
00102     }
00103   else {  /* the general SLOOOOOOOOOW case */
00104     for(k=0; k<len/2; k++) {
00105       tmp = dat[k];
00106       dat[k] = dat[len-1-k];
00107       dat[len-1-k] = tmp;
00108     }
00109   }
00110 }
00111 
00112 /*----------------------------------------------------------------------------*/
00113 /*! \fn void ath_error(char *fmt, ...)
00114  *  \brief Terminate execution and output error message
00115  *  Uses variable-length argument lists provided in <stdarg.h>
00116  */
00117 
00118 void ath_error(char *fmt, ...)
00119 {
00120   va_list ap;
00121    FILE *atherr = atherr_fp();
00122 
00123   fprintf(atherr,"### Fatal error: ");   /* prefix */
00124   va_start(ap, fmt);              /* ap starts with string 'fmt' */
00125   vfprintf(atherr, fmt, ap);      /* print out on atherr */
00126   fflush(atherr);                 /* flush it NOW */
00127   va_end(ap);                     /* end varargs */
00128 
00129 #ifdef MPI_PARALLEL
00130   MPI_Abort(MPI_COMM_WORLD, 1);
00131 #endif
00132 
00133   exit(EXIT_FAILURE);
00134 }
00135 
00136 /*----------------------------------------------------------------------------*/
00137 /*! \fn void minmax1(Real *data, int nx1, Real *dmino, Real *dmaxo)
00138  *  \brief Return the min and max of a 1D array using registers
00139  *  Works on data of type float, not Real.
00140  */
00141 
00142 void minmax1(Real *data, int nx1, Real *dmino, Real *dmaxo)
00143 {
00144   int i;
00145   register Real dmin, dmax;
00146 
00147   dmin = dmax = data[0];
00148   for (i=0; i<nx1; i++) {
00149     dmin = MIN(dmin,data[i]);
00150     dmax = MAX(dmax,data[i]);
00151   }
00152   *dmino = dmin;
00153   *dmaxo = dmax;
00154 }
00155 
00156 /*! \fn void minmax2(Real **data, int nx2, int nx1, Real *dmino, Real *dmaxo)
00157  *  \brief Return the min and max of a 2D array using registers
00158  *  Works on data of type float, not Real.
00159  */
00160 void minmax2(Real **data, int nx2, int nx1, Real *dmino, Real *dmaxo)
00161 {
00162   int i,j;
00163   register Real dmin, dmax;
00164 
00165   dmin = dmax = data[0][0];
00166   for (j=0; j<nx2; j++) {
00167     for (i=0; i<nx1; i++) {
00168       dmin = MIN(dmin,data[j][i]);
00169       dmax = MAX(dmax,data[j][i]);
00170     }
00171   }
00172   *dmino = dmin;
00173   *dmaxo = dmax;
00174 }
00175 
00176 /*! \fn void minmax3(Real ***data, int nx3, int nx2, int nx1, Real *dmino, 
00177  *                   Real *dmaxo)
00178  *  \brief Return the min and max of a 3D array using registers
00179  *  Works on data of type float, not Real.
00180  */
00181 void minmax3(Real ***data, int nx3, int nx2, int nx1, Real *dmino, Real *dmaxo)
00182 {
00183   int i,j,k;
00184   register Real dmin, dmax;
00185 
00186   dmin = dmax = data[0][0][0];
00187   for (k=0; k<nx3; k++) {
00188     for (j=0; j<nx2; j++) {
00189       for (i=0; i<nx1; i++) {
00190         dmin = MIN(dmin,data[k][j][i]);
00191         dmax = MAX(dmax,data[k][j][i]);
00192       }
00193     }
00194   }
00195   *dmino = dmin;
00196   *dmaxo = dmax;
00197 }
00198 
00199 /*----------------------------------------------------------------------------*/
00200 /*! \fn  void do_nothing_bc(GridS *pG)
00201  *
00202  *  \brief DOES ABSOLUTELY NOTHING!  THUS, WHATEVER THE BOUNDARY ARE SET TO 
00203  *  INITIALLY, THEY REMAIN FOR ALL TIME.
00204  */
00205 void do_nothing_bc(GridS *pG)
00206 {
00207 }
00208 
00209 /*============================================================================
00210  * ERROR-ANALYSIS FUNCTIONS
00211  *============================================================================*/
00212 
00213 /*----------------------------------------------------------------------------*/
00214 /*! \fn Real compute_div_b(GridS *pG)
00215  *  \brief COMPUTE THE DIVERGENCE OF THE MAGNETIC FIELD USING FACE-CENTERED 
00216  *  FIELDS OVER THE ENTIRE ACTIVE GRID.  RETURNS THE MAXIMUM OF |DIV B|.
00217  */
00218 Real compute_div_b(GridS *pG)
00219 {
00220 #ifdef MHD
00221   int i,j,k,is,ie,js,je,ks,ke;
00222   Real x1,x2,x3,divB,maxdivB=0.0;
00223   Real lsf=1.0,rsf=1.0,dx2=pG->dx2;
00224 
00225   is = pG->is; ie = pG->ie;
00226   js = pG->js; je = pG->je;
00227   ks = pG->ks; ke = pG->ke;
00228 
00229   for (k=ks; k<=ke; k++) {
00230     for (j=js; j<=je; j++) {
00231       for (i=is; i<=ie; i++) {
00232         cc_pos(pG,i,j,k,&x1,&x2,&x3);
00233 #ifdef CYLINDRICAL
00234         rsf = (x1+0.5*pG->dx1)/x1;  lsf = (x1-0.5*pG->dx1)/x1;
00235         dx2 = x1*pG->dx2;
00236 #endif
00237         divB = (rsf*pG->B1i[k][j][i+1] - lsf*pG->B1i[k][j][i])/pG->dx1;
00238         if (je > js)
00239           divB += (pG->B2i[k][j+1][i] - pG->B2i[k][j][i])/dx2;
00240         if (ke > ks)
00241           divB += (pG->B3i[k+1][j][i] - pG->B3i[k][j][i])/pG->dx3;
00242 
00243         maxdivB = MAX(maxdivB,fabs(divB));
00244       }
00245     }
00246   }
00247 
00248   return maxdivB;
00249 
00250 #else
00251   fprintf(stderr,"[compute_div_b]: This only works for MHD!\n");
00252   exit(EXIT_FAILURE);
00253   return 0.0;
00254 #endif /* MHD */
00255 }
00256 
00257 /*----------------------------------------------------------------------------*/
00258 /*! \fn void compute_l1_error(const char *problem, const MeshS *pM, 
00259  *                            const ConsS ***RootSoln, const int errortype) 
00260  *  \brief COMPUTE THE L1-ERRORS IN ALL VARIABLES AT THE CURRENT 
00261  *  (USUALLY THE FINAL)
00262  *  TIMESTEP USING THE INITIAL SOLUTION.  
00263  *
00264  *  THIS MEANS THAT THE SOLUTION MUST
00265  *  EITHER BE STATIC (STEADY-STATE) OR MUST HAVE COMPLETED A FULL PERIOD OF
00266  *  ROTATION, ETC.  FOR THE ERRORTYPE FLAG, 0 MEANS ABSOLUTE ERROR, AND
00267  *  1 MEANS AVERAGE ERROR PER GRID CELL.
00268  */
00269 void compute_l1_error(const char *problem, const MeshS *pM, 
00270                       const ConsS ***RootSoln, const int errortype)
00271 {
00272   DomainS *pD=&(pM->Domain[0][0]);
00273   GridS   *pG=pM->Domain[0][0].Grid;
00274   int i=0,j=0,k=0;
00275 #if (NSCALARS > 0)
00276    int n;
00277 #endif
00278   int is,ie,js,je,ks,ke;
00279   Real rms_error=0.0;
00280   Real x1,x2,x3,dVol,totVol;
00281   ConsS error,total_error;
00282   FILE *fp;
00283   char *fname, fnamestr[256];
00284   int Nx1,Nx2,Nx3;
00285 #if defined MPI_PARALLEL
00286   double err[8+NSCALARS], tot_err[8+NSCALARS];
00287   int mpi_err;
00288 #endif
00289 
00290   /* Clear out the total_error struct */
00291   memset(&total_error,0.0,sizeof(ConsS));
00292   if (pG == NULL) return;
00293 
00294 /* compute L1 error in each variable, and rms total error */
00295 
00296   is = pG->is; ie = pG->ie;
00297   js = pG->js; je = pG->je;
00298   ks = pG->ks; ke = pG->ke;
00299 
00300   for (k=ks; k<=ke; k++) {
00301     for (j=js; j<=je; j++) {
00302       memset(&error,0.0,sizeof(ConsS));
00303       for (i=is; i<=ie; i++) {
00304         dVol = 1.0;
00305         if (pG->dx1 > 0.0) dVol *= pG->dx1;
00306         if (pG->dx2 > 0.0) dVol *= pG->dx2;
00307         if (pG->dx3 > 0.0) dVol *= pG->dx3;
00308 #ifdef CYLINDRICAL
00309         cc_pos(pG,i,j,k,&x1,&x2,&x3);
00310         dVol *= x1;
00311 #endif
00312 
00313         /* Sum local L1 error for each grid cell I own */
00314         error.d   += dVol*fabs(pG->U[k][j][i].d   - RootSoln[k][j][i].d );
00315         error.M1  += dVol*fabs(pG->U[k][j][i].M1  - RootSoln[k][j][i].M1);
00316         error.M2  += dVol*fabs(pG->U[k][j][i].M2  - RootSoln[k][j][i].M2);
00317         error.M3  += dVol*fabs(pG->U[k][j][i].M3  - RootSoln[k][j][i].M3); 
00318 #ifdef MHD
00319         error.B1c += dVol*fabs(pG->U[k][j][i].B1c - RootSoln[k][j][i].B1c);
00320         error.B2c += dVol*fabs(pG->U[k][j][i].B2c - RootSoln[k][j][i].B2c);
00321         error.B3c += dVol*fabs(pG->U[k][j][i].B3c - RootSoln[k][j][i].B3c);
00322 #endif /* MHD */
00323 #ifndef ISOTHERMAL
00324         error.E   += dVol*fabs(pG->U[k][j][i].E   - RootSoln[k][j][i].E );
00325 #endif /* ISOTHERMAL */
00326 #if (NSCALARS > 0)
00327         for (n=0; n<NSCALARS; n++)
00328           error.s[n] += dVol*fabs(pG->U[k][j][i].s[n] - RootSoln[k][j][i].s[n]);;
00329 #endif
00330       }
00331 
00332       /* total_error is sum of local L1 error */
00333       total_error.d += error.d;
00334       total_error.M1 += error.M1;
00335       total_error.M2 += error.M2;
00336       total_error.M3 += error.M3;
00337 #ifdef MHD
00338       total_error.B1c += error.B1c;
00339       total_error.B2c += error.B2c;
00340       total_error.B3c += error.B3c;
00341 #endif /* MHD */
00342 #ifndef ISOTHERMAL
00343       total_error.E += error.E;
00344 #endif /* ISOTHERMAL */
00345 #if (NSCALARS > 0)
00346       for (n=0; n<NSCALARS; n++) total_error.s[n] += error.s[n];
00347 #endif
00348     }
00349   }
00350 
00351 #ifdef MPI_PARALLEL
00352 /* Now we have to use an All_Reduce to get the total error over all the MPI
00353  * grids.  Begin by copying the error into the err[] array */
00354 
00355   err[0] = total_error.d;
00356   err[1] = total_error.M1;
00357   err[2] = total_error.M2;
00358   err[3] = total_error.M3;
00359 #ifdef MHD
00360   err[4] = total_error.B1c;
00361   err[5] = total_error.B2c;
00362   err[6] = total_error.B3c;
00363 #endif /* MHD */
00364 #ifndef ISOTHERMAL
00365   err[7] = total_error.E;
00366 #endif /* ISOTHERMAL */
00367 #if (NSCALARS > 0)
00368   for (n=0; n<NSCALARS; n++) err[8+n] = total_error.s[n];
00369 #endif
00370 
00371 /* Sum up the Computed Error */
00372   mpi_err = MPI_Reduce(err, tot_err, (8+NSCALARS),
00373                        MPI_DOUBLE, MPI_SUM, 0, pD->Comm_Domain);
00374   if(mpi_err)
00375     ath_error("[compute_l1_error]: MPI_Reduce call returned error = %d\n",
00376               mpi_err);
00377 
00378 /* If I'm the parent, copy the sum back to the total_error variable */
00379   if(pD->DomNumber == 0){ /* I'm the parent */
00380     total_error.d   = tot_err[0];
00381     total_error.M1  = tot_err[1];
00382     total_error.M2  = tot_err[2];
00383     total_error.M3  = tot_err[3];
00384 #ifdef MHD
00385     total_error.B1c = tot_err[4];
00386     total_error.B2c = tot_err[5];
00387     total_error.B3c = tot_err[6];
00388 #endif /* MHD */
00389 #ifndef ISOTHERMAL
00390     total_error.E   = tot_err[7];
00391 #endif /* ISOTHERMAL */
00392 #if (NSCALARS > 0)
00393   for (n=0; n<NSCALARS; n++) total_error.s[n] = err[8+n];
00394 #endif
00395 
00396   }
00397   else return; /* The child grids do not do any of the following code */
00398 #endif /* MPI_PARALLEL */
00399 
00400   /* Compute total number of grid cells */
00401   Nx1 = pD->Nx[0];
00402   Nx2 = pD->Nx[1];
00403   Nx3 = pD->Nx[2];
00404 
00405   totVol = 1.0;
00406   if (errortype == 1) {
00407     if (pD->MaxX[0] > pD->MinX[0]) totVol *= pD->MaxX[0] - pD->MinX[0];
00408     if (pD->MaxX[1] > pD->MinX[1]) totVol *= pD->MaxX[1] - pD->MinX[1];
00409     if (pD->MaxX[2] > pD->MinX[2]) totVol *= pD->MaxX[2] - pD->MinX[2];
00410 #ifdef CYLINDRICAL
00411     totVol *= 0.5*(pD->MinX[0] + pD->MaxX[0]);
00412 #endif
00413   }
00414 
00415 
00416 /* Compute RMS error over all variables, and print out */
00417 
00418   rms_error = SQR(total_error.d) + SQR(total_error.M1) + SQR(total_error.M2)
00419                 + SQR(total_error.M3);
00420 #ifdef MHD
00421   rms_error += SQR(total_error.B1c) + SQR(total_error.B2c) 
00422                + SQR(total_error.B3c);
00423 #endif /* MHD */
00424 #ifndef ISOTHERMAL
00425   rms_error += SQR(total_error.E);
00426 #endif /* ISOTHERMAL */
00427 
00428   rms_error = sqrt(rms_error)/totVol;
00429 
00430 /* Print error to file "BLAH-errors.#.dat"  */
00431    sprintf(fnamestr,"%s-errors",problem);
00432    fname = ath_fname(NULL,fnamestr,NULL,NULL,1,0,NULL,"dat");
00433 
00434 /* The file exists -- reopen the file in append mode */
00435   if((fp=fopen(fname,"r")) != NULL){
00436     if((fp = freopen(fname,"a",fp)) == NULL){
00437       ath_error("[compute_l1_error]: Unable to reopen file.\n");
00438       free(fname);
00439       return;
00440     }
00441   }
00442 /* The file does not exist -- open the file in write mode */
00443   else{
00444     if((fp = fopen(fname,"w")) == NULL){
00445       ath_error("[compute_l1_error]: Unable to open file.\n");
00446       free(fname);
00447       return;
00448     }
00449 /* Now write out some header information */
00450     fprintf(fp,"# Nx1  Nx2  Nx3  RMS-Error  d  M1  M2  M3");
00451 #ifndef ISOTHERMAL
00452     fprintf(fp,"  E");
00453 #endif /* ISOTHERMAL */
00454 #ifdef MHD
00455     fprintf(fp,"  B1c  B2c  B3c");
00456 #endif /* MHD */
00457 #if (NSCALARS > 0)
00458     for (n=0; n<NSCALARS; n++) {
00459       fprintf(fp,"  S[ %d ]",n);
00460     }
00461 #endif
00462     fprintf(fp,"\n#\n");
00463   }
00464 
00465   fprintf(fp,"%d  %d  %d  %e",Nx1,Nx2,Nx3,rms_error);
00466 
00467   fprintf(fp,"  %e  %e  %e  %e",
00468           (total_error.d /totVol),
00469           (total_error.M1/totVol),
00470           (total_error.M2/totVol),
00471           (total_error.M3/totVol));
00472 
00473 #ifndef ISOTHERMAL
00474   fprintf(fp,"  %e",total_error.E/totVol);
00475 #endif /* ISOTHERMAL */
00476 
00477 #ifdef MHD
00478   fprintf(fp,"  %e  %e  %e",
00479           (total_error.B1c/totVol),
00480           (total_error.B2c/totVol),
00481           (total_error.B3c/totVol));
00482 #endif /* MHD */
00483 #if (NSCALARS > 0)
00484     for (n=0; n<NSCALARS; n++) {
00485       fprintf(fp,"  %e",total_error.s[n]/totVol);
00486     }
00487 #endif
00488 
00489   fprintf(fp,"\n");
00490 
00491   fclose(fp);
00492   free(fname);
00493 
00494   return;
00495 }
00496 
00497 /*============================================================================
00498  * ROOT-FINDING FUNCTIONS
00499  *============================================================================*/
00500 
00501 /*----------------------------------------------------------------------------*/
00502 /*! \fn int sign_change(Real (*func)(const Real,const Real), const Real a0, 
00503  *                      const Real b0, const Real x, Real *a, Real *b)
00504  *  \brief SEARCH FOR A SIGN CHANGE.  
00505  *
00506  *  THIS FUNCTION PARTITIONS THE INTERVAL (a0,b0) INTO
00507  *  2^k EQUALLY SPACED GRID POINTS, EVALUATES THE FUNCTION f AT THOSE POINTS,
00508  *  AND THEN SEARCHES FOR A SIGN CHANGE IN f BETWEEN ADJACENT GRID POINTS.  THE
00509  *  FIRST SUCH INTERVAL FOUND, (a,b), IS RETURNED.
00510  */
00511 int sign_change(Real (*func)(const Real,const Real), const Real a0, 
00512                              const Real b0, const Real x, Real *a, Real *b) {
00513   const int kmax=20;
00514   int k, n, i;
00515   Real delta, fk, fkp1;
00516 
00517   for (k=1; k<=kmax; k++) {
00518     n = pow(2,k);
00519     delta = (b0-a0)/(n-1);
00520     *a = a0;
00521     fk = func(x,*a);
00522     for (i=1; i<n; i++) {
00523       *b = *a + delta;
00524       fkp1 = func(x,*b);
00525       if (fkp1*fk < 0)
00526         return 1;
00527       *a = *b;
00528       fk = fkp1;
00529     }
00530   }
00531 //   ath_error("[sign_change]: No sign change was detected in (%f,%f) for x=%f!\n",a0,b0,x);
00532   return 0;
00533 } 
00534 
00535 
00536 /*----------------------------------------------------------------------------*/
00537 /*! \fn int bisection(Real (*func)(const Real,const Real), const Real a0, 
00538  *                    const Real b0, const Real x, Real *root) 
00539  *  \brief THIS FUNCTION IMPLEMENTS THE BISECTION METHOD FOR ROOT FINDING.
00540  */
00541 int bisection(Real (*func)(const Real,const Real), const Real a0, const Real b0,
00542                            const Real x, Real *root) 
00543 {
00544   const Real tol = 1.0E-10;
00545   const int maxiter = 400;
00546   Real a=a0, b=b0, c, fa, fb, fc;
00547   int i;
00548 
00549   fa = func(x,a);
00550   fb = func(x,b);
00551   if (fabs(fa) < tol) {
00552     *root = a;
00553     return 1;
00554   }
00555   if (fabs(fb) < tol) {
00556     *root = b;
00557     return 1;
00558   }
00559 // printf("fa = %f, fb = %f\n", fa, fb);
00560 
00561   for (i = 0; i < maxiter; i++) {
00562     c = 0.5*(a+b);
00563 // printf("x = %f, a = %f, b = %f, c = %f\n", x,a,b,c);
00564 #ifdef MYDEBUG
00565     printf("c = %f\n", c);
00566 #endif
00567     if (fabs((b-a)/c) < tol) {
00568 #ifdef MYDEBUG
00569       printf("Bisection converged within tolerance of %f!\n", eps);
00570 #endif
00571       *root = c;
00572       return 1;
00573     }
00574     fc = func(x,c);
00575     if (fa*fc < 0) {
00576       b = c;
00577       fb = fc;
00578     }
00579     else if (fc*fb < 0) {
00580       a = c;
00581       fa = fc;
00582     }
00583     else if (fc == 0) {
00584       *root = c;
00585       return 1;
00586     }
00587     else {
00588       ath_error("[bisection]:  There is no single root in (%f,%f) for x = %13.10f!!\n", a, b,x);
00589       *root = c;
00590       return 0;
00591     }
00592   }
00593 
00594   ath_error("[bisection]:  Bisection did not converge in %d iterations for x = %13.10f!!\n", maxiter,x);
00595   *root = c;
00596   return 0;
00597 }
00598 
00599 
00600 /*============================================================================
00601  * QUADRATURE FUNCTIONS
00602  *============================================================================*/
00603 
00604 /*----------------------------------------------------------------------------*/
00605 /*! \fn Real trapzd(Real (*func)(Real), const Real a, const Real b, const int n,
00606  *                  const Real s)
00607  *  \brief THIS ROUTINE COMPUTES THE nTH STAGE OF REFINEMENT OF AN EXTENDED 
00608  *  TRAPEZOIDAL RULE.  
00609  *
00610  * func IS INPUT AS A POINTER TO THE FUNCTION TO BE INTEGRATED BETWEEN 
00611  * LIMITS a AND b, ALSO INPUT. WHEN CALLED WITH n=1, THE ROUTINE RETURNS THE 
00612  * CRUDEST ESTIMATE OF \int_a^b f(R) R dR.  SUBSEQUENT CALLS WITH n=2,3,... 
00613  * (IN THAT SEQUENTIAL ORDER) WILL IMPROVE THE ACCURACY BY ADDING 2n-2 
00614  * ADDITIONAL INTERIOR POINTS. 
00615  * ADAPTED FROM NUMERICAL RECIPES BY AARON SKINNER 
00616  */
00617 Real trapzd(Real (*func)(Real), const Real a, const Real b, const int n, 
00618             const Real s)
00619 {
00620   Real x,tnm,sum,dx;
00621   int it,j;
00622 
00623   if (n == 1) {
00624     return 0.5*(b-a)*(func(a)+func(b));
00625   } 
00626   else {
00627     for (it=1,j=1; j<n-1; j++) it <<= 1;  /* it = 2^(n-2) */
00628     tnm = it;
00629     dx = (b-a)/tnm;  /* THIS IS THE SPACING OF THE POINTS TO BE ADDED. */
00630     x = a + 0.5*dx;
00631     for (sum=0.0,j=1; j<=it; j++,x+=dx) sum += func(x);
00632     return 0.5*(s+(b-a)*sum/tnm);  /* THIS REPLACES s BY ITS REFINED VALUE. */
00633   }
00634 }
00635 
00636 /*----------------------------------------------------------------------------*/
00637 /*! \fn Real qsimp(Real (*func)(Real), const Real a, const Real b)
00638  *  \brief RETURNS THE INTEGRAL OF THE FUNCTION func FROM a TO b. 
00639  *
00640  * THE PARAMETER EPS 
00641  * CAN BE SET TO THE DESIRED FRACTIONAL ACCURACY AND JMAX SO THAT 2^(JMAX-1) 
00642  * IS THE MAXIMUM ALLOWED NUMBER OF STEPS.  INTEGRATION IS PERFORMED BY 
00643  * SIMPSON'S RULE.
00644  * ADAPTED FROM NUMERICAL RECIPES BY AARON SKINNER 
00645  */
00646 
00647 #define EPS 1.0e-8
00648 #define JMAX 20
00649 
00650 Real qsimp(Real (*func)(Real), const Real a, const Real b) 
00651 {
00652   int j;
00653   Real s,st,ost,os;
00654 
00655   ost = os = -1.0e30;
00656   for (j=1; j<=JMAX; j++) {
00657     st = trapzd(func,a,b,j,ost);
00658     s = (4.0*st-ost)/3.0;  /* EQUIVALENT TO SIMPSON'S RULE */
00659     if (j > 5)  /* AVOID SPURIOUS EARLY CONVERGENCE. */
00660       if (fabs(s-os) < EPS*fabs(os) || (s == 0.0 && os == 0.0)) return s;
00661     os=s;
00662     ost=st;
00663   }
00664 
00665   ath_error("[qsimp]:  Too many steps!\n");
00666   return 0.0;
00667 }
00668 
00669 
00670 /*----------------------------------------------------------------------------*/
00671 /* FUNCTION avg1d,avg2d,avg3d
00672  *
00673  * RETURNS THE INTEGRAL OF A USER-SUPPLIED FUNCTION func OVER THE ONE-, TWO-,
00674  * OR THREE-DIMENSIONAL GRID CELL (i,j,k).  INTEGRATION IS PERFORMED USING qsimp.
00675  * ADAPTED FROM NUMERICAL RECIPES BY AARON SKINNER 
00676  */
00677 static Real xsav,ysav,zsav,xmin,xmax,ymin,ymax,zmin,zmax;
00678 static Real (*nrfunc)(Real,Real,Real);
00679 
00680 /*! \fn Real avg1d(Real (*func)(Real, Real, Real), const GridS *pG, 
00681  *                 const int i, const int j, const int k)
00682  *  \brief RETURNS THE INTEGRAL OF A USER-SUPPLIED FUNCTION func OVER THE ONE-
00683  * DIMENSIONAL GRID CELL (i,j,k).  
00684  *
00685  * INTEGRATION IS PERFORMED USING qsimp.
00686  * ADAPTED FROM NUMERICAL RECIPES BY AARON SKINNER 
00687  */
00688 Real avg1d(Real (*func)(Real, Real, Real), const GridS *pG, 
00689             const int i, const int j, const int k)
00690 {
00691   Real x1,x2,x3,dvol=pG->dx1;
00692   Real fx(Real x);
00693 
00694   nrfunc=func;
00695   cc_pos(pG,i,j,k,&x1,&x2,&x3);
00696   xmin = x1 - 0.5*pG->dx1;  xmax = x1 + 0.5*pG->dx1;
00697 
00698   ysav = x2;
00699   zsav = x3;
00700 #ifdef CYLINDRICAL
00701   dvol *= x1;
00702 #endif
00703 
00704   return qsimp(fx,xmin,xmax)/dvol;
00705 }
00706 
00707 /*! \fn Real avg2d(Real (*func)(Real, Real, Real), const GridS *pG, 
00708  *                 const int i, const int j, const int k)
00709  *  \brief RETURNS THE INTEGRAL OF A USER-SUPPLIED FUNCTION func OVER THE TWO-
00710  * DIMENSIONAL GRID CELL (i,j,k).  
00711  *
00712  * INTEGRATION IS PERFORMED USING qsimp.
00713  * ADAPTED FROM NUMERICAL RECIPES BY AARON SKINNER 
00714  */
00715 Real avg2d(Real (*func)(Real, Real, Real), const GridS *pG, 
00716             const int i, const int j, const int k)
00717 {
00718   Real x1,x2,x3,dvol=pG->dx1*pG->dx2;
00719   Real fy(Real y);
00720 
00721   nrfunc=func;
00722   cc_pos(pG,i,j,k,&x1,&x2,&x3);
00723   xmin = x1 - 0.5*pG->dx1;  xmax = x1 + 0.5*pG->dx1;
00724   ymin = x2 - 0.5*pG->dx2;  ymax = x2 + 0.5*pG->dx2;
00725 
00726   zsav = x3;
00727 #ifdef CYLINDRICAL
00728   dvol *= x1;
00729 #endif
00730 
00731   return qsimp(fy,ymin,ymax)/dvol;
00732 }
00733 
00734 /*! \fn Real avg3d(Real (*func)(Real, Real, Real), const GridS *pG, 
00735  *                 const int i, const int j, const int k)
00736  *  \brief RETURNS THE INTEGRAL OF A USER-SUPPLIED FUNCTION func OVER THE
00737  * THREE-DIMENSIONAL GRID CELL (i,j,k).  
00738  *
00739  * INTEGRATION IS PERFORMED USING qsimp.
00740  * ADAPTED FROM NUMERICAL RECIPES BY AARON SKINNER
00741  */
00742 Real avg3d(Real (*func)(Real, Real, Real), const GridS *pG, 
00743             const int i, const int j, const int k)
00744 {
00745   Real x1,x2,x3,dvol=pG->dx1*pG->dx2*pG->dx3;
00746   Real fz(Real z);
00747 
00748   nrfunc=func;
00749   cc_pos(pG,i,j,k,&x1,&x2,&x3);
00750   xmin = x1 - 0.5*pG->dx1;  xmax = x1 + 0.5*pG->dx1;
00751   ymin = x2 - 0.5*pG->dx2;  ymax = x2 + 0.5*pG->dx2;
00752   zmin = x3 - 0.5*pG->dx3;  zmax = x3 + 0.5*pG->dx3;
00753 
00754 #ifdef CYLINDRICAL
00755   dvol *= x1;
00756 #endif
00757 
00758   return qsimp(fz,zmin,zmax)/dvol;
00759 }
00760 
00761 Real avgXZ(Real (*func)(Real, Real, Real), const GridS *pG, const int i, const int j, const int k) {
00762   Real x1,x2,x3;
00763 
00764   Real fXZ(Real z);
00765 
00766   nrfunc=func;
00767   cc_pos(pG,i,j,k,&x1,&x2,&x3);
00768   xmin = x1 - 0.5*pG->dx1;  xmax = x1 + 0.5*pG->dx1;
00769   zmin = x3 - 0.5*pG->dx3;  zmax = x3 + 0.5*pG->dx3;
00770 
00771   ysav = x2;
00772   return qsimp(fXZ,zmin,zmax)/(x1*pG->dx1*pG->dx3);
00773 
00774 }
00775 
00776 Real fz(Real z)
00777 {
00778   Real fy(Real y);
00779 
00780   zsav = z;
00781   return qsimp(fy,ymin,ymax);
00782 }
00783 
00784 Real fy(Real y)
00785 {
00786   Real fx(Real x);
00787 
00788   ysav = y;
00789   return qsimp(fx,xmin,xmax);
00790 }
00791 
00792 Real fx(Real x)
00793 {
00794 #ifdef CYLINDRICAL
00795   return x*nrfunc(x,ysav,zsav);
00796 #else
00797   return nrfunc(x,ysav,zsav);
00798 #endif
00799 }
00800 
00801 Real fXZ(Real z) {
00802         Real fx(Real x);
00803 
00804         zsav = z;
00805         return qsimp(fx,xmin,xmax);
00806 
00807 }
00808 
00809 /*----------------------------------------------------------------------------*/
00810 /* FUNCTION vecpot2b1i,vecpot2b2i,vecpot2b3i
00811  *
00812  * THESE FUNCTIONS COMPUTE MAGNETIC FIELD COMPONENTS FROM COMPONENTS OF A
00813  * SPECIFIED VECTOR POTENTIAL USING STOKES' THEOREM AND SIMPSON'S QUADRATURE.
00814  * NOTE:  THIS IS ONLY GUARANTEED TO WORK IF THE POTENTIAL IS OF CLASS C^1.
00815  * WRITTEN BY AARON SKINNER.
00816  */
00817 
00818 static Real (*a1func)(Real,Real,Real);
00819 static Real (*a2func)(Real,Real,Real);
00820 static Real (*a3func)(Real,Real,Real);
00821 
00822 /*! \fn Real vecpot2b1i(Real (*A2)(Real,Real,Real), Real (*A3)(Real,Real,Real),
00823  *                const GridS *pG, const int i, const int j, const int k)
00824  *  \brief Compute B-field components from a vector potential.
00825  *
00826  * THESE FUNCTIONS COMPUTE MAGNETIC FIELD COMPONENTS FROM COMPONENTS OF A
00827  * SPECIFIED VECTOR POTENTIAL USING STOKES' THEOREM AND SIMPSON'S QUADRATURE.
00828  * NOTE:  THIS IS ONLY GUARANTEED TO WORK IF THE POTENTIAL IS OF CLASS C^1.
00829  * WRITTEN BY AARON SKINNER.
00830  */
00831 Real vecpot2b1i(Real (*A2)(Real,Real,Real), Real (*A3)(Real,Real,Real),
00832                 const GridS *pG, const int i, const int j, const int k)
00833 {
00834   Real x1,x2,x3,b1i=0.0,lsf=1.0,rsf=1.0,dx2=pG->dx2;
00835   Real f2(Real y);
00836   Real f3(Real z);
00837 
00838   a2func = A2;
00839   a3func = A3;
00840   cc_pos(pG,i,j,k,&x1,&x2,&x3);
00841   xmin = x1 - 0.5*pG->dx1;  xmax = x1 + 0.5*pG->dx1;
00842   ymin = x2 - 0.5*pG->dx2;  ymax = x2 + 0.5*pG->dx2;
00843   zmin = x3 - 0.5*pG->dx3;  zmax = x3 + 0.5*pG->dx3;
00844 
00845   xsav = xmin;
00846 #ifdef CYLINDRICAL
00847   lsf = xmin;  rsf = xmin;
00848   dx2 = xmin*pG->dx2;
00849 #endif
00850 
00851   if (A2 != NULL) {
00852     if (ymin == ymax)
00853       b1i += rsf*A2(xmin,ymin,zmin) - lsf*A2(xmin,ymin,zmax);
00854     else {
00855       zsav = zmin;
00856       b1i += rsf*qsimp(f2,ymin,ymax);
00857       zsav = zmax;
00858       b1i -= lsf*qsimp(f2,ymin,ymax);
00859     }
00860   }
00861   if (A3 != NULL) {
00862     if (zmin == zmax)
00863       b1i += A3(xmin,ymax,zmin) - A3(xmin,ymin,zmin);
00864     else {
00865       ysav = ymax;
00866       b1i += qsimp(f3,zmin,zmax);
00867       ysav = ymin;
00868       b1i -= qsimp(f3,zmin,zmax);
00869     }
00870   }
00871 
00872   if (pG->dx2 > 0.0) b1i /= dx2;
00873   if (pG->dx3 > 0.0) b1i /= pG->dx3;
00874 
00875   return b1i;
00876 }
00877 
00878 /*! \fn Real vecpot2b2i(Real (*A1)(Real,Real,Real), Real (*A3)(Real,Real,Real),
00879  *                const GridS *pG, const int i, const int j, const int k)
00880  *  \brief Compute B-field components from a vector potential.
00881  *
00882  * THESE FUNCTIONS COMPUTE MAGNETIC FIELD COMPONENTS FROM COMPONENTS OF A
00883  * SPECIFIED VECTOR POTENTIAL USING STOKES' THEOREM AND SIMPSON'S QUADRATURE.
00884  * NOTE:  THIS IS ONLY GUARANTEED TO WORK IF THE POTENTIAL IS OF CLASS C^1.
00885  * WRITTEN BY AARON SKINNER.
00886  */
00887 Real vecpot2b2i(Real (*A1)(Real,Real,Real), Real (*A3)(Real,Real,Real),
00888                 const GridS *pG, const int i, const int j, const int k)
00889 {
00890   Real x1,x2,x3,b2i=0.0;
00891   Real f1(Real x);
00892   Real f3(Real z);
00893 
00894   a1func = A1;
00895   a3func = A3;
00896   cc_pos(pG,i,j,k,&x1,&x2,&x3);
00897   xmin = x1 - 0.5*pG->dx1;  xmax = x1 + 0.5*pG->dx1;
00898   ymin = x2 - 0.5*pG->dx2;  ymax = x2 + 0.5*pG->dx2;
00899   zmin = x3 - 0.5*pG->dx3;  zmax = x3 + 0.5*pG->dx3;
00900 
00901   ysav = ymin;
00902 
00903   if (A1 != NULL) {
00904     if (xmin == xmax)
00905       b2i += A1(xmin,ymin,zmax) - A1(xmin,ymin,zmin);
00906     else {
00907       zsav = zmax;
00908       b2i += qsimp(f1,xmin,xmax);
00909       zsav = zmin;
00910       b2i -= qsimp(f1,xmin,xmax);
00911     }
00912   }
00913   if (A3 != NULL) {
00914     if (zmin == zmax)
00915       b2i += A3(xmin,ymin,zmin) - A3(xmax,ymin,zmin);
00916     else {
00917       xsav = xmin;
00918       b2i += qsimp(f3,zmin,zmax);
00919       xsav = xmax;
00920       b2i -= qsimp(f3,zmin,zmax);
00921     }
00922   }
00923 
00924   if (pG->dx1 > 0.0) b2i /= pG->dx1;
00925   if (pG->dx3 > 0.0) b2i /= pG->dx3;
00926 
00927   return b2i;
00928 }
00929 
00930 /*! \fn Real vecpot2b3i(Real (*A1)(Real,Real,Real), Real (*A2)(Real,Real,Real),
00931  *                const GridS *pG, const int i, const int j, const int k)
00932  *  \brief Compute B-field components from a vector potential.
00933  *
00934  * THESE FUNCTIONS COMPUTE MAGNETIC FIELD COMPONENTS FROM COMPONENTS OF A
00935  * SPECIFIED VECTOR POTENTIAL USING STOKES' THEOREM AND SIMPSON'S QUADRATURE.
00936  * NOTE:  THIS IS ONLY GUARANTEED TO WORK IF THE POTENTIAL IS OF CLASS C^1.
00937  * WRITTEN BY AARON SKINNER.
00938  */
00939 Real vecpot2b3i(Real (*A1)(Real,Real,Real), Real (*A2)(Real,Real,Real),
00940                 const GridS *pG, const int i, const int j, const int k)
00941 {
00942   Real x1,x2,x3,b3i=0.0,lsf=1.0,rsf=1.0,dx2=pG->dx2;
00943   Real f1(Real x);
00944   Real f2(Real y);
00945 
00946   a1func = A1;
00947   a2func = A2;
00948   cc_pos(pG,i,j,k,&x1,&x2,&x3);
00949   xmin = x1 - 0.5*pG->dx1;  xmax = x1 + 0.5*pG->dx1;
00950   ymin = x2 - 0.5*pG->dx2;  ymax = x2 + 0.5*pG->dx2;
00951   zmin = x3 - 0.5*pG->dx3;  zmax = x3 + 0.5*pG->dx3;
00952 
00953   zsav = zmin;
00954 #ifdef CYLINDRICAL
00955   rsf = xmax;  lsf = xmin;
00956   dx2 = x1*pG->dx2;
00957 #endif
00958 
00959   if (A1 != NULL) {
00960     if (xmin == xmax)
00961       b3i += A1(xmin,ymin,zmin) - A1(xmin,ymax,zmin);
00962     else {
00963       ysav = ymin;
00964       b3i += qsimp(f1,xmin,xmax);
00965       ysav = ymax;
00966       b3i -= qsimp(f1,xmin,xmax);
00967     }
00968   }
00969   if (A2 != NULL) {
00970     if (ymin == ymax)
00971       b3i += rsf*A2(xmax,ymin,zmin) - lsf*A2(xmin,ymin,zmin);
00972     else {
00973       xsav = xmax;
00974       b3i += rsf*qsimp(f2,ymin,ymax);
00975       xsav = xmin;
00976       b3i -= lsf*qsimp(f2,ymin,ymax);
00977     }
00978   }
00979 
00980   if (pG->dx1 > 0.0) b3i /= pG->dx1;
00981   if (pG->dx2 > 0.0) b3i /= dx2;
00982 
00983   return b3i;
00984 }
00985 
00986 Real f1(Real x)
00987 {
00988   return a1func(x,ysav,zsav);
00989 }
00990 
00991 Real f2(Real y)
00992 {
00993   return a2func(xsav,y,zsav);
00994 }
00995 
00996 Real f3(Real z)
00997 {
00998   return a3func(xsav,ysav,z);
00999 }
01000 
01001 
01002 #if defined(PARTICLES) || defined(CHEMISTRY)
01003 /*----------------------------------------------------------------------------*/
01004 /*! \fn void ludcmp(Real **a, int n, int *indx, Real *d)
01005  *  \brief LU decomposition from Numerical Recipes
01006  *
01007  * Using Crout's method with partial pivoting
01008  * a is the input matrix, and is returned with LU decomposition readily made,
01009  * n is the matrix size, indx records the history of row permutation,
01010  * whereas d =1(-1) for even(odd) number of permutations.
01011  */
01012 void ludcmp(Real **a, int n, int *indx, Real *d)
01013 {
01014   int i,imax,j,k;
01015   Real big,dum,sum,temp;
01016   Real *rowscale;  /* the implicit scaling of each row */
01017 
01018   rowscale = (Real*)calloc_1d_array(n, sizeof(Real));
01019   *d=1.0;  /* No row interchanges yet */
01020 
01021   for (i=0;i<n;i++)
01022   { /* Loop over rows to get the implicit scaling information */
01023     big=0.0;
01024     for (j=0;j<n;j++)
01025       if ((temp=fabs(a[i][j])) > big) big=temp;
01026     if (big == 0.0) ath_error("[LUdecomp]:Input matrix is singular!");
01027     rowscale[i]=1.0/big;  /* Save the scaling */
01028   }
01029 
01030   for (j=0;j<n;j++) { /* Loop over columns of Crout's method */
01031     /* Calculate the upper block */
01032     for (i=0;i<j;i++) {
01033       sum=a[i][j];
01034       for (k=0;k<i;k++) sum -= a[i][k]*a[k][j];
01035       a[i][j]=sum;
01036     }
01037     /* Calculate the lower block (first step) */
01038     big=0.0;
01039     for (i=j;i<n;i++) {
01040       sum=a[i][j];
01041       for (k=0;k<j;k++)
01042         sum -= a[i][k]*a[k][j];
01043       a[i][j]=sum;
01044       /* search for the largest pivot element */
01045       if ( (dum=rowscale[i]*fabs(sum)) >= big) {
01046         big=dum;
01047         imax=i;
01048       }
01049     }
01050     /* row interchange */
01051     if (j != imax) {
01052       for (k=0;k<n;k++) {
01053         dum=a[imax][k];
01054         a[imax][k]=a[j][k];
01055         a[j][k]=dum;
01056       }
01057       *d = -(*d);
01058       rowscale[imax]=rowscale[j];
01059     }
01060     indx[j]=imax; /* record row interchange history */
01061     /* Calculate the lower block (second step) */
01062     if (a[j][j] == 0.0) a[j][j]=TINY_NUMBER;
01063     dum=1.0/(a[j][j]);
01064     for (i=j+1;i<n;i++) a[i][j] *= dum;
01065   }
01066   free(rowscale);
01067 }
01068 
01069 /*----------------------------------------------------------------------------*/
01070 /*! \fn void lubksb(Real **a, int n, int *indx, Real b[])
01071  *  \brief Backward substitution (from numerical recipies)
01072  *
01073  * a is the input matrix done with LU decomposition, n is the matrix size
01074  * indx id the history of row permutation
01075  * b is the vector on the right (AX=b), and is returned with the solution
01076  */
01077 void lubksb(Real **a, int n, int *indx, Real b[])
01078 {
01079   int i,ii=-1,ip,j;
01080   Real sum;
01081   /* Solve L*y=b */
01082   for (i=0;i<n;i++) {
01083     ip=indx[i];
01084     sum=b[ip];
01085     b[ip]=b[i];
01086     if (ii>=0)
01087       for (j=ii;j<=i-1;j++) sum -= a[i][j]*b[j];
01088     else if (sum) ii=i;
01089     b[i]=sum;
01090   }
01091   /* Solve U*x=y */
01092   for (i=n-1;i>=0;i--) {
01093     sum=b[i];
01094     for (j=i+1;j<n;j++) sum -= a[i][j]*b[j];
01095     b[i]=sum/a[i][i];
01096   }
01097 }
01098 
01099 /*----------------------------------------------------------------------------*/
01100 /*! \fn void InverseMatrix(Real **a, int n, Real **b)
01101  *  \brief Inverse matrix solver
01102  *
01103  * a: input matrix; n: matrix size, b: return matrix
01104  * Note: the input matrix will be DESTROYED
01105  */
01106 void InverseMatrix(Real **a, int n, Real **b)
01107 {
01108   int i,j,*indx;
01109   Real *col,d;
01110 
01111   indx = (int*)calloc_1d_array(n, sizeof(int));
01112   col = (Real*)calloc_1d_array(n, sizeof(Real));
01113 
01114   ludcmp(a,n,indx,&d);
01115 
01116   for (j=0; j<n; j++) {
01117     for (i=0; i<n; i++) col[i]=0.0;
01118     col[j]=1.0;
01119     lubksb(a, n, indx, col);
01120     for (i=0; i<n; i++)    b[i][j] = col[i];
01121   }
01122 
01123   return;
01124 }
01125 
01126 /*----------------------------------------------------------------------------*/
01127 /*! \fn void MatrixMult(Real **a, Real **b, int m, int n, int l, Real **c)
01128  *  \brief Matrix multiplication: a(m*n) * b(n*l) = c(m*l) */
01129 void MatrixMult(Real **a, Real **b, int m, int n, int l, Real **c)
01130 {
01131   int i, j, k;
01132   for (i=0; i<m; i++)
01133     for (j=0; j<l; j++)
01134     {
01135       c[i][j] = 0.0;
01136       for (k=0; k<n; k++) c[i][j] += a[i][k] * b[k][j];
01137     }
01138 }
01139 #endif /* PARTICLES or CHEMISTRY */