function [seqnr,gcd,crs,lonsec,latsec]=...
    gridsec(lat,dlon,c11cmn,lon1lat1,lon2lat2,allnum)
% [seqnr,gcd,crs,lonsec,latsec]=...
%      GRIDSEC(lat,dlon,[c11 cmn],[lon1 lat1],[lon2 lat2],allnum)
%
% Finds the intersections (longitude and latitude) of a given
% great circle with an equal-area (hence: irregular) grid
% generated by the program AUTHALIC.
%
% INPUT:
%
% lat         Latitude boundaries of the grid [degrees]
% dlon        One longitudinal cell size for every latitude cell [degrees]
% c11, cmn    Begin- and endpoint of the grid, in degrees
% lon1lat1    Beginpoint of the ray, in radians
% lon2lat2    Endpoint of the ray, in radians
% allnum      Number of columns per row, returned by 'authalic'.
%
% OUTPUT:
%
% seqnr       Sequential cell number, looping over 'columns' first
% gcd         Great-circle distance, in km
% crs         True course, in degrees
% lonsec      Longitudes of the intersection points
% latsec      Latitudes of the intersection points
%
% EXAMPLE
%
% gridsec('demo')
%
% SEE ALSO: SURFTRACE
%
% PROBLEMS:
%
% RECINDEX doesn't want two arrays both to have zeros
% Seems like starting ON a node line gets you in trouble
%
% Last modified by fjsimons-at-alum.mit.edu, May 11th, 2004

