;# the only input parameter is `date', checks for multiple files per day
date = '19981005'
gpts = 64
dir  = '/home/sundowner/Data/diffrot/'
SPAWN, 'ls -1 ' + dir + 'par' + date + '*.sav', sav_files
SPAWN, 'ls -1 ' + dir + 'lct' + date + '*.bin', lct_files
SPAWN, 'ls -1 ' + dir + 'img' + date + '*.bin', img_files
n_files = SIZE ( sav_files, /N_ELEMENTS )



    ;# main event loop over multiple files per day
ii = 0
sum_lct = FLTARR ( 2, gpts, gpts )
SET_PLOT, 'z'
DEVICE, SET_RESOLUTION = [ 1280, 1024 ]
FOR j = 0, n_files - 1 DO BEGIN

    ;# restore variables created by `diffrot.pro'
    RESTORE, sav_files ( j )
    IF j EQ 0 THEN start_obs = time_obs0
    IF j EQ ( n_files - 1 ) THEN end_obs = time_obs1
    
    ;# do some automatic error checking
    n = SIZE ( dc_int, /N_ELEMENTS )
    flag = BYTARR ( n ) + 1
    
    ;# handle bad images interactively
;    IF j EQ 0 THEN flag ( 429 : 430 ) = 0        ; 19980803 
;    IF j EQ 0 THEN flag ( 14 : 15 ) = 0          ; 19980805 
;    IF j EQ 0 THEN flag ( 83 : 96 ) = 0
;    IF j EQ 1 THEN flag ( 15 : 16 ) = 0 
;    IF j EQ 1 THEN flag ( 27 : 28 ) = 0 
;    IF j EQ 1 THEN flag ( 31 : 32 ) = 0 
;    IF j EQ 1 THEN flag ( 40 : 41 ) = 0 
;    IF j EQ 1 THEN flag ( 56 : 58 ) = 0 
;    IF j EQ 1 THEN flag ( 71 : 75 ) = 0 
;    IF j EQ 1 THEN flag ( 79 : 81 ) = 0 
;    IF j EQ 1 THEN flag ( 86 : 87 ) = 0
;    IF j EQ 0 THEN flag ( 84 : 85 ) = 0          ; 19980813 
;    IF j EQ 0 THEN flag ( 103 : * ) = 0
;    IF j EQ 0 THEN flag ( 12 : 16 ) = 0          ; 19980814
;    IF j EQ 0 THEN flag ( 138 : 139 ) = 0        ; 19980815
;    IF j EQ 0 THEN flag ( 182 : 183 ) = 0
;    IF j EQ 0 THEN flag ( 234 : 235 ) = 0        ; 19980818
;    IF j EQ 0 THEN flag ( 246 : 247 ) = 0
;    IF j EQ 0 THEN flag ( 606 : 608 ) = 0        ; 19980819
;    IF j EQ 0 THEN flag ( 614 : 615 ) = 0
;    IF j EQ 0 THEN flag ( 622 : 623 ) = 0
;    IF j EQ 0 THEN flag ( 627 : 628 ) = 0
;    IF j EQ 0 THEN flag ( 213 : 214 ) = 0        ; 19980825
;    IF j EQ 0 THEN flag ( 231 : 232 ) = 0
;    IF j EQ 0 THEN flag ( 297 : 299 ) = 0
;    IF j EQ 0 THEN flag ( 301 : 303 ) = 0
;    IF j EQ 0 THEN flag ( 686 : 687 ) = 0 
;    IF j EQ 0 THEN flag ( 472 : 473 ) = 0        ; 19980826
;    IF j EQ 0 THEN flag ( 561 : 562 ) = 0
;    IF j EQ 0 THEN flag ( 101 : 102 ) = 0        ; 19980910
;    IF j EQ 0 THEN flag ( 71 : 72 ) = 0          ; 19980914
;    IF j EQ 0 THEN flag ( 180 : 181 ) = 0
;    IF j EQ 0 THEN flag ( 359 : 364 ) = 0
;    IF j EQ 0 THEN flag ( 376 : 384 ) = 0
;    IF j EQ 0 THEN flag ( 390 : 400 ) = 0
;    IF j EQ 0 THEN flag ( 404 : 410 ) = 0
;    IF j EQ 0 THEN flag ( 421 : 422 ) = 0
;    IF j EQ 0 THEN flag ( 19 : 20 ) = 0          ; 19980922
;    IF j EQ 0 THEN flag ( 36 : 37 ) = 0          ; 19980923
;    IF j EQ 0 THEN flag ( 11 : 19 ) = 0          ; 19980928
;    IF j EQ 0 THEN flag ( 12 : 13 ) = 0          ; 19981002
;    IF j EQ 0 THEN flag ( 206 : 207 ) = 0 
;    IF j EQ 0 THEN flag ( 34 : 35 ) = 0          ; 19981004      
;    IF j EQ 0 THEN flag ( 226 : 227 ) = 0        ; 19981005       
      
    index = WHERE ( exposure GE 50, count )
    IF count GT 0 THEN flag ( index ) = 0
    index = ( index + 1 ) < ( n - 1 )
    IF count GT 0 THEN flag ( index ) = 0
    
    index = WHERE ( delta_t GE 100, count )
    IF count GT 0 THEN flag ( index ) = 0
    index = ( index + 1 ) < ( n - 1 )
    IF count GT 0 THEN flag ( index ) = 0
    
    index = WHERE ( ABS ( disp_x ) GE 50, count )
    IF count GT 0 THEN flag ( index ) = 0
    index = ( index + 1 ) < ( n - 1 )
    IF count GT 0 THEN flag ( index ) = 0
    
    index = WHERE ( ABS ( disp_y ) GE 50, count )
    IF count GT 0 THEN flag ( index ) = 0
    index = ( index + 1 ) < ( n - 1 )
    IF count GT 0 THEN flag ( index ) = 0
    
    index = WHERE ( dc_int LE 100, count )
    IF count GT 0 THEN flag ( index ) = 0
    index = ( index + 1 ) < ( n - 1 )
    IF count GT 0 THEN flag ( index ) = 0
    
    ;# compute average flow map
    OPENR, 1, lct_files ( j )
    OPENR, 2, img_files ( j )
    lct = ASSOC ( 1, FLTARR ( 2, gpts, gpts ) )
    img = ASSOC ( 2, BYTARR ( 2 * gpts, 2 * gpts ) )
    FOR i = 1, n - 1 DO BEGIN
        IF flag ( i - 1 ) EQ 1 THEN BEGIN
    
    ;# average flow map
            IF ostype EQ 'linux' THEN sum_lct = sum_lct + lct ( i ) / $
                delta_t ( i - 1 ) ELSE sum_lct = sum_lct + $
                SWAP_ENDIAN ( lct ( i ) ) / delta_t ( i - 1 )
    
    ;# make JPEG images of the time series
            num = ii MOD 80
            xx = ( num MOD 10 ) * 128
            yy = 896 - ( num / 10 ) * 128            
            IF ostype EQ 'linux' THEN TVSCL, img ( i ), xx, yy $
                ELSE TVSCL, SWAP_ENDIAN ( img ( i ) ), xx, yy
            XYOUTS, xx + 10, yy + 10, STRING ( FORMAT = '(I4.4)', i ), $
                COLOR = 255, /DEVICE  
            IF ( num EQ 79 ) THEN BEGIN
                pic = TVRD ( )
                WRITE_JPEG, dir + 'pic' + date + STRING ( FORMAT = '( "_", I2.2)', $
                    ii / 80 ) + '.jpg', pic
                ERASE
            ENDIF
            ii = ii + 1
        ENDIF
    ENDFOR
    CLOSE, 1, 2
