;+
; NAME:
;	EarthTransit3DLoc
; PURPOSE:
;	Provides 3D heliographic coordinates for all lines of sight used for a 'transit sweep' across the sky
; CATEGORY:
;	WWW: skymap
; CALLING SEQUENCE:
	FUNCTION EarthTransit3DLoc, UT, $
		cv2carrington=cv2carrington, $
		nra 		= nra		, $
		nde 		= nde		, $
		nrr 		= nrr		, $
		drr 		= drr		, $
		equator 	= equator	, $
		degrees 	= degrees	, $
		band		= band		, $
		solar		= solar 	, $
		geolng		= geolng	, $
		ra			= ra		, $
		dec 		= dec		, $
		zero_phase	= zero_phase, $
		zero_point	= zero_point, $
		rr_earth	= rr_earth	, $
		ut_earth	= ut_earth	, $
		pa_earth 	= pa_earth 	, $
		elo_earth	= elo_earth	, $
		elo_sun 	= elo_sun	, $
		_extra		= _extra
; INPUTS:
;	UT		scalar; type: time structure
;				universal time, UT
;				(usually UT corresponding to local noon at geographic longitude 'geolng')
; OPTIONAL INPUT PARAMETERS:
;	/degrees	if set then all in- and ouput angles are in degrees (default: radians)
;
;	Four parameters define the grid used to calculate the sky map (see PROCEDURE)
;
;	nra		scalar; type: int; default: 36
;				# points in skymap in longitudinal direction (either right ascensions
;				of ecliptic longitude). Note that nRA may be modified if keyword
;				band is used.
;	nde		scalar; type: int; default: 18
;				# points in skymap in latitudinal direction(either declination or
;				ecliptic latitude)
;	nrr		scalar; type: int; default: 20
;				# 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)
;
;	geolng=geolng
;			scalar; type: any; default: 0
;				geographic longitude of observing location
;	/equator
;			by default the local meridian is set up as a strip of sky centered on the
;				intersection with the ecliptic. This results in a sweep with the Sun
;				centered, and the ecliptic along the horizontal axis (but see RESTRICTIONS).
;				If /equator is set then the strip of sky is centered on the equator, i.e.
;				the strip covers the local meridian from equatorial north to equatorial south.\
;				This is a more strict interpretation of the concept of a 'transit sweep', but
;				has the disadvantage that the Sun is not at the center of the sky map.
;	band	scalar; type: float
;				Width of strip to be processed (in hours)
;				If band is not set then meridian strips covering 24 hours centered
;				around UT are selected. Use this keyword to extract only a partial strip
;				of sky covering less than 24 hour.
;	/solar	by default the meridian strips are centered around the meridian strip at time UT.
;				Setting /solar rearranges the strips, putting the meridian strip at noon in
;				the center (see PROCEDURE)
;	/cv2carrington
; OUTPUTS:
;	R		array[3,nra,nde,nrr]; type: float
;				3D locations (in heliographic coordinates) where values need to obtained
;				by interpolation on the 3D arrays.
;	nra		scalar; type: integer
;				if keyword band is set then nra is set to the number of meridians
;				inside the RA range of width 'band'.
; OPTIONAL OUTPUT PARAMETERS:
;	ra		array[nra]; type: float
;				Right ascension for nra UT times covering 24 hours centered on UT
;	dec		array[nde]; type: float
;				dec = -pi/2+pi/nde*[0.5,1.5,...,nde-0.5]
;				ra and dec define a regular grid in sky coordinates for the centers of
;				the boxes for the lines of sight used to sweep of the sky
;	zero_point=zero_point
;			scalar; type: float
;				/solar NOT set: RA angle for local meridian at 'geoloc' on time UT
;				/solar set    : RA angle of Sun at time UT
;	zero_phase=zero_phase
;			/solar NOT set:
;				scalar, type: float
;					RA angle for local meridian at 'geoloc' on time UT
;					(i.e. same as zero_point)
;			/solar set    :
;				array[nra]; type: float
;					RA angles of Sun for nRA UT times covering 24 hours centered
;					on UT (same times as used for output array RA).
;					(i.e. 'zero_point' is approximately centered in 'zero_phase')
;	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)
;	rr_earth=rr_earth
;			array[3,nra]; type: float
;				heliocentric locations Earth in heliographic coordinates
; INCLUDE:
	@compile_opt.pro				; On error, return to caller
