'SOLRAS - simple solar radio spectrograph


GOSUB Setups
CLS
offon = 0
DO
mins% = INT(TIMER / 60)
IF offon = 0 AND mins% >= 450 AND mins% <= 1050 THEN
  offon = 1
  GOSUB NowSpec
  GOSUB InitDisplay
  GOSUB Openfile
END IF
IF offon = 1 AND (INT(TIMER) MOD 6 = 0) THEN GOSUB CollectData
IF offon = 1 AND (mins% < 450 OR mins% > 1050) THEN
  offon = 0
  CLOSE #1
  TIMER OFF
END IF
LOOP WHILE INKEY$ <> esc$
END


Setups:
 pi = 3.14159
 dr = pi / 180
 cr$ = CHR$(13)
 esc$ = CHR$(27)
 site$ = "Learmonth"         'site name
 siteabrev$ = "LM"
 LO = 114                     'longitude / degrees  East +   West -
 LT = -22                     'latitude / degrees North +   South -
 tz = 8                       'local time offset / hours East of UT
RETURN


CollectData:
 t0 = TIMER
 smin% = 255
 smax% = 0
 rec$ = ""
 GOSUB TimeLabel
 tr% = sec& / 6
 tr$ = MKI$(tr%)
 xt% = 37 + (sec& MOD 3600) / 6
 GOSUB LoBand
 GOSUB NowSpec
 GOSUB SolPos
 COLOR 7
 GOSUB GetData
 GOSUB HiBand
 LOCATE 2, 60
 PRINT USING "HA ###.#"; HAh;
 LOCATE 3, 60
 PRINT USING "DE ###.#"; SDd;
 LOCATE 4, 60
 PRINT USING "AZ ###.#"; AZ;
 LOCATE 5, 60
 PRINT USING "EL ###.#"; EL;
 GOSUB GetData
 GOSUB SaveRecord
 COLOR 15
 LOCATE 2, 77
 PRINT USING "#.#"; TIMER - t0;
 LOCATE 3, 77
 PRINT USING "###"; smin%;
 LOCATE 4, 77
 PRINT USING "###"; smax%;
 LOCATE 5, 75
 PRINT USING "#####"; tr%;
 LOCATE
RETURN


GetData:
 FOR i% = 0 TO 160 STEP 2
   OUT 888, i%
   indx% = i% + LoHi%
   rec% = STICK(0)
   IF rec% < smin% THEN smin% = rec%
   IF rec% > smax% THEN smax% = rec%
   rec$ = rec$ + CHR$(rec%)
   clr% = rec% / 16
   PSET (xt%, 423 - indx%), clr%
 NEXT i%
RETURN


Delay:
 FOR zz = 1 TO 3000
 NEXT zz
RETURN


LoBand:
OUT 888, 0
OUT 890, 1
LoHi% = 1
GOSUB Delay
RETURN


HiBand:
OUT 888, 0
OUT 890, 6
LoHi% = 175
GOSUB Delay
RETURN


Palet:
SCREEN 12
'set up palette
r = 1
G = 256
B = 65536
YEL = r + G
WHT = r + G + B
FOR i% = 0 TO 3
  PALETTE i%, (15 * i%) * B
NEXT i%
FOR i% = 0 TO 3
  PALETTE i% + 4, (7 * i% + 20) * YEL
NEXT i%
FOR i% = 0 TO 3
  PALETTE i% + 8, (8 * i% + 30) * r
NEXT i%
FOR i% = 0 TO 3
  PALETTE i% + 12, (8 * i% + 39) * WHT
NEXT i%
RETURN


InitDisplay:
GOSUB Palet
LINE (1, 1)-(639, 479), 8, BF
LINE (36, 86)-(637, 249), 15, B
LINE (37, 87)-(636, 248), 0, BF
LINE (36, 261)-(637, 423), 15, B
LINE (37, 262)-(636, 422), 0, BF
LOCATE 2, 2
COLOR 15
PRINT " SOLAR RADIOSPECTROGRAPH DISPLAY ";
FOR x% = 37 TO 637 STEP 20                      'minute markers
  LINE (x%, 250)-(x%, 253), 15
  LINE (x%, 424)-(x%, 426), 15
NEXT x%
FOR x% = 37 TO 637 STEP 100                     '5 minute markers
  LINE (x%, 253)-(x%, 260), 15
  LINE (x%, 427)-(x%, 429), 15
NEXT x%
LOCATE 29, 35
PRINT " TIME / minutes ";
LOCATE 6, 2
PRINT "220";
LOCATE 16, 2
PRINT "120";
LOCATE 11, 2
PRINT "MHz";
LOCATE 27, 3
PRINT "50";
LOCATE 22, 2
PRINT "MHz";
LOCATE 17, 2
PRINT "100";
LOCATE 4, 3
PRINT " Site: "; site$; " ";
'show palette
FOR i% = 0 TO 15
  FOR j% = 0 TO 9
    x% = 10 * i% + j% + 300
    LINE (x%, 15)-(x%, 35), i%
  NEXT j%
