'Auroral Oval Visualisation

DECLARE SUB mag2geo (mlon!, mlat!, GLON!, GLAT!)
DECLARE SUB GEO2MAG (GLON!, GLAT!, mlon!, mlat!)
DECLARE FUNCTION ASIN! (x!)
DECLARE FUNCTION ACOS! (x!)

DIM clat(100), clon(100)        'arrays for country data

GOSUB Setups
GOSUB Getdata
GOSUB DrawEarth
GOSUB Grid
GOSUB Label
GOSUB Sun
GOSUB GSP                       'indicate geomagnetic south pole
GOSUB Aurora
DO                              'wait loop for exit command
LOOP WHILE INKEY$ <> esc$
END


Setups:
 pi = 3.14159
 SCREEN 9
 ar = .75               'aspect ratio
 xc = 320                'screen center in x co-ordinate
 yc = 175                'screen center in y co-ordinate
 rad = 180               'Earth radius in y pixels
 dr = pi / 180           'degrees to radians
 esc$ = CHR$(27)         'escape to exit the program
RETURN


Getdata:
 CLS
 PRINT "AURORAL OVAL PLOTTING PROGRAM"
 PRINT
 INPUT "Month "; month
 INPUT "UT "; UT
 INPUT "Kp "; Kp
 CLS
RETURN


DrawEarth:
 GOSUB Equator
 GOSUB Australia
 GOSUB SouthAmerica
 GOSUB Africa
 GOSUB NZ
 GOSUB Antarctica
RETURN


Equator:
 CIRCLE (xc, yc), rad, 1, , , ar
 PAINT (xc, yc), 1
RETURN


Australia:
 'plot Tasmania
 DATA -41, 146,-41,149,-43,147
 np = 3
 ccol = 2
 GOSUB PlotCountry
 lat = -41.5                       'paint Tasmania
 lon = 147.5
 GOSUB Convert
 PAINT (x, y), ccol
 'plot mainland
 DATA -22,114,-20,120,-14,127,-13,136,-17,135,-18,140,-11,142
 DATA -25,154,-39,147,-33,130,-35,117
 np = 11
 GOSUB PlotCountry
 lat = -20
 lon = 122
 GOSUB Convert
 PAINT (x, y), ccol
RETURN


PlotCountry:
 FOR i = 1 TO np
   READ clat(i), clon(i)
 NEXT i
 clat(np + 1) = clat(1)
 clon(np + 1) = clon(1)
 FOR i = 1 TO np
   lat = clat(i)
   lon = clon(i)
   GOSUB Convert
   x1 = x
   y1 = y
   lat = clat(i + 1)
   lon = clon(i + 1)
   GOSUB Convert
   x2 = x
   y2 = y
   LINE (x1, y1)-(x2, y2), ccol
 NEXT i
RETURN


Convert:
 r = (90 + lat) * rad / 90
 theta = lon * dr
 x = r * COS(theta) + xc
 y = yc + r * SIN(theta) * ar
RETURN


SouthAmerica:
 DATA -55,-70,-18,-70,-4,-81,0,-80,0,-75,0,-70,0,-65,0,-60,0,-55,0,-50,-7,-38
 np = 11
 ccol = 2
 GOSUB PlotCountry
 lat = -5
 lon = -70
 GOSUB Convert
 PAINT (x, y), ccol
RETURN


Africa:
 DATA 0,10,-34,20,-30,30,-4,39,0,42,0,40,0,35,0,30,0,25
 DATA 0,20,0,15
 np = 11
 ccol = 2
 GOSUB PlotCountry
 lat = -10
 lon = 20
 GOSUB Convert
 PAINT (x, y), ccol
RETURN


NZ:
 'north island
 DATA -35,173,-38,178,-41.8,175,-39,174
 np = 4
 ccol = 2
 GOSUB PlotCountry
 lat = -39
 lon = 176
 GOSUB Convert
 PAINT (x, y), ccol
 'south island
 DATA -41,172.5,-42,174,-46.5,169,-46,166.5
 np = 4
 GOSUB PlotCountry
 lat = -42
 lon = 173
 GOSUB Convert
 PAINT (x, y), ccol
RETURN


Antarctica:
 DATA -62,-60,-75,-60,-70,-20,-68,0,-68,30,-65,50,-66,90,-65,120
 DATA -65,150,-70,170,-80,170,-75,-140,-72,-110,-72,-80
 np = 14
 ccol = 7
 GOSUB PlotCountry
 lat = -89
 lon = -10
 GOSUB Convert
 PAINT (x, y), ccol
