<SYSTEM_TASK:>
Solve the following problem using Python, implementing the functions described below, one line at a time
<END_TASK>
<USER_TASK:>
Description:
 def VincentyDistance(lon1_in,lat1_in,lon2_in,lat2_in,ELLIPSOID='GRS-80'):

 """
 Calculate the geodesic distance between two points using the formula
 devised by Thaddeus Vincenty, with an accurate ellipsoidal model of the
 earth.
 
 The class attribute `ELLIPSOID` indicates which ellipsoidal model of the
 earth to use. If it is a string, it is looked up in the `ELLIPSOIDS`
 dictionary to obtain the major and minor semiaxes and the flattening.
 Otherwise, it should be a tuple with those values. The most globally
 accurate model is WGS-84. See the comments above the `ELLIPSOIDS`
 dictionary for more information.
 
 """

    lng1 = np.deg2rad(lon1_in)
    lat1 = np.deg2rad(lat1_in)
    lng2 = np.deg2rad(lon2_in)
    lat2 = np.deg2rad(lat2_in)

    major, minor, f = ELLIPSOIDS[ELLIPSOID]

    delta_lng = lng2 - lng1

    reduced_lat1 = np.arctan((1 - f) * np.tan(lat1))
    reduced_lat2 = np.arctan((1 - f) * np.tan(lat2))

    sin_reduced1, cos_reduced1 = np.sin(reduced_lat1), np.cos(reduced_lat1)
    sin_reduced2, cos_reduced2 = np.sin(reduced_lat2), np.cos(reduced_lat2)

    lambda_lng = delta_lng
    lambda_prime = 2 * np.pi

    iter_limit = 20

    while abs(lambda_lng - lambda_prime) > 10e-12 and iter_limit > 0:
        sin_lambda_lng, cos_lambda_lng = np.sin(lambda_lng), np.cos(lambda_lng)

        sin_sigma = np.sqrt( (cos_reduced2 * sin_lambda_lng) ** 2 +
                             (cos_reduced1 * sin_reduced2 -
                              sin_reduced1 * cos_reduced2 * cos_lambda_lng) ** 2
                            )

        if sin_sigma == 0:
            return 0 # Coincident points

        cos_sigma = (
            sin_reduced1 * sin_reduced2 +
            cos_reduced1 * cos_reduced2 * cos_lambda_lng
        )

        sigma = np.arctan2(sin_sigma, cos_sigma)

        sin_alpha = (
            cos_reduced1 * cos_reduced2 * sin_lambda_lng / sin_sigma
        )
        cos_sq_alpha = 1 - sin_alpha ** 2

        if cos_sq_alpha != 0:
            cos2_sigma_m = cos_sigma - 2 * (
                sin_reduced1 * sin_reduced2 / cos_sq_alpha
            )
        else:
            cos2_sigma_m = 0.0 # Equatorial line

        C = f / 16. * cos_sq_alpha * (4 + f * (4 - 3 * cos_sq_alpha))

        lambda_prime = lambda_lng
        lambda_lng = (
            delta_lng + (1 - C) * f * sin_alpha * (
                sigma + C * sin_sigma * (
                    cos2_sigma_m + C * cos_sigma * (
                        -1 + 2 * cos2_sigma_m ** 2
                    )
                )
            )
        )
        iter_limit -= 1

    if iter_limit == 0:
        raise ValueError("Vincenty formula failed to converge!")

    u_sq = cos_sq_alpha * (major ** 2 - minor ** 2) / minor ** 2

    A = 1 + u_sq / 16384. * (
        4096 + u_sq * (-768 + u_sq * (320 - 175 * u_sq))
    )

    B = u_sq / 1024. * (256 + u_sq * (-128 + u_sq * (74 - 47 * u_sq)))

    delta_sigma = (
        B * sin_sigma * (
            cos2_sigma_m + B / 4. * (
                cos_sigma * (
                    -1 + 2 * cos2_sigma_m ** 2
                ) - B / 6. * cos2_sigma_m * (
                    -3 + 4 * sin_sigma ** 2
                ) * (
                    -3 + 4 * cos2_sigma_m ** 2
                )
            )
        )
    )

    s = minor * A * (sigma - delta_sigma)
    return s