subroutine car2kep ( car, kep, mu )
C
C	Converts Cartesian elements to Keplerian Elements
C	Algorithm Given by Bate, Mueller, and White, "Fundamentals of
C	Astrodynamics", Dover
C
C	B. Shapiro 11/24/92 - created to remove dependence on TOPEX libraries
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       parameter name   type description & units
C       --------------   ---- -------------------------------------------
C	input:  car (6)  r*8  x, y, z, vx, vy, vz, km and km/sec
C		mu       r*8  earth gravity, km-3/sec-2
C
C       output: kep (6)  r*8  a (km), e, i,Node,Perigee, M (degrees)
C
C*******************************************************************************

      double precision kep(6), car(6), mu

      double precision eVec(3), r, v, rdotv, first, second, h(3), e, 
     &      N(3), Zhat(3), p, sma, INCL, RAAN, AOP, Ehat(3), Nhat(3), 
     &      Hhat(3), U, RHAT(3), THETA, EANOM, MANOM, SINTHETA,
     &	    COSTHETA, cosE, sinE, RhatDotNhat, NhatDotEhat
      data Zhat /0.0d0, 0.0d0, 1.0d0 /

      double precision dot
      external dot

      double precision pi
      parameter ( pi=3.14159 26535 89793 23846)

      r = sqrt( dot ( car(1), car(1) ) )
      v = sqrt( dot ( car(4), car(4) ) )
      RdotV = dot ( car(1), car(4) )

      first = ( v*v - MU/r ) / MU
      second = -RdotV / MU

      do i = 1, 3
         evec(i) = first * car(i) + second * car(i+3)
      end do
   
      e = sqrt ( dot ( evec, evec ) )

      call UNIT ( EVEC , EHAT)

      call cross ( car(1), car(4), h )
      call UNIT ( H, HHAT) 

      p = dot ( H, H ) / MU
      SMA = P/(1-e*e)

      if (HHAT(3) .ge. 1.0) then
          INCL = 0.0d0
      else if (hhat(3) .le. -1.0) then
          INCL = 180.0d0
      else
          INCL = ACOS ( HHAT(3) ) * 180.0d0 / PI
      end if

      call cross ( zhat, h, n)
      call UNIT ( N, NHAT)

      if (NHAT(1) .ge. 1.0) then
         raan = 0.0d0
      else if (NHAT(2) .le. -1.0) then
         raan = 180.0d0
      else
         RAAN = ACOS ( NHAT(1) ) * 180.0d0 / PI
      end if

      IF ( N(2) .LT. 0.0D0 ) then
         RAAN = 360.0d0 - RAAN
      end if

      if (RAAN .LT. 0.0D0) RAAN = RAAN + 360.0D0

      NhatDotEhat = DOT( NHAT, EHAT )

      if ( nhatdotehat .ge. 1.0d0 ) then
	 aop = 0.0d0
      else if ( nhatdotehat .le. -1.0d0 ) then
	 aop = 180.0d0
      else
         AOP = ACOS ( NhatDotEhat ) * 180.0d0 / PI
      end if

      if (EVEC(3) .LT. 0.0D0 ) then
         AOP = 360.0D0 - AOP
      end if

      IF (AOP .LT. 0.0D0) AOP = AOP + 360.0D0

      CALL UNIT ( CAR(1), RHAT )

C      write(8,*) 'h=',h
C      write(8,*) 'n=',n
C      write(8,*) 'car=',car
C      write(8,*) 'rhat=',rhat
C      write(8,*) 'nhat=',nhat
C      write(8,*) 'dot(rhat,nhat)=',dot(rhat,nhat)

     
      RhatDotNhat = DOT(RHAT,NHAT)
      if (RhatDotNhat .ge. 1.0) then
         U = 0.0d0
      else if (RhatDotNhat .le. -1.0) then
         U = 180.0D0
      else
         U = ACOS ( RhatDotNhat ) * 180.0d0 / PI
      end if

      if (CAR(3) .LT. 0 ) then
         U = 360.0d0 - U
      end if
      THETA = U - AOP

      sinTHETA = SIN ( THETA * PI/180.0d0 )
      cosTHETA = COS ( THETA * PI/180.0d0 )

      cosE = (e + cosTHETA)/(1+e*cosTHETA)
      sinE = (sqrt(1-e*e)*sinTHETA)/(1+e*cosTHETA)
      EANOM = ATAN2 ( sinE, cosE )
      MANOM = ( EANOM - e * sinE ) * 180.0d0 / PI

      if (MANOM .LT. 0.0D0 ) MANOM = MANOM + 360.0D0

      KEP(1) = SMA
      KEP(2) = e
      KEP(3) = INCL
      KEP(4) = RAAN
      KEP(5) = AOP
      KEP(6) = MANOM

      return 
      end


C       program test
C       double precision mu, car(6), kep(6)
C       data mu/398600.436/
C       data car/ 5485.65525126522754817, 423.93438555786292454,
C     &		-5408.29999696842162393, 3.81409491479695123,
C     &		4.39940058458674010, 4.21417533546375223 /
C
C       call car2kep(car, kep, mu)
C       write(6,*) 'car=',car
C       write(6,*) 'kep=',kep
C       end
C