RETURN


Grid:
 'do latitude circles
 FOR lc = 20 TO 160 STEP 20
   CIRCLE (xc, yc), lc, 0, , , ar
 NEXT lc
 'do longitude lines
 lat = 0
 FOR lon = 0 TO 330 STEP 30
   GOSUB Convert
   LINE (xc, yc)-(x, y), 0
 NEXT lon
RETURN


Label:
 LINE (1, 7)-(68, 60), 14, B
 LOCATE 2, 2
 COLOR 12
 PRINT "AURORAL";
 LOCATE 3, 2
 PRINT " OVAL  ";
 LOCATE 4, 2
 PRINT " PLOT  ";
 COLOR 9
 LOCATE 2, 26
 PRINT "Southern Terrestrial Hemisphere";
 COLOR 7
 'label longitudinal lines
 LOCATE 13, 15
 PRINT "180";
 LOCATE 13, 64
 PRINT "0";
 LOCATE 3, 40
 PRINT "90W";
 LOCATE 23, 40
 PRINT "90E";
 'add data info
 LOCATE 25, 5
 PRINT USING "UT (hours) = ##.#"; UT;
 LOCATE 25, 55
 PRINT USING "Magnetic index Kp = #"; Kp;
 'add auroral legend
 COLOR 12
 LINE (620, 250)-(630, 300), 12, BF
 LOCATE 19, 70
 PRINT "=";
 LOCATE 20, 70
 PRINT "Area";
 LOCATE 21, 70
 PRINT "of";
 LOCATE 22, 70
 PRINT "Auroral";
 LOCATE 23, 70
 PRINT "Activity";
 'add visibility legend
 LINE (5, 250)-(5, 300), 4
 COLOR 4
 LOCATE 19, 3
 PRINT "=";
 LOCATE 20, 3
 PRINT "Northern";
 LOCATE 21, 3
 PRINT "Limit of";
 LOCATE 22, 3
 PRINT "Auroral";
 LOCATE 23, 3
 PRINT "Visibility";
RETURN


Sun:
 'compute centre position for sun
 lon = (12 - UT) * 15
 lat = 20
 GOSUB Convert
 CIRCLE (x, y), 10, 14
 PAINT (x, y), 14
 LINE (x, y - 15 * ar)-(x, y + 15 * ar), 14
 LINE (x - 15, y)-(x + 15, y), 14
 'draw solar legend and label
 x = 600: y = 20
 CIRCLE (x, y), 10, 14
 PAINT (x, y), 14
 LINE (x, y - 15 * ar)-(x, y + 15 * ar), 14
 LINE (x - 15, y)-(x + 15, y), 14
 LOCATE 2, 70
 COLOR 14
 PRINT "=";
 LOCATE 3, 70
 PRINT "Solar";
 LOCATE 4, 70
 PRINT "Position";
RETURN


GSP:
 'indicate geomagnetic south pole
 CALL mag2geo(0, -90, lon, lat)
 GSPlon = lon
 GSPlat = lat
 GOSUB Convert
 CIRCLE (x, y), 2, 15
 PAINT (x, y), 15
 'add legend for GSP
 CIRCLE (5, 76), 2, 15
 PAINT (5, 76), 15
 COLOR 15
 LOCATE 6, 3
 PRINT "=";
 LOCATE 7, 1
 PRINT "Geomagnetic";
 LOCATE 8, 1
 PRINT "South";
 LOCATE 9, 1
 PRINT "Pole";
RETURN


Aurora:
 'find subsolar point in geomagnetic coords
 YD = (month - 1) * 30                  'approximate day of the year
 solon = 15 * (12 - UT)                 'approx solar geolongitude
 solat = 23.5 * SIN((YD - 82) * dr)     'approx solar geolatitude
 CALL GEO2MAG(solon, solat, mslon, mslat)
 zhltmlon = mslon - 180                 'zero hr local time mag longitude
 'GOSUB SolMagMer                       'draw solar magnetic meridian if req
 lolat = -90
 hilat = 90
 Q = 1.33 * Kp + EXP(Kp - 8)
 'plot equatorward boundary
 ovalcolour = 12
 va = 0                                 'implies zenithal observation
 magcoef = .9
 timecoef = 5.1
 GOSUB OvalPlot
 'plot northern visibility boundary
 ovalcolour = 4
 va = 12.3                              'implies horizon for 150km altitude
 GOSUB OvalPlot
 'plot poleward boundary
 ovalcolour = 12
 va = 0
 magcoef = .3
 timecoef = timecoef * EXP(-Q / 4)
 GOSUB OvalPlot
 'paint auroral oval
 CALL mag2geo(lolon, (lolat + hilat) / 2, lon, lat)
 GOSUB Convert
 PAINT (x, y), 12
