c write_energy.f c Write output file with total energy (in barycentric frame). c Miroslav Broz (miroslav.broz@email.cz), Dec 11th 2008 subroutine write_energy(n, t, Gm, r, v, outfile) implicit none integer n, i, j, ierr double precision t, Gm(n), r(n), v(n), E character*80 outfile open(unit=10, file=outfile, form="formatted", access="append", : iostat=ierr) if (ierr.ne.0) then write(*,*) "Error opening output file ", outfile stop endif E = 0.d0 do i = 0, n-1 c kinetic energy E = E + 0.5d0 * Gm(i+1) * ( : v(3*i+1)**2 + : v(3*i+2)**2 + : v(3*i+3)**2 ) c gravitational potential energy do j = i+1, n-1 E = E - Gm(i+1)*Gm(j+1)/sqrt( : (r(3*i+1)-r(3*j+1))**2 + : (r(3*i+2)-r(3*j+2))**2 + : (r(3*i+3)-r(3*j+3))**2 ) enddo enddo write(10,10,err=990) t/365.25d6, E 10 format(2(1x,e22.16)) close(10) return 990 write(*,*) "Error writing output to ", outfile stop end