C+
C NAME:
C	ThomsonFarY
C PURPOSE:
C	(internal use only)
C CALLING SEQUENCE:
	subroutine ThomsonFarY(Y,Elo,I,It,Itr)	! Elo in radians
C INPUTS:
C OUTPUTS:
C INCLUDE:
	include		'math.h'
	include		'thomson.h'
C MODIFICATION HISTORY:
C	1996, Paul Hick (UCSD/CASS; pphick@ucsd.edu)
C-
	    real	Y
	    real	Elo
	    real	I
	    real	It
	    real	Itr

	Elo = max(Elo,MATH__TINY)

	S  = sin(Elo)
	S2 = S*S
	C  = cos(Elo)

	if (Y .ge. LOS__PINF) Y = MATH__PINF

	if (Y .eq. -C) then			! S = 0: start of line of sight
	    I   = 0.
	    It  = 0.
	    Itr = 0.
						! Elongations close to 0 or 180 degrees
	else if (S .lt. sin(LOS__TINYELO*MATH__RPD)) then
	    Y0 = 0.
	    if (Y .gt. 0 .and. Elo .lt. MATH__PI/2.) Y0 = MATH__PI/(S2*S)

	    P = 1.
	    if (Elo .gt. MATH__PI/2) P = -1.

	    if (Y .eq. MATH__PINF) then		! S = Infinity
		Y3 = 0.
		Y5 = 0
	    else				! S > 0
		Y3 = 1./Y**3.
		Y5 = 1./Y**5.
	    end if

	    It  = (1./2.)*Y0-1./3.*(Y3+P)+1./5.*(2*Y5-1./2.*P)*S2
	    Itr = (3./8.)*Y0		 -1./5.*(  Y5+      P)*S2
	    I   = (5./8.)*Y0-2./3.*(Y3+P)+	   Y5	      *S2

	else
	    if (Y .eq. MATH__PINF) then		! S = Infinity
		E   = (MATH__PI-Elo)/S
		It  =               C+E
		Itr = (1.+2./3.*S2)*C+E
		I   = (1.-2./5.*S2)*C+E
	    else				! S > 0
		Y2 = Y*Y
		YS = S2+Y2
		E   = (atan2(Y,S)+MATH__PI/2-Elo)/S
		It  = Y		     / YS	+		 C + E
		Itr = Y*(Y2+5./3.*S2)/(YS*YS)	+ (1+2./3.*S2)	*C + E
		I   = Y*(Y2+3./5.*S2)/(YS*YS)	+ (1-2./5.*S2)	*C + E
	    end if

	    It  = (1./2.)/S2*It
	    Itr = (3./8.)/S2*Itr
	    I   = (5./8.)/S2*I

	end if

	return
	end