PROGRAM nbody

USE Mod_Orbits

IMPLICIT NONE

INTEGER, PARAMETER :: nbodies = 3
REAL*8, DIMENSION (nbodies) :: masses, xs, ys, vx, vy
REAL*8 :: dt, t, tmax, tdump, tolerance
INTEGER :: tdumps, n
TYPE(Orbits) :: orbit

!### put in winter values rotated to the x-axis ###
DATA masses /1.d0, 3.003d-6, 9.548d-4/
DATA xs /-4.97d-3, -1.02d0, 5.47638d0/
DATA ys /0.d0, 0.d0, 0.d0/
DATA vx /0.d0, 0.d0, 0.d0/

!### solve for the vy's ###
vy(1) = 0.d0
DO n = 2, nbodies

	vy(n) = DSIGN(xs(n),1.d0)*DSQRT(3.94784176044d1/DABS(xs(n)))

END DO

!### make an initial request at dt ###
DATA dt /1.d-8/
DATA tmax /12.9/
DATA tdumps /300/

!### define a tolerance for position ###
DATA tolerance /1.d-6/

!### construct the orbit simulator ###
CALL OrbitsConstruct(orbit, nbodies, dt, tolerance)

!### initialize the simulation ###
CALL OrbitsInit(orbit, masses, xs, ys, vx, vy)

!### we'll start at t = 0 ###
t     = 0.d0
tdump = tmax / tdumps

OPEN(UNIT=12,FILE="test.dat")

!### print out a startup message ###
PRINT '("Orbital simulation started at T = ",ES20.8," with a time step size DT = ",ES20.8)', t, dt
PRINT '("#######################")'

!### run the simulation ###
DO WHILE (t.LT.tmax)

	!### advance the system ###
	CALL OrbitAdaptStep(orbit)
	!PRINT '("The time is now T = ",ES20.8," with DT = ",ES20.8," and error = ",ES20.8)', t, orbit%dt, orbit%error

	!### show the results of the advance ###
	IF (t.GE.tdump) THEN
		PRINT '("The time is now T = ",ES20.8", with DT = ",ES20.8)', t, orbit%dt
		PRINT '("  ",6(ES16.4))', orbit%bodies(1)%x, orbit%bodies(1)%y, orbit%bodies(2)%x, orbit%bodies(2)%y, orbit%bodies(3)%x, orbit%bodies(3)%y
		WRITE(12,'("  ",6(ES16.4))') orbit%bodies(1)%x, orbit%bodies(1)%y, orbit%bodies(2)%x, orbit%bodies(2)%y, orbit%bodies(3)%x, orbit%bodies(3)%y

		tdump = tdump + tmax / tdumps
	END IF

	t = orbit%dt + t

END DO

CLOSE(12)

END PROGRAM