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