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 % 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=latx1min(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