def _gccalc(lon, lat, azimuth, maxdist=None):
    """
    Original javascript on http://williams.best.vwh.net/gccalc.htm
    Translated into python by Thomas Lecocq
    This function is a black box, because trigonometry is difficult
    
    """
    glat1 = lat * np.pi / 180.
    glon1 = lon * np.pi / 180.
    s = maxdist / 1.852243
    faz = azimuth * np.pi / 180.
 
    EPS = 0.00000000005
    if ((np.abs(np.cos(glat1)) < EPS) and not (np.abs(np.sin(faz)) < EPS)):
        raise CourseException("Only North-South courses are meaningful")

    a = 6378.137 / 1.852243
    f = 1 / 298.257223563
    r = 1 - f
    tu = r * np.tan(glat1)
    sf = np.sin(faz)
    cf = np.cos(faz)
    if (cf == 0):
        b = 0.
    else:
        b = 2. * np.arctan2 (tu, cf)

    cu = 1. / np.sqrt(1 + tu * tu)
    su = tu * cu
    sa = cu * sf
    c2a = 1 - sa * sa
    x = 1. + np.sqrt(1. + c2a * (1. / (r * r) - 1.))
    x = (x - 2.) / x
    c = 1. - x
    c = (x * x / 4. + 1.) / c
    d = (0.375 * x * x - 1.) * x
    tu = s / (r * a * c)
    y = tu
    c = y + 1
    while (np.abs (y - c) > EPS):
        sy = np.sin(y)
        cy = np.cos(y)
        cz = np.cos(b + y)
        e = 2. * cz * cz - 1.
        c = y
        x = e * cy
        y = e + e - 1.
        y = (((sy * sy * 4. - 3.) * y * cz * d / 6. + x) *
            d / 4. - cz) * sy * d + tu

    b = cu * cy * cf - su * sy
    c = r * np.sqrt(sa * sa + b * b)
    d = su * cy + cu * sy * cf
    glat2 = (np.arctan2(d, c) + np.pi) % (2*np.pi) - np.pi
    c = cu * cy - su * sy * cf
    x = np.arctan2(sy * sf, c)
    c = ((-3. * c2a + 4.) * f + 4.) * c2a * f / 16.
    d = ((e * cy * c + cz) * sy * c + y) * sa
    glon2 = ((glon1 + x - (1. - c) * d * f + np.pi) % (2*np.pi)) - np.pi

    baz = (np.arctan2(sa, b) + np.pi) % (2 * np.pi)

    glon2 *= 180./np.pi
    glat2 *= 180./np.pi
    baz *= 180./np.pi
    return (glon2, glat2, baz)
def _gccalc ( lon , lat , azimuth , maxdist = None ) : glat1 = lat * np . pi / 180. glon1 = lon * np . pi / 180. s = maxdist / 1.852243 faz = azimuth * np . pi / 180. EPS = 0.00000000005 if ( ( np . abs ( np . cos ( glat1 ) ) < EPS ) and not ( np . abs ( np . sin ( faz ) ) < EPS ) ) : raise CourseException ( "Only North-South courses are meaningful" ) a = 6378.137 / 1.852243 f = 1 / 298.257223563 r = 1 - f tu = r * np . tan ( glat1 ) sf = np . sin ( faz ) cf = np . cos ( faz ) if ( cf == 0 ) : b = 0. else : b = 2. * np . arctan2 ( tu , cf ) cu = 1. / np . sqrt ( 1 + tu * tu ) su = tu * cu sa = cu * sf c2a = 1 - sa * sa x = 1. + np . sqrt ( 1. + c2a * ( 1. / ( r * r ) - 1. ) ) x = ( x - 2. ) / x c = 1. - x c = ( x * x / 4. + 1. ) / c d = ( 0.375 * x * x - 1. ) * x tu = s / ( r * a * c ) y = tu c = y + 1 while ( np . abs ( y - c ) > EPS ) : sy = np . sin ( y ) cy = np . cos ( y ) cz = np . cos ( b + y ) e = 2. * cz * cz - 1. c = y x = e * cy y = e + e - 1. y = ( ( ( sy * sy * 4. - 3. ) * y * cz * d / 6. + x ) * d / 4. - cz ) * sy * d + tu b = cu * cy * cf - su * sy c = r * np . sqrt ( sa * sa + b * b ) d = su * cy + cu * sy * cf glat2 = ( np . arctan2 ( d , c ) + np . pi ) % ( 2 * np . pi ) - np . pi c = cu * cy - su * sy * cf x = np . arctan2 ( sy * sf , c ) c = ( ( - 3. * c2a + 4. ) * f + 4. ) * c2a * f / 16. d = ( ( e * cy * c + cz ) * sy * c + y ) * sa glon2 = ( ( glon1 + x - ( 1. - c ) * d * f + np . pi ) % ( 2 * np . pi ) ) - np . pi baz = ( np . arctan2 ( sa , b ) + np . pi ) % ( 2 * np . pi ) glon2 *= 180. / np . pi glat2 *= 180. / np . pi baz *= 180. / np . pi return ( glon2 , glat2 , baz )
Original javascript on http://williams.best.vwh.net/gccalc.htm
    Translated into python by Thomas Lecocq
    This function is a black box, because trigonometry is difficult