; Assumes cen = mrdfits('obs_acen1.fits',1) has been done
; and that slots = [0,1,2] has been done
;   1   time                 s            Real8             74785401.31880015:   74785427.96880114 label for field   1
;   2   slot                              Byte                                label for field   2
;   3   alg                               Byte                                label for field   3
;   4   cent_i               pixel        Real4          -Inf:+Inf            label for field   4
;   5   cent_j               pixel        Real4          -Inf:+Inf            label for field   5
;   6   ang_y                deg          Real4          -Inf:+Inf            label for field   6
;   7   ang_z                deg          Real4          -Inf:+Inf            label for field   7
;   8   ang_y_sm             deg          Real4          -Inf:+Inf            label for field   8
;   9   ang_z_sm             deg          Real4          -Inf:+Inf            label for field   9
;  10   counts               count        Real4          -Inf:+Inf            label for field  10
;  11   bkg                  count        Real4          -Inf:+Inf            label for field  11
;  12   chisq                             Real4          -Inf:+Inf            label for field  12
;  13   status                            Byte                                label for field  13

M0 = 10.32                      ; ACA magnitude corresponding to C0 electrons/second
C0 = 5263.0                     ; Electrons/second
t_int = 1.696                   ; Assume fixed readout time of 1.696 seconds

obsid = sxpar(hdr, 'OBS_ID')

for islt = 0, n_elements(slots)-1 do begin
    ok = where(cen.slot eq slots[islt] and cen.alg eq 8 and cen.status eq 0, nok)
    if (nok gt 3) then begin
        cok = cen[ok]
        ind_y = sort(cok.ang_y_sm)
        ind_z = sort(cok.ang_z_sm)
        cts_med = median(cok.counts)
        mag_med = M0 - 2.5 * alog10(cts_med/t_int / C0)

        ; If possible, skip the first 4 ksec of data to avoid thermal transients
        if (nok gt 3000) then ok0 = 2000 else ok0 = 0

        ang_y = cok[ind_y].ang_y_sm*3600.0
        ang_z = cok[ind_z].ang_z_sm*3600.0

        print, 'slot         = ', slots[islt]
        print, 'tstart       = ', cok[0].time
        print, 'obsid        = ', obsid
        print, 'date_obs     = ', sxpar(hdr, 'DATE-OBS')
        print, 'ang_y_med    = ', median(cok[ok0:*].ang_y_sm)*3600
        print, 'ang_z_med    = ', median(cok[ok0:*].ang_z_sm)*3600
        print, 'ang_y_5th    = ', ang_y[long(nok*0.05)]
        print, 'ang_y_95th   = ', ang_y[long(nok*0.95)]
        print, 'ang_z_5th    = ', ang_z[long(nok*0.05)]
        print, 'ang_z_95th   = ', ang_z[long(nok*0.95)]
        print, 'mag_med      = ', mag_med
        print, 'exp_time     = ', (sxpar(hdr, 'TSTOP') - sxpar(hdr, 'TSTART'))/1000.
        print, 'END_RECORD'

    end
end

end