<< dt_dia << ",\n  and snapshot output interval dt_out = "	   << dt_out << "." << endl;    real tcpu = MPI_Wtime();           // check CPU usage    real epot = 0;                     // potential energy of the n-body system    real coll_time = VERY_LARGE_NUMBER;// collision (close encounter) time scale    get_acc_jerk_pot_coll(p, n, epot, coll_time, pipe);    int nsteps = 0;               // number of integration time steps completed    real einit;                   // initial total energy of the system    write_diagnostics(p, n, t, epot, nsteps, einit,		      true, x_flag, tcpu);    if (init_out)                                    // flag for initial output      put_snapshot(p, n, t, particletype);    real t_dia = t + dt_dia;           // next time for diagnostics output    real t_out = t + dt_out;           // next time for snapshot output    real t_end = t + dt_tot;           // final time, to finish the integration    while (true){        while (t < t_dia && t < t_out && t < t_end){            real dt_local = dt_param * coll_time;	    real dt;	    MPI_Allreduce(&dt_local, &dt, 1, MPI_DOUBLE, MPI_MIN, MPI_COMM_WORLD); // synchronize time step	    evolve_step(p, n, t, dt, epot, coll_time, pipe);            nsteps++;        }        if (t >= t_dia){            write_diagnostics(p, n, t, epot, nsteps,                              einit, false, x_flag, tcpu);            t_dia += dt_dia;        }        if (t >= t_out){            put_snapshot(p, n, t, particletype);            t_out += dt_out;        }        if (t >= t_end)            break;    }}/*----------------------------------------------------------------------------- *  evolve_step  --  takes one integration step for an N-body system, using the *                   Hermite algorithm. *----------------------------------------------------------------------------- */void evolve_step(Particle p[], int n, real & t,                 real dt, real & epot, real & coll_time,		 void *pipe){  Particle *po = new Particle[n];    for (int i = 0; i < n ; i++)        for (int k = 0; k < NDIM ; k++){          po[i].pos[k] = p[i].pos[k];          po[i].vel[k] = p[i].vel[k];          po[i].acc[k] = p[i].acc[k];          po[i].jerk[k] = p[i].jerk[k];        }    predict_step(p, n, dt);    get_acc_jerk_pot_coll(p, n, epot, coll_time, pipe);    correct_step(p, po, n, dt);    t += dt;    delete[] po;}/*----------------------------------------------------------------------------- *  predict_step  --  takes the first approximation of one Hermite integration *                    step, advancing the positions and velocities through a *                    Taylor series development up to the order of the jerks. *----------------------------------------------------------------------------- */void predict_step(Particle p[], int n, real dt){    for (int i = 0; i < n ; i++)        for (int k = 0; k < NDIM ; k++){            p[i].pos[k] += p[i].vel[k]*dt + p[i].acc[k]*dt*dt/2                                      + p[i].jerk[k]*dt*dt*dt/6;            p[i].vel[k] += p[i].acc[k]*dt + p[i].jerk[k]*dt*dt/2;        }}/*----------------------------------------------------------------------------- *  correct_step  --  takes one iteration to improve the new values of position *                    and velocities, effectively by using a higher-order *                    Taylor series constructed from the terms up to jerk at *                    the beginning and the end of the time step. *----------------------------------------------------------------------------- */void correct_step(Particle p[], Particle po[], int n, real dt){    for (int i = 0; i < n ; i++)        for (int k = 0; k < NDIM ; k++){            p[i].vel[k] = po[i].vel[k] + (po[i].acc[k] + p[i].acc[k])*dt/2                                       + (po[i].jerk[k] - p[i].jerk[k])*dt*dt/12;            p[i].pos[k] = po[i].pos[k] + (po[i].vel[k] + p[i].vel[k])*dt/2                                       + (po[i].acc[k] - p[i].acc[k])*dt*dt/12;        }}/*----------------------------------------------------------------------------- *  get_acc_jerk_pot_coll  --  calculates accelerations and jerks, and as side *                             effects also calculates potential energy and *                             the time scale coll_time for significant changes *                             in local configurations to occur. *                                                  __                     __ *                                                 |          -->  -->       | *               M                           M     |           r  . v        | *   -->          j    -->       -->          j    | -->        ji   ji -->  | *    a   ==  --------  r    ;    j   ==  -------- |  v   - 3 ---------  r   | *     ji     |-->  |3   ji        ji     |-->  |3 |   ji      |-->  |2   ji | *            | r   |                     | r   |  |           | r   |       | *            |  ji |                     |  ji |  |__         |  ji |     __| *                              *  note: it would be cleaner to calculate potential energy and collision time *        in a separate function.  However, the current function is by far the *        most time consuming part of the whole program, with a double loop *        over all particles that is executed every time step.  Splitting off *        some of the work to another function would significantly increase *        the total computer time (by an amount close to a factor two). * *  We determine the values of all four quantities of interest by walking *  through the system in a double {i,j} loop.  The first three, acceleration, *  jerk, and potential energy, are calculated by adding successive terms; *  the last, the estimate for the collision time, is found by determining the  *  minimum value over all particle pairs and over the two choices of collision *  time, position/velocity and sqrt(position/acceleration), where position and *  velocity indicate their relative values between the two particles, while *  acceleration indicates their pairwise acceleration.  At the start, the *  first three quantities are set to zero, to prepare for accumulation, while *  the last one is set to a very large number, to prepare for minimization. *       The integration loops only over half of the pairs, with j > i, since *  the contributions to the acceleration and jerk of particle j on particle i *  is the same as those of particle i on particle j, apart from a minus sign *  and a different mass factor. *----------------------------------------------------------------------------- */void get_acc_jerk_pot_coll(Particle pl[], int nl, 			   Particle po[], int no, 			   real & epot, real & coll_time){    real coll_time_q = VERY_LARGE_NUMBER;      // collision time to 4th power    real coll_est_q;                           // collision time scale estimate                                               // to 4th power (quartic)    for (int i = 0; i < nl ; i++){        for (int j = 0; j < no ; j++){            // rji[] is the vector from	  if(pl[i].id!=po[j].id) {            real rji[NDIM];                        // particle i to particle j            real vji[NDIM];                        // vji[] = d rji[] / d t            for (int k = 0; k < NDIM ; k++){                rji[k] = po[j].pos[k] - pl[i].pos[k];                vji[k] = po[j].vel[k] - pl[i].vel[k];            }            real r2 = 0;                           // | rji |^2            real v2 = 0;                           // | vji |^2            real rv_r2 = 0;                        // ( rij . vij ) / | rji |^2            for (int k = 0; k < NDIM ; k++){                r2 += rji[k] * rji[k];                v2 += vji[k] * vji[k];                rv_r2 += rji[k] * vji[k];            }            rv_r2 /= r2;            real r = sqrt(r2);                     // | rji |            real r3 = r * r2;                      // | rji |^3// add the {i,j} contribution to the total potential energy for the system:            epot -= pl[i].mass * po[j].mass / r;// add the {j (i)} contribution to the {i (j)} values of acceleration and jerk:            real da[3];                            // main terms in pairwise            real dj[3];                            // acceleration and jerk            for (int k = 0; k < NDIM ; k++){                da[k] = rji[k] / r3;                           // see equations                dj[k] = (vji[k] - 3 * rv_r2 * rji[k]) / r3;    // in the header            }            for (int k = 0; k < NDIM ; k++){                pl[i].acc[k] += po[j].mass * da[k];           // using symmetry		pl[i].jerk[k] += po[j].mass * dj[k];          // acceleration		// in the original version pij = -pji for acc and jerk.		// in this parallel version this is unpractical.                //po[j].acc[k] -= pl[i].mass * da[k];      // find pairwise		//po[j].jerk[k] -= pl[i].mass * dj[k];     // and jerk            }// first collision time estimate, based on unaccelerated linear motion:            coll_est_q = (r2*r2) / (v2*v2);            if (coll_time_q > coll_est_q)                coll_time_q = coll_est_q;// second collision time estimate, based on free fall:            real da2 = 0;                                  // da2 becomes the             for (int k = 0; k < NDIM ; k++)                // square of the                 da2 += da[k] * da[k];                      // pair-wise accel-            double mij = pl[i].mass + po[j].mass;            // eration between            da2 *= mij * mij;                              // particles i and j            coll_est_q = r2/da2;            if (coll_time_q > coll_est_q)                coll_time_q = coll_est_q;	  }        }                                         }        // from q for quartic back to linear collision time and taking the minimum    coll_time = min(coll_time, sqrt(sqrt(coll_time_q))); }                                             /*----------------------------------------------------------------------------- *  get_acc_jerk_pot_coll  --  the distribution of particles over  *                             neighboring processors and subsequent force *                             calculation. *----------------------------------------------------------------------------- */void get_acc_jerk_pot_coll(Particle p[], int n, real &epot, real &coll_time, 			   void *pipe) {  int rank, size;  MPI_Comm_rank( MPI_COMM_WORLD, &rank );  MPI_Comm_size( MPI_COMM_WORLD, &size );    int rlen;    Particle *recvbuf;    for (int i = 0; i < n ; i++)      for (int k = 0; k < NDIM ; k++) {	p[i].acc[k] = p[i].jerk[k] = 0;      }    MPE_Pipe_start( pipe, p, n, 1 ); // load the initial sendbuffer     epot = 0;                        // initialize epot and coll_time    coll_time = VERY_LARGE_NUMBER;    get_acc_jerk_pot_coll(p, n, p, n, epot, coll_time);    for (int step=1; step<size; step++) { // compute forces for other particles      MPE_Pipe_push( pipe, (void**)&recvbuf, &rlen );        // get new data            // Compute forces       get_acc_jerk_pot_coll(p, n, recvbuf, rlen, epot, coll_time);     }}/*----------------------------------------------------------------------------- *                                                                    \\   o *  end of file:  nbody_sh1.C                                         /\\'  O *                                                                   /\     | *============================================================================= */