public LongitudePair getLongitudeOnGreatCircle(Position position, double latitudeDegrees) {
        double lat3 = toRadians(latitudeDegrees);
        double lat1 = toRadians(lat);
        double lon1 = toRadians(lon);
        double lat2 = toRadians(position.getLat());
        double lon2 = toRadians(position.getLon());
        double l12 = lon1 - lon2;
        double sinLat1 = sin(lat1);
        double cosLat2 = cos(lat2);
        double cosLat3 = cos(lat3);
        double cosLat1 = cos(lat1);
        double sinL12 = sin(l12);
        double A = sinLat1 * cosLat2 * cosLat3 * sinL12;
        double B = sinLat1 * cosLat2 * cosLat3 * cos(l12) - cosLat1 * sin(lat2) * cosLat3;
        double C = cosLat1 * cosLat2 * sin(lat3) * sinL12;
        double longitude = atan2(B, A);
        double v = sqrt(sqr(A) + sqr(B));
        if (abs(C) >= v) {
            // not found!
            return null;
        } else {
            double dlon = acos(C / v);
            double lonCandidate1Degrees = to180(FastMath.toDegrees(lon1 + dlon + longitude));
            double lonCandidate2Degrees = to180(FastMath.toDegrees(lon1 - dlon + longitude));
            return new LongitudePair(lonCandidate1Degrees, lonCandidate2Degrees);
        }
    }
public LongitudePair getLongitudeOnGreatCircle ( Position position , double latitudeDegrees ) { double lat3 = toRadians ( latitudeDegrees ) ; double lat1 = toRadians ( lat ) ; double lon1 = toRadians ( lon ) ; double lat2 = toRadians ( position . getLat ( ) ) ; double lon2 = toRadians ( position . getLon ( ) ) ; double l12 = lon1 - lon2 ; double sinLat1 = sin ( lat1 ) ; double cosLat2 = cos ( lat2 ) ; double cosLat3 = cos ( lat3 ) ; double cosLat1 = cos ( lat1 ) ; double sinL12 = sin ( l12 ) ; double A = sinLat1 * cosLat2 * cosLat3 * sinL12 ; double B = sinLat1 * cosLat2 * cosLat3 * cos ( l12 ) - cosLat1 * sin ( lat2 ) * cosLat3 ; double C = cosLat1 * cosLat2 * sin ( lat3 ) * sinL12 ; double longitude = atan2 ( B , A ) ; double v = sqrt ( sqr ( A ) + sqr ( B ) ) ; if ( abs ( C ) >= v ) { return null ; } else { double dlon = acos ( C / v ) ; double lonCandidate1Degrees = to180 ( FastMath . toDegrees ( lon1 + dlon + longitude ) ) ; double lonCandidate2Degrees = to180 ( FastMath . toDegrees ( lon1 - dlon + longitude ) ) ; return new LongitudePair ( lonCandidate1Degrees , lonCandidate2Degrees ) ; } }
TODO add unit test