function msu_xrt_med_dark,index,ddates=ddates,nday=nday,nim=nim,darkff=darkff, $ 
  dfile=dfile,stripdark=stripdark,quiet=quiet
;+
;NAME:
;   xrt_med_dark
;PURPOSE:
;   Search full-frame dark images, make one full-frame medianed dark template, 
;   and then clip an appropriate region for a given FOV.
;
;   If such a medianed dark already exists, program checks to see if better
;   (closer in time to target data image) dark images exist and uses them
;   to make a new medianed dark template.  Otherwise it reuses the previous
;   medianed dark template.
;
;   In the special case of 1x1 images taken after the high-gain antenna
;   failure,  vertical strip (512x2048) darks are used when available
;   and appropriate
;
;   Program exits, setting dark=0., if it cannot derive a median dark
;   given the restrictions of nim and nday
;
;CALLING:
; dark =xrt_med_dark(index [,ddates=ddates,darkff=darkff,nday=nday, $
;                   dfile=dfile,stripdark=stripdark,quiet=quiet])
;
;INPUT:
;   index   - index of (a single) image.
;
;OPTIONAL INPUTS:
;   ddates=ddates - string array of dates/times of the dark files
;              in an anytim() readable format corresponding to previous 
;              medianed dark below 
;   darkff=darkff - the previous radiation-free full-frame dark template 
;               generated by dark files with dates = ddates.
;   nday=nday - define a half-period to search dark images.
;               (default day=14; <7 is not advised)
;   nim=nim - define number of dark images to median
;               (default nim=5; odd numbers are best)
;  quiet=quiet - run without printing [default - not quiet]
;
;OUTPUT:
;   dark    - a medianed dark image for the CCD FOV defined by index.
;
;OPTIONAL OUTPUTS:
;   ddates    - string array of dates/times of the dark files
;   darkff    - the radiation-free full-frame dark template 
;                   generated for the FOV defined by index using dark 
;                   filenames = dfile and dates = ddates
;   dfiles    - filename of dark images that were be used.
;   stripdark - flag =1 if 1x1 binned "strip darks"
;                   (512x2048) were used.
;
; EXAMPLES:
;     generate a single medianed dark (typically for use in dark subtraction
;       or baseline calibration) corresponding to data described by index:
;
;      dark =xrt_med_dark(index)
;
;     generate a single medianed dark corresponding to data described by 
;     together with index, median over 7 darks but decrease the search range
;     to 10 days, save the full-frame median dark, and dark files used:
;
;      dark =xrt_med_dark(index,nday=10,nim=7,darkff=darkff,dfile=dfile )
;
;    when correcting multiple images close in time, it is advisable 
;    to output the full-frame median dark, and the file dates.  Then
;    on the subsequent calls of the program,  it will check to see if
;    it can reuse the last medianed dark (extracting it as needed from the 
;    full-frame median), which can save significant processing time
;
;      dark =xrt_med_dark(index,ddates=ddates,darkff=darkff)
;
;
;          If the user does NOT have a local XRT data archive, a new program
;          <xrt_search_network.pro> can be used to set proper environment
;          variables to permit copying over the network of the dark files
;          needed for best dark correction, eg,
;
;          IDL> xrt_search_network,/enable
;          IDL> dark =xrt_med_dark(index,ddates=ddates,darkff=darkff)
;
;          see <xrt_search_network.pro> for more details.
;
;
; NOTE:  Output darks have nyquist ringing already removed.  If the
;        user wishes to recover a "raw" medianed dark (ringing still
;        included), they can take the dark generated from xrt_med_dark
;        and
;    
;        IDL> data_nonyq = no_nyquist(data,index)
;        IDL> dark_raw = dark + (data - data_nonyq)
;
;HISTORY:
; 
progver = 'v2009.Jun.12' ; --- (SSaar) Version 1.0 written.  Based in part on 
;                              xrt_dark_img.pro by R. Kano
progver = 'v2009.Jul.09' ; --- (SSaar) Fixed some minor bugs. 
progver = 'v2010.May.20' ; --- (SSaar) Fixed bug when files not found by
;                               xrt_cat; reduced error messages.
progver = 'v2010.May.28' ; --- (SSaar) Permit use of 512x2048 "strip darks" 
;                               for more recent 1x1 data.  Added more 
;                               explanatory notes. 
progver = 'v2010.Aug.19' ; --- (PGrigis,SSaar) Improved to allow web access of
;                               required darks if not available locally.
;                               Notes added explaining how to 
;                               use web access of darks.  
progver = 'v2010.Dec.21' ; --- (SSaar) Improved ability to use "strip darks"
;                               
progver = 'v2011.Jul.22 ; --- (SSaar) Further update to array position keywords
progver = 'v2012-Apr-16' ;--- (SSaar (SAO)) Fixed problem arising when
;                             bad images are mislabeled as darks.
progver = 'v2013-Feb-28' ;--- (SSaar (SAO)) added quiet mode
;                             bad images are mislabeled as darks.
progver = 'v2013-Jun-25' ;--- (SSaar (SAO)) added check against pathological
;                              dark files with many hits or bad offsets.
;                              Added note about nyquist removal.
progver = 'v2013-Aug-30' ;--- (McCauley (SAO)) fixed quiet flag. 
progver = 'v2013-Feb-04' ;--- (SSaar (SAO)) fixed bug in hit removal 
progver = 'v2017-jul-21' ;--- (KY (MSU)) replace xrt_cat.pro with msu_xrt_cat.pro
;                              to get around "fail to find dark data" problem.
; Reference:
;   xrt_cat
;      To use these, you may have to re-define XRT_DATA to the
;      level0 hinode directory 
;
;
;-

