C+
C NAME:
C	POINT_ON_LOS
C PURPOSE:
C	Converts from topocentric to heliocentric coordinates and v.v
C CALLING SEQUENCE:
	subroutine POINT_ON_LOS(XLNG,XLAT,RP,ELOLD,ELNEW,iEorW)
C INPUTS:
C	(either topocentric or heliocentric)
C
C	XLNG	real	longitude (0<=XLNG<=360) and ...
C	XLAT	real	.. latitude (-90<=XLAT<=90) of line of sight (degrees)
C	RP	real	radial distance to point P
C			(in units of the observer-Sun distance)
C OUTPUTS:
C	(either heliocentric or topocentric)
C
C	XLNG	real	longitude (0<=XLNG<360) and ...
C	XLAT	real	.. latitude (-90<=XLAT<=90) of point P (degrees)
C	RP	real	distance of point P in units of observer-Sun distance
C	ELOLD	real	elongation in the old coordinate system (deg); 0<=ELOLD<=180
C	ELNEW	real	elongation in the new coordinate system 0<=ELNEW<=180
C	iEorW	integer	+1 or -1; indicates on which side (East or West) of
C			the Sun-Earth line the point is located
C CALLS:
C	dcosd, dsind, dasind, dacosd, datan2d
C RESTRICTIONS:
C >	There is no check for -90 <= XLAT <= 90
C PROCEDURE:
C >	The topocentric coordinate system (with Earth in the origin)
C	has its X-axis pointing towards the Sun. Y and Z-axis are arbitrary
C	(usually the X-Y plane will be the ecliptic).
C	The topocentric longitude is measured in a positive sense, i.e.
C	counterclockwise as viewed from the positive Z-axis.
C >	The heliocentric coordinate system (with the Sun in the origin
C	has its X-axis pointing towards Earth.
C	The heliocentric longitude is measured in a positive sense.
C >	Spherical coordinates for both systems: longitude (deg), latitude
C	(deg) and radial distance (units of Sun-Earth distance).
C >	Input can be in topocentric or heliocentric coordinates. The output
C	will be in the other coordinate system (the calculation is symmetric).
C >	iEorW is determined from the input elongation XLNG:
C	if 0<=XLNG<=180 then iEorW = 1; if 180<XLNG<360 the iEorW = -1
C >	If input is in topocentric coordinates then
C	- ELOLD is the angle Sun-Earth-P
C	- ELNEW is the angle Earth-Sun-P
C	- iEorW = +1 if P is towards the east of the Sun (viewed from Earth)
C	- iEorW = -1 if P is towards the west of the Sun (viewed from Earth)
C >	If input is in heliocentric coordinates then
C	- ELOLD is the angle Earth-Sun-P
C	- ELNEW is the angle Sun-Earth-P
C	- iEorW = -1 if P is towards the east of the Sun (viewed from Earth)
C	- iEorW = +1 if P is towards the west of the Sun (viewed from Earth)
C >	If RP is negative then the opposite direction (180+XLNG,-XLAT,-RP)
C	is used.
C >	The Sun is located at topocentric longitude 0 deg and latitude 0 deg
C	and radial distance 1.0
C >	The Earth is located at heliocentric longitude 0 deg and latitude
C	0 deg and radial distance 1.0
C >	Internal calculations are done in double precision
C >	The easiest way to check the equations is to work out the relations
C	between the x,y,z components of the vector to P in both coordinates
C	systems.
C MODIFICATION HISTORY:
C	FEB-1990, Paul Hick (UCSD/CASS); adapted from the subroutine SC_ECLIP.FOR
C	JUN-1994, Paul Hick (UCSD/CASS; pphick@ucsd.edu); made the calculation
C	    symmetric so that it is valid also going from heliocentric to topocentric.
C-
	    real	XLNG
	    real	XLAT
	    real	RP
	    real	ELOLD
	    real	ELNEW
	    integer	iEorW

	double precision	DLNG
	double precision	DLAT
	double precision	DRP
	double precision	R
	double precision	cosE
	double precision	sinP
	double precision	cosP

	! Needed for g77 when using goniod.f
	! What happens if these are available as intrinsic functions??

	double precision	dsind
	double precision	dcosd
	double precision	dasind
	double precision	dacosd
	double precision	datan2d

	DLNG = dble(XLNG)		! Convert to double
	DLAT = dble(XLAT)
	DRP  = dble(RP  )
	if (RP .lt. 0) then		! Negative RP:
	    DLNG = 180.0d0+DLNG		! .. change to opposite direction
	    DLAT = -DLAT
	    DRP	 = -DRP
	end if

	cosE	= DLNG
	if (dabs(cosE) .gt. 180.0d0) cosE = cosE-dsign(1d0,cosE)*360.0d0
	iEorW	= dnint( dsign(1d0,cosE) )

	cosE	= dcosd(DLNG)*dcosd(DLAT)! Cos(old elongation)
	ELOLD	= dacosd(cosE)		! Old elongation in [0,180]
	
	R	= 1d0+DRP*(DRP-2d0*cosE)! Cosine rule: (dist Sun-P)**2
	if (R .gt. 0.0d0) then
	    R		= dsqrt(R)			! Dist Sun-P
	    RP		= R				! Real*4 output
	    sinP	= dsind(DLAT)*DRP/R		! Sin(Latitude)
	    sinP	= dmin1(1d0,dmax1(-1d0,sinP))
	    XLAT	= sngl(dasind(sinP))		! Latitude [-90,90]
	    cosP	= 1d0-DRP*cosE
	    cosP	= dmin1(R,dmax1(-R,cosP))	! prop. cos(XLNG)
	    ELNEW	= dacosd(cosP/R)		! New elongation
	    sinP	= -DRP*dcosd(DLAT)*dsind(DLNG)	! prop. sin(XLNG)
	    XLNG	= sngl(datan2d(sinP,cosP))
	    XLNG	= amod(amod(XLNG,360.0)+360.0,360.0)! 0<=XLNG<360

	else
	    RP		= 0.0		! (just to prevent 'divide by zero'
	    XLAT 	= 0.0		! condition)
	    XLNG	= 0.0
	    ELOLD	= 0.0
	    ELNEW	= 0.0

	end if

	return
	end

	subroutine PointOnLos(XLNG,XLAT,RP,ELOLD,ELNEW,iEorW)
	    real	XLNG
	    real	XLAT
	    real	RP
	    real	ELOLD
	    real	ELNEW
	    integer	iEorW

	DLNG = XLNG
	DLAT = XLAT
	DRP  = RP
	if (RP .lt. 0) then		! Negative RP:
	    DLNG = 180+DLNG		! .. change to opposite direction
	    DLAT = -DLAT
	    DRP	 = -DRP
	end if

	cosE	= DLNG
	if (abs(cosE) .gt. 180.) cosE = cosE-sign(1.,cosE)*360.
	iEorW	= nint( sign(1.,cosE) )

	cosE	= cosd(DLNG)*cosd(DLAT)	! Cos(old elongation)
	ELOLD	= acosd(cosE)		! Old elongation in [0,180]

	R	= 1+DRP*(DRP-2.*cosE)	! Cosine rule: (dist Sun-P)**2
	if (R .gt. 0) then
	    R		= sqrt(R)			! Dist Sun-P
	    RP		= R				! Real*4 output
	    sinP	= sind(DLAT)*DRP/R		! Sin(Latitude)
	    sinP	= min(1.,max(-1.,sinP))
	    XLAT	= asind(sinP)			! Latitude [-90,90]
	    cosP	= 1-DRP*cosE
	    cosP	= min(R,max(-R,cosP))		! prop. cos(XLNG)
	    ELNEW	= acosd(cosP/R)			! New elongation
	    sinP	= -DRP*cosd(DLAT)*sind(DLNG)	! prop. sin(XLNG)
	    XLNG	= atan2d(sinP,cosP)
	    XLNG	= amod(amod(XLNG,360.)+360.,360.)! 0<=XLNG<360

	else
	    RP		= 0.		! (just to prevent 'divide by zero'
	    XLAT 	= 0.		! condition)
	    XLNG	= 0.
	    ELOLD	= 0.
	    ELNEW	= 0.

	end if

	return
	end