;+
; NAME:
;
;   CORRELATED_LINE
;
; PURPOSE:
;
;   Monte Carlo ramp fitting with correlated (Poisson) and read noise,
;   comparing a variety of slope estimators.
;
; CATEGORY:
;
;   IRS Ramp fitting
;
; CALLING SEQUENCE:
;
;   correlated_line,[SAMPLES=,SAMPLE_TIME=,READ_NOISE=,N=,RATE=,
;                    /POSTSCRIPT,BIN_COUNT=bc]
;
; KEYWORD PARAMETERS:
; 
;   SAMPLES: The number of samples in the ramp.  Default=16.
;
;   SAMPLE_TIME: The amount of time each sample represents.
;      Default=2s.
;
;   READ_NOISE: The read noise *per sample* in electrons.  Default=150e.
;
;   N: The number of Monte Carlo runs to perform for each rate.
;      Default=50000
;
;   RATE: A vector or scalar giving the target slope in e/s.
;      Default=[200,1000,5000] e/s.
;
;   POSTSCRIPT: If set, the plots will be directed to postscript
;      output files.
;
;   BIN_COUNT: The number of bins to span the standard deviation of
;      the fit.  Default=100.
;
; OUTPUTS:
;
;   <none>
;
; SIDE EFFECTS:
;
;   When POSTSCRIPT is set, histogram comparison plots will be saved
;   to files:
;
;     irs_ramp_fit_rate_samptime_nsamp.eps
;
;   where rate is the slope rate in e/s, samptime is the sample time
;   in s, and nsamp is the number of samples in the ramp.  Otherwise,
;   plots will be made to the current device.
;
; EXAMPLE:
;
;   correlated_line,N=20000,BIN_COUNT=40,RATE=[150,2000],READ_NOISE=100.
;
; MODIFICATION HISTORY:
;
;   July 12, 2005 (JD Smith): Written, based on discussions with Karl
;                             Gordon.
;-

pro correlated_line,SAMPLES=samples,SAMPLE_TIME=sample_time, $
                    READ_NOISE=read_noise,N=n,RATE=rate, $
                    POSTSCRIPT=ps,BIN_COUNT=bc

  tvlct,[255b,0b,0b,255b],[0b,255b,0b,0b],[0b,0b,255b,255b],1
  
  if n_elements(rate) eq 0 then rate=[200,1000,5000]
  if n_elements(samples) eq 0 then samples=16
  if n_elements(sample_time) eq 0 then sample_time=2.
  if n_elements(read_noise) eq 0 then read_noise=150. ;e, per sample
  if n_elements(n) eq 0 then n=50000 ;total M.C. repeats
  if n_elements(bc) eq 0 then bc=25
  
  types=['Weighted Fit','Unweighted Fit','Deltas Mean','Deltas 2!7r!X trim']
  
  fits=fltarr(n_elements(types),n,/NOZERO)

  xoff=0.05 
  x=findgen(samples)
  
  nr=n_elements(rate) 
  
  if ~keyword_set(ps) then begin 
     erase
     m=ceil(sqrt(nr))
     !P.MULTI=[0,m,m]
  endif else !P.MULTI=0
     
  for i=0L,n_elements(rate)-1 do begin
     ;; Build the noiseless ramp, with variance N_electrons=rate*sample_time
     ramp=total(randomu(sd,POISSON=rate[i]*sample_time,samples,n), $
                /CUMULATIVE,1)
     ;; impose the read noise on top of the accumulated ramps
     ramp+=read_noise*randomn(sd,samples,n)
     
     for j=0L,n-1 do begin 
        ;; Line fit method
        f=linfit(x,ramp[*,j],MEASURE_ERRORS=sqrt(read_noise^2+ramp[*,j]))
        fits[0,j]=f[1]/sample_time
        ;; No measurement errors line fit method
        f=linfit(x,ramp[*,j])
        fits[1,j]=f[1]/sample_time
     endfor 
     
     ;; Mean of deltas method
     deltas=shift(ramp,-1,0)-ramp
     fits[2,*]=total(deltas[0:samples-2,*],1)/(samples-1)/sample_time
     
     ;; 2-sigma trimmed mean of deltas method
     diff=deltas-rebin(fits[2,*],samples-1,n,/SAMPLE)
     stddev=sqrt(total(diff^2,1)/(samples-2))
     keep=abs(diff) le 2.*rebin(transpose(stddev),samples-1,n,/SAMPLE)
     fits[3,*]=total(keep*deltas,1)/total(keep,1)/sample_time
     
     ;; Plot results to window or eps file
     yoff=0.95
     if keyword_set(ps) then begin 
        file=string(FORMAT='(%"irs_ramp_fit_test_%0.1f_%0.1f_%0d.eps")', $
                    rate[i],sample_time,samples)
        set_plot,'PS' & device,FILENAME=file,/ENCAPSULATED,/COLOR
     endif 
     
     bs=0.
     for k=0,n_elements(types)-1 do begin 
        h=histogram(fits[k,*],BINSIZE=1.,LOCATIONS=l)
        f=gaussfit(l,h,a,NTERMS=3)      
        if bs eq 0. then bs=a[2]/bc
        h=histogram(fits[k,*],BINSIZE=bs,LOCATIONS=l)
        title=string(FORMAT='(%"N!DSamp!N=%0d; Sample Time=%0.1fs; ' + $
                     'RN=%0de;  Slope=%0de/s")', $
                     samples,sample_time,read_noise,rate[i])
        if k eq 0 then begin 
           plot,[0],[0],XRANGE=rate[i]+[-a[2],a[2]]*5, CHARSIZE=1.1,$
                YRANGE=[0,max(h)*1.25], $
                XTITLE="Measured Slope [e/s]", $
                YTITLE="Number",TITLE=title
           plots,rate[i],!Y.CRANGE,LINESTYLE=1
        endif  
        oplot,l+bs/2,h,COLOR=k+1,THICK=2,PSYM=10
        xp=!X.CRANGE[0]+(!X.CRANGE[1]-!X.CRANGE[0])*xoff
        yp=!Y.CRANGE[0]+(!Y.CRANGE[1]-!Y.CRANGE[0])*(yoff-=.05)
        xyouts,xp,yp,COLOR=k+1, CHARSIZE=1.1,$,
               string(FORMAT='(%"%s (%s=%0.2f)")',types[k],'!7r!X',a[2])
     endfor  
     if keyword_set(ps) then device,/close & set_plot,'X'
  endfor
end