;t11=systime(1)

if n_elements(nday) eq 0 then day=14 else day=nday

if n_elements(nim) eq 0 then nim=5

q_qt=keyword_set(quiet) 

netsearch=0                                    ; set flag for network searching
search_net=getenv('XRT_SEARCH_NETWORK')
if search_net eq 'ENABLE' then netsearch=1

; read nim full dark images around the 1st image.

ind=index[0]
dur=3600.*24*day
st=anytim(anytim(ind.date_obs)-dur,/ECS)
en=anytim(anytim(ind.date_obs)+dur,/ECS)
bin=fix(ind.chip_sum)
na=2048/bin
flagdf = n_elements(darkff) gt 0   ; flag for existence of full-frame dark
flag=1                          ; all is well if flag=1
msu_xrt_cat,st,en,dcat,ofil, $                ; look for full-frame darks
  search=['NAXIS1 = '+strtrim(na,2),'NAXIS2 = '+strtrim(na,2), $
           'EC_IMTYP = 1','CHIP_SUM = '+strtrim(bin,2)],quiet=q_qt

case 1 of               ; set up maximum sig_avg and <rms> limits on median dark
   (bin eq 1): begin
      sav_lim=0.3 
      sig_lim=1.3
   end
   (bin eq 2): begin
      sav_lim=0.6
      sig_lim=1.4
   end
   (bin eq 4): begin
      sav_lim=2.2
      sig_lim=1.8
   end
   else: begin
      sav_lim=2.9
      sig_lim=2.5
   end
endcase

s_per_d = 24.*3600.
t_ant = anytim('2008-01-09T00:00:00.0')    ; time of hi-gain antenna problems 
flagstrip=0                                ; set "strip dark" flag
ncat=n_elements(ofil)
if ncat le 4 then begin                ; if not enough files found first time
   flag=0
   if bin eq 1 and anytim(st) ge t_ant then begin  ; if 1x1 & post-antenna fail
      msu_xrt_cat,st,en,dcat,ofil, $           ; look for 512x2048 strip darks
        search=['NAXIS1 = 512','NAXIS2 = 2048', $
           'EC_IMTYP = 1','CHIP_SUM = 1'],quiet=q_qt
      ncat=n_elements(ofil)
      if ofil(0) ne '' then begin            ; if you find them
         flagstrip=1                         ; set "strip dark" flag
         flag=1                              ; and regular "alls well" flag
      endif 
   endif 
endif 


if flag ne 0 then begin 
   if ncat lt nim then begin                 ; if not enough darks (<nim)
      if q_qt eq 0 then $
        print,'WARNING: Number of full-frame dark images is < '+string(nim) 
      if ncat le 4 then flag=0.           ; if # of images <=4, reject
   endif else begin                       ; if all is well 
      ddates1=dcat.date_obs
      dt=(anytim(ddates1)-anytim(ind.date_obs))/s_per_d
      isdt=sort(abs(dt))                   ; sort in ascending delta time
;      isdt=sort(abs(dt(where(dt ne 0))))   ; exclude dark itself if testing
      isdt=isdt(0:nim-1)
      dcat=dcat[isdt]
      ofil=ofil(isdt)
      ddates1=ddates1(isdt)
      dtmax=max(abs(dt(isdt)))              ; max dt of darks to median
      if q_qt eq 0 then print,strtrim(string(nim),2)+' dark images within ', $ 
        string(dtmax,form='(F6.2)'),' days.'
   endelse
