1 |
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
2 |
%% SUB-FUNCTIONS %% |
3 |
%% USED BY: VOLBET2ISO %% |
4 |
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
5 |
% USE DIRECTLY AT YOUR OWN RISK ! % |
6 |
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
7 |
|
8 |
|
9 |
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
10 |
% Master sub-function: |
11 |
function varargout = subfct_getvol(CHP,Z,Y,X,LIMITS) |
12 |
|
13 |
% Limits: |
14 |
%disp(strcat('Limits: ',num2str(LIMITS))); |
15 |
O = LIMITS(1); |
16 |
MZ = LIMITS(2); |
17 |
MY = sort( LIMITS(3:4) ); |
18 |
MX = sort( LIMITS(5:6) ); |
19 |
|
20 |
% Compute the volume: |
21 |
[V Vmat dV] = getvol(Z,Y,X,O,MZ,MY,MX,CHP); |
22 |
|
23 |
% Outputs: |
24 |
switch nargout |
25 |
case 1 |
26 |
varargout(1) = {V}; |
27 |
case 2 |
28 |
varargout(1) = {V}; |
29 |
varargout(2) = {Vmat}; |
30 |
case 3 |
31 |
varargout(1) = {V}; |
32 |
varargout(2) = {Vmat}; |
33 |
varargout(3) = {dV}; |
34 |
end %switch nargout |
35 |
|
36 |
|
37 |
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
38 |
% This function computes the volume limited southward by |
39 |
% MY(1), northward by iso-O (or MY(2) if iso-O reaches it), |
40 |
% eastward by MX(1), westward by MX(2) and downward by iso-O |
41 |
% (or MZ if iso-O reaches it). |
42 |
% |
43 |
% Updates: |
44 |
% 20060615: add a default meridional gradient (negative) |
45 |
% test on x,y,z limits |
46 |
% |
47 |
function varargout = getvol(Z,Y,X,O,MZ,MY,MX,CHP) |
48 |
|
49 |
%% Dim: |
50 |
nz = length(Z); |
51 |
ny = length(Y); |
52 |
nx = length(X); |
53 |
%disp(num2str([nz ny nx])); |
54 |
|
55 |
%% Indices: |
56 |
izmax = min( find( Z>MZ ) ); if isempty(izmax),izmax=nz;end; |
57 |
iymin = max( find( Y<MY(1) ) ); if isempty(iymin),iymin=1; end; |
58 |
iymax = min( find( Y>MY(2) ) ); if isempty(iymax),iymax=ny;end; |
59 |
ixmin = max( find( X<MX(1) ) ); if isempty(ixmin),ixmin=1; end; |
60 |
ixmax = min( find( X>MX(2) ) ); if isempty(ixmax),ixmax=nx;end; |
61 |
%disp(num2str([1 izmax iymin iymax ixmin ixmax])); |
62 |
|
63 |
%% 1- determine the 3D matrix of volume elements defined by |
64 |
% the grid: |
65 |
dV = getdV(Z,Y,X); |
66 |
|
67 |
%% 2- compute the 3D volume matrix where 1 means dV must be |
68 |
% counted and 0 must not: |
69 |
|
70 |
V = ones(nz,ny,nx); % initialy keep all points |
71 |
|
72 |
% Exclude northward iso-O limits: |
73 |
% NB: here the test depends on the meridional gradient of CHP |
74 |
% if CHP increase (resp. decreases) with LAT then we must |
75 |
% keep lower (resp. higher) values than O limit |
76 |
% a: determine test type: |
77 |
N = iymax - iymin + 1; % Number of Y points in the domain |
78 |
CHPsouth = nanmean(nanmean(squeeze(CHP(1,iymin:iymin+fix(N/2),ixmin:ixmax)))); |
79 |
CHPnorth = nanmean(nanmean(squeeze(CHP(1,iymin+fix(N/2):iymax,ixmin:ixmax)))); |
80 |
SNgrad = (CHPnorth - CHPsouth)./abs(CHPnorth - CHPsouth); |
81 |
if isnan(SNgrad), SNgrad=-1; end % Assume negative gradient by default |
82 |
%disp(strcat('Northward gradient sign is:',num2str(SNgrad))); |
83 |
switch SNgrad |
84 |
case 1, testype = 'le'; % Less than or equal |
85 |
case -1, testype = 'ge'; % Greater than or equal |
86 |
end %switch |
87 |
% b: exclude points |
88 |
for iz=1:izmax |
89 |
chpZ = squeeze(CHP(iz,:,:)); |
90 |
V(iz,:,:)=double(feval(testype,chpZ,O)); |
91 |
end %for iz |
92 |
|
93 |
% Exclude southward limit: |
94 |
V(:,1:iymin,:) = zeros(nz,iymin,nx); |
95 |
|
96 |
% Exclude northward limit: |
97 |
V(:,iymax:ny,:) = zeros(nz,(ny-iymax)+1,nx); |
98 |
|
99 |
% Exclude westward limit: |
100 |
V(:,:,1:ixmin) = zeros(nz,ny,ixmin); |
101 |
|
102 |
% Exclude eastward limit: |
103 |
V(:,:,ixmax:nx) = zeros(nz,ny,(nx-ixmax)+1); |
104 |
|
105 |
% Exclude downward limit: |
106 |
V(izmax:nz,:,:) = zeros((nz-izmax)+1,ny,nx); |
107 |
|
108 |
|
109 |
%% 3- Then compute the volume by summing dV elements |
110 |
% for non 0 V points |
111 |
Vkeep = V.*dV; |
112 |
Vkeep = sum(sum(sum(Vkeep))); |
113 |
|
114 |
%% 4- Outputs: |
115 |
switch nargout |
116 |
case 1 |
117 |
varargout(1) = {Vkeep}; % Volume single value |
118 |
case 2 |
119 |
varargout(1) = {Vkeep}; |
120 |
varargout(2) = {V}; % Logical V matrix with included/excluded points |
121 |
case 3 |
122 |
varargout(1) = {Vkeep}; |
123 |
varargout(2) = {V}; |
124 |
varargout(3) = {dV}; % Volume elements matrix |
125 |
end %switch nargout |
126 |
|
127 |
|
128 |
|
129 |
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
130 |
% This function computes the 3D dV volume elements. |
131 |
function dv = getdV(Z,Y,X); |
132 |
|
133 |
nz = length(Z); |
134 |
ny = length(Y); |
135 |
nx = length(X); |
136 |
|
137 |
%%% Compute the DY: |
138 |
% Assuming Y is independant of ix: |
139 |
d = m_lldist([1 1]*X(1),Y); |
140 |
dy = [d(1)/2 (d(2:length(d))+d(1:length(d)-1))/2 d(length(d))/2]; |
141 |
dy = meshgrid(dy,X)'; |
142 |
|
143 |
%%% Compute the DX: |
144 |
clear d |
145 |
for iy = 1 : ny |
146 |
d(:,iy) = m_lldist(X,Y([iy iy])); |
147 |
end |
148 |
dx = [d(1,:)/2 ; ( d(2:size(d,1),:) + d(1:size(d,1)-1,:) )./2 ; d(size(d,1),:)/2]; |
149 |
dx = dx'; |
150 |
|
151 |
%% Compute the horizontal DS surface element: |
152 |
ds = dx.*dy; |
153 |
|
154 |
%% Compute the DZ: |
155 |
d = diff(Z); |
156 |
dz = [d(1)/2 ( d(2:length(d)) + d(1:length(d)-1) )./2 d(length(d))/2]; |
157 |
|
158 |
%% Then compute the DV volume elements: |
159 |
dv = ones(nz,ny,nx)*NaN; |
160 |
for iz=1:nz |
161 |
dv(iz,:,:) = dz(iz).*ds; |
162 |
end |
163 |
|