ENDFOR
sum_lct = sum_lct / ( ii + 1 )
IF ( num NE 79 ) THEN BEGIN
    pic = TVRD ( )
    WRITE_JPEG, dir + 'pic' + date + STRING ( FORMAT = '( "_", I2.2)', $
        ii / 80 ) + '.jpg', pic
    ERASE
ENDIF
SET_PLOT, 'x'



    ;# get B0
ccyy = STRMID ( date, 0, 4 )
mm = STRMID ( date, 4, 2 )
dd = STRMID ( date, 6, 2 )
start_obs = FIX ( STRMID ( start_obs, 0, 2 ) ) * 60. + $
    FIX ( STRMID ( start_obs, 3, 2 ) )
end_obs = FIX ( STRMID ( end_obs, 0, 2 ) ) * 60. + $
    FIX ( STRMID ( end_obs, 3, 2 ) )
et = ( end_obs - start_obs ) / 60.
IF et LT 0 THEN et = et + 24.
et = et / 2.
et = start_obs / 60. + et
IF et GE 24 THEN et = 23.99
SUN, ccyy, mm, dd, et, LAT0 = b0
b0 = b0 / !RADEG



    ;# remap flow map
cx = ( nx - 1. ) / 2.
cy = ( ny - 1. ) / 2.
nxg = LONG ( gpts )
nyg = LONG ( gpts )
ngw = FIX ( 2.0 * ( nx / nxg > ny / nyg ) )
nw = FIX ( 1.25 * ngw )
wx = FLOAT ( nx ) / nxg
wy = FLOAT ( ny ) / nyg
gx = FINDGEN ( nxg, nyg ) MOD nxg
gy = TRANSPOSE ( gx ) * wy + wy / 2. - 1.
gx = gx * wx + wx / 2. - 1.
eta = ( gx - cx ) / rsun
zeta = ( gy - cy ) / rsun
rho = SQRT ( eta^2 + zeta^2 )
mask = rho LT 0.95
REMAP_LCT, b0, gx, gy, nx, ny, rsun, mask, sum_lct, theta, phi, dtheta, dphi



    ;# display the results
