Accretion Disk

[OldMoses]

My goal with this one was to simulate a black hole's accretion disk using my standard vector library tools. It will eventually settle down to the entropic disk that I was after, but I wasn't expecting the initial splash patterns. I like how the particles will arrange themselves into discrete orbit shells. That said, this is NOT a gravity simulation, as attraction does not vary with distance.

Use the space bar to display a green "event horizon", which shows how close the orbits skim it.

It's an amusement to change some of the parameters and see how it affects the pattern and/or number of survivors.

The number in the upper right corner is the number of surviving particles, the rest are lost, forever lost!

If you stare at it long enough the world will start spinning too. Maybe you might also begin to feel sleepy and have an irresistible urge to send me a sizeable check.

Code: (Select All)

OPTION _EXPLICIT

TYPE V2 '                                                       R2 vector format: used with R2_ subs & functions
    x AS SINGLE
    y AS SINGLE
END TYPE

TYPE particle
    ex AS INTEGER '                                             particle exists -1=yes 0=no
    p AS V2 '                                                   position
    d AS V2 '                                                   displacement aka movement
END TYPE

DIM c&, bs&, rs&
DIM vl%, vecs%, attrad%, wd%, ht%, eh%
DIM AS V2 tmp, op, cn
DIM speed!, attstr!

wd% = _DESKTOPWIDTH
ht% = _DESKTOPHEIGHT

SCREEN _NEWIMAGE(wd%, ht%, 32)
DO UNTIL _SCREENEXISTS: LOOP
_SCREENMOVE 0, 0
_TITLE "Accretion disk  {esc to quit}"

vecs% = 10000 '                                                 total starting particles
attrad% = 50 '                                                  radius of the attractor
attstr! = .1 '                                                  pull of the attractor
speed! = 1 '                                                    particle orbit velocity
bs& = &HFFAFAFFF '                                              particle color blue shift
rs& = &HFFFFAFAF '                                              particle color red shift

cn.x = _SHR(wd%, 1): cn.y = _SHR(ht%, 1) '                      screen center

'set initial particle data
RANDOMIZE TIMER
DIM SHARED V(vecs%) AS particle
FOR vl% = 0 TO UBOUND(V)
    V(vl%).ex = -1 '                                            particle exists at start
    DO
        V(vl%).p.x = INT(RND * wd%) '                           random position start
        V(vl%).p.y = INT(RND * wd%) - _SHR(wd% - ht%, 1)
        R2_P2P tmp, V(vl%).p, cn '                              get particle position vector <tmp> relative to screen center
    LOOP UNTIL R2_Mag(tmp) < _SHR(wd%, 1) '                     redo until particle within a circular area
    R2_Orth tmp, tmp, -1 '                                      convert position vector to an orthogonal rotation vector  - CCW  + CW
    R2_Norm V(vl%).d, tmp, speed! '                             set magnitude of rotation vector to speed! and assign it to particle displacement
NEXT vl%

'Main program loop
DO
    IF _KEYDOWN(32) THEN eh% = NOT eh%
    CLS
    _PRINTSTRING (0, 0), STR$(vecs%) '                          number of surviving particles
    IF eh% THEN CIRCLE (cn.x, cn.y), attrad%, &HFF00FF00 '      space bar to toggle the boundary of the attractor in green if desired
    FOR vl% = 0 TO UBOUND(V)
        IF NOT V(vl%).ex THEN _CONTINUE '                       skip particles that have been swallowed by attractor
        op = V(vl%).p
        R2_P2P tmp, op, cn '                                    obtain the inward pointing vector between particle and attractor
        c& = _IIF(R2_Dot(tmp, V(vl%).d) > 0, rs&, bs&) '        blue shift = approaching apogee, red shift = approaching perigee
        IF R2_Mag(tmp) < attrad% THEN '                         if particle is within "event horizon"
            V(vl%).ex = 0: vecs% = vecs% - 1 '                  particle has been absorbed by attractor
        ELSE
            R2_Norm tmp, tmp, attstr! '                         attractor force vector
            R2_Add V(vl%).d, tmp, 1 '                           add that to V(vl%).d
            R2_Add V(vl%).p, V(vl%).d, 1 '                      add displacement to position aka move it
            LINE (op.x, op.y)-(V(vl%).p.x, V(vl%).p.y), c& '    show the vector
        END IF
    NEXT vl%
    _LIMIT 30
    _DISPLAY
LOOP UNTIL _KEYDOWN(27)
END

'------------------------------------------------------------------------------
'SUBROUTINES
'------------------------------------------------------------------------------
SUB R2_P2P (re AS V2, st AS V2, nd AS V2)
    'ÞßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßÝ
    'Þ Return a 2D position vector <re> from point (st) to point (nd)          Ý
    'ÞÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÝ
    re.x = nd.x - st.x
    re.y = nd.y - st.y
END SUB 'R2_P2P

'------------------------------------------------------------------------------
SUB R2_Add (re AS V2, se AS V2, m AS SINGLE)
    'ÞßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßÝ
    'Þ Adds (or subtracts) a scalar multiple (m) of vector <se> to vector <re> Ý
    'Þ if <re> is to be preserved in its prior state then a copy should be     Ý
    'Þ made for the subroutine call. Send a negative scalar (m) to subtract    Ý
    'Þ Returns: vector <re>                                                    Ý
    'ÞÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÝ
    re.x = re.x + se.x * m
    re.y = re.y + se.y * m
END SUB 'R2_Add

'------------------------------------------------------------------------------
FUNCTION R2_Dot (v AS V2, v2 AS V2)
    'ÞßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßÝ
    'Þ Return scalar dot product between two V2 vectors <v> & <v2>            Ý
    'ÞÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÝ
    R2_Dot = v.x * v2.x + v.y * v2.y
END FUNCTION 'R3_Dot

'------------------------------------------------------------------------------
FUNCTION R2_Mag! (v AS V2)
    'ÞßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßÝ
    'Þ Return the scalar magnitude of 2D vector (v)                           Ý
    'ÞÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÝ
    R2_Mag! = _HYPOT(v.x, v.y)
END FUNCTION 'R2_Mag!

'------------------------------------------------------------------------------
SUB R2_Orth (re AS V2, v AS V2, d AS INTEGER)
    'ÞßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßÝ
    'Þ Returns vector <re> as a 90ø orthogonal of vector <v>. Sending a       Ý
    'Þ positive value for d rotates clockwise, while a negative value for     Ý
    'Þ d rotates counterclockwise. <re> & <v> can be the same.                Ý
    'ÞÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÝ
    DIM x!, y!
    x! = v.x: y! = v.y
    re.x = y! * SGN(d)
    re.y = -x! * SGN(d)
END SUB 'R2_Orth

'------------------------------------------------------------------------------
SUB R2_Norm (re AS V2, v AS V2, scalar AS SINGLE)
    'ÞßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßÝ
    'Þ Grow/Shrink V2 vector <v> to (scalar) length. <re> & <v> can be the    Ý
    'Þ same vectors, overwriting the original. Vector direction is preserved. Ý
    'Þ Returns: vector <re> as new length of <v>                              Ý
    'ÞÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÝ
    DIM m!
    DIM t AS V2 '                                               temporary vector
    t = v '                                                     V vector to temporary container
    m! = R2_Mag!(t) '                                           magnitude of temporary
    IF m! = 0 THEN
        re.x = 0: re.y = 0 '                                    zero vector- no direction information available
    ELSE
        re.x = (t.x / m!) * scalar '                            reduce to unit and expand to scalar
        re.y = (t.y / m!) * scalar '                            and assign to return vector
    END IF
END SUB 'R2_Norm