if ~isstr(lat)
  
  % Deal
  [c11 cmn]=deal(c11cmn(1:2),c11cmn(3:4));
  [lon1,lat1]=deal(lon1lat1(:,1),lon1lat1(:,2));
  [lon2,lat2]=deal(lon2lat2(:,1),lon2lat2(:,2));
  nmr=allnum;
  
  % Check the great circle is well within the grid
  gridcheck([lon1 lat1 lon2 lat2],c11,cmn)
  
  %disp(sprintf('%8.3f %8.3f %8.3f %8.3f',[lon1 lat1 lon2 lat2]))
  
  % Find longitudes of intersection with the given grid latitudes
  [lonx1,latx1]=lat2lon([lon1 lat1],[lon2 lat2],lat,2);
  
  % This part only if there are such crossings
  if ~isempty(lonx1) & ~isempty(latx1)
    
    % If the only crossings are the lat1 point then switch 1 and 2
    if all(latx1==lat1*ones(size(latx1)))
      [lon1,lat1,lon2,lat2]=deal(lon2,lat2,lon1,lat1);
      flago=0;
    end
    
    % Let's make sure longitudes go from left to right 0<lon<360
    [lonx1,sord]=sort(lonx1);
    latx1=latx1(sord);
    
    % The latx1 will be a subset of lat (save rounding errors)
    prec=10;
    lat=fround(lat,prec);
    latx1=fround(latx1,prec);
    lonx1=fround(lonx1,prec);
    
    % a  is the index number into the latitudes of the latitude crossings
    [a,b]=recindex(lat,latx1);
    relvantlat=a;
    % Find the corresponding dlon at those latitudes
    % Make grid so the rays never cross the two most extreme latitudes
    
    % This part only if not all lonx1 are equal (meridional path)
    % but the condition does not hold if it's only one point
    if sum(diff(lonx1))^length(diff(lonx1))
      
      if prod(size(latx1))~=1
	%----------------------------------------
	% Possible longitudinal crossings
	difflat=sign(diff(latx1));
	% Think the first and last element won't change course
	% There can only be one zero - the highest or lowest point
	% Make the zero a plus or a minus
	a=[a(1)+difflat(1) a]'; % Add a first point
	difflat=[difflat(1) ;  difflat ; difflat(end)];
	
	if  ~prod(difflat) % If there is a turning point
	  if ~difflat(1) % In case first element is zero
	    difflat(1)=NaN;
	  end
	  % In latest Matlab, ~NaN is not allowed anymore, it used to be 0
	  % Next three lines added by fjsimons-at-alum.mit.edu, May 11th, 2004
	  convi=difflat(1:end-1);
	  if ~any(isnan(convi))
	    convo=find(~convi);
	  else
	    convi(isnan(convi))=-987654321;
	    convo=find(~convi);
	  end
	  difflat(convo)=-(latx1(convo)==min(latx1));
	  if ~difflat(end)  % In case last element is zero    
	    difflat(end)=-difflat(end-1);      
	  end 
	  if isnan(difflat(1))
	    difflat(1)=-difflat(2);
	    if ~isnan(find(~difflat))
	      difflat(find(~difflat))=(-(latx1(find(~difflat))==min(latx1)));
	    end
	  end 
	end
	difflat=(-difflat-1)/2;
	% This is either 0 or -1, one per segment of the curve between lat x ->
	% so take the higher or lower latitude - depending on direction
	% 'relvantlon' is the row in the grid
	%----------------------------------------
      else
	a=[a a];
	difflat=[-1 0]+[1 -1]*([lat1>lat2]==[lon1>lon2]);
	%      disp(sprintf('%8.3f %8.3f %8.3f %8.3f',[lon1 lat1 lon2 lat2]))
      end  
      
      relvantlon=a+difflat;    
      dlonx2=dlon(relvantlon);
      
      % The number of cells in the relevant row
      numor=allnum(relvantlon);
      allnum=cumsum(allnum);
      
      % The longitudes of the begin and end points and those that
      % belong to the latitude crossings in between - referenced to c11(1)
      lonanc=([min(lon1,lon2) ; lonx1 ; max(lon1,lon2)]-c11(1));
      
      % The number of times the interval of the possible crossings
      % So this is the cell number in the particular row
      lonxpos=[ceil(lonanc(1:end-1)./dlonx2) floor(lonanc(2:end)./dlonx2)];
      
      % Find the possible longitudinal crossings
      istda=lonxpos(:,1)<=lonxpos(:,2);
      
      % Of course it's possible that no additional crossings are needed
      if sum(istda)
	dlonx2=dlonx2(istda);
	numor=numor(istda);
	relvantlon=relvantlon(istda);
	lonxpos=lonxpos(istda,:)'; % Important transpose
	
	folded=gamini(dlonx2',lonxpos(2,:)-lonxpos(1,:)+1);
	
	foldnumor=gamini(allnum(relvantlon-1)',...
			 lonxpos(2,:)-lonxpos(1,:)+1);
	
	% lonnum is the sequence number of the longitude-crossings
	lonnum=matranges(lonxpos(:)')+foldnumor;
	lonxpos=c11(1)+matranges(lonxpos(:)').*folded;
	
	% Find the corresponding latitudes
	[lonx2,latx2]=lon2lat([lon1 lat1],[lon2 lat2],lonxpos,1);
      else
	[lonx2,latx2,lonnum]=deal([]);
      end
      
      % Still need to calculate the sequence number of the latitude crossings
      % based on relvantlat
      latnum=ceil((lonx1-c11(1))./dlon(relvantlat'-(1+difflat(2:end))));
      latnum=allnum(relvantlat'-(2+difflat(2:end)))+latnum;
      
    else
      lonx2=[];
      latx2=[];
      lonnum=[];
      welk=latx1<max(lat1,lat2) & latx1>min(lat1,lat2) ;
      latx1=latx1(welk);
      lonx1=lonx1(welk);
      relvantlat=relvantlat(welk);
      latnum=ceil((lonx1-c11(1))./dlon(relvantlat));
      allnum=cumsum(allnum);
      latnum=allnum(relvantlat-1)+latnum;
    end
  else  
    relvantlon=find(lat==min(lat(lat>lat1 & lat>lat2)));
    relvantlat=relvantlon;
    lonx1=sort([lon1 lon2])-c11(1);
    dlonx2=dlon(relvantlon);
    lonnum=ceil(lonx1(1:end-1)./dlonx2):floor(lonx1(2:end)./dlonx2);
    lonx1=[c11(1)+dlonx2*lonnum]';
    [lonx2,latx2]=lon2lat([lon1 lat1],[lon2 lat2],lonx1,1);
    allnum=cumsum(allnum);
    lonnum=allnum(relvantlon-1)+lonnum;
    lonx1=[];
    latnum=[];
  end
  
  % This is the end of the calculation proper
  [lonsec,latsec,seqnr]=...
      deal([lonx1 ; lonx2],[latx1 ; latx2],[latnum ; lonnum']);
  
  % seqnr should be unique - just a check we recalculated them later
  if ~all(size(seqnr)==size(unique(seqnr)))
    error('Cell sequence numbers are not unique')
  end
  
  % Calculation of sequence numbers------------------------------
  % Add first and last points to the sequence; one is already in there
  % depending on the sorting in the program 
  lonsec=[lon1 ; lonsec ; lon2];
  latsec=[lat1 ; latsec ; lat2];
  
  % Get distance from FIRST point
  gcd=grcdist([lon1 lat1],[lonsec latsec]);
  outsort=sortrows([gcd lonsec latsec]);
  [lonsec,latsec,gcd]=deal(outsort(:,2),outsort(:,3),outsort(:,1));
  
  if gcd(1)==gcd(2)
    lonsec=lonsec(2:end);
    latsec=latsec(2:end);
    gcd=gcd(2:end);
    %  disp(sprintf('%8.3f %8.3f %8.3f %8.3f',[lon1 lat1 lon2 lat2]))
  end
  
  if gcd(end)==gcd(end-1)
    lonsec=lonsec(1:end-1);
    latsec=latsec(1:end-1);
    gcd=gcd(1:end-1);
  end
  
  seqnr=acor2ind(lat,dlon,nmr,c11,...
		 ([lonsec(1:end-1) latsec(1:end-1)]+...
		  [lonsec(2:end) latsec(2:end)])/2);
  
  % Great-circle distance and azimuth of the ray
  gcd=diff([gcd]);
  crs=truecourse([lonsec(1:end-1) latsec(1:end-1)],...
		 [lonsec(2:end) latsec(2:end)]);
  
  % Check if sum of distances equals total distance up to the meter
  delta=grcdist([lon1 lat1],[lon2 lat2]);
  akku=1/1000;
  if abs(delta-sum(gcd))>akku
    error('Distances are not accurate enough'); 
  end
  
  if exist('flago')
    [seqnr,gcd,crs,lonsec,latsec]=...
	deal(flipud(seqnr),flipud(gcd),flipud(crs),...
	     flipud(lonsec),flipud(latsec));
  end
elseif strcmp(lat,'demo')
  % Runs the example from Simons et al., GJI (2002) - formerly GRIDTEST
  more off
  % Supply the data file - update the path!
  load('/u/fjsimons/MyPapers/2002/GJI-2002/DATA/gji2002lonlat')

  lonlat=gji2002lonlat;

  %lonlat=[120 -89.9 150 -75];

  % Input in degrees
  % First of all, adhere to the convention 0<=lon<=360
  lonlat(:,1)=lonlat(:,1)+(lonlat(:,1)<0)*360;
  lonlat(:,3)=lonlat(:,3)+(lonlat(:,3)<0)*360;

  % First of all, make a grid
  c11=[0 90];
  cmn=[360 -90];
  glat=1;
  glon=1;

  % Check the boundedness of the grid
  gridcheck(lonlat,c11,cmn)

  % Construct the grid
  [lat,dlon,refarea,numonrow]=authalic(c11,cmn,glat,glon);

  % For checking purposes only, submit lonlat in random order
  [sjuf,unsjuf]=shuffle(lonlat);
  %sjuf=sjuf(unsjuf,:);

  % Loop over paths
  [ms,ns]=size(sjuf);
  counter=0;
  for index=1:ms
    counter=counter+1;
    %  disp(sprintf('%4i %s',[counter 'events processed']))
    lon1=sjuf(index,1);
    lat1=sjuf(index,2);
    lon2=sjuf(index,3);
    lat2=sjuf(index,4);

    %  clf
    %  plot([lon1 lon2],[lat1 lat2],'s'); hold on
    %  bb=eqplot(lat,dlon,[c11 cmn]);

    [seqnr,gcd,crs,lonsec,latsec]=...
	gridsec(lat,dlon,[c11 cmn],[lon1 lat1],[lon2 lat2],numonrow);

    reversetest=1;
    if reversetest
      [seqnr,igcd,crs,lonsec,latsec]=...
	  gridsec(lat,dlon,[c11 cmn],[lon2 lat2],[lon1 lat1],numonrow); 
      if ~all(size(gcd)==size(igcd))
	error
      else
	if sum(gcd-flipud(igcd))>0.1
	  error
	end
      end
    end
    
    gcd=igcd;

    % For plotting only----------------------------------------   
    plo=1;
    if plo
      colormap(gray)
      clf
      [longclatgc,delta]=grcircle([lon1 lat1]*pi/180,[lon2 lat2]*pi/180,100);
      %  fillauth(lat,dlon,numonrow,c11,seqnr,gcd); hold on
      fillauth(lat,dlon,numonrow,c11,seqnr,1); hold on
      %fillauth(lat,dlon,numonrow,c11,seqnr,cumsum(gcd)); hold on
      plot(lonsec,latsec,'y+','LineWidth',2); hold on
      plot(longclatgc(:,1)*180/pi,longclatgc(:,2)*180/pi,...
	   'Color','y','LineWidth',0.5); 
      plot([lon1 lon2],[lat1 lat2],'s')
      title(num2str(counter))
      % plotcont([90 10],[180 -60])
      %  bb=eqplot(lat,dlon,[c11 cmn]);
      drawnow
      pause
    end
    % For plotting only----------------------------------------   

  end
  more on
end