Gravity

GRAVITY.BAS

'A gravity simulation of the solar system written in QuickBASIC
'Similar to something I wrote way back when I was still in school
'Units are km, kg and s

'Type declarations
TYPE vector
    x AS DOUBLE
    y AS DOUBLE
    z AS DOUBLE
END TYPE

TYPE object
    col AS INTEGER 'color
    m AS DOUBLE 'mass
    r AS vector 'position
    v AS vector 'velocity
END TYPE

'Subroutine declarations
DECLARE SUB advance (t AS DOUBLE)
DECLARE SUB drawsystem ()
DECLARE SUB vabs (s AS DOUBLE, v AS vector)
DECLARE SUB vadd (v1 AS vector, v2 AS vector, v3 AS vector)
DECLARE SUB vdiv (v1 AS vector, v2 AS vector, s AS DOUBLE)
DECLARE SUB vmul (v1 AS vector, v2 AS vector, s AS DOUBLE)
DECLARE SUB vsub (v1 AS vector, v2 AS vector, v3 AS vector)

'Solar system data
DIM SHARED count AS LONG
DIM SHARED objects(1 TO 1000) AS object
'Sun
count = count + 1
objects(count).col = 14
objects(count).m = 1.9885D+30
objects(count).r.x = -1067706.805380953#
objects(count).r.y = -418275.271819446#
objects(count).r.z = 30861.8172547681#
objects(count).v.x = 9.31257192652046D-03
objects(count).v.y = -1.282475570794162D-02
objects(count).v.z = -1.633507186350417D-04
'Mercury
count = count + 1
objects(count).col = 7
objects(count).m = 3.302D+23
objects(count).r.x = -20529433.1612346#
objects(count).r.y = -67331550.5353431#
objects(count).r.z = -3648992.526494758#
objects(count).v.x = 37.00430442920562#
objects(count).v.y = -11.17724068132644#
objects(count).v.z = -4.307791469376847#
'Venus
count = count + 1
objects(count).col = 15
objects(count).m = 4.868499999999993D+24
objects(count).r.x = -108524200.8575712#
objects(count).r.y = -5303290.24769197#
objects(count).r.z = 6166496.116973154#
objects(count).v.x = 1.391218601189967#
objects(count).v.y = -35.15311993215458#
objects(count).v.z = -.5602056890007152#
'Earth
count = count + 1
objects(count).col = 9
objects(count).m = 5.972189999999991D+24
objects(count).r.x = -27566740.48281139#
objects(count).r.y = 144279021.5207299#
objects(count).r.z = 30250.66782881313#
objects(count).v.x = -29.78494749851082#
objects(count).v.y = -5.482119695478537#
objects(count).v.z = 1.843295986780902D-05
'Mars
count = count + 1
objects(count).col = 4
objects(count).m = 6.4171D+23
objects(count).r.x = 206980433.8364605#
objects(count).r.y = -2425327.899844665#
objects(count).r.z = -5125427.142013245#
objects(count).v.x = 1.171984975692042#
objects(count).v.y = 26.28323978975659#
objects(count).v.z = .5221336722760512#
'Jupiter
count = count + 1
objects(count).col = 6
objects(count).m = 1.89818722D+27
objects(count).r.x = 597499917.8516816#
objects(count).r.y = 439186404.6763522#
objects(count).r.z = -15195999.85573265#
objects(count).v.x = -7.900547720245474#
objects(count).v.y = 11.14339277065934#
objects(count).v.z = .1307023308637314#
'Saturn
count = count + 1
objects(count).col = 6
objects(count).m = 5.683399999999991D+26
objects(count).r.x = 957317652.1103378#
objects(count).r.y = 982438007.6875058#
objects(count).r.z = -55182117.8815002#
objects(count).v.x = -7.421900386838106#
objects(count).v.y = 6.723930997200821#
objects(count).v.z = .1775749426205731#
'Uranus
count = count + 1
objects(count).col = 9
objects(count).m = 8.681299999999986D+25
objects(count).r.x = 2157907331.374323#
objects(count).r.y = -2055043522.898539#
objects(count).r.z = -35594601.96113003#
objects(count).v.x = 4.64658636945848#
objects(count).v.y = 4.614774391559485#
objects(count).v.z = -4.308124107692005D-02
'Neptune
count = count + 1
objects(count).col = 9
objects(count).m = 1.02409D+26
objects(count).r.x = 2513978490.096805#
objects(count).r.y = -3739132780.868785#
objects(count).r.z = 19063301.32971376#
objects(count).v.x = 4.47485846545944#
objects(count).v.y = 3.06388160579694#
objects(count).v.z = -.1659044011083587#
'Pluto (was still a planet back then)
count = count + 1
objects(count).col = 7
objects(count).m = 1.307D+22
objects(count).r.x = -1478398655.39348#
objects(count).r.y = -4182993264.891326#
objects(count).r.z = 875246398.9486204#
objects(count).v.x = 5.269124016486745#
objects(count).v.y = -2.669250607340872#
objects(count).v.z = -1.250716402210748#

'Start the universe ;-)
SCREEN 12
FOR i = 1 TO 100000000
    drawsystem
    advance 60
NEXT i
PRINT "Press any key to continue"
SLEEP

SUB advance (t AS DOUBLE)
    DIM g AS DOUBLE 'constant of gravitation
    DIM a(1 TO count) AS vector 'accelerations
    g = 6.67429999999999D-20
    FOR i = 1 TO count
        FOR j = i + 1 TO count
            DIM r AS vector
            DIM rabs AS DOUBLE
            DIM f AS vector
            vsub r, objects(j).r, objects(i).r
            vabs rabs, r
            vmul f, r, g * objects(i).m * objects(j).m / rabs ^ 3
            vadd a(i), a(i), f
            vsub a(j), a(j), f
        NEXT j
        DIM temp AS vector
        vdiv a(i), a(i), objects(i).m
        vmul temp, a(i), .5 * t ^ 2
        vadd objects(i).r, objects(i).r, temp
        vmul temp, objects(i).v, t
        vadd objects(i).r, objects(i).r, temp
        vmul temp, a(i), t
        vadd objects(i).v, objects(i).v, temp
    NEXT i
END SUB

SUB drawsystem
    DIM scale AS DOUBLE
    scale = 2.999999999999999D-08
    FOR i = 1 TO count
        DIM x AS DOUBLE, y AS DOUBLE
        x = 320 + scale * objects(i).r.x
        y = 240 - scale * objects(i).r.y
        PSET (x, y), objects(i).col
    NEXT i
END SUB

SUB vabs (s AS DOUBLE, v AS vector)
    s = SQR(v.x ^ 2 + v.y ^ 2 + v.z ^ 2)
END SUB

SUB vadd (v1 AS vector, v2 AS vector, v3 AS vector)
    v1.x = v2.x + v3.x
    v1.y = v2.y + v3.y
    v1.z = v2.z + v3.z
END SUB

SUB vdiv (v1 AS vector, v2 AS vector, s AS DOUBLE)
    v1.x = v2.x / s
    v1.y = v2.y / s
    v1.z = v2.z / s
END SUB

SUB vmul (v1 AS vector, v2 AS vector, s AS DOUBLE)
    v1.x = v2.x * s
    v1.y = v2.y * s
    v1.z = v2.z * s
END SUB

SUB vsub (v1 AS vector, v2 AS vector, v3 AS vector)
    v1.x = v2.x - v3.x
    v1.y = v2.y - v3.y
    v1.z = v2.z - v3.z
END SUB