function [lat,dlon,refarea,nmr]=authalic(c11,cmn,dlat,dlon) % [lat,dlon,refarea,nmr]=AUTHALIC(c11,cmn,dlat,dlon) % % Calculates an equal-area grid for a regional tomography experiment % % INPUT: % % c11 [x,y]/[lon,lat]-coordinates of the upper left grid corner [degrees] % cmn [x,y]/[lon,lat]-coordinates of the bottom right corner [degrees] % dlat latitude interval, maintained throughout [degrees] % dlon longitude interval at the equator [degrees] % % OUTPUT: % % lat latitudes describing the grid, column vector [degrees] % dlon longitudes describing the grid, column vector [degrees] % refarea reference area, size of one grid cell at the equator % nmr the number of cells per row in the grid % % Only needs to calculate DLON at every LAT to have a complete % representation of the grid. Out comes a column vector, in degrees. For % a given latitude, gives the longitude interval so the area spanned by % that latitude and the next one in the row is conserved. Because we give % C11 and CMN, the 'next' latitude is always the more meridional one. So % the outcoming LAT should always have LAT(1)>LAT(end). Also calculates % the number of columns in each row. Even the not finished cells count as % one, ie those at the end of the row that cannot have the full width % % EXAMPLE I: % % [lat,dlon,c,nmr]=authalic([110 0],[180 -50],1,0.5); % latsup=lat(1:end-1); latsdwn=lat(2:end); % lonslft=zeros(size(latsup)); lonsrgt=dlon; % % This should be the reference area again % difer(spharea([lonslft latsup],[lonsrgt latsdwn])*4*pi*6371^2-c,9); % % EXAMPLE II: % %% A grid for isotropic inversions, one parameter per cell % [lat,dlon,c,nmr1]=authalic([60 15],[200 -65],2/sqrt(3),2/sqrt(3)-0.0047); %% A coarser grid for azimuthal anisotropy, three parameters per cell % [lat,dlon,c,nmr2]=authalic([60 15],[200 -65],2,2); %% Check that both grids yield the same number of parameters % difer(sum(nmr1)-sum(nmr2)*3) % % Last modified by fjsimons-at-alum.mit.edu, 06/12/2007 % Compute reference area on Earth (in km^2) refarea=spharea([0 dlat/2],[dlon -dlat/2])*4*pi*(fralmanac('Radius')/1000)^2; disp(['Reference area is ',num2str(fround(refarea,0)),' square km']) % Adjust longitudinal intervals lat=c11(2):-abs(dlat):cmn(2); lat=lat'; % Conversion to radians lat=lat*pi/180; dlat=dlat*pi/180; % Output (in degrees) dlon=dlon.*abs(2*sin(dlat/2))./abs(sin(lat(1:end-1))-sin(lat(2:end))); lat=lat/pi*180; % Number of cells per row nmr=ceil(abs([cmn(1)-c11(1)])./dlon);