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.
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
DO: LOOP: DO: LOOP
sha_na_na_na_na_na_na_na_na_na:
sha_na_na_na_na_na_na_na_na_na:
