C$Procedure                         ORBBP
C
      DOUBLE PRECISION FUNCTION ORBBP( ORBIT , SUNORB )
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
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
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 )
C
C$ Declarations_of_Local_Variables
C
      DOUBLE PRECISION     AX
      DOUBLE PRECISION     AY
      DOUBLE PRECISION     AZ
      DOUBLE PRECISION     B
      DOUBLE PRECISION     COSB
      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     WX
      DOUBLE PRECISION     WY
      DOUBLE PRECISION     WZ
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 )

      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 )

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 )

C      write ( 6, *)  'M = ',MSUN/D2R
C      write ( 6, * ) 'E = ',ECANOM/D2R
C      write ( 6, * ) 'V = ',VSUN/D2R
C      write ( 6, * ) 'U = ', USUN/D2R
C      write ( 6, * ) 'SUN = ',Ax,' ',ay,' ',az
C      write ( 6, * ) 'S/C Normal = ',Wx,' ',wy,' ',wz
C1    Compute beta angle via dot product of orbit normal and sun position.

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

      B    = DACOS( COSB ) / D2R

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

      ORBBP = 90 - B

      RETURN
      END