Subroutine ORB2LatLong (ORB, DATE, MU, FLAT, re, 
     *		TOLERANCE, GDLAT, GCLAT, LONGITUDE)
C
C
C	Computes ground track latitude & longitude from a state
C		vector using Escobal transformation (3),
C		"Methods of Orbit Determination", Page 398-399
C
C INPUT:
C		ORB = GTARG STATE VECTOR in standard GTARG format
C			a,e,i,raan,aop,M (deg & km)
C               data = epoch of orb in timetrans format
C               MU = earth GM 
C               FLAT = earth flattening
C               RE = earth equatorial radius in km
C               TOLERANCE = tolerance in degrees of final answer
C OUTPUT:
C               GDLAT = geodetic latitude in degrees
C               GCLAT = geocentric latitude in degrees
C               longitude = longitude in degrees

        double precision orb(6), mu, xyz(6), rag, sidang,
     *     r, alpha, lat, longitude, phi, delta, flat, rc,
     *     H, phiprime, deltaprime, newphiprime, error, 
     *     tolerance, gclat, gdlat, re
  
        external sidang
        character*25 date
        double precision pi
        parameter  (pi = 3.14159 26535 89793 23846 )
        integer maxsteps, istep
        data maxsteps /25/
C
C	convert to cartesian state vector
C        
        call kep2car ( orb, xyz, mu) 
        r = sqrt( xyz(1)**2 + xyz(2)**2 + xyz(3)**2 )
C
C	determine longitude
C
        alpha = atan2 ( xyz(2), xyz(1) ) * ( 180.0d0 / PI )
        rag = sidang ( date, 0.0d0 )
        longitude = alpha - rag

C       write(8,*) 'LONGITUDE = ',longitude
C
C	iterate for latitude
C
        delta = atan2 ( xyz(3), sqrt( xyz(2)**2 + xyz(1)**2 ) )
        newphiprime = delta
        ERROR = 10.0d0
        istep = 0
        do while( ( error .gt. tolerance ) .and.
     &		  ( istep .lt. maxsteps )  )
           istep = istep + 1
           phiprime = newphiprime
           rc = ( 1.0 - (2*flat-flat*flat) ) /
     &		(1.0 - (2*flat-flat*flat)*cos(phiprime)**2 )
           rc = re * sqrt(rc) 

           phi = atan ( tan(phiprime) / ((1-flat)**2) )
           h = sqrt ( r*r - rc*rc*( sin(phi-phiprime)**2 ) ) -
     &		rc * cos(phi-phiprime)
           deltaprime = asin ( (h/r)*sin(phi-phiprime) )

           newphiprime = delta - deltaprime
	   error = dabs ( newphiprime - phiprime ) * (180.0d0/PI)

C	   write(8,*) 'GD = ', newphiprime*180.0d0/pi
        end do
        phiprime = newphiprime
        gclat = phiprime*180.0d0/PI
        gdlat = phi*180.0d0/PI
        return
        end