FUNCTION no_nyquist,image_in,index_in,fit

; =========================================================================
;+
; PROJECT:
;       Solar-B / XRT
;
; NAME:
;       NO_NYQUIST
;
; CATEGORY:
;       Data calibration
;
; PURPOSE:
;       Remove prevalent Nyquist frequency ringing in x
;
; CALLING SEQUENCE:
;       image_out = NO_NYQUIST (image_in , index [,fit]) 
;
; INPUTS:
;       IMAGE_IN  - [Mandatory] (2-dim float array, [Nx,Ny])
;                   The image containing read-out signals.
;       INDEX_IN   - [Mandatory] Index structure coresponding to image_in
;
; OUTPUTS:
;       IMAGE_OUT - [Mandatory] (2-dim float array, [Nx,Ny])
;                   The image with the fit to nyquist ringing removed.
;       FIT       - [Optional] (2-dim float array, [Nx,Ny])
;                   The fit to the nyquist ringing which is subtracted
;                    from the data
;
; EXAMPLES:
;
;       Basic usage:
;       IDL> image_out = no_nyquist(image_in, index_in)
;
;       You want the fit for some reason:
;       IDL> image_out = no_nyquist(image_in, index_in, fit)
;
;
;  COMMON BLOCKS:
;       None.
;
; NOTES:
;
;       1) the average rms of the y-average (i.e., avg(dark,1)) of 57 1x1 darks 
;          went from 1.347 to 0.165 DN after processing with this routine 
;       2) Does no error checking! (this is done by xrt_clean_ro 
;          which calls it)
;
; MODIFICATION HISTORY:
;
progver = 'v2007.Mar.06' ; --- (S. Saar, M. Weber) Written by SS. Cleaned up 
;                                 by MW.
progver = 'v2008.Oct.08' ; --- (S. Saar) Corrected phase shifts for binned 
;                                 data, dropped correction for chip_sum=8
progver = 'v2013.Feb.17' ; --- (S. Saar) Corrected for phase shifts when
;                                 opposite readport is used.
;
;
;-
; =========================================================================

;===== Get dimensions & chip_sum of image.
  n1 = index_in.naxis1
  n2 = index_in.naxis2
  chipsum = index_in.chip_sum 
  x1 = findgen(n1)

;===== The (apparently constant) phase offset.
  case 1 of
    (chipsum eq 1): phi = -0.5
    (chipsum eq 2): phi = 0.5
    (chipsum eq 4): phi = 0.5
     else:
  endcase

;===== Switch phase 180 deg of using the other readport.
  rp = index_in.readport
  if rp eq 'L' then phi = -phi

;=== nyquist negligible for chipsum=8      

  if chipsum lt 8 then begin

;=== Find and subtract fitted signals =================================

;===== Find and prepare to exclude areas of large abs(deriv(image)),
;===== where image has first been averaged over all y.
     dimy = deriv(avg(image_in,1))
     ik = where(abs(dimy) le median(abs(dimy)))

;===== Compute amplitude from median of difference of odd x vs.
;===== even x (avg over all y) with high deriv areas excluded.
     diff = avg(image_in,1) - shift(avg(image_in,1),-1)
     amp = median(abs(diff(ik)))/2.


;===== Set up nyquist sin wave, period = 2 pixels.
     fit0 = amp*sin(2*!pi*(x1+phi)/2.0)

;===== Make 2D.
     fit = rebin(fit0,n1,n2)

;===== no correction for chipsum=8
   endif else fit = image_in*0 

;===== Subtract, and you're done!
  image_out = image_in - fit 

  return, image_out
end
;
;
;