The program is written in Fortran 90.
The real numbers are all 8-byte reals, and this is expected to be necessary in a simulation where the gravitational fields differ by many orders of magnitude. Distances are measured in Astronomical Units (AU), the mean distance of the Earth from the Sun.
! ss_9b.f90 16-Sep-13 John Collins ! 16-Sep-13 John Collins PROGRAM ss_9_body IMPLICIT NONE INTEGER,PARAMETER :: kr = KIND(1.0D0) ! Real kind INTEGER,PARAMETER :: frames_per_day = 32 INTEGER,PARAMETER :: sample_interval = 30*frames_per_day ! 30 days INTEGER,PARAMETER :: end_frame = (10*365+2)*frames_per_day ! Note: We measure time in days, distance in AU, mass in kg REAL(KIND=kr),PARAMETER :: au = 149597870700.0D0 REAL(KIND=kr),PARAMETER :: gu = -6.67300D-11 ! m**3/kg/s**2 REAL(KIND=kr),PARAMETER :: g = ((60.0D0*60.0D0*24.0D0)**2 / au**3) *gu ! au**3/kg/day**2 REAL(KIND=kr) & time & , p_time & ! Previous time , dt & , end_time & , mass(9) & , s(1:3,1:9) & ! Position (axis,body) , ds(1:3,1:9,1:9) & ! Difference in position, (axis,body1,body2) , s2(1:3) & ! Square of difference in position, body pair , r2(1:9,1:9) & ! Square of radial distance, body to body , r(1:9,1:9) & ! Radial distances, body to body , gf(1:9,1:9) & ! Total gravitational field, body to body , proj(1:3,1:9,1:9) & ! Cosine projections on axes (axis,body1,body2) , vf(1:3,1:9,1:9) & ! Vector field components (axis,body1,body2) , a(1:3,1:9) & ! Accelerations (axis,body) , pv(1:3,1:9) & ! Previous velocities (axis,body) , v(1:3,1:9) & ! Velocities (axis,body) , pa(1:3,1:9) ! Previous accelerations (axis,body) INTEGER & axis & , body1 & , body2 & , frame LOGICAL end_run DATA & mass / 1.988544D+30 & ! Sun , 3.302D+23 & ! Mercury , 48.685D+23 & ! Venus , 59.7219D+23 & ! Earth , 6.4185D+23 & ! Mars , 18981.3D+23 & ! Jupiter , 5683.19D+23 & ! Saturn , 868.103D+23 & ! Uranus , 1024.1D+23 / & ! Neptune ! Initialisation data from JPL, California Institute of Technology, http://ssd.jpl.nasa.gov/horizons.cgi#top ! Positions in AU , s / -3.747147775505799D-03, 2.926587035303590D-03, 4.446807296978120D-06 & ! Sun , 4.715831111175716D-02, 3.056015767171350D-01, 2.006272440161816D-02 & ! Mercury , 4.959853421725741D-02,-7.222408775972484D-01,-1.300479358383397D-02 & ! Venus ,-1.797649390463135D-01, 9.703470886653921D-01,-1.743826995752289D-05 & ! Earth ,-7.334184274376777D-01, 1.457212408029754D+00, 4.839432410544958D-02 & ! Mars , 4.505316309281244D+00,-2.163555523180046D+00,-9.190499900147728D-02 & ! Jupiter ,-9.468003745121990D+00, 2.616109582176896D-01, 3.722149420390660D-01 & ! Saturn , 2.003287446208430D+01,-1.529803159589830D+00,-2.652256724942619D-01 & ! Uranus , 2.481340360454546D+01,-1.689465149553225D+01,-2.239358924631250D-01/ & ! Neptune ! Velocities in AU per day , v / -2.990258994506073D-06,-5.530878190433255D-06, 6.981801439481427D-08 & ! Sun ,-3.338934398882634D-02, 5.702520017360675D-03, 3.529925953990557D-03 & ! Mercury , 2.003227503747324D-02, 1.405250683890825D-03,-1.136887572032341D-03 & ! Venus ,-1.720173260060367D-02,-3.148295197112937D-03, 8.961124784916037D-07 & ! Earth ,-1.198027375283752D-02,-5.092533276958321D-03, 1.875893458494313D-04 & ! Mars , 3.174273016845083D-03, 7.161028006264430D-03,-1.007937605937214D-04 & ! Jupiter ,-4.520953258457782D-04,-5.588245332485017D-03, 1.148915369554184D-04 & ! Saturn , 2.706976434750568D-04, 3.738322243478300D-03, 1.038310207754998D-05 & ! Uranus , 1.746288049471184D-03, 2.613667805270097D-03,-9.397429335776478D-05/ ! Naptune frame = 0 time = 0.0D0 p_time = 0.0D0 1000 CONTINUE IF (mod(frame,sample_interval) == 0) THEN WRITE (*,'(F8.2,18E21.12)')time,((s(axis,body1),axis=1,2),body1=1,9) ENDIF frame = frame+1 ! Note that time will jitter with finite word length. Therefore: time = REAL(frame)/REAL(frames_per_day) dt = time - p_time p_time = time DO body1=1,9 DO body2=body1+1,9 DO axis=1,3 ! Difference in position ds(axis,body1,body2) = s(axis,body1)-s(axis,body2) ! Square of difference in position s2(axis) = ds(axis,body1,body2) * ds(axis,body1,body2) ENDDO ! Square of radial distance r2(body1,body2) = s2(1) + s2(2) + s2(3) ! Radial distance r(body1,body2) = sqrt(r2(body1,body2)) ! Total gravitation field from body2 gf(body1,body2) = g * mass(body2) / r2(body1,body2) ! Total gravitational field from body1 gf(body2,body1) = g * mass(body1) / r2(body1,body2) ! Cosines and vector fields to and from each body DO axis=1,3 proj(axis,body1,body2) = ds(axis,body1,body2) / r(body1,body2) proj(axis,body2,body1) = -proj(axis,body1,body2) vf(axis,body1,body2) = proj(axis,body1,body2) * gf(body1,body2) vf(axis,body2,body1) = proj(axis,body2,body1) * gf(body2,body1) ENDDO ENDDO ENDDO ! Sum accelerations and integrate DO axis=1,3 ! Total accelerations of each body a(axis,1)= vf(axis,1,2)+vf(axis,1,3)+vf(axis,1,4)+vf(axis,1,5)+vf(axis,1,6)+vf(axis,1,7)+vf(axis,1,8)+vf(axis,1,9) a(axis,2)=vf(axis,2,1) +vf(axis,2,3)+vf(axis,2,4)+vf(axis,2,5)+vf(axis,2,6)+vf(axis,2,7)+vf(axis,2,8)+vf(axis,2,9) a(axis,3)=vf(axis,3,1)+vf(axis,3,2) +vf(axis,3,4)+vf(axis,3,5)+vf(axis,3,6)+vf(axis,3,7)+vf(axis,3,8)+vf(axis,3,9) a(axis,4)=vf(axis,4,1)+vf(axis,4,2)+vf(axis,4,3) +vf(axis,4,5)+vf(axis,4,6)+vf(axis,4,7)+vf(axis,4,8)+vf(axis,4,9) a(axis,5)=vf(axis,5,1)+vf(axis,5,2)+vf(axis,5,3)+vf(axis,5,4) +vf(axis,5,6)+vf(axis,5,7)+vf(axis,5,8)+vf(axis,5,9) a(axis,6)=vf(axis,6,1)+vf(axis,6,2)+vf(axis,6,3)+vf(axis,6,4)+vf(axis,6,5) +vf(axis,6,7)+vf(axis,6,8)+vf(axis,6,9) a(axis,7)=vf(axis,7,1)+vf(axis,7,2)+vf(axis,7,3)+vf(axis,7,4)+vf(axis,7,5)+vf(axis,7,6) +vf(axis,7,8)+vf(axis,7,9) a(axis,8)=vf(axis,8,1)+vf(axis,8,2)+vf(axis,8,3)+vf(axis,8,4)+vf(axis,8,5)+vf(axis,8,6)+vf(axis,8,7) +vf(axis,8,9) a(axis,9)=vf(axis,9,1)+vf(axis,9,2)+vf(axis,9,3)+vf(axis,9,4)+vf(axis,9,5)+vf(axis,9,6)+vf(axis,9,7)+vf(axis,9,8) ! Integrate for each body DO body1=1,9 IF (frame == 1) THEN ! For the first frame we make the previous acceleration equal to the current. This causes an Euler step for ! the velocity integration in this frame. pa(axis,body1) = a(axis,body1) ENDIF ! Keep previous velocities for Tustin's method pv(axis,body1)=v(axis,body1) ! Trapezoidal integration for velocities v(axis,body1) = v(axis,body1) + ( 3*a(axis,body1) - pa(axis,body1) ) * dt * 0.5D0 ! Tustin's method for positions s(axis,body1) = s(axis,body1) + 0.5D0 * ( v(axis,body1) + pv(axis,body1) ) * dt ! Keep previous accelerations for next frame pa(axis,body1) = a(axis,body1) ENDDO ENDDO end_run = frame >= end_frame IF (.NOT. end_run) GOTO 1000 WRITE (*,'(/,F8.2,18E21.12)')time,((s(axis,body1),axis=1,2),body1=1,9) END PROGRAM