;+
;
; Project   : s3drs reconstruction software
;                   
; Name      : p_simplenoise
;               
; Purpose   : part of Pixon, part of the s3drs GUI
;               
; Explanation: part of the Pixon wrappers used to access and
; manipulate the various arcane Pixon structures and objects.
;               
; returns a quick sqrt(photon) noise for a given pxwrapper set
; Incoming dataset(s) can be a single array or a set of arrays
; (and can be pointers or actual arrays)
;
; Use       : typically called from 's3drs.pro' 
;    
; Inputs    : (various)
;               
; Outputs   : (various)
;
; Keywords  : (various)
;               
; Common    : none
;               
; Restrictions: none
;               
; Side effects: none
;               
; Category    : reconstruction, pixon
;               
; Prev. Hist. : None.
;
; Written     : Sandy Antunes, NRC, March-April 2005-2009
;               
;-            
; test case:
; quickcme,64,'20071231_022220','20071231_012220',/cor2,/setup,/checksn,/debug

; returns a quick sqrt(photon) noise for a given pxwrapper set
; Incoming dataset(s) can be a single array or a set of arrays
; (and can be pointers or actual arrays)
;
; requires either 'pxwrapper' entity or MSBperDN and DNperPhoton,
; otherwise assumes incident data is already in photons.
;
; If given one dataset (single, array of data, ptr, or ptr array),
;    returns its noise.
; If given two datasets (each a single, array of data, ptr or ptr array),
; returns their total noise presuming they were subtracted
;   (e.g. a difference image or bckg subtraction).
; Assumed units are MSB (and uses MSBperDN or pxwrapper.sr.MSBperDN to enforce)
; If data units are DN and MSBperDN is 1, then also works.
;
; If you give pxwrapper to get the constants but the incoming data
; is only a singlet and not the same as pxwrapper.sr.Nd, then you
; really should give the 'iele' element, which otherwise defaults
; to using the 0th element of pxwrapper.sr.  'iele' should only
; be used when processing single data frames, not sets of data.
;
; Optional flag 'photons' returns noise in photons, not in the
; same units as the original input data.
;
; Also adds the given biasdev, if any (as input or from pxwrapper).
;  Alternately, it will default to a Cor2 readout bias of 0.019DN
;
; As an alternative mode, it accepts input as DN (as per the
; lab's pre-made total brightness images in DN) and does the
; appropriate calculation into MSB->photons, then returns a
; noise in the same input DN units.
;

FUNCTION p_simplenoise,dataM,pxwrapper,dataN

  mcheck=max(dataM)
  ;print,'datamax is ',mcheck
  noise=dataM*0.0
  for i=0,n_elements(dataM[0,0,*])-1 do begin
      
      dat=dataM[*,*,i]
      ie=where (dat le 0)
      if (ie[0] ne -1) then dat[ie]=0

      n1= (1.0/pxwrapper.sr.MSBperDN[i]) * (1.0/pxwrapper.sr.DNperPhoton[i])
      n1 = n1 * pxwrapper.sr.exptime[i]
      if (pxwrapper.sr.polar[i] eq 'b') then $
        n1 = n1 * 3.0           ; for tB
      ;n1 = dataM[*,*,i] * n1
      n1 = dat * n1
      n1=sqrt(n1) * pxwrapper.sr.DNperPhoton[i]
      n1 = sqrt (n1*n1 + pxwrapper.sr.biasdev[i]*pxwrapper.sr.biasdev[i])
      n1 = n1 * (1.0/pxwrapper.sr.exptime[i]) * pxwrapper.sr.MSBperDN[i]
      if (pxwrapper.sr.polar[i] eq 'b') then $
        n1 = n1 * (1.0/3.0)     ; for tB

;      print,'scaling is ',pxwrapper.sr.scaling[i]
;      print,mean(n1)
      n1=n1 * sqrt(1.0/pxwrapper.sr.scaling[i]);

;      print,mean(n1)
      n1 = sqrt(n1^2 + 2.0* (0.027*pxwrapper.sr.MSBperDN[i])^2)
;      print,mean(n1)

      noise[*,*,i]=n1

  end

  if (n_elements(dataN) ne 0) then begin
      n2=p_simplenoise(dataN,pxwrapper)
      noise = sqrt (noise*noise + n2*n2)
  end

  return,noise

END