function [slon,slat] = GreatCircle(lon,lat,azim,dist);
%   calculates latitudes and longitudes along a great circle from a reference location and
%   azimuth
% usage: [slon,slat] = GreatCircle(lon,lat,azm,dist);
%
% input parameters:
%
%   lat = latitude on small circle
%   lon = longitude on small circle
%   azim = azimuth of perpendicular to small circle
%   dist = (optional) leave it out and a great circle path is  calculated
%           out to 150 degrees.  include it (as an angle in degrees) and points
%           will be calculate out to that distance.
%
% outout parameters:
%
%   slat = latitudes along great circle
%   slon = glongitudes along great circle
%
% calls scrot and euler_trans
% see also lld2xyz and xyz2lld
if nargin==4,
    dist=0:dist;
else
    dist=[0:1:150];
end
dist=dist(:);
f  = pi/180;
e2 = (lat-90)*f;
e3 =-(lon+90)*f;
phi=zeros(size(dist));

%  need to loop over each distinct azimuth because the rotation matrix can
%  handle multiple distances, but must be reevaluated for each new azimuth.

n=length(azim);
if n==1,
  e1 = (azim-90)*f;
  [th2,ph2,az0]=scrot(euler_trans(e1,e2,e3), dist*f, phi*f);
else
  th2=zeros(size(dist));ph2=th2;az0=th2;
  temp=sort(azim);
  nnn=sum(diff(temp)~=0)+1;            % number of distinct values of azim
  if nnn/n<.8,                         % if more than 20% of azimuths are repeated then
    AZIM=temp([1;find(diff(temp))+1]); % make a list of distinct values of azim
    for i=1:nnn                        % and calculate only one rotation matrix for 
      az1=AZIM(i);                     % each distinct azimuth
      e1 = (az1-90)*f;
      ii=find(azim==AZIM(i));
      [th2(ii),ph2(ii),az0(ii)]=scrot(euler_trans(e1,e2,e3), dist(ii)*f, phi(ii)*f);
    end
  else                                 % else loop over all azimuths
    for ii=1:n
      az1=azim(ii);
      e1 = (az1-90)*f;
      [th2(ii),ph2(ii),az0(ii)]=scrot(euler_trans(e1,e2,e3), dist(ii)*f, phi(ii)*f);
    end
  end
end
slat=90-th2/f;
slon=ph2/f;


function [th2,ph2,az0]=scrot(cxfmat,th1,ph1);
%   scrot         rotation of spherical coordinates
% usage: [th2,ph2,az0]=scrot(cxfmat,th1,ph1);
%
%  input parameters:
%    cxfmat (3x3) transformation matrix (computed using 'euler_trans')
%    th1,ph1      initial coordinates (scalars or column vectors)
%
%  output parameters:
%    th2,ph2  coordinates in rotated system 
%    az0      azimuth in rotated system of northward-directed meridian 
%             in initial system
%   
%  dimensions of th2,ph2,az0 match those of th1,ph1
%  see also ROT, EULER_TRANS, and XYZ2LLD

cth1 = cos(th1);
sth1 = sin(th1);
sph1 = sin(ph1);
cph1 = cos(ph1);
x    = [sth1.*cph1, sth1.*sph1, cth1];
dxdt = [cth1.*cph1, cth1.*sph1,-sth1];
xp   = (cxfmat*x')';
dxpdt= (cxfmat*dxdt')';
temp = sqrt(xp(:,1).*xp(:,1)+xp(:,2).*xp(:,2));
th2  = atan2(temp,xp(:,3));
ph2  = atan2(xp(:,2),xp(:,1));
ph2  = ph2 - 2*pi*sign(ph2).*(abs(ph2)>pi);  % set range to -pi -> pi
az0  = atan2(-xp(:,1).*dxpdt(:,2)+xp(:,2).*dxpdt(:,1) , -dxpdt(:,3));


function a=euler_trans(phi,theta,psi);
%   euler_trans   Euler Transform rotation matrix
% usage: a=euler_trans(phi,theta,psi);
%
%  euler_trans - transformation matrix for coordinate rotation
%  transformation matrix for an arbitrary coordinate rotation specified
%  in terms of the eulerian angles.  conventions used are those of
%  goldstein, classical mechanics, 1959, chapter 4.
%
%  phi,theta,psi must be scalars, output is a 3x3 rotation matrix 
%  see also rot and scrot


cf = cos(phi);
sf = sin(phi);
ct = cos(theta);
st = sin(theta);
cs = cos(psi);
ss = sin(psi);
ctsf = ct*sf;
ctcf = ct*cf;
a = [ cs*cf-ctsf*ss   ,  cs*sf+ctcf*ss  ,  ss*st;
     -ss*cf-ctsf*cs   , -ss*sf+ctcf*cs  ,  cs*st;
      st*sf           , -st*cf          ,  ct    ];

