;+
; NAME:
;	Positivity
; PURPOSE:
;	For use in nonlinear functional fitting or minimization
;	that require a positivity constraint on solution parameters.
;	Take unconstrained parameter array x (usually an image),
;	and map it uniquely and smoothly into positive values.
;	Negative values of x get mapped to interval ( 0, sqrt( EPSILON )/2 ],
;	positive values go to ( sqrt( EPSILON )/2, oo )
;	with derivative approaching unity. Derivative is always 1/2 at x=0.
;	Derivative is used by the MRL deconvolution algorithm.
;	If EPSILON=0 then mapping reduces to truncation > 0.
;	If EPSILON LT 0 then mapping reduces to identity (no change).
; CALLING:
;	px = positivity( x [, EPS=1, /DERIV ] )
; INPUTS:
;	x = any numerical scalar or array.
; KEYWORDS:
;	EPSILON = scalar parameter, small positive number (default = 1.e-8).
;	/DERIVATIVE : causes derivative at x to be computed and returned.
; OUTPUTS:
;	Function returns variable of same type and size as input.
; PROCEDURE:
;	Use function ( x + sqrt( x^2 + EPSILON ) )/2, a rotated hyperbola.
; HISTORY:
;	written: F.Varosi NASA/GSFC 1992, as suggested by R.Pina UCSD.
;-

function positivity, x, DERIVATIVE=deriv, EPSILON=epsilon

	if N_elements( epsilon ) NE 1 then epsilon = 1.e-8

	if keyword_set( deriv ) then begin
		if (epsilon GT 0) then return,(1 + x/sqrt( x^2 + epsilon ))/2 $
				  else if (epsilon EQ 0) then return,( x GT 0 )
		return, replicate( 1.0, N_elements( x ) )
	  endif else begin
		if (epsilon GT 0) then return,( x + sqrt( x^2 + epsilon ) )/2 $
				  else if (epsilon LT 0) then return, x       $
				  else return,( x > 0 )
	   endelse
end