matlab_class/gcmfaces_calc/gcmfaces_section.m

function [secX,secY,secFLD]=gcmfaces_section(lons,lats,fld,varargin);
%purpose: extract a great circle section (defined by two points) from a field
%         or a latitude circle (defined by one latitude)
%
%inputs:        lons/lats are the longitude/latitude vector
%               fld is the gcmfaces field (can incl. depth/time dimensions)
%optional:      sortByLon to sort point by longitude (default = 0 -> sort by latgitude) 
%outputs:       secX/secY is the vector of grid points longitude/latitude
%               secFLD is the vector/matrix of grid point values (from fld)

if nargin>3; sortByLon=varargin{1}; else; sortByLon=0; end;

global mygrid;

if length(lats)==2;
  line_cur=gcmfaces_lines_transp(lons,lats,{'tmp'});
elseif length(lats)==1;
  tmp1=abs(mygrid.LATS-lats);
  tmp2=find(tmp1==min(tmp1));
  tmp2=tmp2(1);
  line_cur=mygrid.LATS_MASKS(tmp2);
else;
  error('wrong specification of lons,lats');
end;
secP=find(line_cur.mskCedge==1);
secN=length(secP);

%lon/lat vectors:
secX=zeros(secN,1); secY=zeros(secN,1); 
%sections:
n3=max(size(fld{1},3),1); n4=max(size(fld{1},4),1); secFLD=zeros(secN,n3,n4);
%counter:
ii0=0; 
for ff=1:secP.nFaces;
  tmp0=secP{ff}; [tmpI,tmpJ]=ind2sub(size(mygrid.XC{ff}),tmp0);
  tmp1=mygrid.XC{ff}; for ii=1:length(tmpI); secX(ii+ii0)=tmp1(tmpI(ii),tmpJ(ii)); end;
  tmp1=mygrid.YC{ff}; for ii=1:length(tmpI); secY(ii+ii0)=tmp1(tmpI(ii),tmpJ(ii)); end;
  tmp1=fld{ff}; for ii=1:length(tmpI); secFLD(ii+ii0,:,:)=squeeze(tmp1(tmpI(ii),tmpJ(ii),:,:)); end;
  ii0=ii0+length(tmpI);
end;

%sort according to increasing latitude or longitude:
if sortByLon; 
  [tmp1,ii]=sort(secX); %sort according to increasing longitude
else;
  [tmp1,ii]=sort(secY); %sort according to increasing latitude
end;
secX=secX(ii); secY=secY(ii); secFLD=secFLD(ii,:,:);

1	gforget	1.2	function [secX,secY,secFLD]=gcmfaces_section(lons,lats,fld,varargin);
2	gforget	1.1	%purpose: extract a great circle section (defined by two points) from a field
3	gforget	1.2	% or a latitude circle (defined by one latitude)
4	gforget	1.1	%
5			%inputs: lons/lats are the longitude/latitude vector
6			% fld is the gcmfaces field (can incl. depth/time dimensions)
7	gforget	1.2	%optional: sortByLon to sort point by longitude (default = 0 -> sort by latgitude)
8	gforget	1.1	%outputs: secX/secY is the vector of grid points longitude/latitude
9			% secFLD is the vector/matrix of grid point values (from fld)
10
11	gforget	1.2	if nargin>3; sortByLon=varargin{1}; else; sortByLon=0; end;
12
13	gforget	1.1	global mygrid;
14
15	gforget	1.2	if length(lats)==2;
16			line_cur=gcmfaces_lines_transp(lons,lats,{'tmp'});
17			elseif length(lats)==1;
18			tmp1=abs(mygrid.LATS-lats);
19			tmp2=find(tmp1==min(tmp1));
20			tmp2=tmp2(1);
21			line_cur=mygrid.LATS_MASKS(tmp2);
22			else;
23			error('wrong specification of lons,lats');
24			end;
25			secP=find(line_cur.mskCedge==1);
26			secN=length(secP);
27	gforget	1.1
28			%lon/lat vectors:
29			secX=zeros(secN,1); secY=zeros(secN,1);
30			%sections:
31	gforget	1.3	n3=max(size(fld{1},3),1); n4=max(size(fld{1},4),1); secFLD=zeros(secN,n3,n4);
32	gforget	1.1	%counter:
33			ii0=0;
34			for ff=1:secP.nFaces;
35			tmp0=secP{ff}; [tmpI,tmpJ]=ind2sub(size(mygrid.XC{ff}),tmp0);
36			tmp1=mygrid.XC{ff}; for ii=1:length(tmpI); secX(ii+ii0)=tmp1(tmpI(ii),tmpJ(ii)); end;
37			tmp1=mygrid.YC{ff}; for ii=1:length(tmpI); secY(ii+ii0)=tmp1(tmpI(ii),tmpJ(ii)); end;
38			tmp1=fld{ff}; for ii=1:length(tmpI); secFLD(ii+ii0,:,:)=squeeze(tmp1(tmpI(ii),tmpJ(ii),:,:)); end;
39			ii0=ii0+length(tmpI);
40			end;
41
42	gforget	1.2	%sort according to increasing latitude or longitude:
43			if sortByLon;
44			[tmp1,ii]=sort(secX); %sort according to increasing longitude
45			else;
46			[tmp1,ii]=sort(secY); %sort according to increasing latitude
47			end;
48	gforget	1.1	secX=secX(ii); secY=secY(ii); secFLD=secFLD(ii,:,:);
49