pro get_t7_ap, date, ft7, ft7a, fap
;
; Retrieve f10.7, f10.7 81-day smooth, and Fredericksburg Ap indices
; for input date(s) of the form yyyyddd >= 1947091 (1 Apr 1947).  The
; f10.7 indices returned (but not the 81-day smooth) are for the days
; one day prior to the input dates, so that all of the inputs for the
; MSIS-90 model atmosphere may be retrieved with a single call to
; this procedure.
;
; B. Knapp, 1993-04-14
;
; For any input date later than the most recent data available, an
; appropriate extrapolation is made, and a warning message is issued.
; The first time this procedure is called there will be a noticable
; delay while recent time series for f10.7 and Ap are retrieved,
; and the f10.7 81-day smooth is constructed.
;
; Change log
;    1995-08-16  B. Knapp: Changed first_day from 1980001 to 1947091.
;    1998-06-10  B. Knapp: IDL v. 5 compliance
;    2000-02-24  B. Knapp: Change smoothing window from 91 to 81 days (per
;                          e-mail of A. Hedin to S. Solomon, 1999-Nov-01)
;
  common t7ap_save, f, fs, a, last_day, t7_stop, ap_stop
;
; Set up smooth width
  swidth = 81L
  hwidth = swidth/2
;
  dsz = size( date )
  date_type = dsz[0]
  if date_type gt 1 then begin
     print, " DATE must be a scalar or 1-dimensional array"
     return
  endif
;
  case date_type of
     0: last_date = long( date )
     1: last_date = long( date[dsz[1]-1] )
  endcase
  t7_start = 1947045L
  first_day = yd2( t7_start, hwidth+1 )
  new_last_day = ( last_date > long( vms2yd( vmstime() ) ) ) > first_day
  if n_elements( last_day ) eq 0 then last_day = -1L
;
  if new_last_day ne last_day then begin
     last_day = new_last_day
     print,' Initializing F10.7, Ap arrays...'
;
;    First, get the entire F10.7 dataset and extend it into the future
;    so that we can do a 81-day smooth up to the present date.
     get_f10_7,t7_start,last_day,f,t,/fix_missing
     n = n_elements(f)
     t7_stop = yd2(1947045L,n-1)
     print,'    Latest F10.7 date: ',t7_stop
;
;    Determine the double-solar cycle period (about 22 years), so that
;    the present f10.7 data set may be extended by patching on data
;    from two 11-year cycles ago.
;    recent = f[n-5844:n-1]         ;16 years (more than one full cycle)
;    rmin = 1e38
;    for j=n-8766,n-7305 do begin   ;24 to 20 years ago
;    for j=n-8400,n-7670 do begin   ;23 to 21 years ago
;       r = total( (f[j-5843:j]-recent)^2 )   ;j aligns with n-1
;       if r lt rmin then begin
;          rmin=r
;          jmin=j
;       endif
;    endfor
;    SolarCycleYrs = (n-jmin-1)/365.25
     SolarCycleYrs = 22.0397d0
;    print,' RMS Residual =',sqrt(rmin/n), $
;          '; SolarCycleYrs =',SolarCycleYrs,' years'
;    print,' '
;
     jmin = long( SolarCycleYrs*365.25d0 )  ;double solar cycle
     x = ( dyd( t7_stop, last_day ) + hwidth ) > hwidth
     fx = [f,f[jmin+1:jmin+x]]
;
;    Now, the 81-day smooth
     fs = smooth(fx,swidth)
     j = dyd( t7_start, first_day )
     n = n_elements( fx )
     f = fx[j:n-hwidth-1]
     fs = fs[j:n-hwidth-1]
;
;    Now get the corresponding Ap dataset
     get_ap,first_day,last_day,a,t,/fix_missing
     n = n_elements( a )
     ap_stop = yd2( first_day, n-1 )
     print,'    Latest Ap date:    ',ap_stop
;
;    If last date is later than ap_stop, we need to extrapolate
;    (use 81-day mean from Ap stop date to last date).
     if last_day gt ap_stop then begin
        x = dyd( ap_stop, last_day )
        ap_mean = total(a[n-swidth:n-1])/float(swidth)
        ax = fltarr(x)+ap_mean
        a = [a,ax]
     endif
;
  endif
;
  if last_date gt t7_stop then print,string( t7_stop, $
     "(' *** Warning: F10.7 data later than',I8,' are extrapolated ***')")
;
  if last_date gt ap_stop then print,string( ap_stop, $
     "(' *** Warning: Ap data later than',I8,' are extrapolated ***')")
;
;
  t7a_pick = dyd( first_day, date )
  ap_pick = t7a_pick
  t7_pick = t7a_pick-1
;
  ft7a = fs[t7a_pick]
  ft7 = f[t7_pick]
  fap = a[ap_pick]
;
  return
  end