*ESC OFF
      *FLOAT 64
      MODE 8 : OFF
      ORIGIN @vdu%!208, @vdu%!212

      ON ERROR OSCLI "REFRESH ON":ON:CLS:PRINT'" ";:REPORT:PRINT" at line ";ERL:VDU7:END

      MinRadials% = 3
      MaxRadials% = 15
      QuadsPerRadial% = 10
      MaxPts% = 4 * QuadsPerRadial% * MaxRadials%
      DIM p(MaxPts%-1,2), q(MaxPts%-1,2), m(2,2)
      DIM col{(QuadsPerRadial%-1) r%, g%, b%}

      rho% = 700 : REM added to all rotated Z co-ordinates

      REM Time-dependent scale function parameters:
      MaxScale% = 1000
      t1% = 50  : REM object reaches full scale at time = t1%
      t2% = 750 : REM object remains at full scale while t1% <= time < t2%
      t3% = 800 : REM scale diminished to zero when t2% <= time < t3%
      c = MaxScale% - (MaxScale%/(t3%-t2%))*t2%

      newObj% = TRUE : REM Create a new object when this flag is TRUE

      *REFRESH OFF

      R% = RND(-TIME)

      time0% = TIME

      REPEAT
  
        CLS
  
        IF newObj% THEN
          numPts% = FNBuildNewObject
          newObj% = FALSE
        ENDIF
  
        T% = TIME
  
        dt% = T% - time0%
  
        REM A time-dependent compound function to control the scale of the object:
        IF dt% < t1% THEN scale = (MaxScale%/t1%)*dt%
        IF dt% >= t1% AND dt% < t2% THEN scale = MaxScale%
        IF dt% >= t2% AND dt% < t3% THEN scale = MaxScale% - (dt%-t2%)*MaxScale%/(t3%-t2%)
        IF dt% >= t3% THEN
          time0% = T%
          scale = 0
          newObj% = TRUE
        ENDIF
  
        REM Define time-dependent rotation angles:
        A = T%/122
        B = T%/176
        C = T%/100
        D = 2*PI*SIN(T%/520)
        E = 3*PI*COS(T%/720)
        sA = SIN(A)
        cA = COS(A)
        sB = SIN(B)
        cB = COS(B)
        sC = SIN(C)
        cC = COS(C)
  
        REM Rotate all of the object's points:
        m() = (cB*cC), -cB*sC, sB,  cA*sC+sA*sB*cC, cA*cC-sA*sB*sC, -sA*cB, sA*sC-cA*sB*cC, sA*cC+cA*sB*sC, cA*cB
        q() = p().m()
  
        REM Render all the polygons (quads) in the current object:
        colIndex% = 0
        n% = 0
        FOR I% = 0 TO numPts%-1 STEP 4
          IF n% = 0 THEN
            COLOUR 7, col{(colIndex%)}.r%, col{(colIndex%)}.g%, col{(colIndex%)}.b%
          ENDIF
          z0 = scale / (rho%+q(I%,2))
          z1 = scale / (rho%+q(I%+1,2))
          z2 = scale / (rho%+q(I%+2,2))
          z3 = scale / (rho%+q(I%+3,2))
          MOVE z0*q(I%,0), z0*q(I%,1)
          MOVE z1*q(I%+1,0), z1*q(I%+1,1)
          PLOT 85, z2*q(I%+2,0), z2*q(I%+2,1)
          PLOT 85, z3*q(I%+3,0), z3*q(I%+3,1)
          n% += 1
          IF n% = NumRadials% THEN
            colIndex% += 1
            n% = 0
          ENDIF
        NEXT
  
        *REFRESH
        WAIT 1
  
      UNTIL FALSE
      END

      DEF FNBuildNewObject
      LOCAL I%, J%, K%
      LOCAL r1, r2, r3, a, da, x1, y1, x2, y2, x3, y3, x4, y4
      NumRadials% = MinRadials% + RND(MaxRadials% - MinRadials% )
      FOR I% = 0 TO QuadsPerRadial%-1
        col{( I% )}.r% = 150 + RND(105)
        col{( I% )}.g% = col{(I%)}.r% + 1 - RND(10)
        col{( I% )}.b% = 220 + RND(35)
      NEXT I%
      FOR J% = 1 TO QuadsPerRadial%
        r1 = 256 * RND(1)
        r2 = r1 + 30 + 50*RND(1)
        r3 = r2 + 40 + 100*RND(1)
        da = 3 + 12*RND(1)
        FOR I% = 0 TO NumRadials%-1
          a = I% * (360/NumRadials%)
          x1 = r1*SINRAD(a)
          y1 = r1*COSRAD(a)
          x2 = r2*SINRAD(a-da)
          y2 = r2*COSRAD(a-da)
          x3 = r2*SINRAD(a+da)
          y3 = r2*COSRAD(a+da)
          x4 = r3*SINRAD(a)
          y4 = r3*COSRAD(a)
          p(K%,0) = x1
          p(K%,1) = y1
          p(K%+1,0) = x2
          p(K%+1,1) = y2
          p(K%+2,0) = x3
          p(K%+2,1) = y3
          p(K%+3,0) = x4
          p(K%+3,1) = y4
          K% += 4
        NEXT I%
      NEXT J%
      = K%