pro stisnoise,file,freq,magnitude,exten=exten,outfile=outfile,db=db,$
dc=dc,digfilter=digfilter,median=median,boxcar=boxcar,lee=lee,wipe=wipe,$
window=window

 if n_params(0) lt 3 then begin
  print,'stisnoise, file,freq,magnitude,exten=exten,outfile=outfile,db=db,'
  print,'dc=dc,digfilter=digfilter,median=median,boxcar=boxcar,lee=lee,'
  print,'wipe=wipe,window=window'
  return
 endif

 ; Description: computes an FFT on STIS CCD frames to evaluate fixed
 ; pattern noise.  Fixed pattern noise is most obvious in FFT of bias frames.
 ; Optional filtering to correct the fixed pattern noise is provided
 ; through keywords median, boxcar, lee, digfilter, wipe, & window.  Filtered
 ; data can be saved as an output file.

 ; Author: Thomas M. Brown (STScI)

 ; Input:
 ; file= STIS FITS file (or database number if db=1)

 ; Output:
 ; freq = frequency in power spectrum (hz)
 ; magnitude = magnitude in power spectrum

 ; Optional input (keywords):
 ; exten = fits extension to be read (or readout if db=1)
 ; dc:  when 0 (default), the power in the first freq bin is set to zero
 ;      to allow easier plotting of the power spectrum. 
 ; digfilter = 4-element array specifying the digital filter parameters:
 ;          flow, fhigh, power relative to Gibbs (50 is good), nterms.
 ;          Note: flow>fhigh for band stop, flow<fhigh for bandpass.
 ;          See IDL help on digital_filter for more information.
 ; median = width of median filter applied.
 ; boxcar = width of boxcar smoothing to be applied.
 ; lee = 2-element array specifying width and sigma in the Lee filter
 ;        (note that actual width is set to 2*width+1).  See
 ;        IDL help on leefilt.  Note that leefilt has a bug
 ;        in IDL as of now (version 5.5), such that the loop variables
 ;        need to be changed to type long.  Thus you need a corrected version
 ;        of leefilt.pro for this option to work. 
 ; wipe = 3-element array, specifying how to modify the data in frequency 
 ;        space. If set, the image is converted to a 1-D time
 ;        series, fourier transformed to frequency space, modified,
 ;        inverse transformed back to time space, and converted back
 ;        to a 2-D image.  The first and second elements specify the range
 ;        in frequencies to be scaled (in hz), and the third
 ;        element specifies the scaling factor (should be 0-1).
 ; window = a 3 element array, specifying how to modify the data in frequency 
 ;        space.  The first element is the center of the window (in hz).  
 ;        The second element is the width of the window (in hz).  The
 ;        third element controls the tapering of the window - it is the
 ;        scale (in hz) of the tapering width.  Specifically, a square 
 ;        bandstop is convolved with a gaussian having the FWHM 
 ;        given by the third parameter.

 ; db = set to 1 if reading from STIS IDT Database instead of file.  Useful
 ;        for members of the STIS IDT (GSFC) only.

 ; Optional output (keyword):
 ; outimage = filtered image

 ; history:
 ; 11/5/2001  TMB - version 1.  Basic idea comes from ACS analysis software
 ;                  used for analyzing read noise 
 ;                  (dino.pro; W.J. McCann & G. Hartig)
 ; 11/6/2001  TMB - version 2 added other amps, error checking
 ; 11/6/2001  TMB - version 2.1 added check on sci ext
 ; 11/6/2001  TMB - version 3 added various filter options
 ; 11/6/2001  TMB - version 3.1 added ability to read from STIS IDT DB
 ; 11/7/2001  TMB - version 4 added wipe filter and output images
 ; 11/9/2001  TMB - version 4.1 added new parameter to wipe routine, changed
 ;                  output to a file with header preservation.
 ; 11/20/2001 TMB - version 5.0 added window routine, which does the
 ;                  filtering of "wipe" with a more gradual scaling
 ; 11/26/2001 TMB - version 5.1 cleaned up the code comments.

 if not keyword_set(exten) then exten=1
 if not keyword_set(dc) then dc=0
 if not keyword_set(db) then db=0

 ; DEFINE PHYSICAL CHARACTERISTICS OF STIS CCD
 pst = 640.0 ; parallel shift time (us)
 sst = 22.0 ; serial shift period (us)
 nc0 = 1062 ; number of columns in raw data
 nr0 = 1044 ; number of rows in raw data
 fltxy = 1024 ; number of columns and rows in calibrated data
 nos = 19 ; number of physical overscan columns
 pps = pst / sst  ; number of serial shift intervals in parallel interval

 ; RETRIEVE EXPOSURE INFORMATION FROM HEADER
 if db eq 0 then begin
  inimage=readfits(file,extheader,exten=exten,/silent)
  himage=readfits(file,header,exten=0,/silent)
  extname=strtrim(sxpar(extheader,'EXTNAME'))
 endif else begin
  stis_read,file,header,inimage,readout=exten,/no_overscan
  !p.title=''
  extname="SCI"
 endelse
 amp=strtrim(sxpar(header,'CCDAMP'))
 gain=strtrim(sxpar(header,'CCDGAIN'))
 targ=strtrim(sxpar(header,'TARGNAME'))
 print,'Amp = ',amp,', target = ',targ,', gain = ',gain,f='(a,a,a,a,a,i1)'

 ; CHECK TO ENSURE THE SCI EXTENSION IS BEING USED
 if extname ne "SCI" then begin
  print,'You should only run this on a science extension.'
  print,'Extension read: ',extname
  print,'Aborting.'
  return
 endif

 image_sz = size( inimage )
 nc = image_sz[1]
 nr = image_sz[2]
 if (nc eq nc0 and nr eq nr0) then type='raw' else $
 if (nc eq fltxy and nr eq fltxy) then type='flt' else begin
  print,'This program should be run on 1062x1044 or 1024x1024 data only.'
  return
 endelse

 ; PAD DATA WITH FAKE "OVERSCAN" IF DATA HAVE BEEN OVERSCAN TRIMMED
 if type eq 'flt' then begin
  temp=fltarr(nc0,fltxy)
  for row=0,fltxy-1 do begin
   temp[*,row]=median(inimage[*,row])
  endfor
  temp[nos:(nc0-nos-1),*]=inimage
  nc=nc0
 endif else begin
  temp=inimage
 endelse

 ; TRANSLATE FRAME SO THAT IT IS IN READOUT ORDER
 ; amp D data -> rotate by 180 degrees
 if amp eq 'D' then image=rotate(temp,2) else $
 ; amp A data -> leave as is
 if amp eq 'A' then image=temp else $
 ; amp B data -> flip left<->right
 if amp eq 'B' then image=reverse(temp,1) else $
 ; amp C data -> flip top<->bottom 
 if amp eq 'C' then image=reverse(temp,2) else begin
 print,'No amplifier given in header.  Aborting.'
 return
 endelse