endif 

if flag eq 0 and q_qt eq 0 then $
  print,'Insufficient darks found to median within '+ $ 
  strtrim(string(fix(day)),2)+' days. Using unadjusted model dark.'

ndates=n_elements(ddates)                 ; number of input dates

if flag eq 1 and ndates ne 0 then begin   ; if ddates (dates of last med dark)
   tnow=anytim(ddates1)                 ;  exist, get current dark list times
   told=anytim(ddates)                  ; and get ddate times
   nmatch=0
   for j=0,nim-1 do begin                
      imatch=where(tnow eq told(j),nok) ; check for time matches 
      if nok ne 0 then nmatch=nmatch+1
   endfor
endif else nmatch = 0                   ; if ddates not there, make dark median 

if flag eq 1 then begin
   if nmatch ne nim or flagdf eq 0 then begin ; if # matches<nim, or no input
      qr=0                                    ; set remote read flag = OK
                                              ; darkff, make median dark
      dfile = ofil                             ; convert catalog to file list
      pp=strpos(dfile[0],'d.fits')
      if strmid(dfile[0],pp-2,1) ne '.' then $  ; only for older cat database
         dfile=strmid(dfile,0,pp-2)+'d.fits'
      fexist=file_exist(dfile)
      irem=where(fexist eq 0,nrem)         ; check for non-local (remote) files 
      iloc=where(fexist eq 1,nloc)         ; check for local files 
      if nrem gt 0 then begin               ; if not in local archive...
         if netsearch eq 1 then begin       ; and net searching is enabled...
            if q_qt eq 0 then print,'Not in local archive. Getting data ' + $
                'remotely (will be slower).'
            if flagstrip eq 0 then begin           ; and "normal" situation
               msu_xrt_cat,st,en,dcatr,ofilr, $        ; get full-frame files w URL 
                  search=['NAXIS1 = '+strtrim(na,2),'NAXIS2 = '+strtrim(na,2), $
                   'EC_IMTYP = 1','CHIP_SUM = '+strtrim(bin,2)],/url,quiet=q_qt 
                                        ; to get data from remote site e.g.,SAO
            endif else begin                       ; if you need strip darks..                
               msu_xrt_cat,st,en,dcatr,ofilr, $       ; get 512x2048 strip dark URLs
                 search=['NAXIS1 = 512','NAXIS2 = 2048', $ ; get remote data 
                    'EC_IMTYP = 1','CHIP_SUM = 1'],/url,quiet=q_qt  
            endelse
            dcatr=dcatr[isdt(irem)]                 ; select closest in time
            ofilr=ofilr(isdt(irem))                 ; to data file being prepped
            for j=0,nrem-1 do begin
               break_file,ofil[j],disk,dir,file,suffix ; parse (local) file name
               if file_exist(dir) EQ 0 THEN file_mkdir,dir  ; make dir if DNE 
               if q_qt eq 0 then begin
                  print,'Now retrieving file   '+ofilr[j]      ; tell the user
                  print,'and saving locally as '+ofil[irem(j)] ; what's going on
               endif
               sock_copy,ofilr[j],out_dir=dir          ; copy remote files 
            endfor 
            fileok=file_exist(ofil(irem))          ; check if all files are read
            if total(fileok) eq nrem then qr=0 else qr=1  ; set remote read flag
            if nloc ne 0 then begin                ; now get local files, if any
               dfile_loc=dfile(isdt(iloc))             ; select closest in time
               read_xrt,dfile_loc,dindex_loc,ddata_loc,qabort=q,quiet=q_qt  
                                                       ; read locals &
               read_xrt,ofil(irem),dindex_rem,ddata_rem,qabort=q,quiet=q_qt 
                                                             ; (former) remotes
               dindex=[dindex_rem,dindex_loc]         ; combine local & remote
               ddata=[ddata_rem,ddata_loc]
               isort=sort([irem,iloc])                ; sort index & data back 
               dindex=dindex(isort)                   ; into original (ascending
               ddata=ddata(isort)                     ; delta t) order
            endif else read_xrt,ofil(irem),dindex,ddata,qabort=q,quiet=q_qt 
                                                       ; read all
                                                       ; as (former) remotes
         endif else begin
            if q_qt eq 0 then $ 
              print,'Not in local archive.  Remote search not enabled.' 
                                                       ; else print caution and
            q=1                                        ; set read fail flag
         endelse
      endif else begin
         read_xrt,dfile,dindex,ddata,qabort=q,quiet=q_qt ; get all files locally
      endelse
      if q eq 0 and qr eq 0 then begin ; if read ok, make fullframe drk template
         sd=size(ddata)              ; get dimensions
         con=3.                      ; set histogram trim constant = con
         n1=sd(1)
         n2=sd(2)
         if sd(0) eq 3 then nd=sd(3) else nd=1
         darkc=ddata
         j=0
         while j lt nd and flag eq 1 do begin   ; loop thru darks to clean them
            datj=reform(ddata(*,*,j))            
            indj=dindex(j)
            bin=indj.chip_sum
            datc=no_nyquist(datj,indj)          ; remove nyquist
            datc=lsback_away(datc,indj,modj)    ; temporarily remove ramp
            h=histogram(datc,min=fix(min(datc)))   ; histogram the dark
            nh=n_elements(h)
            xh=findgen(nh)+fix(min(datc))
            hmx=max(h)                     ; find histogram peak
            imx=where(h eq hmx)            ; and where it is 
            imx=imx[0]
            rat=3.3e-5*bin^2. - 1.6e-5     ; magic scaling for hist base~rat*hmx
            ik=where(h le con*hmx*rat,nk)  ; find where hist < con*base; these
                                           ; are outliers <> main core of dark
                                           ; points -reset these ->global median
            mdatc=median(datc)             ; find global median
            if nk gt 0 then begin          ; if there are counts > base level
               ilo=where(ik lt imx,nlo)     ; find bad points with DN<hist_max
               if nlo ne 0 then begin       ; if there are some bad points....
                  ilo=ik(ilo)
                  xlomx=max(xh(ilo))
                  ilo2=where(datc le xlomx,nlo2)
                  if nlo2 ne 0 then datc(ilo2)=mdatc ; set them = current median
               endif
               ihi=where(ik gt imx,nhi)      ; find bad points with DN>hist_max
               if nhi ne 0 then begin        ;  if there are some....
                  ihi=ik(ihi)
                  xhimx=min(xh(ihi))
                  ihi2=where(datc ge xhimx,nhi2)
                  if nhi2 ne 0 then datc(ihi2)=mdatc ; set them = current median
               endif
            endif
            if stdev(datc) gt 5*sig_lim or avg(datc) gt 3*sav_lim then begin 
               iok=where(indgen(nd) ne j,nd0)   ; test for, identify & 
               darkc=darkc(*,*,iok)            ; remove bad darks
               if nd0 lt 4 then begin          ; if < 4 darks remain, use model 
                  flag=0.
                  if q_qt eq 0 then $
                    print,'Too many flawed darks.  Defaulting to model dark.'
               endif
               nd=nd0
            endif else darkc(*,*,j)=datc + modj  ; "cleaned" dark (ramp readded)
            j=j+1
         endwhile 
