'Short Wave Fade Plot  [SWF.BAS]
DIM sn(360), cs(360)

GOSUB Setups                    'constants and data
GOSUB Indata                    'get model input data
GOSUB ScreenInfo                'set up screen for plot
GOSUB Plot                      'and plot each point
GOSUB Map                       'overlay map
DO
LOOP WHILE INKEY$ <> esc$       'exit
END

Setups:
  'constants
  cr$ = CHR$(13)                                'carriage return
  esc$ = CHR$(27)                               'escape
  pi = 3.141592653#
  tpi = pi * 2                                  '2 pi
  dr = pi / 180                                 'degrees to radians
  incr = 2 / 3
RETURN

ScreenInfo:
SCREEN 12                       '640 x 480 pixels
CLS
LOCATE 1, 2
PRINT "SHORT WAVE FADE MAP"
LOCATE 3, 12
PRINT USING "Day ##  Month ##  UT ##.##      X-ray Flux =##.#^^^^ W/m2"; day; month; UT; Xflux;
FOR i% = 0 TO 15                'plot colour bar
  LOCATE i% + 7, 77
  COLOR i%
  PRINT CHR$(219);
  PRINT USING "##"; i% * 2;     '& label with respective frequency
NEXT i%
COLOR 7
LOCATE 24, 77
PRINT "MHz";                    'print frequency units
RETURN


Indata:                                 'get model input data
 PRINT "SHORTWAVE FADE PLOT"
 INPUT "Month = (10) "; month           'month for plot
 IF month = 0 THEN month = 10
 INPUT "Day = (23) "; day               'day for plot
 IF day = 0 THEN day = 23
 YD = (month - 1) * 30 + day            'day of year
 INPUT "Universal Time = (6.3)"; UT     'time for plot (UT)
 INPUT "Solar X-ray flux (2E-4)"; Xflux 'solar X-ray flux
 IF Xflux = 0 THEN Xflux = .0002        'default flux is X2
RETURN


Plot:                                   'Main model calc and output
  soldec = 23.5 * SIN((YD - 82) * dr)   'solar declination
  solon = 15 * (12 - UT)                'solar longitude
  FOR lat = -90 TO 90 STEP .5
    FOR lon = 0 TO 360 STEP incr
      cchi = SIN(lat * dr) * SIN(soldec * dr) + COS(lat * dr) * COS(soldec * dr) * COS((solon - lon) * dr)
      IF cchi > 0 THEN          'cchi is cos of solar zenith angle
        fL = 174 * (Xflux + .0000018) ^ .25 * cchi ^ .75        'ALF
        clr = INT(fL / 2)
        IF clr > 15 THEN clr = 15
        xp = 50 + lon / incr
        yp = 240 - lat * 2
        PSET (xp, yp), clr
      END IF
    NEXT lon
  NEXT lat
RETURN

Map:
'overlay world map

LINE (50, 60)-(590, 420), 2, B          'define map edges
ryl = SQR(2)                             'y axis conversion constant
'draw lines of latitude every 30 degrees
FOR y = 120 TO 360 STEP 60
  LINE (50, y)-(590, y), 2, , &H1111
NEXT y
'draw lines of longitude every 30 degrees
FOR x = 95 TO 545 STEP 45
  LINE (x, 60)-(x, 420), 2, , &H1111
NEXT x
'label lattitude lines
LOCATE 5, 3
PRINT "+90";
LOCATE 8, 3
PRINT "+60";
LOCATE 12, 3
PRINT "+30";
LOCATE 15, 3
PRINT "Eqr";
LOCATE 19, 3
PRINT "-30";
LOCATE 23, 3
PRINT "-60";
LOCATE 26, 3
PRINT "-90";
'label longitude lines
LOCATE 27, 7
PRINT "0"
LOCATE 27, 23
PRINT "90E"
LOCATE 27, 39
PRINT "180"
LOCATE 27, 56
PRINT "90W"
LOCATE 27, 74
PRINT "0"
LOCATE 29, 35
PRINT "LONGITUDE";
LINE (50, 240)-(590, 240), 2            'equator

'now read and plot map
OPEN "ws.map" FOR INPUT AS #1
DO WHILE NOT EOF(1)
  INPUT #1, xp%, yp%
  PSET (xp%, yp%), 15
LOOP
CLOSE #1
RETURN