Skip to content
Permalink
Branch: master
Find file Copy path
Fetching contributors…
Cannot retrieve contributors at this time
134 lines (109 sloc) 3.79 KB
function [longclatgc,delta]=grcircle(lon1lat1,lon2lat2,longr,latgr)
% [longclatgc,delta]=GRCIRCLE([lon1 lat1],[lon2 lat2],longr,latgr)
%
% Calculates great circle points through a pair of points (maybe on a grid)
%
% INPUT:
%
% [lon1 lat1] longitude and latitude of the start point [radians]
% [lon2 lat2] longitude and latitude of the end point [radians]
% longr regular-grid longitudes OR integer number of points to be calculated
% latgr regular-grid latitudes [optional]
%
% OUTPUT:
%
% longclatgc longitudes/latitudes of the great circle [radians]
% delta incremental distance along the great circle
%
% If the longitudinal difference is 2pi inserts a NaN.
%
% EXAMPLE:
%
% LaJolla=[242.7340 32.8660]*pi/180; StLouis=[269.8160 38.6280]*pi/180;
% Harvard=[288.4420 42.5060]*pi/180; LosAngeles=[241.6330 34.0830]*pi/180;
%% Calculate great circle on a grid
% longr=linspace(0,2*pi,100); latgr=linspace(-pi/2,pi/2,100);
% [lola,delta]=grcircle(Harvard,LosAngeles,longr,latgr);
% plot(lola(:,1)*180/pi,lola(:,2)*180/pi,'b+'); hold on
% [lo,la]=meshgrid(longr,latgr); fridplot(lo*180/pi,la*180/pi)
%% or calculate just a number of points
% [lola,delta]=grcircle(LaJolla,StLouis,100);
% plot(lola(:,1)*180/pi,lola(:,2)*180/pi,'ro');
% plotcont([220 60],[320 20]); hold off; axis image; axis([220 320 20
% 60])
%
% EXAMPLE:
%
%%Compare with:
% [LAT,LON]=track2(LaJolla(2)*180/pi,LaJolla(1)*180/pi,...
% StLouis(2)*180/pi,StLouis(1)*180/pi,[],'degrees',100);
% norm(norm([LON+360 LAT]-[ lola*180/pi]))
%
% Last modified by fjsimons-at-alum.mit.edu, 02/24/2017
defval('longr',100);
[lon1,lat1]=deal(lon1lat1(1),lon1lat1(2));
[lon2,lat2]=deal(lon2lat2(1),lon2lat2(2));
lon1=lon1+(lon1<0)*2*pi;
lon2=lon2+(lon2<0)*2*pi;
[latmax,lonmax]=apex(lat1,lon1,lat2,lon2);
if nargin<=3 & round(longr)==longr
% Fixed number of points
% Calculate total distance
totdis=grcdist([lon1 lat1]*180/pi,[lon2 lat2]*180/pi)/6371;
% 'longr' points along total distance
delta=linspace(0,totdis,longr)';
[lon,lat]=distgrc([lon1 lat1],[lon2 lat2],delta);
lon=lon+2*pi*(lon<0);
longclatgc=[lon lat];
else % Calculation on grid
% Make sure 'latgr' and 'longr' are column vectors
latgr=latgr(:); longr=longr(:);
%-------------------------------------------------
if abs(lon1-lon2)<pi
l=find(longr>min(lon1,lon2) & longr < max(lon1,lon2));
else
l=find(longr<min(lon1,lon2) | longr > max(lon1,lon2));
end
lonmin=rem(latmax,2*pi);
e=tan(latmax);
if sin(lon1-lonmax)*sin(lon2-lonmax)<0
if cos(lon1-lonmax)+cos(lon2-lonmax)>0
w1=find(lat1<latgr & latmax>latgr);
w2=find(lat2<latgr & latmax>latgr);
w=[w1;w2];
if sin(lon1-lonmax)>0
minus=[zeros(size(w1)); ones(size(w2))];
else
minus=[ones(size(w1)); zeros(size(w2))];
end
else
w1=find(lat1>latgr & -latmax<latgr);
w2=find(lat2>latgr & -latmax<latgr);
w=[w1;w2];
if sin(lon1-lonmax)>0
minus=[zeros(size(w1)); ones(size(w2))];
else
minus=[ones(size(w1)); zeros(size(w2))];
end
end
else
w=find(min(lat1,lat2)<latgr & max(lat1,lat2)>latgr);
if sin(lon1-lonmax)>0
minus=zeros(size(w));
else
minus=ones(size(w));
end
end
latca=atan(e.*cos(longr(l)-lonmax));
lonca=(1-2*minus).*acos(tan(latgr(w))./e)+lonmax;
latgc=[latgr(w); latca; lat1; lat2];
longc=[rem(lonca, 2*pi)+2*pi*(lonca<0); longr(l); lon1; lon2];
delta=real(acos(sin(lat1).*sin(latgc)+cos(lat1).*cos(latgc).*cos(lon1-longc)));
[delta,j]=sort(delta);
longclatgc=[longc(j) latgc(j)];
end
% Should only be one number
lopro=find(abs(diff(longclatgc(:,1))-2*pi)<(2*pi/length(longclatgc)));
if ~isempty(lopro)
longclatgc=[longclatgc(1:lopro,:) ; NaN NaN ; longclatgc(lopro+1:end,:)];
end
You can’t perform that action at this time.
You signed in with another tab or window. Reload to refresh your session. You signed out in another tab or window. Reload to refresh your session.