'IPSHOK/STOA-INTERPLANETARY SHOCKWAVE PREDICTION PROGRAM
'
  DIM M$(12), M(12)     'space for month names and # days in month
'
' Input and initialisation
'
 CLS
 PRINT "Interplanetary Shockwave Prediction"
 PRINT "Solar Event - FLARE"
 INPUT "Longitude (+/- deg)"; LO
 INPUT "Latitude  (+/- deg)"; LA
 INPUT "Flare Type / Description  (optional)"; EVNT$
 INPUT "(Driver) Duration / mins"; TS: T2 = TS * 60
 PRINT "Solar Event - Type II sweep"
 INPUT "Shock Velocity (km/s)"; VK: VS = VK * 1000
 INPUT "  Year (eg 1987)"; TY
 INPUT "  Month (eg 8)"; TM
 INPUT "  Day"; TD
 INPUT "  Hour"; TH
 INPUT "  Minute"; TN
 INPUT "Solar Wind Speed (default is 400 km/s)"; VWK
'
 RS = 6.96E+08                  ' solar radius
 AU = 214.94 * RS               ' 1 astronomical unit
 SS = 45000!                    ' IMF default sound speed at earth
 IF VWK = 0 THEN VWK = 400      ' default solar wind speed at earth
 VW = VWK * 1000                ' put in metres/sec
'
' Compute shock shape factor
'
 DR = 3.1416 / 180: RO = LO * DR: RA = LA * DR   ' to radians
 CS = COS(RO) * COS(RA)
 SM = (CS + 1) / 2
'
' Track blast wave out from sun to earth
'
'    - start at 1.5 solar radii
 R = 1.5 * RS
'    - move out shock vel * driver time
 V0 = VS * SM
 R = R + V0 * T2
 R0 = R: RB = R
 T = T2
 '    - now integrate to earth orbit
 DT = 600                        ' 10 minutes
 DO
   SV = V0 * SQR(R0 / RB)          ' current shock velocity
   DR = (SV + VW) * DT             ' distance increment in time dt
   R = R + DR: T = T + DT          ' new time, new position
   RB = RB + SV * DT               ' blast effective travel
 LOOP WHILE R < AU                  ' earth orbit yet?
'    - we have now reached earth
'
' Test for shock and print
'
 MN = SV / SS                    ' mach number ( >1 for shock )
'
  DATA "JAN","FEB","MAR","APR","MAY","JUN"
  DATA "JUL","AUG","SEP","OCT","NOV","DEC"
  FOR I = 1 TO 12: READ M$(I): NEXT I
  DATA 31,28,31,30,31,30,31,31,30,31,30,31
  FOR I = 1 TO 12: READ M(I): NEXT I
  IF (TY - 4 * INT(TY / 4)) = 0 THEN M(2) = 29
  T = T / 60                      ' transit time in minutes
  DH = INT(T / 60)                ' adjust for day, month, year
  DM = T - 60 * DH                ' as appropriate
  DD = INT(DH / 24)
  DH = DH - 24 * DD
  AN = TN + DM
  IF AN > 59 THEN AN = AN - 60: DH = DH + 1
  AH = TH + DH
  IF AH > 23 THEN AH = AH - 24: DD = DD + 1
  AD = TD + DD: AM = TM: AY = TY
  IF AD > M(TM) THEN AD = AD - M(TM): AM = AM + 1
  IF AM > 12 THEN AM = 1: AY = TY + 1
'
' Print
'
 CLS
 PRINT "INTERPLANETARY SHOCKWAVE PREDICTION PROGRAM"
 PRINT "*******************************************"
 PRINT : PRINT "SOLAR EVENT  "; EVNT$
 PRINT USING "     ## "; TD;
 PRINT M$(TM); TY
 TI = 100 * TH + TN
 PRINT USING "     #### UT"; TI
 PRINT "     HELIOGRAPHIC";
 PRINT USING " ##"; ABS(LA);
 IF LA < 0 THEN PRINT "S";  ELSE PRINT "N";
 PRINT USING " ##"; ABS(LO);
 IF LO < 0 THEN PRINT "E" ELSE PRINT "W"
 PRINT USING "     DURATION ### MINS"; TS
 PRINT : PRINT "TYPE II SHOCK VELOCITY"
 PRINT USING "     #### KM/S"; VK
 PRINT : PRINT "SOLAR WIND SPEED"
 PRINT USING "     #### KM/S"; VWK
 PRINT
'
' Check to see if disturbance is really a shock
'
 IF MN >= 1 THEN
   PRINT "SHOCK ARRIVAL AT EARTH"
   PRINT USING "     EXPECTED  ## "; AD;
   PRINT M$(AM); " "; AY
   PRINT USING "               ##:## UT"; AH; AN
   PRINT
   PRINT USING "     TRAVEL TIME  ####.# HRS"; T / 60
   PRINT USING "     MACH NUMBER  ##.#"; MN
 ELSE
   PRINT "NO SHOCKWAVE EXPECTED"
   PRINT "AT EARTH ORBIT. MINOR"
   PRINT "GEOMAGNETIC FLUCTUATION"
   PRINT "MAY OCCUR AROUND"
   PRINT USING "   ## "; AD;
   PRINT M$(AM); " "; AY
   PRINT USING "   ##:00 UT"; AH
 END IF
END