;+
;NAME: ivm_noise_calc
;PURPOSE: calculate a rough signal,noise for the IVM spectral data
;CATEGORY: image processing
;CALLING SEQUENCE: noise_array=ivm_noise_calc(im,hstr)
;INPUTS: im = fltarr(nx,ny,4,nw) spectral data, as output from RIVM.pro
;OPTIONAL INPUT PARAMETERS:
;OPTIONAL KEYWORD PARAMETERS: 
;	do_spectral_clean: do the prefilter, continuum slope clean-up, 
;		telluric-line removal.  Newer versions of RIVM should
;		have this as part of the data-reduction.
;	avint: average intensity calculated and used for normalizing 
;		noise levels from intensity back to fractional.
;	signal: if named, returns the normalized fractional "signal".
;		such that signalout/noise of same [i,q,u,v]
;		gives a nominal, estimated, ROUGH!!! S/N level.
;OUTPUTS:
;	noise: fltarr(nx,ny,4) of fractional noise levels.
;COMMON BLOCKS:
;SIDE EFFECTS:
;RESTRICTIONS:
;PROCEDURE:
;MODIFICATION HISTORY:
;
;	ELW-20100104 - Introduced into New Rivm
;  			taken from ~leka/IDL/KDL_code_test/NEW_RIVM/
;-

FUNCTION ivm_noise_calc,im, hstr,noise,$
   do_spectral_clean=do_spectral_clean,avint=avint,signal=signal,outfile=outfile

; Handle inputs
nx = hstr.spix
ny = hstr.ppix 
np = 4
nw = hstr.n_image_sets

; Wavelength scale
w = hstr.fp_set(0:nw-1) * IVM_FP_STEP( hstr.f_wavelength )
w = w - MIN(w, iwmin)
; wmax = MAX(w, iwmax)				; 20100204-ELW Unused
w = REBIN( w, nw, np )
FOR i=1,np-1 DO w(*,i) = w(*,i) + i
; wfit = REFORM( w(0:nw-2,*), (nw-1)*np )	; 20100204-ELW Unused
; dw = MEDIAN( DERIV( w(*,0) ) )		; 20100204-ELW Unused
; icon = nw - 1					; 20100204-ELW Unused

; since this is lifted from ivm_triplet_fit, which puts
; polarization into intensity units (rather than fractional polarization),
; here we go

IF keyword_set(do_spectral_clean) THEN BEGIN 

; Pre-clean the Stokes spectra  =============================================
; Prefilter correction
   tmp = IVM_PREFILTER( REFORM( im(*,*,0,*) ), hstr )
; Telluric line correction
   tmp = IVM_TELLURIC( tmp, hstr )
; Trim up continuum slope
   lsp = FLTARR(nw)
   FOR i=0,nw-1 DO lsp(i) = MEDIAN( tmp(*,*,i) )
   t2 = WHERE( lsp GT MEDIAN(lsp) )
   res = POLY_FIT( w(t2,0), lsp(t2), 1 )
   tc = w(*,0)*res(0,1) + res(0,0)
   tc = tc / AVG(tc)
   FOR k=0,nw-1 DO BEGIN
           tmp(*,*,k) = tmp(*,*,k) / tc(k)
   ENDFOR
   tim = FLTARR(nx, ny, np, nw)
   tim(*,*,0,*) = tmp
   FOR i=1,np-1 DO tim(*,*,i,*) = im(*,*,i,*)*tmp
   avint = REBIN( tmp, nx, ny )
ENDIF ELSE BEGIN
   tim=im
   FOR i=0,nx-1 DO BEGIN
      FOR j=0,ny-1 DO BEGIN
         FOR k=1,3 DO BEGIN
	    tim[i,j,k,*]=tim[i,j,k,*]*tim[i,j,0,nw-1]
	 ENDFOR
      ENDFOR
   ENDFOR
   avint=REBIN( tim[*,*,0,*],nx,ny)
ENDELSE

signal = FLTARR(nx,ny,np)
noise = FLTARR(nx,ny,np)
ilc = FLTARR(nx,ny)
wlc = FLTARR(nx,ny)
wlcerr=FLTARR(nx,ny)

; Generate the maps of polarized signal, noise, line center intensity, position
FOR j=0,ny-1 DO BEGIN
        FOR i=0,nx-1 DO BEGIN
               y = TRANSPOSE( REFORM( tim(i,j,*,*) ) )
               FOR k=0,np-1 DO BEGIN
                   s = SMOOTH( y(*,k), 5 )
                   signal(i,j,k) = STDEV( s(2:nw-3) )
                   noise(i,j,k) = STDEV( y(2:nw-3,k) - s(2:nw-3) )
               ENDFOR
               ilc(i,j) = MIN( y(nw/4:3*nw/4,0), iwm )
               i0 = iwm + nw/4
               res = POLY_FIT( w(i0-3:i0+3,0)-w(i0,0), $
                  y(i0-3:i0+3,0), 2 ,yfit,yband,sigma, $
                  STATUS=status, /DOUBLE)
	       cc=res(0) & bb=res(1) & aa=res(2)
	       wtmp= -1*bb/(2.*aa)
	       wlc(i,j) = w(i0,0) + wtmp
; Get +/- x-axis error in wavelength
               wwtmp=cc + bb*wtmp + aa*wtmp^2 + sigma
               wlcerr(i,j) = SQRT((wwtmp-cc+(bb^2/(4*aa)))/aa)
        ENDFOR
ENDFOR
;snr = signal/noise
;mns = MEDIAN( snr(*,*,0) )
;mni = MEDIAN( noise(*,*,0) )
;mnq = MEDIAN( noise(*,*,1) )
;mnu = MEDIAN( noise(*,*,2) )
;mnv = MEDIAN( noise(*,*,3) )
FOR i=0,3 DO BEGIN
   noise[*,*,i]=noise[*,*,i]/avint
   signal[*,*,i]=signal[*,*,i]/avint
ENDFOR

if keyword_set(outfile) then begin
    ;print, '    Noise file:  ', outfile
    openw,  lun, outfile, /get_lun	; ELW20100304 -removed /F77
    writeu, lun, noise			; Write noise to a file
    close,  lun & free_lun, lun
endif

return,noise
END