; CONVERT 2-D ARRAY TO 1-D TIME SERIES
 nx = nc + pps
 time_series = fltarr(nx*nr)
 ds = fltarr( pps )
 for i=0, nr-1 do begin
  time_series[i*nx] = image[*,i] 
  time_series[i*nx+nc]     = ds + median( image[*,i] ) ; pad dead-time
; (note that non-integer nx prevents phase wandering)
 endfor

; BEGIN FILTERING OPTIONS **************

 ; frequency filtering using digital_filter algorithm
 if keyword_set(digfilter) then begin
  nyqf=1.0/(2*sst*1.0e-6)   ; Nyquist frequency
  flow=digfilter[0]/nyqf       ; define frequency boundaries 
  fhigh=digfilter[1]/nyqf
  powgibbs=digfilter[2]
  nterms=digfilter[3]
  print,"Filtering with flow = ",flow," and fhigh = ",fhigh,$
  " (in terms of Nyquist frequency)",f='(a,f5.3,a,f5.3,a)'
  coeff=digital_filter(flow,fhigh,powgibbs,nterms)  
  time_series=convol(time_series,coeff)
 endif

 ; median filtering
 if keyword_set(median) then begin
  time_series=median(time_series,median)
 endif

 ; boxcar smoothing
 if keyword_set(boxcar) then begin
  time_series=smooth(time_series,boxcar)
 endif

 ; Lee filter
 if keyword_set(lee) then begin
  time_series=leefilt(time_series,lee[0],lee[1],/exact)
 endif

 ; wipe frequencies in fft and reverse transform
 if keyword_set(wipe) then begin
  ntime=n_elements(time_series)
  ; if ntime is an odd number the fft will take forever so make
  ; it even with a small prime factor sum (not as quick as power of 2 but 
  ; still much quicker).
  ; Note that padding data out to next power of two is too much padding
  ; (number of elements is just a bit over 2^20 for STIS data)
  if type eq 'raw' then ntimep=ntime+13 else ntimep=ntime+7
  t2=fltarr(ntimep)
  t2[0]=time_series
  freq=findgen(ntimep)/(ntimep * sst * 1.0e-6)
  freq[(ntimep/2+1):(ntimep-1)]=reverse(freq[1:(ntimep/2-1)])
  tran=fft(t2)
  ind=where(freq gt wipe[0] and freq lt wipe[1])
  tran[ind]=tran[ind]*wipe[2]
  inv=fft(tran,/inverse)
  time_series=(float(inv))[0:ntime]
 endif

 ; window filter
 ; gradually wipe frequencies in fft and reverse transform
 if keyword_set(window) then begin
  ntime=n_elements(time_series)
  ; if ntime is an odd number the fft will take forever so make
  ; it even with a small prime factor sum (not as quick as power of 2 but 
  ; still much quicker).
  ; Note that padding data out to next power of two is too much padding
  ; (number of elements is just a bit over 2^20 for STIS data)
  if type eq 'raw' then ntimep=ntime+13 else ntimep=ntime+7
  t2=fltarr(ntimep)
  t2[0]=time_series
  freq=findgen(ntimep)/(ntimep * sst * 1.0e-6)
  freq[(ntimep/2+1):(ntimep-1)]=reverse(freq[1:(ntimep/2-1)])
  tran=fft(t2)
  filter=fltarr(ntimep)+1.0
  ind=where(freq gt (window[0]-window[1]/2.0) and $
      freq lt (window[0]+window[1]/2.0))
  filter[ind]=0.0
  freqstep=1.0/(ntimep * sst * 1.0e-6)
  width=window[2]/freqstep ; specify window width in freq steps
  sigma=width/2.354820044  ; convert fwhm to sigma
  kernw=fix(5*sigma)      ; make kernel have width of 5 sigma
  if (kernw mod 2) eq 0 then kernw=kernw+1 ; make kernel odd
  kernx=findgen(kernw)   
  gauss,kernx,kernw/2,sigma,1.0,kerny  ; gaussian kernel
  kerny=kerny/total(kerny)    
  filterc=convol(filter,kerny,/edge_truncate,/center)
  tran=tran*filterc
  inv=fft(tran,/inverse)
  time_series=(float(inv))[0:ntime]
 endif
     
