PRO rdrifts,date,times,vasun,vdusk,v_azi,v_hor,fillgap=fillgap

ds=STRING(date,FORMAT='(I6)')

; Read in drift velocities derived from individual patch
; analysis, and average over half-hour intervals
; The following code was largely chopped out of mkmtg.pro.
; ---------------------------------------------------------------------
path='/jasper/cnsr3_data1/pgm/'+ds+'02.mtg'
mtg=qgetfile(path)
n=N_ELEMENTS(mtg)
READS,mtg(n-1),lastn
npix2=FIX(lastn+1)
PRINT,npix2,' images analysed',FORMAT='(I0,A)'
info0=n-npix2
times=fltarr(npix2)
vdusk=times
vasun=times
p1=''
FOR i=0,npix2-1 DO BEGIN
    j=i+info0
    READS,mtg(j),p0,p1,p2,p3,p4,p5,p6,p7,p8 $
         ,FORMAT='(I3,2X,A12,2X,F8.2,2X,I3,2X,F8.4,2(2X,I3),2(2X,F8.5))'
    times(i)=p2
    vasun(i)=p7
    vdusk(i)=p8
    ENDFOR

IF KEYWORD_SET(fillgap) THEN BEGIN
    novels=BYTARR(npix2)
    gap=WHERE((vasun EQ 0.) and (vdusk EQ 0.),ngap)
    IF ngap GT 0 THEN BEGIN
        novels(gap)=1
        good=REPLICATE(1B,N_ELEMENTS(novels))-novels
        FOR i=0,ngap-1 DO BEGIN
	        tmp=extrac(vasun,gap(i)-3,7)
	        vasun(gap(i))=TOTAL(tmp)/TOTAL(tmp GT 0)
	        tmp=extrac(vdusk,gap(i)-3,7)
	        vdusk(gap(i))=TOTAL(tmp)/TOTAL(tmp GT 0)
            ENDFOR
        ENDIF
    ENDIF

vasun=vasun*111.		; convert from deg/s to km/s;
vdusk=vdusk*111.		; assume 111 km per deg of latitude

ymd=date
dst=ymd2date(1900+ymd/10000,(ymd/100) MOD 100,ymd MOD 100 $
	    ,format='d$ n$ Y$')
PRINT,dst
mint=MIN(times/3600.,MAX=maxt)
PRINT,'Patch velocities span ',mint,' to ',maxt,' hours UT' $
     ,FORMAT='(2(A,F5.2),A)'

IF NOT KEYWORD_SET(fillgap) $
    THEN vgood=WHERE((vasun NE 0.) OR (vdusk NE 0.)) $
    ELSE vgood=INDGEN(npix2)

v_hor=FLTARR(npix2)
v_azi=v_hor
v_hor(vgood)=SQRT(vasun(vgood)^2+vdusk(vgood)^2)
v_azi(vgood)=(ATAN(-vdusk(vgood),-vasun(vgood))*!RADEG+360.) MOD 360.

RETURN
END