'PROTON PREDICTION MODEL ' - based on PPS76 by Smart and Shea ' - with additions DIM EL(7), M$(12), ndm(12), PGE(7) GOSUB Setups 'initialisation GOSUB ManualInput 'get input data GOSUB PreCompute 'compute initial values for model GOSUB Headings 'printout headings GOSUB Prediction 'compute & display prediction values GOSUB TimeHistory 'compute & display event time history PRINT "PPM output sent to ppm.dat" PRINT "Program finished" END Setups: 'Initialisation section PI = 3.14159: DR = PI / 180 'Degrees to Radians 'set up energy levels in MeV in array EL EL(0) = 1000: EL(1) = 100: EL(2) = 50: EL(3) = 20 EL(4) = 10: EL(5) = 5: EL(6) = 2: EL(7) = 1 'load M$ array with days of month DATA "JAN","FEB","MAR","APR","MAY","JUN" DATA "JUL","AUG","SEP","OCT","NOV","DEC" FOR I = 1 TO 12: READ M$(I): NEXT I 'load number of days in each month array ndm DATA 31,28,31,30,31,30,31,31,30,31,30,31 FOR I = 1 TO 12: READ ndm(I): NEXT I 'open file for output OPEN "ppm.dat" FOR OUTPUT AS #1 RETURN ManualInput: 'Input required data manually from keyboard PRINT "Proton Prediction Program" PRINT " -default values are shown in parentheses" PRINT "FLARE LOCATION" INPUT " Latitude-degrees N (0) "; FLA INPUT " Longitude-degrees W (0) "; FLO INPUT "Solar Wind Speed (400 km/s) "; VSW: IF VSW = 0 THEN VSW = 400 INPUT "Semi-integrated flux on 8800 MHz in SFU-sec (1.0E5) "; SFLUX IF SFLUX = 0 THEN SFLUX = 100000! '10^5 is minimum for good proton event INPUT "Castelli-U Omega ratio W3/W2 (5) "; OMEGA: IF OMEGA = 0 THEN OMEGA = 5 PRINT "CASTELLI-U PEAK BURST TIME" INPUT " Year (1990) "; DY: IF DY = 0 THEN DY = 1990 IF ((DY - 4 * INT(DY / 4)) = 0) AND (DY <> 2000) THEN ndm(2) = 29 INPUT " Month (1) "; DM: IF DM = 0 THEN DM = 1 INPUT " Day (1) "; DD: IF DD = 0 THEN DD = 1 INPUT " Hour (0) "; DH INPUT " Minute (0)"; DHM: DH = DH + DHM / 60 RETURN PreCompute: 'Compute initial values for model PHIA = 400 / VSW 'solar CMD of Archimedian spiral to Earth ' compute great circle distance from flare to foot of spiral [THETA] CTHETA = COS(PHIA - FLO * DR) * COS(FLA * DR) THETA = ATN(SQR(1 - CTHETA * CTHETA) / CTHETA) 'inverse cos function ' compute distance along spiral to Earth (in A.U.) [D] D = SQR(PHIA * PHIA + 1) D = .5 * (D + (LOG(PHIA + D)) / PHIA) ' compute proton attenuation factor - due to coronal diffusion [ATF1] ATF1 = EXP(-3 * THETA) ' compute initial proton flux spectral index (Initial Gamma Zero) IG0 = -2.73 / LOG(OMEGA) 'integral G0 = IG0 - 1 'differential ' compute proton flux over 5 MeV - J(>5 MeV) J5 = .0017 * (SFLUX) ^ .812 ' compute coefficient K for diferential proton spectra FTOT = 0 'do crude integration FOR E = 5 TO 1000 'integrate to 1000 MeV GOSUB EnergiesAndExponents 'get energy spectral exponent FTOT = FTOT + E ^ GB 'add all contributions above 4 MeV NEXT E K = J5 / FTOT 'coefficient of differential spectrum RETURN Headings: 'Print out headings PRINT #1, : PRINT #1, : PRINT #1, PRINT #1, , "PROTON PREDICTION MODEL" PRINT #1, , "***********************" PRINT #1, : PRINT #1, : PRINT #1, IF FLA < 0 THEN NS$ = "S" ELSE NS$ = "N" IF FLO < 0 THEN EW$ = "E" ELSE EW$ = "W" PRINT #1, "Solar event at "; ABS(FLA); NS$; " "; ABS(FLO); EW$ PRINT #1, "Event reference time is "; DY; M$(DM); DD; " at "; DH; "Hours UT" PRINT #1, PRINT #1, "Solar wind speed = "; VSW; " km/sec" PRINT #1, PRINT #1, "8800 MHz Semi-Integrated Flux (Start to Peak) = "; SFLUX PRINT #1, "Castelli-U Omega Ratio (W3/W2) = "; OMEGA PRINT #1, : PRINT #1, : PRINT #1, PRINT #1, , "*** Fast Prediction Parameters ***" PRINT #1, , " (At Geosynchronous Orbit)" PRINT #1, PRINT #1, "Proton Energy", "Onset", "Maximum", "Peak flux >E", "Decay Constant" PRINT #1, " E (MeV)", "Time(hrs)", "Time(hrs)", " (PFU)", " (hours)" PRINT #1, RETURN Prediction: 'Compute and Display Prediction Parameters F$ = " ### ###.# ###.# ###### ###.#" PF = 0 FOR I = 1 TO 7 'step through energy levels E = EL(I) 'get current energy GOSUB TimesAndFluxes 'compute times and fluxes IF I = 1 THEN TMIN = INT(TS) PF = J5 * (E / 5) ^ IGB * ATF1 * ATF2 PRINT #1, USING F$; E; TS; TM; PF; TD NEXT I RETURN TimeHistory: 'Compute & display event time history 'print headings FOR I = 1 TO 3: PRINT #1, : NEXT I PRINT #1, , "*** Proton Event History ***": PRINT #1, PRINT #1, "Time after ----------------Proton Flux (PFU)---------------- 30MHz Riometer" PRINT #1, "event(hrs) >100MeV >50MeV >20MeV >10Mev >5MeV >2Mev >1MeV A(dB)/d A(dB)/n" PRINT #1, H$ = " ### ###### ###### ###### ###### ###### ###### ###### ##.# ##.#" FOR T = TMIN TO 50 PGE(0) = 0 FOR I = 1 TO 7 PGE(I) = PGE(I - 1) FOR E = EL(I - 1) - 1 TO EL(I) STEP -1 GOSUB TimesAndFluxes IF T >= TS THEN IF T < TM THEN PGE(I) = PGE(I) + K * E ^ GB * ATF1 * ATF2 * EXP((T - TM) / (T - TS)) ELSE PGE(I) = PGE(I) + K * E ^ GB * ATF1 * EXP(-T / TD) END IF END IF NEXT E NEXT I IF PGE(7) >= 1 THEN AD = .6 * SQR(PGE(5)): AN = .25 * SQR(PGE(6)) PRINT #1, USING H$; T; PGE(1); PGE(2); PGE(3); PGE(4); PGE(5); PGE(6); PGE(7); AD; AN END IF NEXT T RETURN TimesAndFluxes: 'Compute times, fluxes, etc. GOSUB EnergiesAndExponents 'get particle speed & spectral exponents 'compute onset (start) time - hours TS = .25 + .133 * D / B + 4 * THETA * THETA 'get max flux time - hours TM = 2 * (.133 * D * D / B) + 8 * THETA * THETA 'compute decay time TD = 3 * D / (4 * VSW * (ABS(GB) + 1)) * 1.5E+08 / 3600 'compute peak attenuation factor ATF2 = EXP(-TM / TD) RETURN EnergiesAndExponents: 'Compute particle energies and exponents ' compute particle speed B=v/c B = SQR(1 - 1 / (1 + E / 938) ^ 2) '938 MeV is proton rest energy ' spectral exponents GB = G0 * (1 + THETA / 2) + .15 - B 'differential IGB = GB + 1 'integral RETURN