c=========================================================== c PHY/CAM 381C: Homework 3, Problem 3 c c History: sode c c nbody integrates the equations of motion for c 'np' point particles interacting gravitationally. c c See Homework handout for description of arguments c standard input, standard output etc. c c Note: One 'trick' used here is to take advantage c of FORTRAN storage convention to have nominally 1-d c arrays 'y', 'ydot' (used by the subroutine c 'fcn' which LSODA calls) interpreted as 3-d c arrays for convenience in implementing the equations c of motion and computation of auxiliary quantities c (energy, linear and angular momentum, etc.). See c source code for 'fcn1.f' for more details. c=========================================================== program nbody implicit none real*8 dvvdot character*5 cdnm parameter ( cdnm = 'nbody' ) integer iargc, indlnb, i4arg real*8 r8arg real*8 r8_never parameter ( r8_never = -1.0d-60 ) c----------------------------------------------------------- c tfinal: Final integration time c dtout: Output interval c ntout: # of output times (computed) c trace: Enables tracing of "conserved quantities" c----------------------------------------------------------- real*8 tfinal, dtout integer ntout logical trace c----------------------------------------------------------- c Common communication with routine 'fcn' in 'fcn.f' ... c c Includes defn of maximum # of particles, storage c for particle masses, Newton's gravitational c constant ... c----------------------------------------------------------- include 'fcn.inc' integer maxneq, neq parameter ( maxneq = 6 * maxnp ) c----------------------------------------------------------- c Storage for energy, momentum, center of mass ... c----------------------------------------------------------- real*8 ke(maxnp), pe(maxnp) real*8 etot, ketot, petot real*8 pmom(d,maxnp), jmom(d,maxnp), & ptot(d), jtot(d) real*8 com(d) c----------------------------------------------------------- c LSODA Variables. c----------------------------------------------------------- external fcn, jac real*8 y(d,2,maxnp) real*8 tbgn, tend integer itol real*8 rtol, atol integer itask, istate, iopt integer lrw parameter ( lrw = 22 + maxneq * 16 ) real*8 rwork(lrw) integer liw parameter ( liw = 20 + maxneq ) integer iwork(liw) integer jt real*8 tol real*8 default_tol parameter ( default_tol = 1.0d-6 ) c----------------------------------------------------------- c Other locals c----------------------------------------------------------- integer a, i, j, & itout c----------------------------------------------------------- c Fix gravitational constant ... c----------------------------------------------------------- G = 1.0d0 c----------------------------------------------------------- c Parse command line arguments (initial values) ... c----------------------------------------------------------- if( iargc() .lt. 2 ) go to 900 tfinal = r8arg(1,r8_never) dtout = r8arg(2,r8_never) tol = r8arg(3,default_tol) trace = iargc() .ge. 4 if( tfinal .eq. r8_never .or. dtout .eq. r8_never & .or. dtout .le. 0.0d0 ) go to 900 c----------------------------------------------------------- c Compute number of uniform time steps requested. c----------------------------------------------------------- ntout = tfinal / dtout + 1.5d0 c----------------------------------------------------------- c Get particle masses, initial positions and c velocities. c----------------------------------------------------------- call getivs('-',m,y,d,maxnp,np) if( np .eq. 0 ) then write(0,*) cdnm//': No particles' stop end if neq = 6 * np c----------------------------------------------------------- c Dump # of particles, particle masses, initial c time and initial particle positions to standard c output. c----------------------------------------------------------- write(*,*) np do i = 1 , np write(*,*) m(i) end do tbgn = 0.0d0 write(*,*) tbgn do i = 1 , np write(*,*) y(1,1,i), y(2,1,i), y(3,1,i) end do c----------------------------------------------------------- c Compute initial energy, center of mass, linear mom c and ang mom about center of mass and output if c standard error tracing is enabled. c----------------------------------------------------------- if( trace ) then call clce(y,m,ke,pe,ketot,petot,etot,G,d,np) call clccom(y,m,com,d,np) call clcmom(y,m,pmom,jmom,ptot,jtot,com,d,np) write(0,5000) tbgn, & com(1), com(2), com(3), & etot, ketot, petot, & sqrt(dvvdot(ptot,ptot,d)), & sqrt(dvvdot(jtot,jtot,d)) end if c----------------------------------------------------------- c Set LSODA parameters ... c----------------------------------------------------------- itol = 1 rtol = tol atol = tol itask = 1 iopt = 0 jt = 2 istate = 1 c----------------------------------------------------------- c Do the integration ... c----------------------------------------------------------- do itout = 2 , ntout c----------------------------------------------------------- c Set next output time. c----------------------------------------------------------- tend = tbgn + dtout c----------------------------------------------------------- c Integrate to the next output time. c----------------------------------------------------------- call lsoda(fcn,neq,y,tbgn,tend, & itol,rtol,atol,itask, & istate,iopt,rwork,lrw,iwork,liw,jac,jt) c----------------------------------------------------------- c Check for LSODA error and exit if one occured. c----------------------------------------------------------- if( istate .lt. 0 ) then write(0,1500) cdnm, istate, itout, ntout, & tbgn, tend 1500 format(/' ',a,': Error return ',i2, & ' from integrator LSODA.'/ & ' At output time ',i5,' of ',i5/ & ' Interval ',1pe11.3,' .. ', & 1pe11.3/) go to 950 end if c----------------------------------------------------------- c Compute energy, COM, linear momentum and angular c momentum about the COM and output if tracing c is enabled. c----------------------------------------------------------- if( trace ) then call clce(y,m,ke,pe,ketot,petot,etot,G,d,np) call clccom(y,m,com,d,np) call clcmom(y,m,pmom,jmom,ptot,jtot,com,d,np) write(0,5000) tend, & com(1), com(2), com(3), & etot, ketot, petot, & sqrt(dvvdot(ptot,ptot,d)), & sqrt(dvvdot(jtot,jtot,d)) 5000 format(1p,10E25.16) end if c----------------------------------------------------------- c Output time and particle coordinates. c----------------------------------------------------------- write(*,*) tend do i = 1 , np write(*,*) y(1,1,i), y(2,1,i), y(3,1,i) end do c----------------------------------------------------------- c End of integration loop. c----------------------------------------------------------- end do 500 continue stop c----------------------------------------------------------- c Usage exit. c----------------------------------------------------------- 900 continue write(0,*) 'usage: '//cdnm// & '
[ ]' write(0,*) ' ' write(0,*) ' Masses, initial positions and' write(0,*) ' velocities of particles read from' write(0,*) ' standard input' stop c----------------------------------------------------------- c LSODA error exit. c----------------------------------------------------------- 950 continue write(0,*) 'nbody: Error return from LSODA' stop stop end c=========================================================== c Routine called by LSODA. Simply calls auxiliary c routine 'fcn1' which interprets 'y' and 'ydot' c as three dimensional arrays for coding convenience, c and which actually implements equations of motion. c=========================================================== subroutine fcn(neq,t,y,ydot) implicit none integer neq real*8 t, y(neq), ydot(neq) call fcn1(y,ydot,neq/6) return end c=========================================================== c Implements differential equations for gravitational c nbody problem. Note use of FORTRAN storage-order c convention to "re-interpret" arrays 'y' and 'ydot' c as 3-dimensional arrays: c c First index: dimension (coordinate) c Second index: positions (1) or velocities (2) c Third index: particle number c=========================================================== subroutine fcn1(y,ydot,npart) implicit none include 'fcn.inc' integer neq integer npart real*8 y(3,2,npart), ydot(3,2,npart) integer a, i, j real*8 rsq, c1 c----------------------------------------------------------- c Implement equations which define velocities, c and initialize accelerations to 0. c----------------------------------------------------------- do i = 1 , npart do a = 1 , d ydot(a,1,i) = y(a,2,i) ydot(a,2,i) = 0.0d0 end do end do c----------------------------------------------------------- c Loop over pairs of particles, compute pair-wise c forces and update accelerations of *both* c particles affected by each force. c----------------------------------------------------------- do i = 1 , npart do j = 1 , i - 1 c----------------------------------------------------------- c Compute square of separation. c----------------------------------------------------------- rsq = 0.0d0 do a = 1 , d rsq = rsq + (y(a,1,j) - y(a,1,i))**2 end do c1 = G / (rsq ** 1.5d0) c----------------------------------------------------------- c Update accns, making use of F_ji = - F_ij c----------------------------------------------------------- do a = 1 , d ydot(a,2,i) = ydot(a,2,i) + & c1 * m(j) * (y(a,1,j) - y(a,1,i)) ydot(a,2,j) = ydot(a,2,j) - & c1 * m(i) * (y(a,1,j) - y(a,1,i)) end do end do end do return end c=========================================================== c Dummy Jacobian routine. c=========================================================== subroutine jac implicit none include 'fcn.inc' return end c=========================================================== c Reads particle masses and initial positions and c velocities from file (including standard input). c c This version will work only for d = 3. c=========================================================== subroutine getivs(fname,m,y0,d,maxnp,np) implicit none integer indlnb, getu character*(*) fname integer maxnp, np, d real*8 m(maxnp), y0(d,2,maxnp) integer ustdin parameter ( ustdin = 5 ) real*8 mi, x(3), v(3) integer ufrom, rc, i np = 0 if( fname .eq. '-' ) then ufrom = ustdin else ufrom = getu() open(ufrom,file=fname(1:indlnb(fname)), & form='formatted',status='old',iostat=rc) if( rc .ne. 0 ) then write(0,*) 'dvfrom: Error opening ', & fname(1:indlnb(fname)) return end if end if 100 continue read(ufrom,*,iostat=rc,end=200) mi, x, v if( rc .eq. 0 ) then np = np + 1 if( np .gt. maxnp ) then write(0,*) 'getivs: Read maximum of ', & maxnp, ' from ', & fname(1:indlnb(fname)) np = maxnp go to 200 end if m(np) = mi do i = 1 , d y0(i,1,np) = x(i) y0(i,2,np) = v(i) end do end if go to 100 200 continue c----------------------------------------------------------- c If we are reading from a file, close the file. c This releases the unit number for subsequent use. c----------------------------------------------------------- if( ufrom .ne. ustdin ) then close(ufrom) end if return end c=========================================================== c Computes energy quantities: individual kinetic c and gravitational potential energies, and total KE c and PE and total mechanical energy. c=========================================================== subroutine clce(y,m,ke,pe,ketot,petot,etot,G,d,np) implicit none real*8 dvsum integer d, np real*8 y(d,2,np), m(np), ke(np), pe(np), & ketot, petot, etot, G real*8 vsq, rsq, c1 integer a, i, j do i =1 , np vsq = 0.0d0 do a = 1 , d vsq = vsq + y(a,2,i)**2 end do ke(i) = 0.5d0 * m(i) * vsq pe(i) = 0.0d0 c1 = -G * m(i) do j = 1 , np if( j .ne. i ) then rsq = 0.0d0 do a = 1 , d rsq = rsq + (y(a,1,j) - y(a,1,i))**2 end do c----------------------------------------------------------- c Associate 1/2 the potential energy of an c interaction with each particle. c----------------------------------------------------------- pe(i) = pe(i) + 0.5d0 * c1 * m(j) / sqrt(rsq) end if end do end do ketot = dvsum(ke,np) petot = dvsum(pe,np) etot = ketot + petot return end c=========================================================== c Computes individual linear and angular momentum about c specified origin and total linear and angular momentum c about said origin. Currently implemented only for c d = 3. c=========================================================== subroutine clcmom(y,m,pmom,jmom,ptot,jtot,o,d,np) implicit none integer d, np real*8 y(d,2,np), pmom(d,np), jmom(d,np), & m(np), ptot(d), jtot(d), & o(d) integer a, i call dvls(ptot,0.0d0,d) call dvls(jtot,0.0d0,d) if( d .ne. 3 ) then write(0,*) 'clcmom: Not implemented for d .ne. 3 ' write(0,*) 'clcmom: Invoked with d = ', d return end if c----------------------------------------------------------- c Linear momentum c----------------------------------------------------------- do i = 1 , np do a = 1 , d pmom(a,i) = m(i) * y(a,2,i) ptot(a) = ptot(a) + pmom(a,i) end do end do c----------------------------------------------------------- c Angular momentum c----------------------------------------------------------- do i = 1 , np jmom(1,i) = (y(2,1,i) - o(2)) * pmom(3,i) - & (y(3,1,i) - o(3)) * pmom(2,i) jmom(2,i) = (y(3,1,i) - o(3)) * pmom(1,i) - & (y(1,1,i) - o(1)) * pmom(3,i) jmom(3,i) = (y(1,1,i) - o(1)) * pmom(2,i) - & (y(2,1,i) - o(2)) * pmom(1,i) do a = 1 , d jtot(a) = jtot(a) + jmom(a,i) end do end do return end c=========================================================== c Computes center-of-mass of distribution of particles. c=========================================================== subroutine clccom(y,m,com,d,np) implicit none real*8 dvsum integer d, np real*8 y(d,2,np), m(np), com(d) real*8 mtotm1 integer a, i call dvls(com,0.0d0,d) do i = 1 , np do a = 1 , d com(a) = com(a) + m(i) * y(a,1,i) end do end do call dvsm(com,1.0d0/dvsum(m,np),com,np) return end c=========================================================== c Vector dot product. c=========================================================== double precision function dvvdot(v1,v2,n) implicit none real*8 v1(*), v2(*) integer i, n dvvdot = 0.0d0 do i = 1 , n dvvdot = dvvdot + v1(i) * v2(i) end do return end c=========================================================== c Vector sum. c=========================================================== double precision function dvsum(v,n) implicit none real*8 v(*) integer i, n if( n .gt. 0 ) then dvsum = v(1) do i = 2 , n dvsum = dvsum + v(i) end do end if return end c=========================================================== c Load vector with scalar. c=========================================================== subroutine dvls(v1,s1,n) implicit none real*8 v1(*) real*8 s1 integer i, n do i = 1 , n v1(i) = s1 end do return end c=========================================================== c Vector-scalar multiply c=========================================================== subroutine dvsm(v1,s1,v2,n) implicit none real*8 v1(*), v2(*) real*8 s1 integer i, n do i = 1 , n v2(i) = s1 * v1(i) end do return end