;+
; NAME:
;	smei_star_standard
; PURPOSE:
;	Read in the standard star map
; CATEGORY:
;	camera/idl/star
; CALLING SEQUENCE:
	FUNCTION smei_star_standard, camera 		, $
		bkgnd_origin		= bkgnd_origin		, $
		bkgnd_radius		= bkgnd_radius		, $
		wing_origin 		= wing_origin		, $
		wing_radius 		= wing_radius		, $
		wing_delta			= wing_delta 		, $
		min_wing_radius		= min_wing_radius	, $
		psfstretch			= psfstretch		, $
		degrees 			= degrees			, $

		get_origin			= get_origin		, $
		get_ndim			= get_ndim			, $
		get_fov 			= get_fov			, $
		get_scale			= get_scale 		, $

		get_nomansland_id	= get_nomansland_id , $
		get_background_id	= get_background_id , $
		get_wayoutside_id	= get_wayoutside_id , $

		use_filled_psf		= use_filled_psf	, $
		in_psf				= in_psf			, $
		show				= show				, $
		scale				= scale 			, $
		origin				= origin			, $
		extend_corepsf		= extend_corepsf	, $

		silent				= silent
; INPUTS:
;	camera			scalar; type: integer; default: 2
;						camera number (1,2,3)
; OPTIONAL INPUT PARAMETERS:
;	!!! DO NOT MIX the bkgnd_* and wing_* keywords: results will be unpredictable.
;
;	/degrees		all angles are in degrees (default is radians)
;					This includes bkgnd_radius, wing_radius, psfstretch 
;					bkgnd_origin, wing_origin, wing_delta and min_wing_radius
;	bkgnd_origin=bkgnd_origin
;					array[2]; type: float; default: [0,0]
;						only used if bkgnd_radius is set.
;						Shifts the origin of the ellipsoidal background area
;						The units are the same as bkgnd_radius (i.e. radians or
;						degrees, depending on setting of /degrees).
;	bkgnd_radius=bkgnd_radius
;					scalar or array[2]; type: float; default: none
;						Radius of circular or ellipsoidal area around centroid
;						of standard star used to define the background area around the PSF
;						If two radii are specified than the ring in between
;						the two circles/ellipses is used to define the background area.
;						(an ellipse results from the use of 'psfstretch' see below.)
;	wing_origin=wing_origin
;					array[2]; type: float; default: [0,0]
;						only used if wing_radius is set.
;						Shifts the origin of the ellipsoidal wing area
;						The units are the same as wing_radius (i.e. radians or
;						degrees, depending on setting of /degrees).
;	wing_radius=wing_radius
;					scalar or array[2]; type: float; default: none
;						Axes of ellipse defining a wing area around the PSF
;						If scalar then the ellipse becomes a circle.
;						Inside this area the PSF is extended using an
;						exponential dropoff in directions along and
;						perpendicular to the PSF symmetry axis.
;	wing_delta=wing_delta
;					scalar or array[2]; type: float; default: none
;						correction to wing_radius
;	min_wing_radius=min_wing_radius
;					scalar; type: float; default: 0.5 degrees
;						sets a lower value on the effective wing radius
;						(wing_radius+wing_delta) to avoid using too low
;						values (wing_delta is used to scale the effective
;						wing radius with magnitude; this could lead to
;						small, or even negative values for very faint stars).
;	psfstretch=psfstretch
;					array[2]; type: float: default: [0.0,0.0]
;						Only used if bkgnd_radius or wing_radius is set
;						to non-zero value
;						Setting psfstretch to a value other than one will turn the
;						background (if bkgnd_* keywords are uses) or wing (if wing_*
;						keywords are used) area into an ellipse with a semi-major
;						axis of bkgnd_radius/psfstretch, i.e. larger values of
;						psfstretch shrink the axes.
;						(usually psfstretch=cos(fovangle) where fovangle is the
;						direction cosine angle in the long fov dimension)
;	/get_origin 	if set return origin (array indices for RA=0, dec=0) of
;						standard star map
;	/get_ndim		if set return dimensions of standard star map
;	/get_fov		if set return width of standard map in angular units
;						This is the product of the return values for the /get_ndim
;						and /get_scale keywords.
;	/get_scale		if set return scale (angular units per bin) of standard
;						star map
;	/show			if set display standard star map
;
;	/use_filled_psf if set, use the alternative PSF map
; OUTPUTS:
;	See /get_origin, /get_scale, /get_ndim
;	for return value if one of these keywords is set. If no keyword is set:
;
;	std_map			array[n,n]; type: float
;						2D standard star map of 200x200 bins
;						map >  0 inside PSF (brightness of standard star)
;						map = -1 inside inner radius, but outside PSF
;								(not present if bkgnd_radius is a scalar (outer radius only),
;								or if wing_radius is used)
;						map = -2 inside outer radius, but outside PSF and outside
;								inner radius (if specified);
;								or inside wing_radius but outside PSF.
;						map = -3 outside outer radius, and outside PSF
;	in_psf			array[n,n]; type: byte
;						1 inside the non-zero area of the PSF; 0 outside
;						This is the effective non-zero area of the PSF, or 'core PSF'
;						!!! If keyword wing_radius is set this EXCLUDES the circular
;						!!! area out to radius 'wing_radius', but only covers the part
;						!!! that is non-zero in the standard_star file.
;
;	origin=origin	array[2]; float
;						array indices for RA=0,dec=0) of the standard star map
;	scale=scale		scalar; type: float
;						scale (angular units per bin) of standard star map
; INCLUDE:
	@compile_opt.pro		; On error, return to caller
