function tc=truecourse(lon1lat1,lon2lat2) % tc=TRUECOURSE([lon1 lat1],[lon2 lat2]) % % Calculates the true course at the first point between great-circle points % % INPUT: % % [lon1 lat1] longitude and latitude of the start points [degrees] % [lon2 lat2] longitude and latitude of the end points [degrees] % % OUTPUT: % % tc the requested true course % % (True course is the angle between the local meridian and the course line % measured anticlockwise.) Everything in degrees. Input may be two arrays % of (mx2) coordinates. From http://www.best.com/~williams/avform.htm % % EXAMPLE: % % plot([110 100],[-4 -10]); truecourse([110 -4],[100 -10]) % truecourse([100 -10],[110 -4]) % % Last modified by fjsimons-at-alum.mit.edu, 06/13/2007 % For longitudes close together, get minus signs for % both the cosine and the sine of the different angles. % Once you take the double angle, the cos and sine % are the same (if the longitudes are very close together) % Conversion to radians lon1lat1=lon1lat1*pi/180; lon2lat2=lon2lat2*pi/180; [lon1,lat1]=deal(lon1lat1(:,1),lon1lat1(:,2)); [lon2,lat2]=deal(lon2lat2(:,1),lon2lat2(:,2)); tc=mod(atan2(sin(lon1-lon2).*cos(lat2),... cos(lat1).*sin(lat2)-sin(lat1).*cos(lat2).*cos(lon1-lon2)), 2*pi); % Special case: course should be % pi starting from N-pole and 0 starting from S-pole checkmat=cos(lat1) < eps; tc(checkmat)=0; tc(checkmat & lat1>0)=pi; % Conversion to degrees tc=tc*180/pi;