Accretion Disk - Printable Version

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

' QB64 - RELATIVISTIC SINGULARITY (FINAL CINEMATIC VERSION) ' --------------------------------------------------------------------------- wd = 800: ht = 600 SCREEN _NEWIMAGE(wd, ht, 32) vecs = 18000 ' Max density for beauty DIM Px(vecs), Py(vecs), Vx(vecs), Vy(vecs) G = 2.5: M = 950.0: RS = 24.0 FOR i = 0 TO vecs a = RND * 6.283: d = RS + 40 + RND * 380 Px(i) = COS(a) * d: Py(i) = SIN(a) * d v_orb = SQR((G * M) / d) Vx(i) = -SIN(a) * v_orb * 0.99: Vy(i) = COS(a) * v_orb * 0.99 NEXT i DO CLS ' --- THE EVENT HORIZON --- FOR r_loop = RS TO 0 STEP -1: CIRCLE (wd / 2, ht / 2), r_loop, _RGB32(0, 0, 0): NEXT r_loop FOR i = 0 TO vecs r = _HYPOT(Px(i), Py(i)) ' --- CONSUMPTION & RESET --- resetted = 0 IF r < RS OR r > 580 THEN d = 400 + RND * 50: a = RND * 6.283 Px(i) = COS(a) * d: Py(i) = SIN(a) * d v = SQR((G * M) / d): Vx(i) = -SIN(a) * v: Vy(i) = COS(a) * v resetted = 1 END IF ' --- PHYSICS --- mag = (G * M) / (r * r) Vx(i) = Vx(i) - (Px(i) / r) * mag: Vy(i) = Vy(i) - (Py(i) / r) * mag Px(i) = Px(i) + Vx(i): Py(i) = Py(i) + Vy(i) IF resetted = 0 THEN ' --- DOPPLER / RELATIVISTIC BEAMING --- ' Particles moving LEFT (Vx < 0) are moving TOWARD the observer. ' We boost their brightness and add a "Blue" tint. doppler = 1.0 - (Vx(i) / 15.0) ' Boost if Vx is negative (moving left) ' --- HEAT GRADIENT --- heat_val = (RS * 18) / (r - RS + 5) * 15 * doppler IF heat_val > 255 THEN heat_val = 255 r_c = heat_val * 1.6: IF r_c > 255 THEN r_c = 255 g_c = heat_val * 0.9: IF g_c > 255 THEN g_c = 255 b_c = heat_val * (0.3 + (doppler * 0.2)): IF b_c > 255 THEN b_c = 255 clr = _RGB32(r_c, g_c, b_c) ' --- PROJECTIONS --- screenX = (wd / 2) + Px(i) screenY = (ht / 2) + (Py(i) * 0.22) ' Main Disk Render IF screenX > 0 AND screenX < wd AND screenY > 0 AND screenY < ht THEN PSET (screenX, screenY), clr END IF ' --- THE EINSTEIN RING (Lensed Light) --- ' Masked to prevent the "horizontal lines" glitch IF Py(i) < -2 AND ABS(Px(i)) < RS * 3.8 THEN distort = (RS * 30) / (ABS(Px(i)) + 20) warpY1 = (ht / 2) - (RS + distort) warpY2 = (ht / 2) + (RS + distort) alpha = 1.0 - (ABS(Px(i)) / (RS * 3.8)) IF alpha > 0 THEN l_clr = _RGB32(r_c * alpha, g_c * alpha, b_c * alpha) IF screenX > 0 AND screenX < wd THEN IF warpY1 > 0 AND warpY1 < ht THEN PSET (screenX, warpY1), l_clr IF warpY2 > 0 AND warpY2 < ht THEN PSET (screenX, warpY2), l_clr END IF END IF END IF END IF NEXT i CIRCLE (wd / 2, ht / 2), RS, _RGB32(255, 140, 50) _DISPLAY: _LIMIT 60 LOOP UNTIL _KEYDOWN(27)