;+
; NAME:
;	CvPointOnLos
; PURPOSE:
;	Converts from topocentric to heliocentric coordinates and v.v
; CALLING SEQUENCE:
	FUNCTION CvPointOnLos, RP,$
		oldelo	= ELOLD 	, $
		oldphase= PHOLD 	, $
		newelo	= ELNEW 	, $
		eorw	= EorW		, $
		degrees = degrees
; INPUTS: (either topocentric or heliocentric)
;	RP					array[3,*]; type: float
;							RP[0,*]	= longitude (0<=LNG<=360) and ...
;							RP[1,*] = latitude (-90<=LAT<=90) of line of sight
;							RP[2,*] = radial distance to point on los
;								(in units of the observer-Sun distance)
; OPTIONAL INPUT PARAMETERS:
;	/degrees			if set all in- and output angles are in degrees
;						(default is radians)
; OUTPUTS: (either heliocentric or topocentric)
;	Result				array[3,*]; type: float
;							RP[0,*]	= longitude (0<=XLNG<=360) and ...
;							RP[1,*] = latitude (-90<=XLAT<=90) of line of sight
;							RP[2,*] = radial distance to point on los
;								(in units of the observer-Sun distance)
;	oldelo=OldElo		array[*]; type: float
;							elongation in the old coordinate system; 0<=ELOLD<=180
;	oldphase=OldPhase	array[*]; type: float
;							phase angle measured counterclockwise from the north
;	newelo=NewElo		array[*]; type: float
;							elongation in the new coordinate system 0<=ELNEW<=180
;	EorW				array[*]; type: integer
;							+1 or -1; indicates on which side (East or West) of
;							the Sun-Earth line the point is located
; INCLUDE:
	@compile_opt.pro		; On error, return to caller
; CALLS:
;	ToRadians, SyncDims
; RESTRICTIONS:
; >	There is no check for -90 <= XLAT <= 90
; PROCEDURE:
; >	The topocentric coordinate system (with Earth in the origin)
;	has its X-axis pointing towards the Sun. Y and Z-axis are arbitrary
;	(usually the X-Y plane will be the ecliptic).
;	The topocentric longitude is measured in a positive sense, i.e.
;	counterclockwise as viewed from the positive Z-axis.
; >	The heliocentric coordinate system (with the Sun in the origin
;	has its X-axis pointing towards Earth.
;	The heliocentric longitude is measured in a positive sense.
; >	Spherical coordinates for both systems: longitude (deg), latitude
;	(deg) and radial distance (units of Sun-Earth distance).
; >	Input can be in topocentric or heliocentric coordinates. The output
;	will be in the other coordinate system (the calculation is symmetric).
; >	EorW is determined from the input elongation XLNG:
;	if 0<=XLNG<=180 then EorW = 1; if 180<XLNG<360 the EorW = -1
; >	If input is in topocentric coordinates then
;	- OldElo is the angle Sun-Earth-P
;	- NewElo is the angle Earth-Sun-P
;	- EorW = +1 if P is towards the east of the Sun (viewed from Earth)
;	- EorW = -1 if P is towards the west of the Sun (viewed from Earth)
; >	If input is in heliocentric coordinates then
;	- OldElo is the angle Earth-Sun-P
;	- NewElo is the angle Sun-Earth-P
;	- EorW = -1 if P is towards the east of the Sun (viewed from Earth)
;	- EorW = +1 if P is towards the west of the Sun (viewed from Earth)
; >	If RP is negative then the opposite direction (180+XLNG,-XLAT,-RP)
;	is used.
; >	The Sun is located at topocentric longitude 0 deg and latitude 0 deg
;	and radial distance 1.0
; >	The Earth is located at heliocentric longitude 0 deg and latitude
;	0 deg and radial distance 1.0
; >	Internal calculations are done in double precision
; >	The easiest way to check the equations is to work out the relations
;	between the x,y,z components of the vector to P in both coordinates
;	systems.
; MODIFICATION HISTORY:
;	FEB-1990, Paul Hick (UCSD)
;		Adapted from the subroutine SC_ECLIP.FOR
;	JUN-1994, Paul Hick (UCSD)
;		Made the calculation symmetric so that it is valid also going
;		from heliocentric to topocentric.
;-

rpm = ToRadians(degrees=degrees)

sz = size(RP)

IF sz[1] EQ 2 THEN BEGIN

	; If only two coordinates (lng,lat) are given use point-P distance

	IF IsType(RP,/generic_int) THEN RP = float(RP)
	boost, RP, pushdim=0, plus=1
	RP = SubArray(RP,element=2,add=cos(SubArray(RP,element=0)*rpm)*cos(SubArray(RP,element=1)*rpm))
	sz = size(RP)

ENDIF

sz = size(RP)
nR = sz[sz[0]+2]/sz[1]					; # position vectors
IF nR EQ 1 THEN sz1 = size(RP[0]) ELSE sz1 = size( reform(RP[0,*] ) )

RP = reform(RP,sz[1],sz[sz[0]+2]/sz[1])

I		= RP[2,*] LT 0					; For negative RP, change to opposite direction
DLNG	= double(RP[0,*]*rpm)+I*!dpi	; .. add !pi to longitude
I		= 1-2*I
DLAT	= I*double(RP[1,*]*rpm)			; .. opposite latitude
DRP		= I*double(RP[2,*])				; .. make radius positive

COSELO	= cos(DLNG)*cos(DLAT)			; Cos(old elongation)
ELOLD	= float(acos(COSELO))			; Old elongation in [0,!pi]
PHOLD	= float(atan(sin(DLNG)*cos(DLAT),sin(DLAT)))
										; Map longitude to [-!pi,!pi]
EorW	= DLNG-(abs(DLNG) GT !dpi)*(2*(DLNG GE 0)-1)*(2*!dpi)
EorW	= 2*(EorW GE 0)-1				; sign(1.,EorW)

RP = 0.*RP								; Initialize output to zero
ELNEW = reform(RP[0,*])

R		= 1+DRP*DRP-2*DRP*COSELO		; Cosine rule: (dist Sun-P)^2
I		= where(R GT 0)
IF I[0] NE -1 THEN BEGIN
	DLNG	= DLNG[I]
	DLAT	= DLAT[I]
	DRP		= DRP [I]
	R		= sqrt(R[I])				; Dist Sun-P

	RP[2,I]	= R
	SINP	= sin(DLAT)*DRP/R			; Sin(Latitude)
	SINP	= (SINP > (-1.D0)) < 1.D0
	RP[1,I] = asin(SINP)				; Latitude [-!pi/2,!pi/2]
	COSP	= 1.d0-DRP*COSELO
	COSP	= (COSP > (-R)) < R			; prop. cos(LNG)
	ELNEW[I]= acos(COSP/R)				; New elongation
	SINP	= -DRP*cos(DLAT)*sin(DLNG)	; prop. sin(LNG)
	R		= atan(SINP,COSP)			; Longitude in [-!pi,!pi]
	RP[0,I]	= R+(R LT 0)*(2*!dpi)		; 0<=LNG<2*!pi
ENDIF

RP[0:1,*] = RP[0:1,*]/rpm
ELOLD = ELOLD/rpm
PHOLD = PHOLD/rpm
ELNEW = ELNEW/rpm

SyncDims, RP   , size=sz
SyncDims, ELOLD, size=sz1
SyncDims, ELNEW, size=sz1
SyncDims, EORW , size=sz1

RETURN, RP  &  END