WINDOW, XSIZE = 800, YSIZE = 800, XPOS = 380, YPOS = 224

    ;# full disk flow map
pic = BYTARR ( gpts, gpts )
x = FINDGEN ( gpts ) / gpts * 2032
y = FINDGEN ( gpts ) / gpts * 2028
vx = REFORM ( sum_lct ( 0, *, * ) ) * mask
vy = REFORM ( sum_lct ( 1, *, * ) ) * mask
SHOW_IMAGE_VECT, pic, vx, vy, x, y, LENGTH = 2
WAIT, 2

    ;# remaped velocities and legendre polynomial fit
index = WHERE ( mask EQ 1 )
theta =  COS ( theta ( index ) )
dphi = dphi ( index ) * 24. * 60. * 60. * 180. / !PI
index = SORT ( theta )
theta = theta ( index )
dphi = dphi ( index ) 
PLOT, theta, dphi, PSYM = 3, XRANGE = [ -1, 1 ], YRANGE = [ 0, 30 ]
coef = SVDFIT ( theta, dphi, 10, CHISQ = chisq, COVAR = covar, $
    SIGMA = sigma, VARIANCE = var, YFIT = yfit, /LEGENDRE )
OPLOT, theta, yfit

    ;# set the output device to PostScript
SET_PLOT, 'ps'

    ;# full disk flow map 
DEVICE, XSIZE = 18, YSIZE = 18, XOFF = 1.75, YOFF = 5, $
    FILENAME = dir + 'map' + date + '.ps'
MY_SHOW_IMAGE_VECT, pic, vx, vy, x, y, LENGTH = 2
XYOUTS, 0.0, 1.1, date, FONT = 0, CHARSIZE = 2, /NORM
DEVICE, /CLOSE

    ;# remaped velocities and legendre polynomial fit
DEVICE, XSIZE = 18, YSIZE = 18, XOFF = 1.75, YOFF = 6, $
    FILENAME = dir + 'fit' + date + '.ps'
PLOT, theta, dphi, FONT =0, PSYM = 3, THICK = 2, XRANGE = [ -1, 1 ], $
    XTHICK = 2,  YRANGE = [ 0, 30 ], YTHICK = 2
OPLOT, theta, yfit, THICK = 2
XYOUTS, 0.08, 1.05, date, FONT = 0, CHARSIZE = 2, /NORM
XYOUTS, 0.08, -0.02, STRING ( FORMAT = '( "Legendre Coefficients: ", 10(f8.2) )', $
    coef ), CHARSIZE = 0.8, FONT = 0, /NORM
XYOUTS, 0.08, -0.05, STRING ( FORMAT = '( "Standard Deviation: ", 10(f8.2) )', $
    sigma ), CHARSIZE = 0.8, FONT = 0, /NORM
XYOUTS, 0.08, -0.08, STRING ( FORMAT = '( "Variance:", 10(f8.2) )', $
    var ), CHARSIZE = 0.8, FONT = 0, /NORM
XYOUTS, 0.08, -0.11, 'Sum of Squared Errors: ' + STRING ( chisq ), $
    CHARSIZE = 0.8, FONT = 0, /NORM
XYOUTS, 0.08, -0.14, 'Number of Individual Flow Maps: ' + STRING ( ii ), $
    CHARSIZE = 0.8, FONT = 0, /NORM
XYOUTS, 0.08, -0.17, 'B0: ' + STRING ( FORMAT = '( f8.2, " deg at ",  a8 )', $
    b0 * !RADEG, STRSEC ( et * 3600. ) ), CHARSIZE = 0.8, FONT = 0, /NORM
DEVICE, /CLOSE    



    ;# reset the output device
SET_PLOT, 'x'

END