void _geoAzDistanceRads()

in src/h3lib/lib/geoCoord.c [205:259]
50 lines of code
15 McCabe index (conditional complexity)

void _geoAzDistanceRads(const GeoCoord* p1, double az, double distance,
                        GeoCoord* p2) {
    if (distance < EPSILON) {
        *p2 = *p1;
        return;
    }

    double sinlat, sinlon, coslon;

    az = _posAngleRads(az);

    // check for due north/south azimuth
    if (az < EPSILON || fabs(az - M_PI) < EPSILON) {
        if (az < EPSILON)  // due north
            p2->lat = p1->lat + distance;
        else  // due south
            p2->lat = p1->lat - distance;

        if (fabs(p2->lat - M_PI_2) < EPSILON)  // north pole
        {
            p2->lat = M_PI_2;
            p2->lon = 0.0;
        } else if (fabs(p2->lat + M_PI_2) < EPSILON)  // south pole
        {
            p2->lat = -M_PI_2;
            p2->lon = 0.0;
        } else
            p2->lon = constrainLng(p1->lon);
    } else  // not due north or south
    {
        sinlat = sin(p1->lat) * cos(distance) +
                 cos(p1->lat) * sin(distance) * cos(az);
        if (sinlat > 1.0) sinlat = 1.0;
        if (sinlat < -1.0) sinlat = -1.0;
        p2->lat = asin(sinlat);
        if (fabs(p2->lat - M_PI_2) < EPSILON)  // north pole
        {
            p2->lat = M_PI_2;
            p2->lon = 0.0;
        } else if (fabs(p2->lat + M_PI_2) < EPSILON)  // south pole
        {
            p2->lat = -M_PI_2;
            p2->lon = 0.0;
        } else {
            sinlon = sin(az) * sin(distance) / cos(p2->lat);
            coslon = (cos(distance) - sin(p1->lat) * sin(p2->lat)) /
                     cos(p1->lat) / cos(p2->lat);
            if (sinlon > 1.0) sinlon = 1.0;
            if (sinlon < -1.0) sinlon = -1.0;
            if (coslon > 1.0) coslon = 1.0;
            if (coslon < -1.0) coslon = -1.0;
            p2->lon = constrainLng(p1->lon + atan2(sinlon, coslon));
        }
    }
}