| 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 |
% 20060616: bug fixed in CHPsouth, CHPnorth nanmean |
| 43 |
% 20060615: add a default meridional gradient (negative) |
| 44 |
% test on x,y limits |
| 45 |
% |
| 46 |
function varargout = getsurf(Y,X,O,MY,MX,CHP) |
| 47 |
|
| 48 |
%% Dim: |
| 49 |
ny = length(Y); |
| 50 |
nx = length(X); |
| 51 |
%disp(num2str([ny nx])); |
| 52 |
|
| 53 |
%% Indices: |
| 54 |
iymin = max( find( Y<MY(1) ) ); if isempty(iymin),iymin=1; end; |
| 55 |
iymax = min( find( Y>MY(2) ) ); if isempty(iymax),iymax=ny;end; |
| 56 |
ixmin = max( find( X<MX(1) ) ); if isempty(ixmin),ixmin=1; end; |
| 57 |
ixmax = min( find( X>MX(2) ) ); if isempty(ixmax),ixmax=nx;end; |
| 58 |
%disp(num2str([iymin iymax ixmin ixmax])); |
| 59 |
|
| 60 |
%% 1- determine the 2D matrix of surface elements defined by |
| 61 |
% the grid: |
| 62 |
dS = getdS(Y,X); |
| 63 |
|
| 64 |
%% 2- compute the 2D surface matrix where 1 means dS must be |
| 65 |
% counted and 0 must not: |
| 66 |
|
| 67 |
S = ones(ny,nx); % initialy keep all points |
| 68 |
|
| 69 |
% Exclude northward iso-O limits: |
| 70 |
% NB: here the test depends on the meridional gradient of CHP |
| 71 |
% if CHP increase (resp. decreases) with LAT then we must |
| 72 |
% keep lower (resp. higher) values than iso-O limit |
| 73 |
% a: determine test type: |
| 74 |
N = iymax - iymin + 1; % Number of Y points in the domain |
| 75 |
CHPsouth = nanmean(nanmean(squeeze(CHP(iymin:iymin+fix(N/2),ixmin:ixmax)),1),2); |
| 76 |
CHPnorth = nanmean(nanmean(squeeze(CHP(iymin+fix(N/2):iymax,ixmin:ixmax)),1),2); |
| 77 |
SNgrad = (CHPnorth - CHPsouth)./abs(CHPnorth - CHPsouth); |
| 78 |
if isnan(SNgrad), SNgrad=-1; end % Assume negative gradient by default |
| 79 |
%disp(strcat('Northward gradient sign is:',num2str(SNgrad))); |
| 80 |
switch SNgrad |
| 81 |
case 1, testype = 'le'; % Less than or equal |
| 82 |
case -1, testype = 'ge'; % Greater than or equal |
| 83 |
end %switch |
| 84 |
% b: exclude points |
| 85 |
S = double(feval(testype,CHP,O)); |
| 86 |
|
| 87 |
% Exclude southward limit: |
| 88 |
S(1:iymin,:) = zeros(iymin,nx); |
| 89 |
|
| 90 |
% Exclude northward limit: |
| 91 |
S(iymax:ny,:) = zeros((ny-iymax)+1,nx); |
| 92 |
|
| 93 |
% Exclude westward limit: |
| 94 |
S(:,1:ixmin) = zeros(ny,ixmin); |
| 95 |
|
| 96 |
% Exclude eastward limit: |
| 97 |
S(:,ixmax:nx) = zeros(ny,(nx-ixmax)+1); |
| 98 |
|
| 99 |
|
| 100 |
%% 3- Then compute the surface by summing dS elements |
| 101 |
% for non nul S points |
| 102 |
Skeep = S.*dS; |
| 103 |
Skeep = sum(sum(sum(Skeep))); |
| 104 |
|
| 105 |
%% 4- Outputs: |
| 106 |
switch nargout |
| 107 |
case 1 |
| 108 |
varargout(1) = {Skeep}; % Surface single value |
| 109 |
case 2 |
| 110 |
varargout(1) = {Skeep}; |
| 111 |
varargout(2) = {S}; % Logical S matrix with included/excluded points |
| 112 |
case 3 |
| 113 |
varargout(1) = {Skeep}; |
| 114 |
varargout(2) = {S}; |
| 115 |
varargout(3) = {dS}; % Surface elements matrix |
| 116 |
end %switch nargout |
| 117 |
|
| 118 |
|
| 119 |
|
| 120 |
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
| 121 |
% This function computes the 2D dS surface elements. |
| 122 |
function ds = getdS(Y,X); |
| 123 |
|
| 124 |
ny = length(Y); |
| 125 |
nx = length(X); |
| 126 |
|
| 127 |
%%% Compute the DY: |
| 128 |
% Assuming Y is independant of ix: |
| 129 |
d = m_lldist([1 1]*X(1),Y); |
| 130 |
dy = [d(1)/2 (d(2:length(d))+d(1:length(d)-1))/2 d(length(d))/2]; |
| 131 |
dy = meshgrid(dy,X)'; |
| 132 |
|
| 133 |
%%% Compute the DX: |
| 134 |
clear d |
| 135 |
for iy = 1 : ny |
| 136 |
d(:,iy) = m_lldist(X,Y([iy iy])); |
| 137 |
end |
| 138 |
dx = [d(1,:)/2 ; ( d(2:size(d,1),:) + d(1:size(d,1)-1,:) )./2 ; d(size(d,1),:)/2]; |
| 139 |
dx = dx'; |
| 140 |
|
| 141 |
%% Compute the horizontal DS surface element: |
| 142 |
ds = dx.*dy; |
| 143 |
|
Parent Directory
| 
Revision Log
| 
Revision Graph