FUNCTION ffilter_err,image_in,index_in

; =========================================================================
;+
; PROJECT:
;       Solar-B / XRT
;
; NAME:
;       FFILTER_ERR
;
; CATEGORY:
;       Data calibration
;
; PURPOSE:
;       Estimates maximum errors due to uncertainties in the Fourier
;       filter in <xrt_prep> (if used) in DN/s 
;
; CALLING SEQUENCE:
;       error_image_out = FFILTER_ERR (image_in , index_in) 
;
; INPUTS:
;       IMAGE_IN  - [Mandatory] (2-dim float array, [Nx,Ny])
;                   The prepped, normed image containing read-out signals.
;       INDEX_IN  - [Mandatory] (Index structure corresponding to image_in)
;
; OUTPUTS:
;       ERROR_IMAGE_OUT - [Mandatory] (2-dim float array, [Nx,Ny])
;                   The image of the Fourier filter errors (in DN/s).
;
; EXAMPLES:
;
;       Basic usage:
;       IDL> error_image_out = ffilter_err(image_in, index_in)
;
;
;  COMMON BLOCKS:
;       None.
;
; NOTES:
;
;       1) a set of 56 2048^2 1x1 images before contamination spots
;          and 56 images in each spot epoch were studied
;       2) Does no error checking! (this is done by xrt_clean_ro 
;          which calls it)
;
; MODIFICATION HISTORY:
;
progver = 'v2012.Jul.03' ; --- (SSaar) Written by SSaar.  
;
progver = 'v2013.Jan.20' ; --- (SSaar) Modified to include 2 epoch spot
;                                      contamination effects; streamlined
progver = 'v2013.Mar.05' ; --- (SSaar) Modified adjust smoothing for smaller
;                                      images
;
;-
; =========================================================================

;===== Get dimensions and parameters of image.
  nx = index_in.naxis1
  ny = index_in.naxis2
  chip_sum = index_in.chip_sum   
;===== correct exposure time if exposure occurs near clock rollover
  t_exp = index_in.exccdex/1e6  
  tmx = 1677721500L
  if t_exp lt 0 then t_exp =(tmx -index_in.e_sh_ope +index_in.e_sh_clo)/1e6


; after 23-JUL-2007 09:10 is contam epoch 1, 01-feb-2008 00:00 is epoch 2
  t0=anytim2jd(index_in.date_obs)
  tspot1=anytim2jd('23-jul-2007 09:10')
  tspot2=anytim2jd('01-feb-2008 00:00')
  t0a=t0.int+t0.frac
  tspot1a=tspot1.int+tspot1.frac
  tspot2a=tspot2.int+tspot2.frac
 
  ader = avg(abs(deriv(image_in)))    ; average amplitude of gradients
  av=avg(image_in)                    ; average DN/s
  case 1 of
     (t0a lt tspot1a): begin          ; if before contam spot
        ccon = [-0.61881,1.22440]     ; power law coeffs relating const level  
                                      ; w/ <|deriv|> of data (sig = 0.224 dex)
        cnorm=  [1.42217,-3.39947]    ; norm w/ <|deriv|>/<I>^0.5 (sig = 0.463)
        csmoo = [1.60371,-0.52988]    ; powerlaw relating <deriv>/<data> 
                                      ; to nsmoo (sig = 0.263 dex)
        xcon =  ader                  ; con dependence
        xnorm = ader/sqrt(av)         ; norm dependence
        xsmoo = ader/av               ; smoo dependence 
      end
     (t0a ge tspot2a): begin          ; contam spot #2
        ccon = [-0.58548,1.17931]     ; power law coeffs relating const level
                                      ; w/ <|deriv|> of data (sig = 0.160 dex)
        cnorm = [1.89582,0.58744]     ; norm w/ <I>  (sig = 0.475)
        csmoo = [1.45293,-0.32464]    ; smoo w/ <|deriv|>/<I>  (sig = 0.447)
        xcon =  ader                  ; con dependence
        xnorm = av                    ; norm dependence
        xsmoo = ader/av^1.5           ; smoo dependence 
     end
     else: begin                      ; contam spot #1 *
        ccon = [-0.59384,1.19208]     ; power law coeffs relating const level  
                                      ; w/ <|deriv|> of data (sig = 0.183 dex)
        cnorm = [1.88666,0.545614]    ; norm w/ <I>  (sig = 0.368)
        csmoo = [1.40854,-0.536169]   ; smoo w/ <|deriv|>/<I>  (sig = 0.438)
        xcon =  ader                  ; con dependence
        xnorm = av                    ; norm dependence
        xsmoo = ader/av
     end
  endcase
  con = 10.^poly(alog10(xcon),ccon)                 ; power law for const
  norm = 10^poly(alog10(xnorm),cnorm)               ; power law for norm
  nsmoo = round(10^poly(alog10(xsmoo),csmoo))       ; power law for nsmoo, 
  nsmoo=nsmoo/chip_sum>3<512                        ; adjust for binning
                                                    ;   with limits
  base = 50.                                        ; threshold level
  nmin=min([nx,ny])                                 ; min dimension
  if nsmoo le nmin/4. then begin                    ; test size of nsmoo, if ok
     ds = smooth(image_in>base,nsmoo,/edge_trunc)
     for j=0,2 do ds=smooth(ds, nsmoo,/edge_trunc)  ; smooth 
  endif else begin                                  ; otherwise...
     big=fltarr(2048,2048)                          ; make a full frame image
     mask=image_in*0 + 1                            ; make 3 pixel wide boundary
     mask(3:nx-4,3:ny-4)=0                          ; mask of image
     mededg=median(image_in(where(mask eq 1)))      ; fill full-frame w/median 
     big=big+mededg                                 ;  of image in boundary mask
     big(1024-nx/2,1024-ny/2)=image_in              ; drop image into full-frame
     ds = smooth(big>base,nsmoo,/edge_trunc)             
     for j=0,2 do ds=smooth(ds, nsmoo,/edge_trunc)  ; smooth as normal 
     ds=ds(1024-nx/2:1024+nx/2-1,1024-ny/2:1024+ny/2-1) ; extract smoothed image
  endelse
  ds = ds - min(ds)                                 ; bring min -> 0

  error_model = con + ds/norm                       ; compute model from
                                                    ; derived parameters
  error_model = error_model/(chip_sum)^1.5          ; binning correction

return, error_model 
end

;
;
;