% calculate intersection coordinates of given coord/bearings
function [lat3, lon3, found] = intersection(lat1, lon1, crs13, lat2, lon2, crs23)
    dst12=2*asin(sqrt((sin((lat1-lat2)/2))^2+cos(lat1)*cos(lat2)*sin((lon1-lon2)/2)^2));
    if sin(lon2-lon1)<0
       crs12=acos((sin(lat2)-sin(lat1)*cos(dst12))/(sin(dst12)*cos(lat1)));
       crs21=2.*pi-acos((sin(lat1)-sin(lat2)*cos(dst12))/(sin(dst12)*cos(lat2)));
    else
       crs12=2.*pi-acos((sin(lat2)-sin(lat1)*cos(dst12))/(sin(dst12)*cos(lat1)));
       crs21=acos((sin(lat1)-sin(lat2)*cos(dst12))/(sin(dst12)*cos(lat2)));
    end

    ang1=mod(crs13-crs12+pi,2.*pi)-pi;
    ang2=mod(crs21-crs23+pi,2.*pi)-pi;

    if (sin(ang1)==0 && sin(ang2)==0)
       %error('infinity of intersections')
       lat3 = 0;
       lon3 = 0;
       found = 0;
    elseif (sin(ang1)*sin(ang2)<0)
       %error('intersection ambiguous')
       lat3 = 0;
       lon3 = 0;
       found = 0;
    else
       ang1=abs(ang1);
       ang2=abs(ang2);
       ang3=acos(-cos(ang1)*cos(ang2)+sin(ang1)*sin(ang2)*cos(dst12));
       dst13=atan2(sin(dst12)*sin(ang1)*sin(ang2),cos(ang2)+cos(ang1)*cos(ang3));
       lat3=asin(sin(lat1)*cos(dst13)+cos(lat1)*sin(dst13)*cos(crs13));
       dlon=atan2(sin(crs13)*sin(dst13)*cos(lat1),cos(dst13)-sin(lat1)*sin(lat3));
       lon3=mod(lon1-dlon+pi,2*pi)-pi;
       found = 1;
    end
end