;+
; NAME: stddevarr
; PURPOSE: compute the standard deviation for each pixel in a cube
;
;NOTES:
;	sd = stddevarr(array,/badval,alsobad=<>)
; 
; 	Just as medarr returns the median at each pixel of an image cube,
; 	this routine returns the standard deviation for each pixel in a
; 	cube. This is useful when combining a stack of images if you want
; 	an estimate of the error at a given pixel. 
;
;   You can specify a badval which is used to mask out all pixels
;   in the cube that have that value and only take the standard deviation
;   of the remaining ones. This is useful if the images do not overlap
;   completely and there are blank regions on some of the layers - just set 
;   the blank regions to BADVAL and set the badval keyword likewise.
;
; INPUTS:
; 	cube	an image cube, M*N*K
; KEYWORDS:
; 	BADVAL	set this equal to the pixel value of pixels to ignore.
; 	ALSOBAD	in case you have a second bad pixel value.
; 	DIM=	which dimension to take the stddev over. Default=3
; 	/DIVSQRT	Divide by the square root of the number of good valuse at
; 				each pixel. In other words, return the standard deviation of
; 				the mean, instead of just the standard deviation of the data
; 				set.
; OUTPUTS:
; 	returns an M*N array containing the stddev at each pixel
;
; HISTORY:
; 	Began 2002-06-21 14:51:12 by Marshall Perrin, based on
; 		errmap_neg.pro by Erik Rosolowsky
; 	2003-07-16  Modifed keyword checking to allow '0' to be specified
; 		as a bad value.
; 	2006-04-17	added dim= keyword; can now work on 2D images.
; 	2006-06-01	Added recursive chunking of large cubes
;-

function stddevarr,cube,badval=badval,alsobad=alsobad,dim=dim,divsqrt=divsqrt
	
sz = size(cube)

	if ~(keyword_set(dim)) then dim=3
	if dim eq 3 then nz = sz[3] else nz=1


		if dim eq 3 and n_elements(cube) gt 3e6 then begin
			; break the problem down recursively for large arrays.
			; This prevents it from using absolutely humongous amounts of
			; memory. Cutoff at 3 million elements is chosen arbitrarily.
			; and will certainly vary depending on the machine in question
			; and properties of the array!
			nx1 = sz[1]/3
;			print,"Breaking down into smaller chunks: nx = "+strc(nx1)
			return, [stddevarr( cube[0:nx1,*,*], badval=badval,alsobad=alsobad,dim=dim,divsqrt=divsqrt ), $
					 stddevarr( cube[nx1+1:2*nx1,*,*], badval=badval,alsobad=alsobad,dim=dim,divsqrt=divsqrt ), $
					 stddevarr( cube[2*nx1+1:*,*,*], badval=badval,alsobad=alsobad,dim=dim,divsqrt=divsqrt )]

		endif

if (n_elements(badval) gt 0) then begin
	if (n_elements(alsobad) gt 0)then begin
		wgood=where((cube ne badval) and (cube ne alsobad) and finite(cube),goodcount)
	endif else begin
		wgood=where(cube ne badval and finite(cube),goodcount)
	endelse
	if goodcount le 0 then message,"No good pixels!"
endif else wgood = where(finite(cube),goodcount)

; are all pixels bad? if so, stddev = 0
if goodcount eq 0 then return, fltarr(sz[1],sz[2])

goodmask = cube*0
goodmask[wgood]=1

if (size(goodmask))[0] lt dim then return,fltarr(sz[1],sz[2])

tg= total(goodmask,DIM,/NaN)

means = total(cube*goodmask, DIM,/NaN)/tg
meancube = rebin(means,sz[1],sz[2],nz)

result = sqrt( total(((cube-meancube)^2)*goodmask, DIM,/NaN)/(tg-1) )
if keyword_set(divsqrt) then result /= sqrt(tg)
return, result

end