C$Procedure                         ORBBP
C
      DOUBLE PRECISION FUNCTION ORBBP( ORBIT , SUNORB , ALPHA )
C
C
C*******************************************************************************
C
C       Copyright (C) 1993, California Institute of Technology.  U.S.
C       Government Sponsorhip under NASA Contract NAS7-918 is
C       acknowledged.
C
C*******************************************************************************
C
C$ Purpose
C
C  ORBBP computes the beta-prime angle in degrees given the classical
C  elements for the satellite and Sun in a central body centered system.
C  For example, for Topex, the Earth is the central body.
C
C$ Input_Arguments
C
C  ORBIT  - central body centered classical elements of satellite 
C           (a,e,i,LAN,w,M) km,deg
C  SUNORB - central body centered classical elements of the Sun 
C           (a,e,i,LAN,w,M) km,deg
C
C$ Output_Arguments
C
C  ORBBP  - beta-prime in degrees
C  ALPHA  - angle in orbit plane from perigee to projection of earth-sun line
C		into orbit plane
C
C$ Log
C
C  24-May-1989 - Eric Cannell     - creation            
C   9-APR-1992 - Bruce Shapiro
C		 Correct 2 errors (1) call to MA2EA must have MA in deg.
C				  (2) correct calculation of true anomaly.
C  11-SEP-1995 - B. S. add alpha output
C
C$ Restrictions
C
C  1- The mean anomaly to eccentric anomaly will not converge for
C     eccentricities close to 1. See source code for MA2EA.
C
C$ Library_Links
C
C  TPXORB
C
C$ Parameters
C
      DOUBLE PRECISION     PI
      PARAMETER          ( PI     = 3. 14159 26535 89793 23846 D0 )

      DOUBLE PRECISION     D2R
      PARAMETER          ( D2R    = PI / 180.0D0 )
C
C$ Declarations_of_Input_and_Output_Arguments
C
      DOUBLE PRECISION     ORBIT  ( 6 )
      DOUBLE PRECISION     SUNORB ( 6 )
      DOUBLE PRECISION     ALPHA
C
C$ Declarations_of_Local_Variables
C
      DOUBLE PRECISION     AX
      DOUBLE PRECISION     AY
      DOUBLE PRECISION     AZ
      DOUBLE PRECISION     B, BP
      DOUBLE PRECISION     COSB, COSBP
      DOUBLE PRECISION     COSV
      DOUBLE PRECISION     ECANOM
      DOUBLE PRECISION     ESUN
      DOUBLE PRECISION     IORB
      DOUBLE PRECISION     ISUN
      DOUBLE PRECISION     MSUN
      DOUBLE PRECISION     OORB
      DOUBLE PRECISION     OSUN
      DOUBLE PRECISION     SINV
      DOUBLE PRECISION     USUN
      DOUBLE PRECISION     VSUN
      DOUBLE PRECISION     WSUN
      double precision     worb
      DOUBLE PRECISION     WX
      DOUBLE PRECISION     WY
      DOUBLE PRECISION     WZ
      DOUBLE PRECISION     PX, PY, PZ, QX, QY, QZ
      DOUBLE PRECISION     SINW, COSW, SINO, COSO, SINI,COSI
      DOUBLE PRECISION     AINX, AINY, AINZ ,COSALPHA, SINALPHA

C
C$ External_Statements
C
      DOUBLE PRECISION     MA2EA
      EXTERNAL             MA2EA
C        
C$ Method
C-&

C1    Place input orbital elements into local storage.

      IORB = D2R * ORBIT( 3 )
      OORB = D2R * ORBIT( 4 )
      WORB = D2R * ORBIT( 5 )

      COSI = DCOS(IORB)
      SINI = DSIN(IORB)

      COSO = DCOS(OORB)
      SINO = DSIN(OORB)
    
      COSW = DCOS(WORB)
      SINW = DSIN(WORB)

      ESUN =       SUNORB( 2 )
      ISUN = D2R * SUNORB( 3 )
      OSUN = D2R * SUNORB( 4 )
      WSUN = D2R * SUNORB( 5 )
      MSUN = D2R * SUNORB( 6 )

C1    Compute orbit normal vector.

      WX =   DSIN( OORB ) * DSIN( IORB )
      WY = - DCOS( OORB ) * DSIN( IORB )
      WZ =   DCOS( IORB )

C     Compute Line of Perigee unit vector

      PX = COSW*COSO-SINW*SINO*COSI
      PY = COSW*SINO+SINW*COSO*COSI
      PZ = SINW*SINI

C     COMPLETE TRIAD OF NORMAL VECTORS

      QX = -SINW*COSO-COSW*SINO*COSI
      QY = -SINW*SINO+COSW*COSO*COSI
      QZ = COSW*SINI

C1>   Use MA2EA to convert mean anomaly to eccentric anomaly.
C     NOTE:  MA2EA expects the Mean Anomaly IN DEGREES!!!!

      ECANOM = D2R * MA2EA( ESUN , MSUN/D2R )

C1    Compute true anomaly in radians. Ignore the possibility of
C1    both arguments to atan2 being zero, since this will only
C1    happen if eccentricity is one. See restriction #1.

      cosV =   ( DCOS( ECANOM ) - ESUN ) 
     &       / ( 1.0D0 - ESUN * DCOS( ECANOM ) )

      sinV =   ( DSQRT( 1.0D0 - ESUN**2 ) * DSIN( ECANOM ) )
     &       / ( 1.0D0 - ESUN * DCOS( ECANOM ) )

      VSUN = DATAN2( SINV , COSV )

C1    Sum argument of periapsis and true anomaly in radians.

      USUN = WSUN + VSUN

C1    Compute unit position vector of sun.

      AX =   DCOS( USUN ) * DCOS( OSUN ) 
     &     - DCOS( ISUN ) * DSIN( USUN ) * DSIN( OSUN )
      AY =   DCOS( USUN ) * DSIN( OSUN ) 
     &     + DCOS( ISUN ) * DSIN( USUN ) * DCOS( OSUN )
      AZ =   DSIN( USUN ) * DSIN( ISUN )

C1    Compute beta angle via dot product of orbit normal and sun position.

      COSB = WX * AX + WY * AY + WZ * AZ

      B    = DACOS( COSB ) / D2R

      BP   = (90 - B)*D2R
      COSBP = DCOS(BP)

C     COMPUTE PROJECTION OF SUN LINE INTO ORBIT PLANE
	
	AINX = AX - WX*COSBP
	AINY = AY - WY*COSBP
	AINZ = AZ - WZ*COSBP

C     COMPUTE DIRECTION COSINES OF SUN IN ORBIT PLANE

	COSALPHA = AINX*PX + AINY*PY + AINZ*PZ
	SINALPHA = AINX*QX + AINY*QY + AINZ*QZ 

	ALPHA = ATAN2(SINALPHA, COSALPHA)/D2R

C1    Return beta-prime, the complement of beta angle.

      ORBBP = 90 - B

      RETURN
      END