#!/usr/bin/env gbs3

Class Body

  Static Private SOLAR_MASS As Float = 4 * Pi * Pi
  Private Const DAYS_PER_YEAR As Float = 365.24
  
  Public X As Float
  Public Y As Float
  Public Z As Float
  Public VX As Float
  Public VY As Float
  Public VZ As Float
  Public Mass As Float

  Static Public Sub Jupiter() As Body

    Dim P As New Body
    p.x = 4.84143144246472090e+00
    p.y = -1.16032004402742839e+00
    p.z = -1.03622044471123109e-01
    p.vx = 1.66007664274403694e-03 * DAYS_PER_YEAR
    p.vy = 7.69901118419740425e-03 * DAYS_PER_YEAR
    p.vz = -6.90460016972063023e-05 * DAYS_PER_YEAR
    p.mass = 9.54791938424326609e-04 * SOLAR_MASS
    return p

  End

  Static Public Sub Saturn() As Body

    Dim P As New Body
    p.x = 8.34336671824457987e+00
    p.y = 4.12479856412430479e+00
    p.z = -4.03523417114321381e-01
    p.vx = -2.76742510726862411e-03 * DAYS_PER_YEAR
    p.vy = 4.99852801234917238e-03 * DAYS_PER_YEAR
    p.vz = 2.30417297573763929e-05 * DAYS_PER_YEAR
    p.mass = 2.85885980666130812e-04 * SOLAR_MASS
    return p

  End

  Static Public Sub Uranus() As Body

    Dim P As New Body

    p.x = 1.28943695621391310e+01
    p.y = -1.51111514016986312e+01
    p.z = -2.23307578892655734e-01
    p.vx = 2.96460137564761618e-03 * DAYS_PER_YEAR
    p.vy = 2.37847173959480950e-03 * DAYS_PER_YEAR
    p.vz = -2.96589568540237556e-05 * DAYS_PER_YEAR
    p.mass = 4.36624404335156298e-05 * SOLAR_MASS
    return p

  End

  Static Public Sub Neptune() As Body

    Dim P As New Body

    p.x = 1.53796971148509165e+01
    p.y = -2.59193146099879641e+01
    p.z = 1.79258772950371181e-01
    p.vx = 2.68067772490389322e-03 * DAYS_PER_YEAR
    p.vy = 1.62824170038242295e-03 * DAYS_PER_YEAR
    p.vz = -9.51592254519715870e-05 * DAYS_PER_YEAR
    p.mass = 5.15138902046611451e-05 * SOLAR_MASS
    return p

  End

  Static Public Sub Sun() As Body

    Dim P As New Body

    p.mass = SOLAR_MASS
    return p

  End

  Public Sub OffsetMomentum(px As Float, py As Float, pz As Float) As Body

    vx = -px / SOLAR_MASS
    vy = -py / SOLAR_MASS
    vz = -pz / SOLAR_MASS

    Return Me

  End

End Class


Class NBodySystem

  Private Bodies As Body[]

  Public Sub _new()

    Dim PX, PY, PZ As Float
    Dim I As Integer

    Bodies = [ Body.Sun(), Body.Jupiter(), Body.Saturn(), Body.Uranus(), Body.Neptune() ]

    For I = 0 To Bodies.Max
      PX += Bodies[I].vx * Bodies[I].mass
      PY += Bodies[I].vy * Bodies[I].mass
      PZ += Bodies[I].vz * Bodies[I].mass
    Next

    Bodies[0].offsetMomentum(PX, PY, PZ)

  End

  Public Sub Advance(dt As Float)

    Dim I, J As Integer
    Dim iBody, jBody As Body
    Dim dx, dy, dz As Float
    Dim dSquared, fMag As Float
    Dim iMass, jMass, iMag, jMag As Float

    For I = 0 To Bodies.Max
      iBody = Bodies[I]
      iMass = iBody.mass

      For J = I + 1 To Bodies.Max
        jBody = Bodies[J]
	jMass = jBody.mass

        dx = iBody.x - jBody.x
        dy = iBody.y - jBody.y
        dz = iBody.z - jBody.z

        dSquared = dx * dx + dy * dy + dz * dz
        fMag = dt / (dSquared * Sqr(dSquared))
	iMag = iMass * fMag
	jMag = jMass * fMag

        iBody.vx -= dx * jMag
	iBody.vy -= dy * jMag
	iBody.vz -= dz * jMag

	jBody.vx += dx * iMag
	jBody.vy += dy * iMag
	jBody.vz += dz * iMag
      Next
    Next

    For Each iBody in Bodies
      iBody.x += dt * iBody.vx
      iBody.y += dt * iBody.vy
      iBody.z += dt * iBody.vz
    Next

  End

  Public Sub Energy() As Float

    Dim dx, dy, dz, distance, E As Float
    Dim iBody, jBody As Body
    Dim I, J As Integer

    For I = 0 To Bodies.Max

      iBody = bodies[i]
      E += 0.5 * iBody.mass * (iBody.vx * iBody.vx + iBody.vy * iBody.vy + iBody.vz * iBody.vz)

      For J = I + 1 To Bodies.Max

	jBody = Bodies[J]
        dx = iBody.x - jBody.x
        dy = iBody.y - jBody.y
        dz = iBody.z - jBody.z

        distance = Sqr(dx*dx + dy*dy + dz*dz)
        E -= (iBody.mass * jBody.mass) / distance

      Next
    Next
    
    return E

  End

End Class


Dim S As New NBodySystem
Dim I, N As Integer

For N = 1 To 10

  Print S.Energy()

  For I = 1 To 100000
    S.Advance(0.01)
  Next

Next

Print S.Energy()