DIM f(2000), s(2000), c(2000), t(12000), V(12000)

dri$ = "\ER2\"
CLS : PRINT "Welcome to ROGUE.BAS"
PRINT "Program creates Fourier component amplitudes & constructs a temp(time) trace"
PRINT "  A total amplitude squared (i.e., power) [K^2] is generated, "
PRINT "  noise is added."
PRINT

PRINT "Will create spectral components & store in f(L) & s(L)"
 RANDOMIZE

 df = .2775  '[cy/day]
 fmax = 500: f = 0

 L = 1: f = df: f(1) = f: s(1) = .7 * 2.24: c(1) = .7 * 2.24
 s(1) = .7 * 2.24: c(1) = -.7 * 2.24

 DO
   L = L + 1:
   f = f + df: f(L) = f
   Mn = 0: Mp = 0
   Mp = .3 + .1 * Log10(f)
   Mn = 1.2 + .5 * Log10(f)
   Asqr0 = 1! * f(L) ^ (-5 / 3)
 
   CALL GAUSS(X)
   IF X < 0 THEN m = Mn ELSE m = Mp
   ratio = 10 ^ (m * X): Asqr = Asqr0 * ratio
   s(L) = SQR(Asqr)
 
   CALL GAUSS(X)
   IF X < 0 THEN m = Mn ELSE m = Mp
   ratio = 10 ^ (m * X): Asqr = Asqr0 * ratio
   c(L) = SQR(Asqr)
 
 LOOP WHILE f(L) < fmax
 Lmax = L - 1

 OPEN "O", 1, "\er2\AB_SPEC.PRN"
   PRINT #1, "Synthetic Spectrum"
   FOR L = 1 TO Lmax
   PRINT #1, L;
   PRINT #1, USING "####.##"; f(L);
   PRINT #1, USING "###.#####"; s(L);
   PRINT #1, USING "###.#####"; c(L);
   PRINT #1, USING "###.#########"; s(L) ^ 2 + c(L) ^ 2
   NEXT L
 CLOSE #1

 Lmax = L: PRINT : PRINT "There are"; Lmax; "terms."
 PRINT "Frequencies are:"; f(1); ","; f(2); " to"; f(Lmax)
 INPUT "Enter wind speed (for calculating wavelengths; def= 24) "; A$
  IF A$ = "" THEN A$ = "24"
  ws = VAL(A$): PRINT "  Wind speed = "; ws; "[m/s]": PRINT

 Wmax = ws * 60 * 1440 / (1000 * f(1)):
 Wmin = ws * 60 * 1440 / (1000 * f(Lmax))
 PRINT "Wavelengths range from"; Wmin; " to"; Wmax; "[km]."

 PRINT : INPUT "Enter beginning time for reconstituted trace"; A$
  IF A$ = "" THEN A$ = "0"
    tbegin = VAL(A$)
  INPUT "Enter ending time for reconstituted trace (def=2798 min)"; A$:
  IF A$ = "" THEN A$ = "2798.0"
    tend = VAL(A$)
  dt = (ws * 60 / 1000) * Wmin / 4
  dt = .964495 * dt     'empirical fudge
  PRINT "  Will reconstitute for time range"; tbegin; "to"; tend; ",";
  PRINT "for dt ="; dt; "minutes.";
  PRINT "  There will be a total of"; INT(tend / dt); "time locations."

 PRINT : INPUT "Enter high-pass cut-off wavelength [km] (def=10)"; A$
   IF A$ = "" THEN A$ = "10"
   HPW = VAL(A$): PRINT "  Will use wavelength high-pass cut-off of"; HPW; "[km]"

 Cfw = ws * 60 * 1440 / 1000 'constant for converting freq to wavelength

 PRINT : PRINT "Wait for reconstitution...": PRINT : L = 0
 t = tbegin - dt: Ndun = 0

NewTime:
 t = t + dt: IF t > tend THEN GOTO DoneWithTimes
 L = 0: V = 0
 phi0 = t * 2 * 3.1415926# / 1440

AddTerm:
 L = L + 1: f = f(L): s = s(L): c = c(L)
 w = Cfw / f:
 IF w < HPW THEN GOTO DoneWithTerms
 phi = phi0 * f
 V = V + s * SIN(phi) + c * COS(phi)
 GOTO AddTerm

DoneWithTerms:
 Ndun = Ndun + 1
 t(Ndun) = t: V(Ndun) = V
 IF Ndun MOD 10 = 0 THEN PRINT USING "########"; Ndun; : PRINT USING "####.###"; V;
 GOTO NewTime

DoneWithTimes:
 Nmax = Ndun
 PRINT : PRINT "Reconstituted trace is ready."
OPEN "O", 2, dri$ + "\" + "AA_RECON.PRN"

 PRINT "Recording reconstituted trace..."
 PRINT #2, "Reconstituted temp(t)"
 FOR j = 1 TO Nmax
  PRINT #2, USING "#####.###"; t(j);
  PRINT #2, USING "#####.###"; V(j)
  IF j MOD 10 = 0 THEN PRINT ".";
 NEXT j
 CLOSE
 BEEP: PRINT : PRINT "DONE!"

SUB GAUSS (X)
  Y = RND
  X = -.46 * LOG(1 / Y - 1)

END SUB

FUNCTION Log10 (X) STATIC
     Log10 = LOG(X) / LOG(10#)
END FUNCTION