;+
; NAME:
;	smei_sky_atlocation
; PURPOSE:
;	Good question
; CATEGORY:
;	Extract function values at specified sky locations from a SMEI
;	skymaps. Put together in kinda adhoc fashion. So sorry.
; CALLING SEQUENCE:
	FUNCTION smei_sky_atlocation, wanted_map, $
		show		= show		, $
		time		= time		, $
		_extra		= _extra
; INPUTS:
;	wanted_map		passed to href=smei_getfile=
; OPTIONAL INPUT PARAMETERS:
;	/show			prints the locations and function values extracted
;					from the skymap on the screen for all extensions
;					in the selected skymap.
;	_extra=_extra	/show SET:
;					One of the from_* keywords to href=cvsmei= is used here to
;					specify one or more skymap locations
;					(just a single location should be specified, I think, otherwise
;					the screen will be flooded with too much output).
;					/show NOT set:
;					keywords are passed to href=smei_sky_field= to determine
;					whether a skymap is requested, and what extension to read,
;					if it is a skymap.
;					Also, one of the from_* keywords to href=cvsmei= is used
;					to specify the sky locations
; OUTPUTS:
;	R				function value 
; OPTIONAL OUTPUT PARAMETERS:
;	time=time		array[1]; type: time structure
;					if /show is set, then the time at which SMEI observed the
;					sky location is returne here.
; INCLUDE:
	@compile_opt.pro		; On error, return to caller
; CALLS:
;	InitVar, smei_getfile, BadValue, smei_sky_field, cvsmei
;	TimeSet, TimeGet, TimeOp
; PROCEDURE:
; MODIFICATION HISTORY:
;	APR-2008, Paul Hick (UCSD/CASS; pphick@ucsd.edu)
;-
InitVar, show, /key

given_map = smei_getfile( $
	wanted_map			, $
	count	= count 	, $
	_extra	= _extra	)

f = BadValue(0.0)

IF count EQ 0 THEN	$
	RETURN, f


CASE show OF

0: BEGIN

	CASE smei_sky_field(_extra=_extra, /frame_map) OF

	0: BEGIN

		exten = smei_sky_field(given_map,_extra=_extra)

		CASE exten NE -1 OF

		0: message, /info, 'requested field not present'

		1: BEGIN

			map = readfits( given_map, hdr, exten=exten, /silent )

			p = cvsmei( $
				_extra=_extra, $; Provides from* keywords
				/to_map 	 , $
				to_mode = smei_sky_hdr2mode(hdr), $
				/silent)

			p = round(p)		; Nearest skybin

			x = reform(p[0,*])
			y = reform(p[1,*])

			n = size(map,/dim)

			i = where( 0 LE x AND x LE n[0]-1 AND	$
					   0 LE y AND y LE n[1]-1 )

			f = replicate(f, n_elements(x))
			IF i[0] NE -1 THEN f[i] = map[x[i],y[i]]

		END
		ENDCASE

	END

	1: message, /info, 'requested extension is not a skymap'

	ENDCASE

END

1: BEGIN

	; Read first extension to determine maptype.
	; If the "maptype" string contains "RA" OR "DEC" then it is a skymap

	exten = 0
	map = readfits( given_map, hdr, exten=exten, /silent )

	maptype = fxpar(hdr,'MAP')
	j = strpos(maptype,'RA' )
	k = strpos(maptype,'DEC')
	skymap = j NE -1 OR k NE -1

	WHILE skymap DO BEGIN		; Process extensions until one is encountered
					; that is not a skymap
		IF j EQ -1 THEN j = k
		maptype = strmid(maptype,0,j-1)

		; Take coordinates specified in from* keyword and convert
		; to array indices into the current skymap extension

		p = cvsmei( $
			_extra=_extra, $; Provides from* keywords
			/to_map 	 , $
			to_mode = smei_sky_field(given_map,exten=exten,/cvsmei), $
			/silent)

		p = round(p)		; Nearest skybin

		x = reform(p[0,*])	; Horizontal array index
		y = reform(p[1,*])	; Vertical array index

		n = size(map,/dim)

		i = where( 0 LE x AND x LE n[0]-1 AND	$
				   0 LE y AND y LE n[1]-1 )

		f = replicate( BadValue(0.0), n_elements(x) )

		; For all locations inside the map (after rounding to the nearest bin
		; extract the skymap value

		IF i[0] NE -1 THEN f[i] = map[x[i],y[i]]

		; Print information extracted from the skymap.
		; NOTE: I think this only makes sense if a single map location is
		;       extracted.

		CASE strpos(maptype,'Time (sec since TORIGIN)') NE -1 OF
		0: print, maptype, x, y, f, format='(A35,2I6,F10.3)'
		1: BEGIN
			t = TimeSet(fxpar(hdr,'TORIGIN'))
			t = TimeOp(/add,t,TimeSet(/diff,f,TimeUnit(/sec)))
			time = t			; Output keyword
			t = TimeGet(/ydoy,t,upto=TimeUnit(/sec),/roundt)
			print, maptype, x, y, f, t, format='(A35,2I6,F10.3,1X,A)'
		END
		ENDCASE

		exten++ 				; Move to next extension
		map = readfits( given_map, hdr, exten=exten, /silent )

		maptype = fxpar(hdr,'MAP')
		j = strpos(maptype,'RA' )
		k = strpos(maptype,'DEC')
		skymap = j NE -1 OR k NE -1

	ENDWHILE

END

ENDCASE

RETURN, f  &  END