#!/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 ix, iy, iz 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
      ix = iBody.x
      iy = iBody.y
      iz = iBody.z

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

        dx = ix - jBody.x
        dy = iy - jBody.y
        dz = iz - 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 I = 0 To Bodies.Max
      iBody = Bodies[i]
      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 5

  Print S.Energy()

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

Next

Print S.Energy()

Error CStr(Jit.Time)