code: def geodetic2geocentric(theta, alt):
      
      ct = np.cos(theta)
      st = np.sin(theta)
      a2 = 40680631.6
      b2 = 40408296.0
      one = a2 * st * st
      two = b2 * ct * ct
      three = one + two
      rho = np.sqrt(three)
      r = np.sqrt(alt * (alt + 2.0 * rho) + (a2 * one + b2 * two) / three)
      cd = (alt + rho) / r
      sd = (a2 - b2) / rho * ct * st / r
      one = ct
      ct = ct * cd - st * sd
      st = st * cd + one * sd
      gccolat = np.arctan2(st, ct)
      d = np.arctan2(sd, cd)
      return gccolat, d, r
docs: Conversion from geodetic to geocentric coordinates by using the WGS84 spheroid.
      :param theta: colatitude (float, rad)
      :param alt: altitude (float, km)
      :return gccolat: geocentric colatitude (float, rad)
              d: gccolat minus theta (float, rad)
              r: geocentric radius (float, km)