SUBROUTINE RADGRAV ( ALT, PHID, G, RE )
C====================================================================C
C                                                                    C
C     SUBROUTINE RADGRAV (ALT, PHID, G, R)                           C
C                                                                    C
C     PURPOSE:                                                       C
C                                                                    C
C     THIS ROUTINE PROVIDES THE ACCELERATION OF GRAVITY AND RADIUS   C
C     OF THE EARTH AT A SPECIFIED GEODETIC LATITUDE AND ALTITUDE     C
C                                                                    C
C     PARAMETERS:                                                    C
C         INPUT                                                      C
C              ALT.....ALTITUDE (KM) ABOVE THE GEOID, MEASURED       C
C                      ALONG THE LOCAL GEOID VERTICAL                C
C              PHID....GEODETIC LATITUDE (DEGREES) OF A SURFACE      C
C                      POINT                                         C
C         OUTPUT                                                     C
C              G.......ACCELERATION OF GRAVITY (M/SEC2)              C
C              RE......RADIUS OF THE ELLIPSOIDAL EARTH               C
C                                                                    C
C                                                                    C
C     SOURCE:        JIM BUGLIA AS MODIFIED BY R. K. SEALS, JR.      C
C                    NASA LANGLEY RESEARCH CENTER                    C
C                    MODIFIED BY EARL THOMPSON                       C
C                    MODIFIED BY B.T. MARSHALL                       C
C                                                                    C
C     DATE STARTED:      APRIL 1989                                  C
C     LATEST REVISION:   MARCH 30,1995                               C
C                                                                    C
C====================================================================C
      REAL*4 MU, J2
      COSD(X) = COS(1.745329252D-2*X)
      SIND(X) = SIN(1.745329252D-2*X)
C      TAND(X) = TAN(1.745329252D-2*X)
      ATAN2D(X,Y) = 57.2957795131D0*ATAN2(X,Y)
C                                      A = EQUATORIAL RADIUS, KM
C                                      B = POLAR RADIUS, KM
C                                      MU = GRAVITATIONAL CONSTANT, KM**3/SEC**2
C                                      J2 = FIRST OBLATENESS CONSTANT
C                                      W = ROTATIONAL RATE OF THE EARTH, RAD/SEC
      A = 6378.140
      B = 6356.755
      MU = 398600.6
      J2 = 1082.28E-6
      W=7.292124E-5
      BA2=(B/A)**2
C                                      X AND Y POSITION CALCULATIONS
      C=A/SQRT((COSD(PHID)**2)+BA2*(SIND(PHID)**2))
      S=BA2*C
      X=(C+ALT)*COSD(PHID)
      Y=(S+ALT)*SIND(PHID)
      R=SQRT(X*X+Y*Y)
C                                      RADIUS OF EARTH
      X1=C*COSD(PHID)
      Y1=S*SIND(PHID)
      RE=SQRT(X1*X1+Y1*Y1)
C                                      GEOCENTRIC LATITUDE CALCULATION
      PHIC=ATAN2D(Y,X)
C                                      CALCULATE GRAVITY VECTOR
  150 G0=MU/(R*R)
      C1=1.0-1.5*((A/R)**2)*J2*(3.*(SIND(PHIC)**2)-1.)
      GR=-G0*C1
      C2=-3.*((A/R)**2)*J2*SIND(PHIC)*COSD(PHIC)
      GP=G0*C2
      GX=GR*COSD(PHIC)-GP*SIND(PHIC)
      GY=GR*SIND(PHIC)+GP*COSD(PHIC)
      GX=GX+W*W*X
      G=SQRT(GX*GX+GY*GY)*1000.
      RETURN
      END