From ioannis@ircamera.as.arizona.edu Tue Nov 14 18:06:22 2000 Return-Path: Received: from astro.as.arizona.edu (IDENT:root@astro.as.arizona.edu [128.196.208.2]) by ngala.as.arizona.edu (8.9.3/8.8.7) with ESMTP id SAA05178 for ; Tue, 14 Nov 2000 18:06:22 -0800 Received: from ircamera.as.arizona.edu (ircamera.as.arizona.edu [128.196.211.20]) by astro.as.arizona.edu (8.9.3/8.9.3) with ESMTP id SAA14365 for ; Tue, 14 Nov 2000 18:25:22 -0700 Received: (from ioannis@localhost) by ircamera.as.arizona.edu (8.9.3/8.8.7) id SAA03351 for dzaritsky@as.arizona.edu; Tue, 14 Nov 2000 18:29:28 -0700 Date: Tue, 14 Nov 2000 18:29:28 -0700 From: John Moustakas Message-Id: <200011150129.SAA03351@ircamera.as.arizona.edu> To: dzaritsky@as.arizona.edu Status: RO nbody0.f SUBROUTINE nbody0 C C S.J. Aarseth's Standard N-Body Program (nbody0) C C Adapted from the source code in Binney & Tremaine's (1987) C book 'Galactic Dynamics', C by Peter Teuben - June '89: C - all variables to be declared and a choice of 'real' or C 'double precision' C - all I/O is done through subroutines in order to allow easy C interface with development different systems, e.g. NEMO89. C - PARAMETERS for NMAX and NDIM via an include file (see Makefile) C Note that NDIM should NEVER be reset from 3, the code does C not handle 2D stuff (yet). C C jun-89 Created C 23-jan-90 Back to double precision - PJT C 5-apr-90 started to declare all variables - PJT C 14-nov-91 finished(!) declaring all variables; C tested using IMPLICIT NONE and compile with -u C 15-nov-91 Connection Machine testing in nbody0.fcm file C 11-feb-98 Fixed array index for f2dot(k) a(16)->a(15) C*********************************************************************** INCLUDE 'fdefs.inc' INCLUDE 'nmax.inc' CSED The next statement can be modified with SED to toggle type DOUBLE PRECISION - time,tnext,dt,s,deltat,tcrit,e,eta,eps2, - t1pr,t2pr,t3pr,dt1,dt2,dt3, - x(NDIM,NMAX),x0(NDIM,NMAX),x0dot(NDIM,NMAX),t0(NMAX), - body(NMAX),step(NMAX), - f(NDIM,NMAX),fdot(NDIM,NMAX), - d1(NDIM,NMAX),d2(NDIM,NMAX),d3(NDIM,NMAX), - t1(NMAX),t2(NMAX),t3(NMAX),a(17), - f1(NDIM),f1dot(NDIM),f2dot(NDIM),f3dot(NDIM) INTEGER nsteps, n, i, j, k DATA time,tnext,nsteps /0.0,0.0,0/ C----------------------------------------------------------------------- CALL INPARS(NMAX,n,eta,deltat,tcrit,eps2) CALL INBODS (n, body, x0, x0dot) C obtain total forces and first derivative for each body DO 20 i=1,n DO 2 k=1,NDIM f(k,i)=0.0 fdot(k,i)=0.0 d2(k,i)=0.0 d3(k,i)=0.0 2 CONTINUE DO 10 j=1,n IF (j.EQ.i) GO TO 10 DO 5 k=1,NDIM a(k) = x0(k,j) - x0(k,i) a(k+3) = x0dot(k,j) - x0dot(k,i) 5 CONTINUE a(7) = 1.0/(a(1)**2 + a(2)**2 + a(3)**2 + eps2) a(8) = body(j)*a(7)*sqrt(a(7)) a(9) = 3.0*(a(1)*a(4)+a(2)*a(5)+a(3)*a(6))*a(7) DO 8 k=1,NDIM f(k,i) = f(k,i) + a(k)*a(8) fdot(k,i) = fdot(k,i) + (a(k+3) - a(k)*a(9))*a(8) 8 CONTINUE 10 CONTINUE 20 CONTINUE C form second and third derivative DO 40 i=1,n DO 30 j=1,n IF (i.EQ.j) GO TO 30 DO 25 k=1,NDIM a(k) = x0(k,j) - x0(k,i) a(k+3) = x0dot(k,j) - x0dot(k,i) a(k+6) = f(k,j) - f(k,i) a(k+9) = fdot(k,j) - fdot(k,i) 25 CONTINUE a(13) = 1.0/(a(1)**2 + a(2)**2 + a(3)**2 + eps2) a(14) = body(j)*a(13)*sqrt (a(13)) a(15) = (a(1)*a(4) + a(2)*a(5) + a(3)*a(6))*a(13) a(16) = (a(4)**2 + a(5)**2 + a(6)**2 + - a(1)*a(7) + a(2)*a(8) + a(3)*a(9))*a(13) + a(15)**2 a(17) = (9.0*(a(4)*a(7) + a(5)*a(8) + a(6)*a(9)) + - 3.0*(a(1)*a(10) + a(2)*a(11) + a(3)*a(12)))*a(13) + - a(15)*(9.0*a(16) - 12.0*a(15)**2) DO 28 k=1,NDIM f1dot(k) = a(k+3) - 3.0*a(15)*a(k) f2dot(k) = (a(k+6)-6*a(15)*f1dot(k)-3*a(16)*a(k))*a(14) f3dot(k) = (a(k+9)-9*a(16)*f1dot(k)-a(17)*a(k))*a(14) d2(k,i) = d2(k,i) + f2dot(k) d3(k,i) = d3(k,i) + f3dot(k) - 9.0*a(15)*f2dot(k) 28 CONTINUE 30 CONTINUE 40 CONTINUE C initialize integration steps and convert to force differences DO 50 i=1,n step(i) = sqrt(eta*sqrt ((f(1,i)**2 + f(2,i)**2 + f(3,i)**2)/ - (d2(1,i)**2 + d2(2,i)**2 + d2(3,i)**2))) t0(i) = time t1(i) = time - step(i) t2(i) = time - 2.0*step(i) t3(i) = time - 3.0*step(i) DO 45 k=1,NDIM d1(k,i) = (d3(k,i)*step(i)/6-0.5*d2(k,i))*step(i)+fdot(k,i) d2(k,i) = 0.5*d2(k,i) - 0.5*d3(k,i)*step(i) d3(k,i) = d3(k,i)/6 f(k,i) = 0.5*f(k,i) fdot(k,i) = fdot(k,i)/6 45 CONTINUE 50 CONTINUE C----------------------------------------------------------------------- C==============> entry point when major output + diagnostics to be done C C energy check and ouput 100 CONTINUE e = 0.0 DO 110 i=1,n dt = tnext - t0(i) DO 101 k=1,NDIM f2dot(k) = d3(k,i)*((t0(i)-t1(i))+(t0(i)-t2(i))) + d2(k,i) x(k,i) = ((((0.05*d3(k,i)*dt + f2dot(k)/12.0)*dt + - fdot(k,i))*dt + f(k,i))*dt + x0dot(k,i))*dt - + x0(k,i) a(k) = (((0.25*d3(k,i)*dt + f2dot(k)/3.0)*dt + - 3.0*fdot(k,i))*dt + 2.0*f(k,i))*dt + x0dot(k,i) 101 CONTINUE CALL OUTBODS (body(i),x(1,i),a,step(i),i) e = e + 0.5*body(i)*(a(1)**2 + a(2)**2 + a(3)**2) 110 CONTINUE DO 130 i=1,n DO 120 j=1,n IF (j.EQ.i) GO TO 120 e = e - 0.5*body(i)*body(j)/sqrt ((x(1,i) - x(1,j))**2 + - (x(2,i) - x(2,j))**2 + (x(3,i) - x(3,j))**2 + eps2) 120 CONTINUE 130 CONTINUE CALL OUTENE(tnext,nsteps,e) IF (time.GT.tcrit) RETURN tnext = tnext + deltat C----------------------------------------------------------------------- C ===============> Normal entry point for next timstep C C C find next body to be advanced and set new time 200 CONTINUE time = 1.0e+10 DO 210 j=1,n IF (time.GT.t0(j)+step(j)) i=j IF (time.GT.t0(j)+step(j)) time = t0(j) + step(j) 210 CONTINUE C predict all coordinates to first order in force derivative DO 220 j=1,n s = time - t0(j) x(1,j) = ((fdot(1,j)*s + f(1,j))*s + x0dot(1,j))*s + x0(1,j) x(2,j) = ((fdot(2,j)*s + f(2,j))*s + x0dot(2,j))*s + x0(2,j) x(3,j) = ((fdot(3,j)*s + f(3,j))*s + x0dot(3,j))*s + x0(3,j) 220 CONTINUE C include 2nd and 3rd order and obtain the velocity dt = time - t0(i) DO 230 k=1,NDIM f2dot(k) = d3(k,i)*((t0(i)-t1(i)) + (t0(i)-t2(i))) + d2(k,i) x(k,i) = (0.05*d3(k,i)*dt + f2dot(k)/12.0)*dt**4 + x(k,i) x0dot(k,i) = (((0.25*d3(k,i)*dt + f2dot(k)/3.0)*dt + - 3.0*fdot(k,i))*dt + 2.0*f(k,i))*dt + x0dot(k,i) f1(k) = 0.0 230 CONTINUE C obtain current force on i-th body DO 240 j=1,n IF (j.EQ.i) GO TO 240 a(1) = x(1,j) - x(1,i) a(2) = x(2,j) - x(2,i) a(3) = x(3,j) - x(3,i) a(4) = 1.0/(a(1)**2 + a(2)**2 + a(3)**2 + eps2) a(5) = body(j)*a(4)*sqrt(a(4)) f1(1) = f1(1) + a(1)*a(5) f1(2) = f1(2) + a(2)*a(5) f1(3) = f1(3) + a(3)*a(5) 240 CONTINUE C set time intervals for new difference and update the times dt1 = time - t1(i) dt2 = time - t2(i) dt3 = time - t3(i) t1pr = t0(i) - t1(i) t2pr = t0(i) - t2(i) t3pr = t0(i) - t3(i) t3(i) = t2(i) t2(i) = t1(i) t1(i) = t0(i) t0(i) = time C form new differences and include 4th order semi-iterative DO 250 k=1,NDIM a(k) = (f1(k) - 2.0*f(k,i))/dt a(k+3) = (a(k) - d1(k,i))/dt1 a(k+6) = (a(k+3) - d2(k,i))/dt2 a(k+9) = (a(k+6) - d3(k,i))/dt3 d1(k,i) = a(k) d2(k,i) = a(k+3) d3(k,i) = a(k+6) f1dot(k) = t1pr*t2pr*t3pr*a(k+9) f2dot(k) = (t1pr*t2pr + t3pr*(t1pr+t2pr))*a(k+9) f3dot(k) = (t1pr + t2pr + t3pr)*a(k+9) x0(k,i) = (((a(k+9)*dt/30.0 + 0.05*f3dot(k))*dt + - f2dot(k)/12.0)*dt + f1dot(k)/6.0)*dt**3 + x(k,i) x0dot(k,i) = (((0.2*a(k+9)*dt + 0.25*f3dot(k))*dt + - f2dot(k)/3.0)*dt + 0.5*f1dot(k))*dt**2 + x0dot(k,i) 250 CONTINUE C scale F and FDOT by factorials and set new integration step DO 260 k=1,NDIM f(k,i) = 0.5*f1(k) fdot(k,i) = ((d3(k,i)*dt1 + d2(k,i))*dt + d1(k,i))/6.0 f2dot(k)= 2.0*(d3(k,i)*(dt+dt1) + d2(k,i)) 260 CONTINUE step(i) = sqrt(eta*sqrt((f1(1)**2 + f1(2)**2 + f1(3)**2)/ - (f2dot(1)**2 + f2dot(2)**2 + f2dot(3)**2))) nsteps = nsteps + 1 IF (time - tnext) 200,100,100 END * NBODYIO.C C-version I/O for Fortran-C interface of Aarseth's nbody0.f * Has a crude F2C interface * 22-nov-91 fixed aix interface */ #include #include #include #include #include #include #include #include local double l_eta; /* local copies of things for Aarseth */ local double l_deltat; /* now double, used to be float */ local double l_tcrit; local double l_eps2; local int nbody; /* actual number of bodies */ local int bits; /* bitmask what's in snapshot */ local bool Qstep; local stream instr; /* - pointer to input file */ local stream outstr; /* - pointer to output file */ local Body *btab = NULL; /* - pointer to body structure */ extern double cputime(void); /*============================================================================*/ /* Fortran-to-C interface stuff: * Since we only handle simple parameters here (i.e. no character * variables), the only thing to be done is to properly define the name * of the C routine, because it will be called as a fortran subroutine. * General rule: avoid character variables and life is fairly simple in F2C. * * In VMS the C names stay the same * in UNICOS the C name is to be in UPPER CASE. * In most (BSD) unix versions C name gets an underscore (_) appended */ #if !defined(VMS) && !defined(aix) # if defined(unicos) # define nbody0 NBODY0 # define inbods INBODS # define inpars INPARS # define outbods OUTBODS # define outene OUTENE # else # define nbody0 nbody0_ # define inbods inbods_ # define inpars inpars_ # define outbods outbods_ # define outene outene_ # endif #endif /*============================================================================*/ /* * pars_2_aarseth: kludge to prepare transport parameters to fortran * and start reading the input snapshot file * */ pars_2_aarseth() { double time; l_eta = (double) getdparam("eta"); l_deltat = (double) getdparam("deltat"); l_tcrit = (double) getdparam("tcrit"); l_eps2 = (double) sqr(getdparam("eps")); if (hasvalue("options")) { Qstep = TRUE; } else Qstep = FALSE; if (NDIM != 3) error("%d: Code not supported for non 3D",NDIM); instr = stropen(getparam("in"),"r"); /* open input file */ get_history(instr); /* get history from input file */ get_snap (instr, &btab, &nbody, &time, &bits); /* read snapshot in */ if ((bits & TimeBit)==0) /* set time to zero if not present */ time=0.0; if ((bits & MassBit)==0) /* check if essentials are present */ error ("No masses in input snapshot"); if ((bits & PhaseSpaceBit)==0) error ("No phase space coordinates in input snapshot"); if (time >= (double) l_tcrit) error("Input snapshot has larger timestamp (%g) than required tcrit(%g)", time, l_tcrit); outstr = stropen(getparam("out"),"w"); /* open output file */ put_history(outstr); /* copy history to output file */ bits = (TimeBit | MassBit | PhaseSpaceBit); /* set things to copy */ if (Qstep) bits |= AuxBit; } /* * call_2_aarseth: the actual routine which does the fortran calling */ call_2_aarseth() { nbody0(); /* call fortran hard core worker routine */ } /* * (fortran) INPARS: read global integration parameters */ inpars (nmax, n, eta, deltat, tcrit, eps2) int *nmax, *n; double *eta, *deltat, *tcrit, *eps2; { if (nbody > *nmax) error("Too many particles (%d) in inputfile, recompile (%d) program", nbody, *nmax); *n = nbody; *eta = l_eta; *deltat = l_deltat; *tcrit = l_tcrit; *eps2 = l_eps2; } /* * (fortran) INBODS: read particle masses and phases into array */ inbods (n, m, x, v) int *n; double *m, *x, *v; { Body *bp; int i; /* WARNING: what about double and DOUBLE here ? */ for (bp=btab, i=0; i<*n; bp++, i++) { *m++ = Mass(bp); /* copy mass & increase pointer */ SETV(x,Pos(bp)); /* position */ SETV(v,Vel(bp)); /* and velocity */ x += NDIM; v += NDIM; /* increase array pointers */ } } /* * (fortran) OUTBODS: outputs characteristics of particle i (i=1,...nbody) * * In C version just prepares data for flushing later on */ outbods (m, x, v, step, i) double *m, *x, *v, *step; int *i; { Body *bp; if (*i < 0 || *i > nbody) error("OUTBODS: counting error particle\n"); bp = (btab + (*i - 1)); /* point to proper position in Body */ Mass(bp) = *m; SETV(Pos(bp), x); SETV(Vel(bp), v); if (Qstep) Aux(bp) = *step; } /* * (fortran) OUTENE: output energy * * In C: also flushes particles buffer to file */ outene (tnext, nsteps, e) double *tnext, *e; int *nsteps; { double time; printf("time = %g steps = %d energy = %g cpu = %10.3g min\n", *tnext, *nsteps, *e, cputime()); time = (double) *tnext; put_snap (outstr, &btab, &nbody, &time, &bits); } C C Input: C nmax max number particles in workspace C Output: (The input integration parameters) C n number of bodies to be read C eta accuracy parameter C deltat output time interval C tcrit length of integration C eps2 square of softening length C INCLUDE 'fdefs.inc' INTEGER nmax,n CSED The next statement can be modified with SED to toggle type DOUBLE PRECISION eta,deltat,tcrit,eps2 WRITE(6,*) 'Enter n,eta,deltat,tcrit,eps2:' READ(5,*) n,eta,deltat,tcrit,eps2 IF (n.LE.nmax) RETURN WRITE(6,*) 'Recompile program with larger value of NMAX' WRITE(6,*) 'Current value of NMAX = ',nmax WRITE(6,*) 'whereas number of bodies to be read, N = ',n STOP END C*********************************************************************** SUBROUTINE INBODS (n, body, x0, x0dot) C input of initial conditions for integrations INCLUDE 'fdefs.inc' INCLUDE 'nmax.inc' INTEGER i,k,n CSED The next statement can be modified with SED to toggle type DOUBLE PRECISION body(*), x0(NDIM,*), x0dot(NDIM,*) DO i=1,n READ(5,*) body(i),(x0(k,i),k=1,NDIM),(x0dot(k,i),k=1,NDIM) ENDDO RETURN END C*********************************************************************** SUBROUTINE OUTBODS (body, x, a, step, i) C output particle stuff INCLUDE 'fdefs.inc' INCLUDE 'nmax.inc' CSED The next statement can be modified with SED to toggle type DOUBLE PRECISION body, x(NDIM), a(NDIM), step INTEGER i,k WRITE(6,105) body,(x(k),k=1,NDIM),(a(k),k=1,NDIM),step, i 105 FORMAT(1h ,f10.2,3x,3f10.2,3x,3f10.2,3x,f12.4,3x,i5) RETURN END C*********************************************************************** SUBROUTINE OUTENE (tnext, nsteps, e) C C Output current timestep, number of steps performed so far and C total energy of system C INCLUDE 'fdefs.inc' CSED The next statement can be modified with SED to toggle type DOUBLE PRECISION tnext, e INTEGER nsteps WRITE(6,140) tnext,nsteps,e 140 FORMAT(1h0,5x,'time =',f7.2,' steps =',i6,' energy =',f10.4,/) RETURN END * MAIN.C : NEMO driver for Aarseth's fortran nbody0 program * * 3-jul-89 V1.1 * 23-jan-90 V1.2 ??? PJT * 17-may-91 V1.2a bit more doc PJT * 6-mar-92 Happy GCC2.0 PJT * 6-aug-97 V1.3 options= PJT * 11-feb-98 V1.3a fixed that bad index bug in nbody0.f/c PJT */ #include string defv[] = { "in=???\n Input (snapshot) file", "out=???\n Output (snapshot) file", "eta=0.02\n Accuracy parameter - determines timesteps", "deltat=0.25\n When to dump major output", "tcrit=2\n When to stop integrating", "eps=0.05\n Softening length", "options=\n Optional output of 'step' into AUX", "VERSION=1.3a\n 11-feb-98 PJT", NULL, }; string usage = "NEMO driver for nbody0"; extern void pars_2_aarseth(void), call_2_aarseth(void); /* * nemo_main is the main() type entry for the program. The functions * called here live in nbodyio.c, which serves as a C-to-Fortran * layer. */ void nemo_main() { pars_2_aarseth(); /* get parameters */ call_2_aarseth(); /* run the code */ } C Fortran driver for Aarseth's NBODY0 program C By picking the fortran MAIN, as opposed to the C C nemo_main(), you can run NBODY0 withouth NEMO's C I/O routines. You need to the fortran I/O stuff C too (nbodyio.f) C C 1-jul-89 created? PJT C 17-may-91 documented PJT C PROGRAM MAIN0 CALL NBODY0 END C IMPLICIT NONE C see Makefile for instructions INTEGER NMAX, NDIM PARAMETER ( NMAX=2048 ) C do not changed the NDIM value, nbody0.f does not support it PARAMETER ( NDIM=3 ) C end of include file * * S.J. Aarseth's Standard N-Body Program (nbody0) * * Adapted from the source code in Binney & Tremaine's (1987) * book 'Galactic Dynamics', * by Peter Teuben - June '89: * - all variables to be declared and a choice of 'real' or * 'double precision' * - all I/O is done through subroutines in order to allow easy * interface with development different systems, e.g. NEMO89. * - PARAMETERS for NMAX and NDIM via an include file (see Makefile) * Note that NDIM should NEVER be reset from 3, the code does * not handle 2D stuff (yet). * * jun-89 Created * 23-jan-90 Back to double precision - PJT * 5-apr-90 started to declare all variables - PJT * 14-nov-91 finished(!) declaring all variables; * tested using IMPLICIT NONE and compile with -u * 15-nov-91 Connection Machine testing in nbody0.fcm file * 10-may-92 f2c conversion - manually optimized * body index (i,j) are 0 based * dim index (k) is 0 based * 11-feb-98 Fixed array index for f2dot(k) a(16)->a(15) */ #define NMAX 512 #define NDIM 3 /* extract X, Y and Z components of a 3-vector */ #define X(a,i) (a[i*NDIM]) #define Y(a,i) (a[i*NDIM+1]) #define Z(a,i) (a[i*NDIM+2]) #define SQR(x) ((x)*(x)) /* Table of constant values */ static int c_nmax = NMAX; /* *********************************************************************** */ int nbody0_() { /* Initialized data */ double time = 0.0, tnext = 0.0; int nsteps = 0; double time0, time1=0.0; /* counters for my timing/debugging */ /* System (f2c) generated locals */ double d__1, d__2, d__3, d__4, d__5, d__6; /* Local variables */ double fdot[NDIM*NMAX], body[NMAX], step[NMAX], f1dot[NDIM], f2dot[NDIM], f3dot[NDIM], x0dot[NDIM*NMAX], a[17], f[NDIM*NMAX], e, s, x[NDIM*NMAX], tcrit, d1[NDIM*NMAX], d2[NDIM*NMAX], d3[NDIM*NMAX], f1[NDIM], t0[NMAX], t1[NMAX], t2[NMAX], t3[NMAX], x0[NDIM*NMAX], dt, deltat, dt1, dt2, dt3, eta, eps2, t1pr, t2pr, t3pr; int i, j, k, n, itmp; /* Builtin functions */ double sqrt(); /* External subroutines/functions */ extern int outbods_(), inbods_(), inpars_(), outene_(); extern double cputime(); inpars_(&c_nmax, &n, &eta, &deltat, &tcrit, &eps2); inbods_(&n, body, x0, x0dot); /* obtain total forces and first derivative for each body */ for (i = 0; i < n; ++i) { for (k = 0; k < NDIM; ++k) { f[k+i*3] = 0.0; fdot[k+i*3] = 0.0; d2[k+i*3] = 0.0; d3[k+i*3] = 0.0; } for (j = 0; j < n; ++j) { if (j == i) continue; for (k = 0; k < NDIM; ++k) { a[k ] = x0[k+j*3] - x0[k+i*3]; a[k + 3] = x0dot[k+j*3] - x0dot[k+i*3]; } d__1 = a[0]; d__2 = a[1]; d__3 = a[2]; a[6] = 1.0 / (SQR(d__1) + SQR(d__2) + SQR(d__3) + eps2); a[7] = body[j] * a[6] * sqrt(a[6]); a[8] = (a[0]*a[3] + a[1]*a[4] + a[2]*a[5]) * 3.0 * a[6]; for (k = 0; k < NDIM; ++k) { f[k+i*3] += a[k] * a[7]; fdot[k+i*3] += (a[k + 3] - a[k] * a[8]) * a[7]; } } } /* for(i) */ /* form second and third derivative */ for (i = 0; i < n; ++i) { for (j = 0; j < n; ++j) { if (i == j) continue; for (k = 0; k < NDIM; ++k) { a[k ] = x0[k+j*3] - x0[k+i*3]; a[k + 3] = x0dot[k+j*3] - x0dot[k+i*3]; a[k + 6] = f[k+j*3] - f[k+i*3]; a[k + 9] = fdot[k+j*3] - fdot[k+i*3]; } d__1 = a[0]; d__2 = a[1]; d__3 = a[2]; a[12] = 1.0 / (SQR(d__1) + SQR(d__2) + SQR(d__3) + eps2); a[13] = body[j] * a[12] * sqrt(a[12]); a[14] = (a[0]*a[3] + a[1]*a[4] + a[2]*a[5]) * a[12]; d__1 = a[3]; d__2 = a[4]; d__3 = a[5]; d__4 = a[14]; a[15] = (SQR(d__1) + SQR(d__2) + SQR(d__3) + a[0]*a[6] + a[1]*a[7] + a[2]*a[8]) * a[12] + SQR(d__4); d__1 = a[14]; a[16] = ((a[3]*a[6] + a[4]*a[7] + a[5]*a[8]) * 9 + (a[0]*a[9] + a[1]*a[10] + a[2]*a[11]) * 3) * a[12] + a[14] * (a[15] * 9.0 - SQR(d__1) * 12.); for (k = 0; k < NDIM; ++k) { f1dot[k] = a[k + 3] - a[14] * 3 * a[k]; f2dot[k] = (a[k + 6] - a[14] * 6 * f1dot[k] - a[15] * 3 * a[k]) * a[13]; f3dot[k] = (a[k + 9] - a[15] * 9 * f1dot[k] - a[16] * a[k]) * a[13]; d2[k+i*3] += f2dot[k]; d3[k+i*3] = d3[k+i*3] + f3dot[k] - a[14] * 9 * f2dot[k]; } } } /* for (i) */ /* initialize integration steps and convert to force diffences */ for (i = 0; i < n; ++i) { d__1 = X(f,i); d__2 = Y(f,i); d__3 = Z(f,i); d__4 = X(d2,i); d__5 = Y(d2,i); d__6 = Z(d2,i); step[i] = sqrt(eta * sqrt( (SQR(d__1)+SQR(d__2)+SQR(d__3)) / (SQR(d__4)+SQR(d__5)+SQR(d__6)) ) ); t0[i] = time; t1[i] = time - step[i]; t2[i] = time - step[i] * 2; t3[i] = time - step[i] * 3; for (k = 0; k < NDIM; ++k) { d1[k+i*3] = (d3[k+i*3] * step[i] / 6 - d2[k+i*3] * 0.5) * step[i] + fdot[k+i*3]; d2[k+i*3] = d2[k+i*3] * 0.5 - d3[k+i*3] * 0.5 * step[i]; d3[k+i*3] /= 6; f[k+i*3] *= 0.5; fdot[k+i*3] /= 6; } } /* for (i) */ /* ----------------------------------------------------------------------- */ /* ==============> entry point when major output + diagnostics to be done */ /* energy check and output */ L100: e = 0.0; for (i = 0; i < n; ++i) { dt = tnext - t0[i]; for (k = 0; k < NDIM; ++k) { f2dot[k] = d3[k+i*3] * (t0[i] - t1[i] + (t0[i] - t2[i])) + d2[k+i*3]; x[k+i*3] = ((((d3[k+i*3] * 0.05 * dt + f2dot[k] / 12.0) * dt + fdot[k+i*3]) * dt + f[k+i*3]) * dt + x0dot[k+i*3]) * dt + x0[k+i*3]; a[k] = (((d3[k+i*3] * 0.25 * dt + f2dot[k] / 3.0) * dt + fdot[k+i*3] * 3.0) * dt + f[k+i*3] * 2.0) * dt + x0dot[k+i*3]; } itmp=i+1; outbods_(&body[i], &x[i*3], a, &step[i], &itmp); d__1 = a[0]; d__2 = a[1]; d__3 = a[2]; e += 0.5 * body[i] * (SQR(d__1) + SQR(d__2) + SQR(d__3)); } for (i = 0; i < n; ++i) { for (j = 0; j < n; ++j) { if (j == i) continue; d__1 = X(x,i) - X(x,j); d__2 = Y(x,i) - Y(x,j); d__3 = Z(x,i) - Z(x,j); e -= body[i] * 0.5 * body[j] / sqrt(SQR(d__1) + SQR(d__2) + SQR(d__3) + eps2); } } outene_(&tnext, &nsteps, &e); if (time > tcrit) { printf("Time spent in searching for next advancement: %g\n", time1*60.0); return 0; } tnext += deltat; /* ----------------------------------------------------------------------- */ /* ===============> Normal entry point for next timstep */ /* find next body (i) to be advanced and set new time * this linear search takes about 5% of the total CPU * see sift() to speed up this process */ L200: time0=cputime(); time = 1e10; for (j = 0; j < n; ++j) { if (time > t0[j] + step[j]) { i = j; time = t0[j] + step[j]; } } time1 += (cputime()-time0); /* predict all coordinates to first order in force derivative */ for (j = 0; j < n; ++j) { s = time - t0[j]; X(x,j) = ((X(fdot,j) * s + X(f,j)) * s + X(x0dot,j)) * s + X(x0,j); Y(x,j) = ((Y(fdot,j) * s + Y(f,j)) * s + Y(x0dot,j)) * s + Y(x0,j); Z(x,j) = ((Z(fdot,j) * s + Z(f,j)) * s + Z(x0dot,j)) * s + Z(x0,j); } /* include 2nd and 3rd order and obtain the velocity */ dt = time - t0[i]; for (k = 0; k < NDIM; ++k) { f2dot[k] = d3[k+i*3] * (t0[i] - t1[i] + (t0[i] - t2[i])) + d2[k+i*3]; d__1 = dt, d__1 *= d__1; x[k+i*3] = (d3[k+i*3] * 0.05 * dt + f2dot[k] / 12.) * SQR(d__1) + x[k+i*3]; x0dot[k+i*3] = (((d3[k+i*3] * 0.25 * dt + f2dot[k] / 3.0) * dt + fdot[k+i*3] * 3.0) * dt + f[k+i*3] * 2.0) * dt + x0dot[k+i*3]; f1[k] = 0.0; } /* obtain current force on i-th body */ for (j = 0; j < n; ++j) { if (j == i) continue; d__1 = a[0] = X(x,j) - X(x,i); d__2 = a[1] = Y(x,j) - Y(x,i); d__3 = a[2] = Z(x,j) - Z(x,i); a[3] = 1.0 / (SQR(d__1) + SQR(d__2) + SQR(d__3) + eps2); a[4] = body[j] * a[3] * sqrt(a[3]); f1[0] += a[0] * a[4]; f1[1] += a[1] * a[4]; f1[2] += a[2] * a[4]; } /* set time intervals for new difference and update the times */ dt1 = time - t1[i]; dt2 = time - t2[i]; dt3 = time - t3[i]; t1pr = t0[i] - t1[i]; t2pr = t0[i] - t2[i]; t3pr = t0[i] - t3[i]; t3[i] = t2[i]; t2[i] = t1[i]; t1[i] = t0[i]; t0[i] = time; /* form new differences and include 4th order semi-iterative */ for (k = 0; k < NDIM; ++k) { a[k ] = (f1[k ] - f[k+i*3 ] * 2.0) / dt; a[k + 3] = (a[k ] - d1[k+i*3]) / dt1; a[k + 6] = (a[k + 3] - d2[k+i*3]) / dt2; a[k + 9] = (a[k + 6] - d3[k+i*3]) / dt3; d1[k+i*3] = a[k ]; d2[k+i*3] = a[k + 3]; d3[k+i*3] = a[k + 6]; f1dot[k] = t1pr * t2pr * t3pr * a[k + 9]; f2dot[k] = (t1pr * t2pr + t3pr * (t1pr + t2pr)) * a[k + 9]; f3dot[k] = (t1pr + t2pr + t3pr) * a[k + 9]; d__1 = dt, d__2 = d__1; x0[k+i*3] = (((a[k + 9] * dt / 30. + f3dot[k] * 0.05) * dt + f2dot[k] / 12.) * dt + f1dot[k] / 6.) * (d__2 * SQR(d__1)) + x[k+i*3]; d__1 = dt; x0dot[k+i*3] = (((a[k + 9] * 0.2 * dt + f3dot[k] * 0.25) * dt + f2dot[k] / 3.0) * dt + f1dot[k] * 0.5) * SQR(d__1) + x0dot[k+i*3]; } /* scale F and FDOT by factorials and set new integration step */ for (k = 0; k < NDIM; ++k) { f[k+i*3] = f1[k] * 0.5; fdot[k+i*3] = ((d3[k+i*3] * dt1 + d2[k+i*3]) * dt + d1[k+i*3]) / 6.0; f2dot[k] = (d3[k+i*3] * (dt + dt1) + d2[k+i*3]) * 2.0; } d__1 = f1[0]; d__2 = f1[1]; d__3 = f1[2]; d__4 = f2dot[0]; d__5 = f2dot[1]; d__6 = f2dot[2]; step[i] = sqrt(eta * sqrt( (SQR(d__1)+SQR(d__2)+SQR(d__3)) / (SQR(d__4)+SQR(d__5)+SQR(d__6)) ) ); ++nsteps; if (time - tnext >= 0.0) { goto L100; } else { goto L200; } } /* nbody0_ */ 0; j < n; ++j) { if (i == j) continue; for (k = 0; k < NDIM; ++k) { a[k ] = x0[k+j*3] - x0[k+i*3]; a[k + 3] = x0dot[k+j*3] - x0dot[k+i*3]; a[k + 6] = f[k+j*3] - f[k+i*3]; a[k + 9] = nbody0.h * include file for nbody0 */ #include #define NMAX 2048 #define NDIM 3 # 21-nov HOSTTYPE # 17-may-91 NEMO V2.x # 16-oct-91 new sun fortran 1.4 compiler - linking via f77 # 10-may-92 handedited f2c version of nbody0.f # 2-apr-97 added Testfile for the new 2.4 release, and # made default flags -O # 9-mar-98 both include files now created via makefile CFLAGS = -O FFLAGS = -O # maximum number of particles hardcoded in fortran and C version NMAX = 2048 # Note the absence of the NDIM macro, as NDIM <> 3 is not supported # in nbody0.f (yet). - july 1989 , PJT L = $(NEMOLIB)/libnemo.a # linking is now done via f77 command, hence no LIBF77 needed anymore #LIBF77 = -lF77 -lI77 -lU77 -lF77 -lI77 -lU77 LIBF77 = SRCDIR = $(NEMO)/src/nbody/evolve/aarseth/nbody0 SRC = nbody0.f nbodyio.c nbodyio.f main.c main.f \ fdefs.inc nmax.inc \ nbody0.c nbody0.h \ Makefile Testfile README BIN = nbody0 nbody00 nbody0_ff MAN1= nbody0.1 MAN = $(MAN1) TAR = $(SRC) $(MAN) nbody0.man help: @echo "Aarseth programs - from Binney and Tremaine" @echo "Various Fortran and C versions available" @echo "target: main: I/O:" @echo "" @echo " ff F F (vanilla I/O)" @echo " cc c c (nemo)" @echo "" @echo "Note:" @echo "User definable macro is the NMAX macro (currently $(NMAX))" @echo "which is used to generate a new include file" @echo "such that this will be the maximum number of particles" @echo "that can be integrated with NBODY0" @echo "Assuming your local fortran compiler FC=$(FC) knows the" @echo "INCLUDE statement." # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # N E M O Install # nemo_lib: @echo No nemo_lib in `pwd` # # nemo_bin: $(BIN) mv $? $(NEMOBIN) rm -f *.o nemo_src: -@for i in $(BIN); do \ echo `pwd` $$i ; done # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # install: old: .install_src .install_man .install_src: $(SRC) cp $? $(SRCDIR) @touch .install_src .install_man: .install_man1 .install_man1: $(MAN1) cp $? $(NEMO)/man/man1 @touch .install_man1 install_bin: $(BIN) clean: rm -f $(BIN) *.o *.a core # # # # # # # # # # # # # # # # # # # # # # # # # nbody0: main.o nbody0.o nbodyio.o $(FC) $(FFLAGS) -o nbody0 main.o nbodyio.o nbody0.o \ $(L) $(LIBF77) -lm nbody00: main.o nbody00.o nbodyio.o $(FC) $(FFLAGS) -o nbody00 main.o nbodyio.o nbody00.o \ $(L) $(LIBF77) -lm nbody00.o: nbody0.c $(CC) $(CFLAGS) -c nbody0.c mv nbody0.o nbody00.o main.o: main.c $(CC) $(CFLAGS) -c main.c nbodyio.o: nbodyio.c $(CC) $(CFLAGS) -c nbodyio.c nbody0.o: nbody0.f nmax.inc $(FC) $(FFLAGS) -c nbody0.f # Vanilla fortran stuff: no NEMO I/O nbody0_ff: main.f nbodyio.f nmax.inc nbody0.o $(FC) $(FFLAGS) -o nbody0_ff main.f nbodyio.f \ nbody0.o -lm double: @echo Tarfile for cray-et-al, all DOUBLEs become REAL @mkdir tmp -@for i in $(SRC); do\ (sed "s/DOUBLE PRECISION/REAL/g" $$i > tmp/$$i);done @(cd tmp;tar cvf ../double.tar *;cd ..;rm -r tmp) @echo All done, file double.tar created tar: nbody0.man @tar cvf nbody0.tar $(TAR) @echo Testing if your man page is not out of date.... @diff nbody0.1 $(NEMO)/man/man1 # this is only for the NEMO maintainer pjt: gzip nbody0.tar; mv nbody0.tar.gz /home/ftp/progs/nemo nbody0.man: nbody0.1 groff -M$(NEMO)/man -Tascii -man nbody0.1 > nbody0.man nmax: .dummy .dummy: @echo "Creating new nmax.inc with NMAX=$(NMAX)" @echo "C include file generated by Makefile - do not edit" > nmax.inc @echo "C see Makefile for instructions" >> nmax.inc @echo " INTEGER NMAX, NDIM" >> nmax.inc @echo " PARAMETER ( NMAX=$(NMAX) )" >> nmax.inc @echo "C do not changed the NDIM value, nbody0.f does not support it" >>nmax.inc @echo " PARAMETER ( NDIM=3 )" >> nmax.inc @echo "C end of include file" >> nmax.inc @echo "A new file nmax.inc with NMAX=$(NMAX) has just been created." @echo "/*" > nbody0.h @echo " * include file for nbody0" >> nbody0.h @echo " */" >> nbody0.h @echo "#include " >> nbody0.h @echo "#define NMAX $(NMAX)" >> nbody0.h @echo "#define NDIM 3" >> nbody0.h @echo "and also a new nbody0.h for the C version." BIN = nbody0 NEED = $(BIN) mkplummer snaptrim help: @echo $(DIR) need: @echo $(NEED) clean: @echo Cleaning $(DIR) @rm -f nbody0.in nbody0.out nbody0.log NBODY = 10 all: $(BIN) nbody0.in: @echo Creating nbody0.in mkplummer nbody0.in $(NBODY) seed=123 nbody0: nbody0.in @echo Running $@ nbody0 nbody0.in nbody0.out tcrit=2 > nbody0.log snaptrim nbody0.out - times=2 | tsf - The modifications are listed in the source code, nbody0.f It also has an optional NEMO I/O interface, as explained in the Makefile. One can also use the vanilla Fortran version, in which case the NEMO library is not needed to complete linking. Hernquist' ascii 205 format (see NEMO programs atos and stoa to convert from/to snapshots) is used to store snapshots. A manual page, nbody0.1, is available which explains the NEMO user interface part. It normally lives in $NEMO/man/man1, but a local copy is kept here, as well as an ascii formatted version of the manual page. A makefile for the Cray is also available, but not exported in the non-NEMO (see http://www.astro.umd.edu/nemo/) version. A C version of the compute engine is also available under the name nbody0.c, though the Makefile has a target name 'nbody00' to make it into an executable. It runs typically 10% slower then the fortran version (on my sparc-1). The C version was generated with the public domain program 'f2c', and subsequently hand edited to optimize and improve readability. Peter Teuben - november 1991, may 1992 .SH NAME nbody0 - direct summation Aarseth N-body integrator .SH SYNOPSIS nbody0 in=snap_in out=snap_out [keyword=value ...] .SH DESCRIPTION \fBnbody0\fP is the NEMO adaptation of version 0 (the \fIMicky Mouse\fP version) of Aarseth's N-body integrator, as published in: Binney & Tremaine (1987). It is a direct N-body integrator, \fIi.e.\fP for each particle it computes the force due to all other N-1 particles, hence the computational time grows approximately as N*N. Although also being referred to as a 'toy version', it is a fully functional N-body integrator. .PP Each particle is followed with its own integration step - an essential feature when the dynamical times of different particles vary a lot. A complete description is given in: S.J. Aarseth, "Multiple Time Scales", ed. U.J. Brackhill & B.I. Cohen, pp377. Orlando: Academic Press. (1985). .PP An important input parameter to the program is the accuracy parameter (called \fIeta\fP below, see also \fBeta=\fP below). The timestep, \fIdt\fP, chosen for a given particle is related to the force, \fIF\fP, and its time derivatives by \fIdt = sqrt(eta * F / (d2Fdt2))\fP, which is a slight modification of the criterion given by Aarseth in the above mentioned reference. .PP A typical value of \fIeta = 0.02\fP usually conserves the energy better than one part in 10000 over one crossing time, in the absence of close encounters. .SH PARAMETERS The following parameters are recognized in any order when the keyword is also given .TP 25 \fBin=\fIin_name\fP Input file in snapshot/205 format. The actual format depends on which I/O module was nbody0 was compile with. See Makefile for details. [no default]. .TP \fBout=\fIout_name\fP Output file in snapshot/205 format [no default]. .TP \fBeta=\fIvalue\fP Accurracy parameter, which determines the integration step [0.02]. .TP \fBdeltat=\fIvalue\fP Time interval of a major output [0.25]. .TP \fBtcrit=\fIvalue\fP Final integration time [2.0]. .TP \fBeps=\fP\fIvalue\fP Softening length [0.05]. .TP \fBoptions=\fP\fBstep\fP Miscellaneous control options, specified as a comma separated list of keywords. Currently implemented are: \fBstep\fP, outputs the current timestep per particle in the \fIAux\fP array of a \fIsnapshot(5NEMO)\fP. Default: none. .SH LIMITATIONS The code has a hardcoded maximum number of particles in the fortran code (File: \fInmax.h\fP), change the relevant PARAMETER statement to whatever is required. Using NEMO install the Makefile macro NMAX can also be used to do this automatically and hence easily generate a new version with a different value for NMAX, \fIe.g.\fP \fBmake nmax NMAX=256 nbody0\fP. The NDIM parameter should not be changed, and remain 3. .PP Close encounters are not treated specially. See any of the modified versions of the Aarseth code (\fInbody2(1NEMO)\fP) or in case regularization is needed (see also newton0_reg in: \fInewton0(1NEMO)\fP). .PP Timesteps are not recomputed until the current timestep has expired. .PP In order to save time, all calculations in the fortran code (nbody0.f) can be done in single precision. A different version of nbody0.f is needed in this case. See code comments ``SED''. .PP A C version is also available, as nbody0.c, though the Makefile needs the target \fBnbody00\fP to install it. The user interface is the same as that of nbody0. .SH BUGS Starting time of initial conditions is set to 0, even if the snapshot had another time. .SH STORAGE A total of order \fI600.NMAX\fP bytes is needed, for a given maximum of \fINMAX\fP particles. This breaks down as follows: .PP The \fBFORTRAN\fP I/O code uses \fI10.NMAX\fP double precision and \fI20.NMAX\fP real storage units for a maximum compiled number of \fINMAX\fP particles. On most machines this adds up to \fI160.NMAX\fP bytes. .PP The \fBC\fP I/O code allocates space dynamically; and through limitations of the fortran code, a maximum of \fI100.NMAX\fP double's are needed, which means \fI400.NMAX\fP bytes. .SH AUTHOR Sverre Aarseth (F77), Peter Teuben (C) .SH GENEALOGY The \fBnbody*\fP series of NBODY integrators were all originally written by Sverre Aarseth, and exist in many forms throughout the scientific community. The list below is a mere zeroth order approximation to the current state of affairs: .nf .ta +1i nbody0 Binney & Tremaine's toy version, with NEMO interface nbody1 variable timestep * email version sent on Fri, 15 Oct 93 17:52 GMT nbody2 with Ahmad-Cohen scheme - (see also Benz et al, 1997) * email version 17-mar-93 nbody3 XXX nbody4 ... for HARP nbody5 with regularization handling triples & binaries nbody6 ... for new parallel machine .fi .SH FILES .nf .ta +2i ~/src/nbody/evolve/aarseth/nbody0 official source code within NEMO ~/usr/aarseth/ SJA's other nbodyX versions (not exported) .fi .SH SEE ALSO hackcode1(1NEMO), runbody2(1NEMO), snapshot(5NEMO), atos(1NEMO), stoa(1NEMO) .PP Binney, J. & Tremaine, S. \fIGalactic Dynamics\fP. Princeton Unversity Press (1987), pp678. .PP S.J.Aarseth, 1972, p.373 in: "Gravitational N-Body Problem", IAU Colloquium #10, M.Lecar (Ed.), Reidel, Dordrecht. .PP S.J.Aarseth, 1985, p.377 in: "Multiple Time Scales", U.J. Brackbill & B.I. Cohen (Eds.), Academic Press, Orlando. .SH HISTORY .nf .ta +1i +4i 30-jun-89 V1.0 created + NEMO interfaces to fortran source PJT 3-jul-89 V1.1 mods to f2c interface, name of keywords PJT 24-jan-90 V1.2 all in double precision PJT 15-nov-91 fixed up pure nbody0_ff version PJT 20-may-92 -- also made the C (f2c) version available PJT 2-apr-97 documentation updated for NEMO 2.4 PJT 6-aug-97 V1.3 added options= PJT 11-feb-98 V1.3a fixed array index bug for higher order term PJT .fi NBODY0(1NEMO) NEMO COMMANDS NBODY0(1NEMO) NNAAMMEE nbody0 - direct summation Aarseth N-body integrator SSYYNNOOPPSSIISS nbody0 in=snap_in out=snap_out [keyword=value ...] DDEESSCCRRIIPPTTIIOONN nnbbooddyy00 is the NEMO adaptation of version 0 (the _M_i_c_k_y _M_o_u_s_e version) of Aarseth's N-body integrator, as published in: Binney & Tremaine (1987). It is a direct N-body integrator, _i_._e_. for each particle it computes the force due to all other N-1 particles, hence the computational time grows approximately as N*N. Although also being referred to as a 'toy version', it is a fully functional N-body integrator. Each particle is followed with its own integration step - an essential feature when the dynamical times of different par- ticles vary a lot. A complete description is given in: S.J. Aarseth, "Multiple Time Scales", ed. U.J. Brackhill & B.I. Cohen, pp377. Orlando: Academic Press. (1985). An important input parameter to the program is the accuracy parameter (called _e_t_a below, see also eettaa== below). The timestep, _d_t, chosen for a given particle is related to the force, _F, and its time derivatives by _d_t _= _s_q_r_t_(_e_t_a _* _F _/ _(_d_2_F_d_t_2_)_), which is a slight modification of the criterion given by Aarseth in the above mentioned reference. A typical value of _e_t_a _= _0_._0_2 usually conserves the energy better than one part in 10000 over one crossing time, in the absence of close encounters. PPAARRAAMMEETTEERRSS The following parameters are recognized in any order when the keyword is also given iinn==_i_n___n_a_m_e Input file in snapshot/205 format. The actual format depends on which I/O module was nbody0 was compile with. See Makefile for details. [no default]. oouutt==_o_u_t___n_a_m_e Output file in snapshot/205 format [no default]. eettaa==_v_a_l_u_e Accurracy parameter, which deter- mines the integration step [0.02]. ddeellttaatt==_v_a_l_u_e Time interval of a major output [0.25]. ttccrriitt==_v_a_l_u_e Final integration time [2.0]. Nemo Release 2.4 Last change: 14 February 1998 1 NBODY0(1NEMO) NEMO COMMANDS NBODY0(1NEMO) eeppss==_v_a_l_u_e Softening length [0.05]. ooppttiioonnss==sstteepp Miscellaneous control options, specified as a comma separated list of keywords. Currently implemented are: sstteepp, outputs the current timestep per particle in the _A_u_x array of a _s_n_a_p_s_h_o_t_(_5_N_E_M_O_). Default: none. LLIIMMIITTAATTIIOONNSS The code has a hardcoded maximum number of particles in the fortran code (File: _n_m_a_x_._h), change the relevant PARAMETER statement to whatever is required. Using NEMO install the Makefile macro NMAX can also be used to do this automati- cally and hence easily generate a new version with a differ- ent value for NMAX, _e_._g_. mmaakkee nnmmaaxx NNMMAAXX==225566 nnbbooddyy00. The NDIM parameter should not be changed, and remain 3. Close encounters are not treated specially. See any of the modified versions of the Aarseth code (_n_b_o_d_y_2_(_1_N_E_M_O_)) or in case regularization is needed (see also newton0_reg in: _n_e_w_- _t_o_n_0_(_1_N_E_M_O_)). Timesteps are not recomputed until the current timestep has expired. In order to save time, all calculations in the fortran code (nbody0.f) can be done in single precision. A different version of nbody0.f is needed in this case. See code com- ments ``SED''. A C version is also available, as nbody0.c, though the Make- file needs the target nnbbooddyy0000 to install it. The user inter- face is the same as that of nbody0. BBUUGGSS Starting time of initial conditions is set to 0, even if the snapshot had another time. SSTTOORRAAGGEE A total of order _6_0_0_._N_M_A_X bytes is needed, for a given maxi- mum of _N_M_A_X particles. This breaks down as follows: The FFOORRTTRRAANN I/O code uses _1_0_._N_M_A_X double precision and _2_0_._N_M_A_X real storage units for a maximum compiled number of _N_M_A_X particles. On most machines this adds up to _1_6_0_._N_M_A_X bytes. The CC I/O code allocates space dynamically; and through lim- itations of the fortran code, a maximum of _1_0_0_._N_M_A_X double's are needed, which means _4_0_0_._N_M_A_X bytes. Nemo Release 2.4 Last change: 14 February 1998 2 NBODY0(1NEMO) NEMO COMMANDS NBODY0(1NEMO) AAUUTTHHOORR Sverre Aarseth (F77), Peter Teuben (C) GGEENNEEAALLOOGGYY The nnbbooddyy** series of NBODY integrators were all originally written by Sverre Aarseth, and exist in many forms through- out the scientific community. The list below is a mere zeroth order approximation to the current state of affairs: nbody0 Binney & Tremaine's toy version, with NEMO interface nbody1 variable timestep * email version sent on Fri, 15 Oct 93 17:52 GMT nbody2 with Ahmad-Cohen scheme - (see also Benz et al, 1997) * email version 17-mar-93 nbody3 XXX nbody4 ... for HARP nbody5 with regularization handling triples & binaries nbody6 ... for new parallel machine FFIILLEESS ~/src/nbody/evolve/aarseth/nbody0official source code within NEMO ~/usr/aarseth/ SJA's other nbodyX versions (not exported) SSEEEE AALLSSOO hackcode1(1NEMO), runbody2(1NEMO), snapshot(5NEMO), atos(1NEMO), stoa(1NEMO) Binney, J. & Tremaine, S. _G_a_l_a_c_t_i_c _D_y_n_a_m_i_c_s. Princeton Unversity Press (1987), pp678. S.J.Aarseth, 1972, p.373 in: "Gravitational N-Body Problem", IAU Colloquium #10, M.Lecar (Ed.), Reidel, Dordrecht. S.J.Aarseth, 1985, p.377 in: "Multiple Time Scales", U.J. Brackbill & B.I. Cohen (Eds.), Academic Press, Orlando. HHIISSTTOORRYY 30-jun-89 V1.0 created + NEMO interfaces to fortran sourcePJT 3-jul-89 V1.1 mods to f2c interface, name of keywords PJT 24-jan-90 V1.2 all in double precision PJT 15-nov-91 fixed up pure nbody0_ff version PJT 20-may-92 -- also made the C (f2c) version available PJT 2-apr-97 documentation updated for NEMO 2.4 PJT 6-aug-97 V1.3 added options= PJT 11-feb-98 V1.3a fixed array index bug for higher order termPJT Nemo Release 2.4 Last change: 14 February 1998 3