void _geoAzDistanceRads()

in h3_latLng.c [227:282]
51 lines of code
15 McCabe index (conditional complexity)

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

    double sinlat, sinlng, coslng;

    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->lng = 0.0;
        } else if (fabs(p2->lat + M_PI_2) < EPSILON)  // south pole
        {
            p2->lat = -M_PI_2;
            p2->lng = 0.0;
        } else
            p2->lng = constrainLng(p1->lng);
    } 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->lng = 0.0;
        } else if (fabs(p2->lat + M_PI_2) < EPSILON)  // south pole
        {
            p2->lat = -M_PI_2;
            p2->lng = 0.0;
        } else {
            double invcosp2lat = 1.0 / cos(p2->lat);
            sinlng = sin(az) * sin(distance) * invcosp2lat;
            coslng = (cos(distance) - sin(p1->lat) * sin(p2->lat)) /
                     cos(p1->lat) * invcosp2lat;
            if (sinlng > 1.0) sinlng = 1.0;
            if (sinlng < -1.0) sinlng = -1.0;
            if (coslng > 1.0) coslng = 1.0;
            if (coslng < -1.0) coslng = -1.0;
            p2->lng = constrainLng(p1->lng + atan2(sinlng, coslng));
        }
    }
}