function INF_POLY_2D, array, p, q, interp, dimx, dimy, $
	nospread = nospread, infval = infval, cubic = cubic, $
	_extra = extra

;+
;
; Wrapper for IDL's POLY_2D function that tries to handle Infs 
; 
; Should be a drop-in replacement for POLY_2D; only difference is that it
; preserves and propagates Infs (POLY_2D turns them into Nans). A little bit
; slower since it checks for Infs at the beginning, but if it doesn't find
; any there should be no additional overhead.
;
; MODIFIED:
;	2015-06-01	: Do not turn NaNs into Infs! (PBoerner)
;
;-

asize = SIZE(array)
if N_ELEMENTS(dimx) eq 0 then dimx = asize[1]
if N_ELEMENTS(dimy) eq 0 then dimy = asize[2]
if N_ELEMENTS(interp) eq 0 then interp = 0
if N_ELEMENTS(cubic) eq 0 then cubic = 0

infs = WHERE(FINITE(array, /infinity), numinf)
if numinf eq N_ELEMENTS(array) then begin
	PRINT, 'INF_POLY_2D: No finite values in array!'
endif

if numinf gt 0 then begin
	; Some Infs present, so handle them 
	atemp = array
	mask = array
	mask[*] = 0.
	mask[infs] = !values.f_nan
	
	tight = KEYWORD_SET(nospread)
	if tight then begin
		if not KEYWORD_SET(infval) then infval = MAX(array[WHERE(FINITE(array), numfin)])
		atemp[infs] = infval
		maskcubic = 0
		maskinterp = 1
	endif else begin
		maskcubic = cubic
		maskinterp = interp
	endelse
	
	result = POLY_2D(atemp, p, q, interp, dimx, dimy, cubic = cubic, _extra = extra)
	infmask = POLY_2D(mask, p, q, maskinterp, dimx, dimy, cubic = maskcubic, _extra = extra)
;	infpix = WHERE(infmask ne infmask, numnewinf)
	infpix = WHERE(FINITE(infmask, /infinity), numnewinf)
	if numnewinf gt 0 then result[infpix] = !values.f_infinity
endif else begin
	; No Infs, so just call poly_2d
	result = POLY_2D(array, p, q, interp, dimx, dimy, cubic = cubic, _extra = extra)
endelse

RETURN, result

end