RETURN


OvalPlot:
 FOR mlon = 0 TO 360 STEP 10
   tCG = (mlon - zhltmlon) / 15
   mlat = -72 + magcoef * Q + timecoef * COS((360 / 24 * tCG - 12) * dr) + va
   CALL mag2geo(mlon, mlat, lon, lat)
   GOSUB Convert
   IF va = 0 AND mlat > lolat THEN
     lolat = mlat
     lolon = mlon
   END IF
   IF mlon = lolon THEN
     hilat = mlat
   END IF
   IF mlon = 0 THEN
     PSET (x, y), ovalcolour
   ELSE
     LINE -(x, y), ovalcolour
   END IF
 NEXT mlon
RETURN


SOLMAGMER:
 'draw solar magnetic meridian
 FOR mlat = -90 TO -20
   CALL mag2geo(mslon, mlat, lon, lat)  'dayside
   GOSUB Convert
   PSET (x, y), 14
   CALL mag2geo(zhltmlon, mlat, lon, lat)'nightside
   GOSUB Convert
   PSET (x, y), 8
 NEXT mlat
RETURN

FUNCTION ACOS (x)
 IF x = 0 THEN ACS = 1.570795 ELSE ACS = ATN(SQR(1 - x * x) / x)
 IF x < 0 THEN ACS = 3.14159 + ACS
 ACOS = ACS
END FUNCTION

FUNCTION ASIN (x)
 IF ABS(x) = 1 THEN
   ASIN = SGN(x) * 1.570795
 ELSE
   ASIN = ATN(x / SQR(1 - x * x))
 END IF
END FUNCTION

SUB GEO2MAG (GLON, GLAT, mlon, mlat)
' calculates geomagnetic coordinates (MLON,MLAT)      ie. outputs
' from geographic coordinates (GLON,GLAT)            ie. inputs
' all angles in degrees : latitude -90 to 90 : longitude 0 to 360 East
 pi = 3.14159: tpi = pi * 2: dr = pi / 180  'degrees to radians constant
 CBG = 10.7 * dr
 CI = COS(CBG)
 SI = SIN(CBG)
 YLG = GLON + 71.5
 CBG = COS(GLAT * dr)
 SBG = SIN(GLAT * dr)
 CLG = COS(YLG * dr)
 SLG = SIN(YLG * dr)
 SBM = SBG * CI + CBG * CLG * SI
 IF ABS(SBM) > 1! THEN SBM = SGN(SBM)
 mlat = ASIN(SBM)
 CBM = COS(mlat)
 SLM = (CBG * SLG) / CBM
 CLM = (-SBG * SI + CBG * CLG * CI) / CBM
 IF ABS(CLM) > 1! THEN CLM = SGN(CLM)
 mlon = ACOS(CLM)
 IF SLM < 0 THEN mlon = tpi - mlon
 mlat = mlat / dr
 mlon = mlon / dr
END SUB

SUB mag2geo (mlon, mlat, GLON, GLAT)
' calculates geographic coordinates (GLON,GLAT)      ie. outputs
' from geomagnetic coordinates (MLON,MLAT)           ie. inputs
' all angles in degrees : latitude -90 to 90 : longitude 0 to 360 East
 pi = 3.14159: tpi = pi * 2: dr = pi / 180      'degrees to radians constant
 CBG = 10.7 * dr                                'angle between mag & geo axes
 CI = COS(CBG)
 SI = SIN(CBG)
 CBM = COS(mlat * dr)
 SBM = SIN(mlat * dr)
 CLM = COS(mlon * dr)
 SLM = SIN(mlon * dr)
 SBG = SBM * CI - CBM * CLM * SI
 IF ABS(SBG) > 1! THEN SBG = SGN(SBG)
 GLAT = ASIN(SBG)
 CBG = COS(GLAT)
 SLG = (CBM * SLM) / CBG
 CLG = (SBM * SI + CBM * CLM * CI) / CBG
 IF ABS(CLG) > 1! THEN CLG = SGN(CLG)
 GLON = ACOS(CLG)
 IF SLG < 0 THEN GLON = tpi - GLON
 GLAT = GLAT / dr
 GLON = GLON / dr
 GLON = GLON - 71.5
 IF GLON < 0 THEN GLON = GLON + 360!
END SUB