; ap_inter_sol: interpolated solar wind hourly averages 
; this auroral power version averages over 4 hours, with declining weights 
; needs time-lagged hourly average IMF data in ascii files (OMNI time-lagged) 
;Patrick Newell, November 2006

pro ap_inter_sol, yearcall, day, sod, mag, bx, by, bz, ni, v,Ec 
common dst_inter_sol,last_year,nrec,maga,bxa,bya,bza,nia,va
nh = 4       ;hours previous to integrate over
wh = 0.65     ;reduce weighting by factor of wh each hour back

hour = fix(sod/3600)
res = (float(sod) - 3600.*hour)/3600.
if( n_elements(last_year) eq 0 )then begin
  last_year = 0                        ;first call to routine
  nrec = 24*366
  maga = fltarr(nrec)
  bxa = fltarr(nrec)
  bya = fltarr(nrec)
  bza = fltarr(nrec)
  nia = fltarr(nrec)
  va = fltarr(nrec)
endif

; input year should be only 2 digits. year, day, hour are integer mag = 0.
bx = 0.
by = 0.
bz = 0.
ni = 0.
v = 0.

r = 0
i = 0
yrs=['1983','1984','1985','1986','1987','1988','1989','1990','1991','1992',$
 '1993','1994','1995','1996','1997','1998','1999','2000','2001','2002','2003',$
'2004','2005']

year = yearcall
if( year gt 200) then year=year-1900
if( (year lt 83) or (year gt 105) ) then begin return endif

if( year ne last_year )then begin
close,23
nrec = 24*365
if ( 4*fix(year/4) eq year )then nrec  = 24*366 
sol_file = '/project/dmsp/solwind/solwind' + yrs(year-83) + '.dat'
openr,23,sol_file


f101 = '(2i4,i3,5f6.1,f6.0)'
yr = 0
dy = 0
hr = 0
for r = 0, nrec-1 do begin
readf,23,format=f101,yr,dy,hr,mag,bx,by,bz,ni,v
testr = hr + 24*(dy-1)
if( r ne testr )then begin
   print,'Discrepancy between expected time and solwind file time:'
   print,'hr,dy,r,testr=',hr,dy,r,testr
endif
maga(r) = mag
bxa(r) = bx
bya(r) = by
bza(r) = bz
nia(r) = ni
va(r) = v
endfor
last_year = year
endif

r1 = hour + 24*(day-1)
;r1 = r + 1
if( r1 ge 8784 ) then begin
  r1 = 8783
endif
if( (4*(year/4)) ne year )then begin
  if( r1 ge 8760 )then begin
    r1 = 8759
  endif
endif

w_arr = fltarr(nh)
mag_arr = fltarr(nh)
bx_arr = fltarr(nh)
by_arr = fltarr(nh)
bz_arr = fltarr(nh)
ni_arr = fltarr(nh)
v_arr = fltarr(nh)
Ec_arr = fltarr(nh)
ng = 0
for i = 0, nh-1 do begin
  ri = r1 - i
  if( ri lt 0 )then begin
     ri = 0
  endif
  mag_arr(i) = maga(ri)
  bx_arr(i) = bxa(ri)
  by_arr(i) = bya(ri)
  bz_arr(i) = bza(ri)
  ni_arr(i) = nia(ri)
  v_arr(i) = va(ri)
  sol_coup,bx_arr(i),by_arr(i),bz_arr(i),v_arr(i),ni_arr(i),Ec
  Ec_arr(i) = Ec
  w_arr(i) = wh^i
  ng = ng + 1
  if( (mag_arr(i) eq 0.) or (v_arr(i) eq 0.) or (ni_arr(i) eq 0.) )then begin
    w_arr(i) = 0.
    ng = ng - 1
  endif
endfor
a = (float(sod) - 3600.*hour)/3600.
w_arr(0) = sqrt(a)*w_arr(0)

if( ng lt 3 )then begin
   mag = 0.
   bx = 0.
   by = 0.
   bz = 0.
   v = 0.
   ni = 0.
   return
endif

wt = total(w_arr)
w_arr = w_arr/wt
mag = total(w_arr*mag_arr)
bx = total(w_arr*bx_arr)
by = total(w_arr*by_arr)
bz = total(w_arr*bz_arr)
v = total(w_arr*v_arr)
ni = total(w_arr*ni_arr)

Ec = total(w_arr*Ec_arr)

return
end