gmaze_pv/subfct/boxcar.m

% PII = boxcar(C3D,H,X,Y,Z,isoC,dC)
% The boxcar function:
%                               {  isoC-dC/2 <= C3D(iZ,iY,iX) < isoC + dC/2
% PII(isoC,C3D(iZ,iY,iX) = 1 if:{  Z(iZ) > H(iY,iX)
%                        = 0 otherwise
%
% Rq:
% H may be a single value
% Z and H should be negative
% Z orientatd downward
%

%function [PII] = boxcar(C3D,H,X,Y,Z,isoC,dC)

function [PII A B C] = boxcar(C3D,H,X,Y,Z,isoC,dC)

nz  = length(Z);
ny  = length(Y);
nx  = length(X);

method = 2;

if length(H) == 1, H = H.*ones(ny,nx); end

switch method
  case 1 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     PII = zeros(nz,ny,nx); 
     warning off
     for ix = 1 : nx
       for iy = 1 : ny
         Cprof = squeeze(C3D(:,iy,ix));
         li = find( isoC-dC/2 <= Cprof   & ...
                        Cprof < isoC+dC/2 & ...
                            Z > H(iy,ix) );
         if ~isempty(li)
           PII(li,iy,ix) = 1;
         end %if
       end %for iy
     end %for ix
     warning on

  case 2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     PII = ones(nz,ny,nx); 
     
     [a b]=meshgrid(Z,H); b=reshape(b,[ny nx nz]);b=permute(b,[3 1 2]);H=b;clear a b
     [a b c]=meshgrid(Z,Y,X);a=permute(a,[2 1 3]);Z=a;clear a b c
     
     PII(find( -dC/2 < C3D-isoC & C3D-isoC <= +dC/2 & H<=Z  ))  = 0;
     PII = 1-PII;


end %switch method


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Also provide the 1/0 matrix of the layer boundaries:
bounds_vert = zeros(nz,ny,nx);
bounds_meri = zeros(nz,ny,nx);

for ix = 1 : nx
piisect   = squeeze(PII(:,:,ix));
boundsect = zeros(nz,ny);
% Determine vertical boundaries of the layer:
for iy = 1 : ny
  li = find(piisect(:,iy)==1);
  if length(li) ~= 0
    boundsect(li(1),iy)  = 1;
    boundsect(li(end),iy) = 1;
  end
end
bounds_vert(:,:,ix) = boundsect;

boundsect = zeros(nz,ny);
% Determine horizontal meridional boundaries of the layer:
for iz = 1 : nz
  li = find(piisect(iz,:)==1);
  if length(li) ~= 0
    boundsect(iz,li(1))   = 1;
    boundsect(iz,li(end)) = 1;
  end
end
bounds_meri(:,:,ix) = boundsect;

end %for ix

bounds_zona = zeros(nz,ny,nx);
for iy = 1 : ny
piisect   = squeeze(PII(:,iy,:));
boundsect = zeros(nz,nx);
% Determine horizontal zonal boundaries of the layer:
for iz = 1 : nz
  li = find(piisect(iz,:)==1);
  if length(li) ~= 0
    boundsect(iz,li(1))   = 1;
    boundsect(iz,li(end)) = 1;
  end
end
bounds_zona(:,iy,:) = boundsect;
end %for iy

A = bounds_vert;
B = bounds_meri;
C = bounds_zona;

1	% PII = boxcar(C3D,H,X,Y,Z,isoC,dC)
2	% The boxcar function:
3	% { isoC-dC/2 <= C3D(iZ,iY,iX) < isoC + dC/2
4	% PII(isoC,C3D(iZ,iY,iX) = 1 if:{ Z(iZ) > H(iY,iX)
5	% = 0 otherwise
6	%
7	% Rq:
8	% H may be a single value
9	% Z and H should be negative
10	% Z orientatd downward
11	%
12
13	%function [PII] = boxcar(C3D,H,X,Y,Z,isoC,dC)
14
15	function [PII A B C] = boxcar(C3D,H,X,Y,Z,isoC,dC)
16
17	nz = length(Z);
18	ny = length(Y);
19	nx = length(X);
20
21	method = 2;
22
23	if length(H) == 1, H = H.*ones(ny,nx); end
24
25	switch method
26	case 1 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
27	PII = zeros(nz,ny,nx);
28	warning off
29	for ix = 1 : nx
30	for iy = 1 : ny
31	Cprof = squeeze(C3D(:,iy,ix));
32	li = find( isoC-dC/2 <= Cprof & ...
33	Cprof < isoC+dC/2 & ...
34	Z > H(iy,ix) );
35	if ~isempty(li)
36	PII(li,iy,ix) = 1;
37	end %if
38	end %for iy
39	end %for ix
40	warning on
41
42	case 2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
43	PII = ones(nz,ny,nx);
44
45	[a b]=meshgrid(Z,H); b=reshape(b,[ny nx nz]);b=permute(b,[3 1 2]);H=b;clear a b
46	[a b c]=meshgrid(Z,Y,X);a=permute(a,[2 1 3]);Z=a;clear a b c
47
48	PII(find( -dC/2 < C3D-isoC & C3D-isoC <= +dC/2 & H<=Z )) = 0;
49	PII = 1-PII;
50
51
52	end %switch method
53
54
55	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
56	%%% Also provide the 1/0 matrix of the layer boundaries:
57	bounds_vert = zeros(nz,ny,nx);
58	bounds_meri = zeros(nz,ny,nx);
59
60	for ix = 1 : nx
61	piisect = squeeze(PII(:,:,ix));
62	boundsect = zeros(nz,ny);
63	% Determine vertical boundaries of the layer:
64	for iy = 1 : ny
65	li = find(piisect(:,iy)==1);
66	if length(li) ~= 0
67	boundsect(li(1),iy) = 1;
68	boundsect(li(end),iy) = 1;
69	end
70	end
71	bounds_vert(:,:,ix) = boundsect;
72
73	boundsect = zeros(nz,ny);
74	% Determine horizontal meridional boundaries of the layer:
75	for iz = 1 : nz
76	li = find(piisect(iz,:)==1);
77	if length(li) ~= 0
78	boundsect(iz,li(1)) = 1;
79	boundsect(iz,li(end)) = 1;
80	end
81	end
82	bounds_meri(:,:,ix) = boundsect;
83
84	end %for ix
85
86	bounds_zona = zeros(nz,ny,nx);
87	for iy = 1 : ny
88	piisect = squeeze(PII(:,iy,:));
89	boundsect = zeros(nz,nx);
90	% Determine horizontal zonal boundaries of the layer:
91	for iz = 1 : nz
92	li = find(piisect(iz,:)==1);
93	if length(li) ~= 0
94	boundsect(iz,li(1)) = 1;
95	boundsect(iz,li(end)) = 1;
96	end
97	end
98	bounds_zona(:,iy,:) = boundsect;
99	end %for iy
100
101	A = bounds_vert;
102	B = bounds_meri;
103	C = bounds_zona;
104