;+
; NAME:
;	EarthSky3DLoc
; PURPOSE:
;	Returns the 3D positions for all segments along a collection of lines
;	of sight. The lines of sight either form a grid for a 2D sky map, or
;	are explicitly input.
; CATEGORY:
;	sat/idl/util
; CALLING SEQUENCE:
	FUNCTION EarthSky3DLoc, UT	, $
		cv2carrington=cv2carrington, $
		nra 		= nra		, $
		nde 		= nde		, $
		nrr 		= nrr		, $
		drr 		= drr		, $
		equator 	= equator	, $
		degrees 	= degrees	, $
		ra			= ra		, $
		dec 		= dec		, $
		zero_point	= zero_point, $
		zero_phase	= zero_phase, $
		zero_shift	= zero_shift, $
		rr_earth	= rr_earth	, $
		ut_earth	= tt		, $
		pa_earth 	= pa_earth 	, $
		elo_earth	= elo_earth , $
		elo_sun 	= elo_sun 	, $
		center		= center	, $
		_extra		= _extra
; INPUTS:
;	UT			scalar; type: standard time structure
;				time when skymap is required
; OPTIONAL INPUT PARAMETERS:
;	/degrees	if set then all in- and ouput angles are in degrees
;				(default: radians)
;
;	/equator	by default the lines of sight are specified in ecliptic coordinates
;				(i.e. in a skymap the horizontal plane is the ecliptic).
;				If /equator is set then equatorial coordinates are assumed (i.e.
;				in a skymap the horizontal plane is Earth's equator)

;	/cv2carrington
;				NOTE: this only makes sense if keyword /to_heliographic is set !!
;
;				By default, the longitudinal (first) dimension in the output
;				R and rr_earth is a phase angle (longitude).
;				If heliographic coordinates are used then instead of an angle
;				a Carrington variable within 0.5 of the Carrington variable
;				corresponding to input time UT.
;
;	zero_shift=zero_shift
;			scalar; type: float; default: 0.0
;				shifts the center of the sky grid. Instead of centering on the
;				Sun the grid is shifted to increasing RA or ecliptic longitude
;				by 'zero_shift'.
;				The units must be consistent with the setting of /degrees.
;
;	nrr		scalar; type: int; default: 20
;				number of steps along lines of sight at each point in the sky
;	drr		scalar; type: float; default: 0.1
;				step size along line of sight (AU)
;
;	Lines of sight are specified as sidereal location in the sky (either
;	RA/dec (/equator SET) or ecliptic longitude/latitude (/equator NOT set).
;
;	If nRA and nDE are set to non-zero values then a grid of lines of sight
;	covering the entire sky is created (see PROCEDURE): 	
;
;	nra		scalar; type: int; default: 72
;				nr of points in skymap in longitudinal direction (either right
;				ascensions of ecliptic longitude)
;				If nra or nde is set to zero then ra and dec are used as input.
;	nde		scalar; type: int; default: 36
;				nr of points in skymap in latitudinal direction (either declination
;				or ecliptic latitude)
;
;	Alternatively, set the product nra*nde=0 and specify lines of sight explicitly
;	through keywords ra,dec:
;
;	ra[nLOS]
;			scalar or array; right ascensions (/equator SET) or
;				ecliptic longitudes (/equator NOT set)
;	dec[nLOS]
;			scalar or array; declinations (/equator SET) or
;				ecliptic latitudes (/equator NOT set).
;
; OUTPUTS:
;	Result	array[3,nra,nde,nrr]; type: float
;			array[3,nlos,nrr]; type: float (if nra or nde is zero)
;				3D locations for all segments along all lines of sight. This is where
;				values need to obtained by interpolation on the 3D arrays.
;				The result is in spherical coordinates in the coordinates specified
;				as one of the /to_* keywords to href=CvSky=.
;
;				R[0,*,*,*]: longitudes
;							If /to_heliographic and /cv2carrington are SET then
;							these are Carrington variables within 0.5 of Carrington(UT)
;				R[1,*,*,*]; latitudes
;				R[2,*,*,*]: heliocentric distances (AU)
; OPTIONAL OUTPUT PARAMETERS:
;	Only if nra*nde = 0 on input:
;
;	ra		array[nra]; type: float
;				right ascensions or ecliptic longitudes across sky
;				(with ra/longitude of Sun in center of array)
;				i.e. ra = (ra Sun at time UT) -pi+2*pi/nra*[0.5,1.5,...,nra-0.5]
;	dec		array[nDE]; type: float
;				declination or ecliptic latitude across sky
;				(with decl/lat zero in center of array)
;				DE = -pi/2+pi/nDE*[0.5,1.5,...,nDE-0.5]
;
;	Depending on setting of /equator, RA and DEC define a regular grid in
;	equatorial or ecliptic sky coordinates for the centers of the boxes for the
;	lines of sight used to build the skymap with the position of the Sun
;	centered in the RA/longitudinal direction.
;
;	zero_phase=zero_phase
;			scalar; type: float
;				right ascension or ecliptic longitude of Sun (plus zero_shift)
;				if specified) at time UT, i.e. this is the value on which ra is centered
;	zero_point=zero_point
;			scalar; type: float
;				right ascension or ecliptic longitude of Sun (plus zero_shift)
;				at time UT (i.e. same as 'zero_phase'; should stay to ensure
;				that the values are consistent when fed e.g. to PlotEarthSkymap)
;	rr_earth
;			array[3]; type: float
;				heliocentric location Earth in heliographic coordinates
;				rr_earth[0]: longitude
;							If /to_heliographic and /cv2carrington are SET then
;							these is a Carrington variable within 0.5 of Carrington(UT)
;				rr_earth[1]; latitude
;				rr_earth[2]: heliocentric distance (AU)
;	pa_earth=pa_earth
;			array[nra,nde,nrr]; type: float
;				line of sight position angle measured counterclockwise
;				from ecliptic north
;	elo_earth=elo_earth
;			array[nra,nde,nrr]; type: float
;				line of sight elongations: angle (Sun-Earth-los segment)
;	elo_sun=elo_sun
;			array[nra,nde,nrr]; type: float
;				angle Earth-Sun-los segment
; INCLUDE:
	@compile_opt.pro							; On error, return to caller
; SEE ALSO:
;	EarthTransit3DLoc, PlotEarthSkymap
; CALLS:
;	InitVar, SuperArray, CvPointOnLos, big_eph, jpl_body, CvSky, gridgen
;	ToRadians, AngleRange, Carrington
; PROCEDURE:
; >	The sky map is build from nRA*nDE line of sight integrations through the
;	F3D array. The locations in the sky of the lines of sight are
;		ra/Lng    = (ra/Lng Sun) -180 + (360./nra)*[0.5, 1.5, .. , nra-0.5 ]
;		Decl./Lat = - 90 + (180./nde)*[0.5, 1.5, .. , nde-0.5 ]
;		Dist along los = drr*[0.5,1.5,...,nrr-1]
;	Note that each line of sight is centered on a 'box' of 360/nra by 180/nde
;	degrees, and that all lines of sight together cover the entire sky.
; MODIFICATION HISTORY:
;	AUG-1999, Paul Hick (UCSD/CASS)
;	OCT-2004, Paul Hick (UCSD/CASS)
;		Added zero_shift keyword.
;	APR-2006, Paul Hick (UCSD/CASS)
;		Added keyword /longitudes (the previous version always returned
;		heliographic longitudes). Add option to explictly specify lines of sight
;		in RA,DEC after setting nRA*nDE=0.
;	FEB-2007, Paul Hick (UCSD/CASS; pphick@ucsd.edu)
;		Replaced keyword /longitudes by keyword /cv2carrington
;		Added _extra keyword to be able to pass one of the /to_* keywords
;		to CvSky. This allows the return value to be in any one of the
;		coordinates supported by CvSky (instead of just heliographic coordinates)
;-
InitVar, equator	  , /key
InitVar, cv2carrington, /key

rpm = ToRadians(degrees=degrees)
pi  = !dpi/rpm
pi2 = pi/2

InitVar, zero_shift, 0.0

; Geocentric ecliptic location of Sun in ecliptic and equatorial coordinates

InitVar, center, jpl_body(/earth,/string)

SunEq = big_eph(UT		, $
	center	= jpl_body(/sun,/string), $
	body	= center	, $
	/to_equatorial		, $
	/to_sphere			, $
	/precess			, $
	degrees = degrees	, $
	/onebody			, $
	/silent 			)

SunEq[0] = AngleRange(SunEq[0]+pi,degrees=degrees)
SunEq[1] = -SunEq[1]

;SunEq = big_eph(UT		, $
;	body    = jpl_body(/sun,/string), $
;	center  = center	, $
;	/to_equatorial		, $
;	/to_sphere			, $
;	/precess			, $
;	degrees = degrees	, $
;	/onebody			, $
;	/silent 			)

SunEcl = CvSky(UT,from_equatorial=SunEq,/to_ecliptic,degrees=degrees,/silent)

InitVar, nra, 360/5 	; # points in ra direction
InitVar, nde, 180/5 	; # points in decl direction

InitVar, nrr, 20		; # steps along line of sight
InitVar, drr, 0.1		; Step size along line of sight (AU)

skygrid = nra*nde NE 0

zero_point = AngleRange(([SunEcl[0],SunEq[0]])[equator]+zero_shift,/degrees)
zero_phase = zero_point	; Longitude/ra Sun

CASE skygrid OF

0: BEGIN
	SyncArgs, ra, dec
	nlos = n_elements(ra)

	;ra = zero_phase+ra							; Center on Sun
	rr = drr*gridgen(nrr, /open)				; Distances along line of sight (AU)

	; Set up collection of lines of sight by combining every ra with every declination
	; and with every distance along line of sight.

	R = fltarr(3,nlos*nrr, /nozero)

	R[0,*] = SuperArray(ra ,nrr ,/trail)
	R[1,*] = SuperArray(dec,nrr ,/trail)
	R[2,*] = SuperArray(rr ,nlos,/lead )
END

1: BEGIN
	nra = long(nra)
	nde = long(nde)
	nrr = long(nrr)

	; ra covers 360 deg in right ascension or ecliptic longitude relative to the
	; Sun, going from 180 deg west of the Sun to 180 deg east of the Sun
	; dec covers declinations or latitudes going from -90 to +90
	; RR covers geocentric distances from 0 to nrr*drr

	ra = gridgen(nra, /open, range=[-pi,pi])	; Right ascensions/ecliptic longitudes [-180,+180]

	; ra are right ascensions or ecliptic longitudes

	ra += zero_phase							; Center on Sun
	dec = gridgen(nde, /open, range=[-pi2,pi2])	; Latitudes/declinations
	rr  = drr*gridgen(nrr, /open)				; Distances along line of sight (AU)

	; Set up collection of lines of sight by combining every RA with every declination
	; and with every distance along line of sight.

	R = fltarr(3,nra*nde*nrr, /nozero)

	R[0,*] = SuperArray(ra,nde*nrr,/trail)
	R[1,*] = SuperArray(SuperArray(dec,nra,/lead),nrr,/trail)
	R[2,*] = SuperArray(rr,nra*nde,/lead)
END

ENDCASE

; Convert ra and declinations to geocentric ecliptic coordinates

IF equator THEN R = CvSky(UT,from_equatorial=R,/to_ecliptic,degrees=degrees,/silent)

; Change from geocentric to heliocentric perspective
; !!! oldelo is angle Sun-Observer-Point on los (i.e. the elongation)
;	newelo is angle Observer-Sun-Point on los
; These angles are needed to compute the angle between the radial direction
;	and the direction perpendicular to the line of sight. This angle in turn
;	is needed in the integration for the IPS velocity.

R[0,*] -= SunEcl[0]							; Ecliptic longitude relative to Sun
R[2,*] /= SunEcl[2]							; Normalize distance to Sun-Earth distance

; Convert from geocentric to heliocentric perspective

R = CvPointOnLos(R,oldelo=elo_earth,oldphase=pa_earth,newelo=elo_sun,degrees=degrees)

R[2,*] *= SunEcl[2]							; Convert distance back to AU
R[0,*] += SunEcl[0]+pi						; Heliocentric ecliptic longitude
											; Convert from heliocentric ecliptic to heliographic coordinates
R = CvSky(UT,from_ecliptic=R,_extra=_extra,degrees=degrees,/silent,euler_info=euler_info)

IF cv2carrington THEN BEGIN
	to_heliographic = strpos(euler_info,'to heliographic') NE -1
	IF NOT to_heliographic THEN message, 'invalid use of "/cv2carrington" keyword: no heliographic coordinates'
	R[0,*] = Carrington(UT, near_longitude=R[0,*], degrees=degrees)
ENDIF

CASE skygrid OF
0: BEGIN
	R = reform(R,3,nlos,nrr, /overwrite)
	elo_earth = reform(elo_earth, nlos, nrr, /overwrite)
	elo_sun   = reform(elo_sun  , nlos, nrr, /overwrite)
	pa_earth  = reform(pa_earth , nlos, nrr, /overwrite)
END
1: BEGIN
	R = reform(R,3,nra,nde,nrr, /overwrite)
	elo_earth = reform(elo_earth, nra,nde,nrr, /overwrite)
	elo_sun   = reform(elo_sun  , nra,nde,nrr, /overwrite)
	pa_earth  = reform(pa_earth , nra,nde,nrr, /overwrite)
END
ENDCASE

; R now contains all 3D locations (in heliographic coordinates) where values
; need to obtained by interpolation on the input 3D arrays.

; Location Earth in heliographic coordinates

rr_earth = [SunEcl[0]+pi,-SunEcl[1],SunEcl[2]]
rr_earth = CvSky(UT,from_ecliptic=rr_earth,_extra=_extra,degrees=degrees,/silent)

IF cv2carrington THEN	$
	rr_earth[0] = Carrington(UT, near_longitude=rr_earth[0], degrees=degrees)

TT = UT

RETURN, R  &  END