; COMMON BLOCKS:
	common smei_star_standard_save, $
		std_psf 		, $
		std_wing		, $
		std_cam 		, $
		std_origin		, $
		std_scale		, $
		std_adus		, $
		std_map 		, $
		std_in_psf
; CALLS:
;	who_am_i, destroyvar, IsType, ToDegrees, gridgen, InitVar
; PROCEDURE:
;	The standard star is stored in the Fits file standard_star.fts.gz
;	in the same directory as this procedure.
;	The standard map is a 200x200 array covering a 5x5 deg area
;	(0.025 deg/bin). for each camera. It contains positive-definite
;	values defining the PSF and is zero everywhere else.
;
;	The center of the map (at array indices [99.5,99.5] is the
;	centroid of the PSF. The standard star represents the shape of a star
;	with its centroid located at RA=0, dec=0 (i.e. on the equator).
;
;	Either the bkgnd_* or the wing_* are used to fill the return
;	array:
;
;	Only outer radius specified in bkgnd_radius:
;
;	The input values bkgnd_radius, bkgnd_origin and psfstretch are
;	used to define an ellipsoidal area around the center of the map.
;	Bins outside the ellipse are set to -3.
;	Bins inside the ellipse (but not in the PSF) are set to -2.
;	The output map contains -3, -2, and positive-definite values.
;
;	Both inner and outer radius specified in bkgnd_radius:
;
;	The input values bkgnd_radius, bkgnd_origin and psfstretch are
;	used to define an ellipsoidal ring around the center of the map. 
;	Bins outside the outer ellipse are set to -3.
;	Bins inside the ring between inner and outer ellipse
;	(but not in the PSF) are set to -2.
;	Bins inside the inner ellipse (but not in the PSF) are set to -1
;	The output map contains -3, -2, -1, and positive-definite values
;
;	Single radius or semi-major axes specified in wing_radius:
;
;	The input values wing_radius (and wing_delta), wing_origin, and
;	psfstretch are used to define an ellipsoidal ring around the
;	center of the map.
;	Bins outside the ellipse (but not in the PSF) are set to -3.
;	The output map contains -3, and positive-definite values
;
;	The bins with value -2 in the output map can be used to define
;	a background area around the PSF.
; MODIFICATION HISTORY:
;	AUG-2006, Paul Hick (UCSD/CASS)
;	AUG-2006, Paul Hick (UCSD/CASS)
;		Added /use_filled_psf keyword to work with alternative PSF map
;	MAR-2007, Paul Hick (UCSD/CASS)
;		Change bkgnd_radius and fovxangle from argument to keyword.
;		Added keywords wing_radius and wing_origin
;	MAY-2007, Paul Hick (UCSD/CASS)
;		Changed scalar fovxangle keyword to 2-element array psfstretch
;		Changed keyword /get_adus to /get_psfadus.
;		Added keyword /get_psfarea
;	JUL-2007, Paul Hick (UCSD/CASS; pphick@ucsd.edu)
;		Added keyword in_psf
;		Removed keywords get_psfarea and get_psfadus.
;-

InitVar, use_filled_psf	, /key
InitVar, camera			, 2
InitVar, show			, /key
InitVar, extend_corepsf	, 0
InitVar, silent			, 0

new_psf  = (['fixd','fill'])[use_filled_psf]
new_wing = IsType(wing_radius,/defined)

dpm = ToDegrees(degrees=degrees)

; Make sure to read proper standard map for specified camera.
; If the wrong camera map is in memory we need read another one.

IF IsType(std_psf, /defined) THEN	$
	IF new_psf NE std_psf OR new_wing NE std_wing OR camera NE std_cam THEN $
		destroyvar, std_psf

IF IsType(std_psf, /undefined) THEN BEGIN

	; Read the 200x200 standard star for the appropriate camera
	; Extract the scaling factors for the equatorial mapping from the header

	; The standard star map std_map has non-zero positive values
	; everywhere inside the PSF

	std_psf 	= new_psf
	std_wing	= IsType(wing_radius,/defined)

	std_cam 	= camera
	std_map 	= readfits(filepath(root=who_am_i(/dir),'standard_star_'+std_psf+'.fts.gz'),exten_no=std_cam-1,hdr, /silent)
	std_origin	= [fxpar(hdr,'CX'), fxpar(hdr,'CY')]-1	; 0-based array index
	std_scale	= 1.0/(fxpar(hdr,'BD'))					; degrees per bin (=0.025)

	std_in_psf	= std_map NE 0							; Defines the 'core PSF'

	IF silent LE 0 THEN message, /info, 'standard_star_'+std_psf+'.fts.gz (camera'+strcompress(std_cam)+')'

	IF NOT use_filled_psf AND std_wing THEN BEGIN

		; Fill the entire 200 x 200 array with non-zero values
		; using the specified fall-off per bin

		base		= 0			; 800
		radial		= 1.259		; Fall-off per bin in radial direction (narrow dim of fov)
		tangential	= 1.585		; Fall-off per bin in tangential direction (long dim of fov)

		ndim = size(std_map,/dim)

		old_map = std_map-base*std_in_psf
		new_map = old_map

		right = ndim[0]/2
		left  = right-1
		band  = 12

		high  = ndim[1]/2
		low   = high-1

		FOR i=left-band,right+band DO BEGIN

			FOR j=low,0,-1 DO				$
				IF new_map[i,j] EQ 0 THEN	$
					new_map[i,j] = new_map[i,j+1]/radial

			FOR j=high,ndim[1]-1,1 DO		$
				IF new_map[i,j] EQ 0 THEN	$
					new_map[i,j] = new_map[i,j-1]/radial

		ENDFOR

		FOR j=0,ndim[1]-1 DO BEGIN

			FOR i=left-band-1,0,-1 DO		$
				IF new_map[i,j] EQ 0 THEN	$
					new_map[i,j] = new_map[i+1,j]/tangential

			FOR i=right+band+1,ndim[0]-1,1 DO 	$
				IF new_map[i,j] EQ 0 THEN		$
					new_map[i,j] = new_map[i-1,j]/tangential

		ENDFOR

		;new_map = smooth(new_map,3)
		;new_map[inside_psf] = old_map[inside_psf]

		std_map = new_map

		;FOR n=0,extend_corepsf-1 DO BEGIN
		;	n = where( (shift(std_in_psf,-1,0) EQ 1 OR shift(std_in_psf,1,0) EQ 1 OR shift(std_in_psf,0,-1) EQ 1 OR shift(std_in_psf,0,1) EQ 1 ) AND std_in_psf EQ 0 )
		;	IF n[0] NE -1 THEN std_in_psf[n] = 1B
		;ENDFOR

	ENDIF

ENDIF

map    = std_map
origin = std_origin
scale  = std_scale/dpm
in_psf = std_in_psf

not_in_psf = 1B-in_psf

nomansland_id = -1.0
background_id = -2.0
wayoutside_id = -3.0

InitVar, get_origin , /key
InitVar, get_scale	, /key
InitVar, get_ndim	, /key
InitVar, get_fov	, /key

InitVar, get_nomansland_id, /key
InitVar, get_background_id, /key
InitVar, get_wayoutside_id, /key

IF get_origin  THEN RETURN, origin
IF get_scale   THEN RETURN, scale
IF get_ndim    THEN RETURN, size(map,/dim)
IF get_fov	   THEN RETURN, size(map,/dim)*scale

IF get_nomansland_id THEN RETURN, nomansland_id
IF get_background_id THEN RETURN, background_id
IF get_wayoutside_id THEN RETURN, wayoutside_id

CASE 1 OF

IsType(bkgnd_radius, /defined): BEGIN

	; Around the PSF create an ellipsoidal area of zeros used
	; to pin down the background. Outside the background area
	; the map is set to -1.0

	p = size(map,/dim)						; Dimensions of standard map
	p = gridgen(p, origin=(p-1.0)/2.0)		; Hor/vert indices across standard
											; map relative to center of map
	; The ellipse has semi-major axes of radius/psfstretch[0] in
	; the horizontal (RA) dimension and radius/psfstretch[1] in
	; the vertical (dec; usually psfstretch[1]=1.0) dimension,
	; i.e. the equation for the ellipse is
	;	(x*psfstretch[0]/radius)^2 + (y*psfstretch[1]/radius)^2 = 1
	; (<1 inside, > 1 outside ellipse).

	InitVar, psfstretch, [1.0,1.0]
	p[0,*] *= psfstretch[0]
	p[1,*] *= psfstretch[1]
	
	p *= scale 								; Convert to angular units

	InitVar, bkgnd_origin, [0.0,0.0]		; Shift origin of ellipse
	p -= bkgnd_origin#(replicate(1,n_elements(p)/2))

	CASE n_elements(bkgnd_radius) OF
	1: BEGIN
		p /= bkgnd_radius					; Outer radius only
		p  = total(p*p,1)

		; Select bins in standard map that are NOT part of the PSF and are
		; outside the ellipse

		p_near = -1 						; No inner radius
		p_back = where(p LE 1.0 AND not_in_psf); Inside ellipse
		p_far  = where(p GT 1.0 AND not_in_psf); Outside outer radius
	END
	2: BEGIN
		q  = p/bkgnd_radius[0]				; Inner radius
		p /=   bkgnd_radius[1]				; Outer radius

		p = total(p*p,1)
		q = total(q*q,1)

		; Select bins in standard map that are NOT part of the PSF and are
		; inside the inner ellipse or outside the outer ellipse

		p_near = where(q LT 1.0 AND 			 not_in_psf) 	; Inside inner radius
		p_back = where(1.0 LE q AND p LE 1.0 AND not_in_psf)	; In ring
		p_far  = where(p GT 1.0 AND 			 not_in_psf)	; Outside outer radius
	END
	ENDCASE

END

IsType(wing_radius,/defined): BEGIN

	; Extend the PSF by an ellipoidal area defined by an exponential
	; dropoff of the PSF intensity.
	
	p = size(map,/dim)						; Dimensions of standard map
	p = gridgen(p, origin=(p-1.0)/2.0)		; Hor/vert indices across standard
											; map relative to center of map
	; The ellipse has semi-major axes of radius/psfstretch[0] in
	; the horizontal (RA) dimension and radius/psfstretch[1]
	; in the vertical (dec; usually psfstretch[1]=1.0) dimension,
	; i.e. the equation for the ellipse is
	;	(x*psfstretch[0]/radius)^2 + (y*psfstretch[1]/radius)^2 = 1
	; (<1 inside, > 1 outside ellipse).

	InitVar, psfstretch, [1.0,1.0]
	p[0,*] *= psfstretch[0]
	p[1,*] *= psfstretch[1]

	p *= scale 								; Convert to angular units

	InitVar, wing_origin, [0.0,0.0]			; Shift origin of ellipse
	p -= wing_origin#(replicate(1,n_elements(p)/2))

	InitVar, wing_delta		, 0.0
	InitVar, min_wing_radius, 0.5/dpm
	p[0,*] /= (wing_radius[0]+wing_delta[0]) > min_wing_radius
	p[1,*] /= (wing_radius[n_elements(wing_radius)-1]+wing_delta[n_elements(wing_delta)-1]) > min_wing_radius

	p = total(p*p,1)

	; Select bins in standard map that are NOT part of the PSF
	; and are outside the ellipse. Set the standard map to
	; zero outside the ellipse.

	q = where(p GT 1.0 AND not_in_psf) 
	IF q[0] NE -1 THEN map[q] = 0.0

	p_near = -1 							; No inner radius
	p_back = -1 							; No outer radius
	p_far  = where(map LE 0.0)				; Outside outer radius

END

ELSE: BEGIN 								; Use the unmodified PSF

	p_near = -1 							; No inner radius
	p_back = -1 							; No outer radius
	p_far  = where(map LE 0.0)				; Outside PSF

END

ENDCASE

IF p_near[0] NE -1 THEN map[p_near] = nomansland_id
IF p_back[0] NE -1 THEN map[p_back] = background_id
IF p_far [0] NE -1 THEN map[p_far ] = wayoutside_id

IF show THEN	$
	view, in=map+(1-map)*(map GT 0)+in_psf, /stretch

RETURN, map  &  END