;+
; NAME: calcphotnoise
; PURPOSE:
; 	Given an image, and a map of the exposure times for that image,
; 	compute the photon noise as a function of position.
;
; INPUTS:
; 	image		an image (normalized to units of DN/1 s)
; 	expmap		an exposure map, in seconds
; KEYWORDS:
; 	gain		Gain in photons/ADU. Default is 10 (for IRCAL)
; 	/readnoise	include read noise too
; OUTPUTS:
; 	Noise map, in the same scaled units as the supplied image.
;
; HISTORY:
; 	Began 2004-04-08 11:55:45 by Marshall Perrin 
; 	2004-05-12 	Added sky level option
;-

function calcphotnoise,image0,expmap,gain=gain,header=header,readnoise=readnoise,$
	ncombine=ncombine,photonnoise=photonnoise,readlevel=readlevel

	if not(keyword_set(gain)) then gain=10. ; photons/ADU
	if not(keyword_set(expmap)) then expmap=1

	image = image0
if keyword_set(header) then begin
	; does the FITS header record a sky level that was subtracted away?
	skylevel = sxpar(header,"SKYLEVEL",count=count)
	if count lt 1 then skylevel=0.0

endif else skylevel = 0.

	image = image+skylevel    ; skylevel is in DN/s - same units as image ought to be in.

 
; 	totalimage = image*expmap ; image after expmap seconds, in DN
; 	photons = totalimage*gain ; total photons after expmap seconds, in photons.
; 	noisephot = sqrt(photons)	
; 	totalnoise = noisephot/gain ; total noise after expmap seconds, in DN. 
; 	SNR = totalimage/totalnoise
; 
; 	;outnoise = image/sqrt(image*expmap*gain)
; 	outnoise2 = image/SNR ; = totalimage/noise*image = 

; This really is the right formula. I've checked it carefully.
;   MDP 2004-05-12
	outnoise = sqrt(image)/sqrt(1.0*expmap*gain)


if keyword_set(readnoise) then begin
	; read noise will be ~ 1.6 counts per 16 CDS reads, at a very rough guess
	noiseperframe = 1.6 ; DN
	noiseperframe = 2.2 ; DN. Maybe better value empirically from 03nov17 darks 1410-1419?
	;noiseperframe=3.40
	if keyword_set(header) and not(keyword_set(ncombine))  then begin
		; does the FITS header record a sky level that was subtracted away?
		ncombine = sxpar(header,"NCOMBINE",count=count)
		if count lt 1 then ncombine=0.0
;		message,/info,'HARD CODED NCOMBINE=80!!!'
;		ncombine=80
	endif

	totalreadnoise = noiseperframe*sqrt(ncombine) ; since they add in quadrature
	;normalizedreadnoise = totalreadnoise/expmap
	normalizedreadnoise=totalreadnoise
	photonnoise=outnoise
	outnoise = sqrt(outnoise^2+(normalizedreadnoise/10)^2)
	readlevel = normalizedreadnoise/1.00/max(expmap,/nan) ; normalize to units of counts/second!
endif
	
	
	;stop
	return,outnoise

end