; CALLS:
;	InitVar, IsType, ToRadians, gridgen, TimeOp, TimeSet, TimeGST, SuperArray
;	CvSky, big_eph, CvPointOnLos, AngleRange, TimeUnit, Carrington
; RESTRICTIONS:
;	If /equator is NOT set then the ecliptic is put on the horizontal axis. This is
;	done using a kludge. Each meridian strip (-90 to +90 degrees in declination with
;	the equator in the center, is shifted vertically to put the ecliptic in the center.
;
;	The sweep across the sky covers 24 hours in UT. Since UT is a solar time, the sweep covers
;	a little more than 360 degree of sky. However, the RA array still covers exactly 360
;	degrees (RA is passed to PlotEarthSkyMap to make a Hammer-Aitoff or fish-eye map).
; PROCEDURE:
; >	UT is used to set up an array of nRA equally spaced times covering UT-12h to UT+12h
;		TT = UT+((0.5+findgen(nra))/nra-0.5)*(24)
;	(if 'band' is specified then only UT +/- band/2 is used).
; >	The geographic longitude, together with GST, is used to calculate the RA
;	of the local meridian (hour angle zero) at times TT
; >	At each ra(TT) value nde equally spaced declinations are taken
;		dec = ((0.5+findgen(nde))/nde-0.5)*90
; >	The resulting strips along the local meridian patched together give
;	a (partial) sweep of the sky. The return array RA is a monotonic increasing
;	array inside [-180,+180].
; >	If /solar is NOT set then the return array zero_phase contains the right ascension
;	of the local meridian at time UT.
;	If /solar is set then the return array zero_phase contains the right ascension of
;	the sun at the times centered on UT used also to calculate RA.
;	The scalar 'zero_point'  and array 'zero_phase' are used to center the proper
;	RA in a skymap by subtracting it from RA.
; >	The output array rr_earth gives the location of Earth at times TT
; MODIFICATION HISTORY:
;	SEP-1999, Paul Hick (UCSD/CASS)
;	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, solar		, /key
InitVar, equator	, /key
InitVar, cv2carrington, /key

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

body_sun   = jpl_body(/sun  ,/string)
body_earth = jpl_body(/earth,/string)

InitVar, geolng, 0.0d0, set=geoloc

InitVar, nra, 360/10						; # points in RA direction
InitVar, nde, 180/10						; # points in decl direction

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

nra = long(nra)
nde = long(nde)
nrr = long(nrr)

ut_earth = gridgen(nra, /open, range=[-12,12]); 24 hours centered on UT

dra = 2*pi/nra								; Bin width for right ascension grid

IF IsType(band, /defined) THEN BEGIN		; Select less than 24 hour
	hr = where(-band/2 LT ut_earth AND ut_earth LE band/2, i)
	IF i GE 1 THEN BEGIN					; Need two columns to determine width in RA
		ut_earth = ut_earth[hr]
		nRA = i
	ENDIF
ENDIF
											; nRA universal times centered on UT
ut_earth = TimeOp(/add,UT,TimeSet(/diff,ut_earth,TimeUnit(/hour)))

; LST = GST+geoloc. At time LST the local meridian (hour angle=0) has RA = LST

ra  = TimeGST(ut_earth, degrees=degrees)+geoloc	; RA for local meridians at ut_earth
dec = gridgen(nde, /open, range=pi2*[-1,1])	; Declinations (radians) covering [-pi/2,+pi/2]
rr  = drr*gridgen(nrr, /open)				; Distances along line of sight (AU)

; Set up the right ascension array. Information about which ra is plotted in the
; center of the skymap is stored in zero_point and zero_phase.

CASE solar OF

0: BEGIN

	zero_point = (TimeGST(UT, degrees=degrees)+geoloc)[0] ; ra meridian at UT
	zero_phase = zero_point					; ra relative to meridian at UT

END

1: BEGIN

	; The skysweep is displayed with the Sun in the center. The ra of the Sun is
	; stored in 'zero_phase'. This will be used to plot the noon meridian in the center.

	; Position Sun at time UT

	zero_point = big_eph(UT	, $
		body   = body_sun	, $
		center = body_earth , $
		degrees= degrees	, $
		/to_equatorial		, $
		/to_sphere			, $
		/precess			, $
		/onebody			, $
		/silent 			)

	zero_point = zero_point[0]

	; Position Sun at times ut_earth

	zero_phase = big_eph(ut_earth, $
		body   = body_sun	, $
		center = body_earth , $
		degrees= degrees	, $
		/to_equatorial		, $
		/to_sphere			, $
		/precess			, $
		/onebody			, $
		/silent 			)

	zero_phase = reform(zero_phase[0,*])	; RA angles of Sun at times ut_earth

END

ENDCASE

CASE equator OF

0: BEGIN

	; If /equator is NOT set then the ecliptic should be put on the horizontal.
	; Currently this is done using a kludge. Each meridian strip ([-90,+90] degrees declination
	; with the equator in the center) is shifted vertically to put the ecliptic in the center.

	B  = 23.439291d0/180d0*!dpi				; Angle equator/ecliptic
	AB = ra*rpm

	; RA is the right ascension of the local meridian at time ut_earth
	; A is the intersection of the local meridian with the equator.
	; B is the vernal equinox (intersection ecliptic and equator)
	; C is the intersection of the local meridian with the ecliptic.
	; We need the declination of C (angle AC).
	; In spherical triangle ABC: angle A=90 degrees, B = 23.4 degrees. AB = RA
	; First calculate BC (ecliptic longitude of C) by applying analogue and cosine formulae:
	;	sin(BC) = cos(AC)*sin(AB)/cos(B)
	;	cos(BC) = cos(AC)*cos(AB)
	; Since cos(AC) > 0, we can drop it in the atan call.

	tmp = sin( atan( sin(AB)/cos(B), cos(AB) ) ); Analogue and cosine formulae

	; Then calculate declination AC using sine rule and analogue formula:
	;	sin(AC) = sin(BC)*sin(B)
	;	cos(AC) = sin(BC)*cos(B)/sin(AB)

	tmp = atan( tmp*sin(B), tmp*cos(B)/sin(AB) ) /rpm

	; Set up declinations for a 180 degree meridian centered on C
	; Combine every right ascension with every declination

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

	R[0,*] = SuperArray(ra , nde, /trail)
	R[1,*] = SuperArray(tmp, nde, /trail) + SuperArray(dec, nra, /lead)

	; R[1,*] now contains some declinations larger than 90 degrees or smaller then -90 degrees.
	; For these points 180 degrees is added to the righ ascension while the declination
	; is put back into the range [-90,+90] degrees

	i = where(R[1,*] GT pi2)				; Declination > 90
	IF i[0] NE -1 THEN BEGIN
		R[0,i] += pi
		R[1,i]  = pi-R[1,i]
	ENDIF

	i = where(R[1,*] LT -pi2)				; Declination < -90
	IF i[0] NE -1 THEN BEGIN
		R[0,i] +=  pi
		R[1,i]  = -pi-R[1,i]
	ENDIF

END

1: BEGIN									; Combine every right ascension

	R = fltarr(3,nra*nde, /nozero)			; .. with every declination

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

END

ENDCASE

; Set up array with a time for all nrr line-of-sight segments on all nra*nde lines of sight

TTall = replicate( !TheTime, nra, nde*nrr )
FOR i=0L,nde*nrr-1 DO TTall[*,i] = ut_earth
TTall = reform(TTall, nra*nde*nrr, /overwrite)

; Convert geocentric equatorial to geocentric ecliptic coordinates

R = CvSky(TTall[0:nra*nde-1], from_equatorial=R, /to_ecliptic, degrees=degrees, /silent)

; Replicate R for all nrr line-of-sight segments

R = SuperArray(R, nrr, /trail)
R = reform(R, 3, nra*nde*nrr, /overwrite)
R[2,*] = SuperArray(RR, nra*nde, /lead)		; Fill in line-of-sight distances

; Change from geocentric to heliocentric perspective
; !!! elold is angle Sun-Observer-Point on los (i.e. the elongation)
;		elnew is angle Observer-Sun-Point on los
; These angles are need 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.

SunEcl = big_eph(ut_earth, $
		body	= body_sun	, $
		center	= body_earth, $
		degrees = degrees	, $
		/to_ecliptic		, $
		/to_sphere			, $
		/precess			, $
		/onebody			, $
		/silent 			)

SunEcl = SuperArray(SunEcl,nde*nrr,/trail)
SunEcl = reform(SunEcl,3,nra*nde*nrr,/overwrite)

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

; Convert from geo- to heliocentric perspective

R = CvPointOnLos(R,oldphase=pa_earth,oldelo=elo_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(TTall, 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(TTall, near_longitude=R[0,*], degrees=degrees)
ENDIF

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)

; 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:nra-1]
rr_earth[0,*] += pi
rr_earth[1,*] *= -1

rr_earth = CvSky(ut_earth, from_ecliptic=rr_earth, _extra=_extra, degrees=degrees, /silent)

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

RETURN, R  &  END