52b25350e1
* NEW: New N-Poly benchmarks. * NEW: All benchmarks times were updated on the wiki. git-svn-id: svn://localhost/gambas/trunk@4582 867c0c6c-44f3-4631-809d-bfa615b0a4ec
205 lines
4.1 KiB
Text
Executable file
205 lines
4.1 KiB
Text
Executable file
#!/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()
|
|
|