; END FILTERING OPTIONS

; RECREATE 2-D IMAGE FROM TIME SERIES
outimage=fltarr(nc,nr)
for i=0,nr-1 do begin
 outimage[*,i]=time_series[i*nx:(i*nx+nc-1)]
endfor
if type eq 'flt' then outimage=outimage[nos:(nc0-nos-1),*]

; RESTORE ORIGINAL IMAGE ORIENTATION
; amp D data -> rotate by 180 degrees
if amp eq 'D' then outimage=rotate(outimage,2) else $
; amp A data -> leave as is
if amp eq 'A' then outimage=outimage else $
; amp B data -> flip left<->right
if amp eq 'B' then outimage=reverse(outimage,1) else $
; amp C data -> flip top<->bottom 
if amp eq 'C' then outimage=reverse(outimage,2) 

; TRIM VECTOR TO POWER OF 2 FOR FFT 
; (this is the fastest fft calculation but it doesn't preserve all
; data, as needed in wipe routine above)
 p2 = fix( alog( nx*nr ) / alog(2) )
 n_ts = 2L ^ p2
 time_series = time_series[0:n_ts-1]

; PERFORM FFT AND RETURN FIRST HALF
 nt = n_ts / 2
 fft_output = fft( time_series )
 magnitude = (abs( fft_output ))[0:nt-1]
 real_part = float( fft_output )
 imag_part = imaginary( fft_output )
 freq = findgen(nt) / (2 * nt * sst * 1.0e-6)
 if dc eq 0 then magnitude[0]=0   ; set first bin in power spectrum 
                                  ; to zero unless dc=0

 if keyword_set (outfile) then begin
  if (db) then begin
   writefits,outfile,outimage,header ; write header & data from IDT stis_read
  endif else begin
   writefits,outfile,0,header       ; write primary header then append ext
   writefits,outfile,outimage,extheader,/append
  endelse
 endif

 return
end