;stop ; *** HERE
         if flag ne 0 then begin
            darkff=median(darkc,dim=3)        ; median the "cleaned" darks
            ddates=ddates1                    ; set ddates to new one
         endif
      endif else begin                     ; if a read error
         if q_qt eq 0 then print,'Read error. Defaulting to model dark'
         flag=0.
      endelse
   endif else begin                       ; otherwise, if full-frame dark exists
      sd=size(darkff)                     ; check size; if Nx <> Ny...
      if sd(1) ne sd(2) then flagstrip=1  ; its really a strip dark! set flag
   endelse 
endif
n1=na                               ; get Nx & Ny for rebinned full-frame dark  
n2=na 
;stop

if flag eq 1 then begin
   if flagstrip eq 1 then dark1=rebin(darkff,2048,2048) $  ; expand strip dark 
     else dark1=darkff                   ; otherwise use input full-frame dark
                            ; NOTE: rebin above is a fudge to use partial 1x1 
                            ; darks to sub for full-frame when the latter 
                            ; are not present. Not optimal, but it works...
   dark1=rebin(dark1,n1,n2)              ; rebin the median full-frame dark 

; make a dark image for the observation from a full frame dark template

   sx=ind.siz_col/bin                   ; find FOV size  rsiz-> siz
   x0=ind.pos_col/bin                   ; find FOV xmin  rpos -> pos
   if x0 gt 2000 then x0=0
   x1=x0+sx-1                            ; find FOV xmax
   sy=ind.siz_row/bin                    ; rsiz-> siz
   y0=0                                  ; always begin with dark ramp bottom!
   y1=y0+sy-1                            ; find FOV ymax 
   dark=dark1[x0:x1,y0:y1]               ; extract dark for FOV

endif else dark=0.

stripdark=flagstrip
;print,systime(1)-t11

return, dark
end

;
;
;