NEXT i%
'draw a little sun
 CIRCLE (570, 40), 20, 7
 PAINT (570, 40), 7
 LINE (540, 40)-(600, 40), 7
 LINE (570, 10)-(570, 70), 7
'add date
 LOCATE 4, 21
 PRINT USING " Date: ####/##/## "; yr; mn; dy
'add time scale
 LOCATE 28, 5
 PRINT "00          10          20           30           40          50          60";
RETURN


Openfile:
 OPEN filnam$ FOR APPEND AS #1
RETURN


TimeLabel:
 LOCATE 4, 40
 tm$ = TIME$
 sec& = TIMER
 min = VAL(MID$(tm$, 4, 2))
 PRINT USING " Time: & LT "; tm$;
RETURN


DisplayRecord:
  FOR i% = 1 TO 256
    xt% = (INT(TIMER) MOD 600)
    clr% = rec(i%) / 16
    PSET (xt% + 36, 423 - i%), clr%
  NEXT i%
RETURN


SaveRecord:
  PRINT #1, tr$ + rec$;
RETURN


NowSpec:
  'Get date and process to year, month and day
  dt$ = DATE$
  yr$ = RIGHT$(dt$, 4)                       'year
  yr2$ = RIGHT$(yr$, 2)
  yr = VAL(yr$)
  IF (yr - 4 * INT(yr / 4)) = 0 THEN MD(2) = 29  'adjust month days for leap year
  mn$ = LEFT$(dt$, 2)                        'month
  mn = VAL(mn$)
  dy$ = MID$(dt$, 4, 2)                      'day
  dy = VAL(dy$)
  filnam$ = siteabrev$ + yr2$ + mn$ + dy$ + ".SSS"
  LT = TIMER / 3600!             'Local time in decimal hours
  hr = INT(LT)
  min = INT((LT - hr) * 60)
  LT$ = TIME$                    'Local time in usual string format
  UT = LT - tz
RETURN


SolPos:
  'this code determines various parameters of the sun
  'and uses formulas given at the end of section C of the
  'Astronomical Almanac
  day = dy + UT / 24            'UT decimal day
  N = INT(yr + FIX((mn - 9) / 7))               'N is days to J2000.0
  N = -INT(3 * (INT(N / 100) + 1) / 4)
  N = N - INT(7 * (INT((mn + 9) / 12) + yr) / 4)
  N = N - 730516.5 + day + 367 * yr + INT(275 * mn / 9)
  L = 280.46 + .9856474 * N             'mean solar longitude
  L = L - INT(L / 360) * 360            'put in range 0 to 360 deg
  G = (357.528 + .9856003 * N) * dr     'mean anomaly
  G = G - INT(G / 2 / pi) * 2 * pi      'put in range 0 to 2*pi
  LA = L + 1.915 * SIN(G) + .02 * SIN(2 * G)    'ecliptic longitude
  EP = (23.439 - .0000004 * N) * dr     'obliquity of ecliptic
  SD = SIN(EP) * SIN(LA * dr)           ' sin(declination)
  SD = ATN(SD / SQR(1 - SD * SD))       ' ASIN function
  AL = ATN(COS(EP) * TAN(LA * dr)) / dr         'right ascension
  IF AL < 0 THEN AL = AL + 180          'put in correct quadrant
  IF LA > 180 THEN AL = AL + 180
  IF LA > 360 THEN AL = AL + 180
  AU = 1.00014 - .01671 * COS(G) - .00014 * COS(2 * G)   'solar distance
  ET = (L - AL) * 4                     'equation of time
  HA = 15 * (UT - 12 + ET / 60 + LO / 15)       'hour angle (deg)
  IF HA > 180 THEN HA = HA - 360        'put in range -180 to 180 deg
  SE = SIN(LT * dr) * SIN(SD) + COS(LT * dr) * COS(SD) * COS(HA * dr)
  IF ABS(SE) = 1 THEN                        'singularity
    EL = 90 * SGN(SE): AZ = 0                'solar elevation and azimuth
  ELSE
    EL = ATN(SE / SQR(1 - SE * SE)) / dr     'solar elevation
    CA = (SIN(SD) - SIN(LT * dr) * SE) / COS(LT * dr) / SQR(1 - SE * SE)
    IF CA = 0 THEN                           'CA is cos(azimuth)
      AZ = 90
    ELSE
      AZ = ATN(SQR(ABS(1 - CA * CA)) / CA) / dr         'ACOS function
      IF AZ < 0 THEN AZ = 180 + AZ              'put in right quadrant
    END IF
    IF HA > 0 THEN AZ = 360 - AZ
  END IF
  SDd = SD / dr                 'solar declination / degrees
  HAh = HA / 360 * 24           'solar hour angle / hours
RETURN