From ioannis@ircamera.as.arizona.edu Tue Nov 14 18:06:29 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 SAA05182 for ; Tue, 14 Nov 2000 18:06:29 -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 SAA14371 for ; Tue, 14 Nov 2000 18:25:29 -0700 Received: (from ioannis@localhost) by ircamera.as.arizona.edu (8.9.3/8.8.7) id SAA03355 for dzaritsky@as.arizona.edu; Tue, 14 Nov 2000 18:29:35 -0700 Date: Tue, 14 Nov 2000 18:29:35 -0700 From: John Moustakas Message-Id: <200011150129.SAA03355@ircamera.as.arizona.edu> To: dzaritsky@as.arizona.edu Status: RO pro twobody, n, ndim, body, x0, x0dot ;+ ; NAME: ; TWOBODY ; ; PURPOSE: ; Default N-body problem to solve using NBODY: two particles of ; equal mass in one plane with a small amount of initial ; velocity. ; ; INPUTS: ; None. ; ; MODIFICATION HISTORY: ; John Moustakas, 2000 November 14 ;- n = 2L ; number of particles ndim = 3L ; number of dimensions (x,y,z) body = [1.0,1.0] ; initial masses x0 = [[5.0,0.0,0.0],$ ; initial positions [-5.0,0.0,0.0]] x0dot = [[0.0,0.05,0.0],$ ; initial velocities [0.0,-0.05,0.0]] return end pro initialize_arrays, t0, t1, t2, t3, step, d1, d2, d3, f, fdot, f1, $ f1dot, f2dot, f3dot, time, tnext, nsteps, x, a, n, ndim ;+ ; NAME: ; INITIALIZE_ARRAYS ; ; PURPOSE: ; Initialize all the arrays that will be filled by NBODY. ; ; MODIFICATION HISTORY: ; John Moustakas, 2000 November 14 ;- t0 = fltarr(n) t1 = fltarr(n) t2 = fltarr(n) t3 = fltarr(n) step = fltarr(n) d1 = fltarr(ndim,n) d2 = fltarr(ndim,n) d3 = fltarr(ndim,n) f = fltarr(ndim,n) f1 = fltarr(ndim) fdot = fltarr(ndim,n) f1dot = fltarr(ndim) f2dot = fltarr(ndim) f3dot = fltarr(ndim) time = 0.0 tnext = 0.0 nsteps = 0L x = fltarr(ndim,n) ; output positions a = fltarr(17) ; a[0:ndim-1L] are output velocities return end pro nbody, eta=eta, deltat=deltat, tcrit=tcrit, eps2=eps2 ;+ ; NAME: ; NBODY ; ; PURPOSE: ; Compute the time evolution of N point masses under the ; influence of their mutual self-gravity. ; ; INPUTS: ; None. ; ; OPTIONAL INPUTS: ; eta : accuracy parameter ; deltat : output interval ; tcrit : length of the integration ; eps2 : square of the softening parameter ; ; OUTPUTS: ; The mass (body), timestep (step), position (x), and velocity ; (a[0:ndim-1L]), as well as the time (tnext), the cumulative ; number of steps (nsteps), and the energy (e) are printed to ; the screen. ; ; COMMON BLOCKS: ; None. ; ; COMMENTS: ; This program was adopted directly from S. J. Aarseth's ; Standard N-Body Program found in Appendix 4.B of Binney & ; Tremaine (1987). The code gives consistent results with the ; Fortran version. ; ; To keep the program simple, no plots have been generated and ; no text files are written. Note that all calculations are ; done with floating point precision. ; ; The default is to solve a simple 2-body problem with simple ; initial conditions. Any procedure can be passed in the place ; of TWOBODY. ; ; PROCEDURES USED: ; INITIALIZE_ARRAYS, TWOBODY ; ; MODIFICATION HISTORY: ; John Moustakas & Andy Marble, 2000 November 14 ;- if not keyword_set(eta) then eta = 0.02 ; accuracy parameter if not keyword_set(deltat) then deltat = 0.5 ; output time interval if not keyword_set(tcrit) then tcrit = 50.0 ; integration length if not keyword_set(eps2) then eps2 = 0.05 ; softening parameter squared twobody, n, ndim, body, x0, x0dot ; two-body problem (default) initialize_arrays, t0, t1, t2, t3, step, d1, d2, d3, $ ; initialize the arrays f, fdot, f1, f1dot, f2dot, f3dot, time, tnext, nsteps, $ x, a, n, ndim ; obtain the total force and the first derivative for each body for i = 0L, n-1L do begin for j = 0L, n-1L do begin if (j eq i) then goto, m10 for k = 0L, ndim-1L do begin a[k] = x0[k,j] - x0[k,i] a[k+3] = x0dot[k,j] - x0dot[k,i] endfor a[6] = 1.0/(a[0]^2 + a[1]^2 + a[2]^2 + eps2) a[7] = body[j]*a[6]*sqrt(a[6]) a[8] = 3.0*(a[0]*a[3] + a[1]*a[4] + a[2]*a[5])*a[6] for k =0L, ndim-1L do begin f[k,i] = f[k,i] + a[k]*a[7] fdot[k,i] = fdot[k,i] + (a[k+3] - a[k]*a[8])*a[7] endfor m10: endfor endfor ; form second and third force derivative for i = 0L, n-1L do begin for j = 0L, n-1L do begin if (j eq i) then goto, m30 for k = 0L, ndim-1L do begin 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] endfor a[12] = 1.0/(a[0]^2 + a[1]^2 + a[2]^2 + 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] a[15] = (a[3]^2 + a[4]^2 + a[5]^2 + a[0]*a[6] + a[1]*a[7] + a[2]*a[8])*a[12] + a[14]^2 a[16] = (9.0*(a[3]*a[6] + a[4]*a[7] + a[5]*a[8]) + 3.0*(a[0]*a[9] + a[1]*a[10] + a[2]*a[11]))*a[12] + $ a[14]*(9.0*a[15] - 12.0*a[14]^2) for k = 0L, ndim-1L do begin f1dot[k] = a[k+3] - 3.0*a[14]*a[k] f2dot[k] = (a[k+6] - 6.0*a[14]*f1dot[k] - 3.0*a[15]*a[k])*a[13] f3dot[k] = (a[k+9] - 9.0*a[15]*f1dot[k] - a[16]*a[k])*a[13] d2[k,i] = d2[k,i] + f2dot[k] d3[k,i] = d3[k,i] + f3dot[k] - 9.0*a[14]*f2dot[k] endfor m30: endfor endfor ; initialize integration steps and convert to force differences for i = 0L, n-1L do begin step[i] = sqrt(eta*sqrt((f[0,i]^2 + f[1,i]^2 + f[2,i]^2)/(d2[0,i]^2 + d2[1,i]^2 + d2[2,i]^2))) t0[i] = time t1[i] = time - step[i] t2[i] = time - 2.0*step[i] t3[i] = time - 3.0*step[i] for k = 0L, ndim-1L do begin d1[k,i] = (d3[k,i]*step[i]/6.0 - 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.0 f[k,i] = 0.5*f[k,i] fdot[k,i] = fdot[k,i]/6.0 endfor endfor ; entry point when major output -- diagnostics to be done. energy check and output m100: e = 0.0 for i = 0L, n-1L do begin dt = tnext - t0[i] for k = 0L, ndim-1L do begin 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] endfor print, i, body[i], step[i], x[*,i], a[0:2], format='(I10,F10.2,F12.4,3x,3F10.2,3x,3F10.2)' e = e + 0.5*body[i]*(a[0]^2 + a[1]^2 + a[2]^2) endfor for i = 0L, n-1L do begin for j = 0L, n-1L do begin if (j eq i) then goto, m120 e = e - 0.5*body[i]*body[j]/sqrt((x[0,i] - x[0,j])^2 + (x[1,i] - x[1,j])^2 + (x[2,i] - x[2,j])^2 + eps2) m120: endfor endfor print print, 'TIME = '+strn(tnext,format='(F7.2)')+', STEPS = '+strn(nsteps,format='(I6)')+$ ', ENERGY = '+strn(e,format='(F7.2)') print if (time gt tcrit) then return tnext = tnext + deltat ; normal entry point for next time step. find next body to be ; advanced and set new time m200: time = 1.0E+10 for j = 0L, n-1L do begin if (time gt t0[j] + step[j]) then i = j if (time gt t0[j] + step[j]) then time = t0[j] + step[j] endfor ; predict all coordinates to first order in force derivative for j = 0L, n-1L do begin s = time - t0[j] x[0,j]=((fdot[0,j]*s+f[0,j])*s + x0dot[0,j])*s+x0[0,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] endfor ; include second and third order and obtain the velocity dt = time - t0[i] for k = 0L, ndim-1L do begin 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 endfor ; obtain the current force on ith body for j = 0L, n-1L do begin if (j eq i) then goto, m240 a[0] = x[0,j] - x[0,i] a[1] = x[1,j] - x[1,i] a[2] = x[2,j] - x[2,i] a[3] = 1.0/(a[0]^2 + a[1]^2 + a[2]^2 + eps2) a[4] = body[j]*a[3]*sqrt(a[3]) f1[0] = f1[0] + a[0]*a[4] f1[1] = f1[1] + a[1]*a[4] f1[2] = f1[2] + a[2]*a[4] m240: endfor ; set time intervals for new differences 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 fourth-order semi-iteration for k = 0L, ndim-1L do begin 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] endfor ; scale f and fdot by factorials and set new integration step for k = 0L, ndim-1L do begin 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]) endfor step[i] = sqrt(eta*sqrt((f1[0]^2+f1[1]^2+f1[2]^2)/(f2dot[0]^2+f2dot[1]^2+f2dot[2]^2))) nsteps = nsteps + 1L if (time-tnext) lt 0.0 then goto, m200 else goto, m100 return end