'Orbital Space Debris - Simplified Model for Skyglow [File OSDMOD-E.BAS]
  'NOTE:  This version allows debris-debris collisions
CLS             'clear screen
OPEN "OSDsky.txt" FOR OUTPUT AS #1      'set print device or file
ufac = 1E+12 / 5  'microdebris factor (from inverse square law)
A1000 = 6.8E+14 'spherical area at 1000 km altitude
isat = 2000     'initial number of active satellites
idebris = 50000 'initial number of debris
t0 = 2000         'start year
A = 100         'annual insertion rate less reentries
chi = 3E-10     'annual collision rate
alpha = 2000    'collision generation coeff
beta = 50       'primary debris coeff (explosions,etc)
epsilon = 1000  'debris-debris collision generation coeff
delta = .0001    'debris removal coefficient
D = .0001       'debris area               
Nsat = isat
Ndebris = idebris
PRINT #1, "SIMPLE ORBITAL SPACE DEBRIS MODEL for SKYGLOW"
PRINT #1, "Date   Sats   Debris   Microdebris  Areal Density"      'heading
PRINT #1, "(yr)   (num)  (num)    (num)        (num/m^2@1000km)"
 f$ = "####  #####  ##.#^^^^   ##.#^^^^       ##.#^^^^ "       'print format
FOR tm = t0 TO t0 + 700      'cycle one year at a time
  IF tm MOD 50 = 0 THEN
    udebris = Ndebris * ufac
    udebdens = udebris / A1000
    PRINT #1, USING f$; tm; Nsat; Ndebris; udebris; udebdens
  END IF
  colrat = chi * Ndebris * Nsat        'collision rate
  dsat = A - colrat             'increase in satellites
  ddebris = beta * A + alpha * colrat - delta * Ndebris  'debris increase
  ddebris = ddebris + epsilon * chi * D * Ndebris ^ 2  'debris-debris collisions added
  Nsat = Nsat + dsat            'new number of satellites
  Ndebris = Ndebris + ddebris   'new debris count
NEXT tm

PRINT #1, "Model Parameters"
PRINT #1, USING "     Initial number of active satellites = #####"; isat
PRINT #1, USING "     Initial number of debris pieces     = ###,###"; idebris
PRINT #1, USING "     Annual insertion rate less reentries = ####"; A
PRINT #1, USING "     Annual collision rate factor = ##.##^^^^"; chi
PRINT #1, USING "     Collision generation coefficient = #####"; alpha
PRINT #1, USING "     Primary debris coefficient = #####"; beta
PRINT #1, "Simulation end"

DO WHILE INKEY$ <> CHR$(27)
LOOP
END