Formula:


θ = atan2(
    sin(Δλ).cos(φ2),
    cos(φ1).sin(φ2) − sin(φ1).cos(φ2).cos(Δλ)
)
Implemented in php-geospatial:

https://github.com/php-geospatial/geospatial/blob/master/geospatial.c#L623