-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathstreamtubeXS.m
144 lines (117 loc) · 5.37 KB
/
streamtubeXS.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
function [XS_Nodes,XS_Tubes,TotalFlow,XS_Mask,XS_ks_Node] = streamtubeXS(XS_X, XS_Y, XS_Vel, XS_Depth, XS_ks, NoHorizTubes, ...
NoVertTubes, MaxDryCellsInTube)
%streamtubeXS Generate streamtubes for a single cross-section
% [Nodes,Tubes,TotalFlow] = streamtubeXS(XS_X, XS_Y, XS_Vel, XS_Depth,...
% XS_ks, NoHorizTubes, NoVertTubes)
%
% Richard Measures, Gu Stecca, NIWA
%
% See also delft3d_streamtubes
%% First, some tidying of the input data
% Remove any negative or nan velocities
XS_Vel(XS_Vel<0|isnan(XS_Vel)) = 0;
% Remove any negative or nan depths
XS_Depth(XS_Depth<0|isnan(XS_Depth)) = 0;
DryCells = XS_Depth<=0;
%% Calcuate some basic cross-section properties
% Calculate cell width
CellWidth = sqrt((XS_X(1:end-1)-XS_X(2:end)).^2 + (XS_Y(1:end-1)-XS_Y(2:end)).^2);
CellWidth(isnan(CellWidth)) = 0;
% Calculate cell flow
CellFlow = XS_Vel .* XS_Depth .* CellWidth;
CumulativeFlow = [0;cumsum(CellFlow)];
TotalFlow = sum(CellFlow);
% Find Banks
FlowingCells = CellFlow > 0;
FlowingEdges = [false;FlowingCells]|[FlowingCells;false];
%% Divide XS horizontally into streamtubes
% 1st make sure CumulativeFlow increases so we can use it in interp1
CumFlowOfFlowingEdges = CumulativeFlow(FlowingEdges);
while sum(CumFlowOfFlowingEdges(2:end) == CumFlowOfFlowingEdges(1:end-1))>0
CumFlowOfFlowingEdges([false;(CumFlowOfFlowingEdges(2:end) == CumFlowOfFlowingEdges(1:end-1))]) = ...
CumFlowOfFlowingEdges([false;(CumFlowOfFlowingEdges(2:end) == CumFlowOfFlowingEdges(1:end-1))]) + 1e-8;
end
TubeEdgePos = interp1(CumFlowOfFlowingEdges,...
find(FlowingEdges),...
[0,(TotalFlow/NoHorizTubes)*(1:NoHorizTubes-1),max(CumFlowOfFlowingEdges)]');
TubeEdgeXYDVks = [interp1([XS_X, XS_Y], TubeEdgePos),...
interp1([XS_Depth(1) , XS_Vel(1) , XS_ks(1) ;...
XS_Depth , XS_Vel , XS_ks ;...
XS_Depth(end), XS_Vel(end), XS_ks(end)],...
(TubeEdgePos+0.5))];
TubeEdgePos_L= floor(TubeEdgePos(1:end-1));
TubeEdgePos_R= ceil(TubeEdgePos(2:end))-1;
MaskTubes = ones(size(TubeEdgeXYDVks, 1)-1, 1);
for TubeNo = 1:size(MaskTubes, 1)
NDry = sum(DryCells(TubeEdgePos_L(TubeNo):TubeEdgePos_R(TubeNo)));
MaskTubes(TubeNo) = 1*(NDry<=MaxDryCellsInTube);
end
%% Divide XS vertically into streamtubes
% Calculate vertical divisions
NoBins = 100; % hardcoded!
TubeEdgeZ = zeros(NoHorizTubes+1,NoVertTubes+1);
for EdgeNo = 2:NoHorizTubes
zBins = [1e-6,(1:NoBins)*(TubeEdgeXYDVks(EdgeNo,3)/NoBins)]';
IntegratedVel = log(11*zBins/TubeEdgeXYDVks(EdgeNo,5)).*zBins;
TubeEdgeZ(EdgeNo,2:end-1) = interp1(IntegratedVel,zBins,(1:NoVertTubes-1)*IntegratedVel(end)/NoVertTubes);
TubeEdgeZ(EdgeNo,1:end-1) = -TubeEdgeXYDVks(EdgeNo,3) + TubeEdgeZ(EdgeNo,1:end-1);
end
%% Convert to nodes and tubes (polygons)
% Nodes
% Node numbering = 1 for Left Bank (EdgeNo=1)
% = (EdgeNo-2)*(NoVertTubes+1) + (2:(NoVertTubes+2)) for intermediate verticals (EdgeNo = 2:NoHorizTubes)
% = (NoHorizTubes-1) * (NoVertTubes+1) + 2 for Right Bank (EdgeNo = NoHorizTubes+1)
XS_Nodes = nan((NoHorizTubes-1) * (NoVertTubes+1) + 2, 3);
XS_Nodes(1,:) = [TubeEdgeXYDVks(1,[1,2]),0]; % Left bank
for EdgeNo = 2:NoHorizTubes;
XS_Nodes((EdgeNo-2)*(NoVertTubes+1) + (2:(NoVertTubes+2)), :) = ...
[repmat(TubeEdgeXYDVks(EdgeNo,[1,2]),NoVertTubes+1,1),TubeEdgeZ(EdgeNo,:)'];
end
XS_Nodes(end,:) = [TubeEdgeXYDVks(end,[1,2]),0]; % Right bank
% Tubes
XS_Tubes = cell(NoVertTubes,NoHorizTubes);
XS_Mask = cell(NoVertTubes,NoHorizTubes);
% Left Bank Tubes
for LayerNo = 1:NoVertTubes
XS_Tubes{LayerNo,1} = [1,LayerNo+(1:2),1];
end
% Middle Tubes
for VertNo = 2:NoHorizTubes-1
for LayerNo = 1:NoVertTubes
XS_Tubes{LayerNo,VertNo} = [(VertNo-2)*(NoVertTubes+1)+LayerNo+1,...
(VertNo-1)*(NoVertTubes+1)+LayerNo+[1,2],...
(VertNo-2)*(NoVertTubes+1)+LayerNo+2];
end
end
% Right Bank Tubes
for LayerNo = 1:NoVertTubes
XS_Tubes{LayerNo,end} = [(NoHorizTubes-2)*(NoVertTubes+1)+LayerNo+1,...
(NoHorizTubes-1)*(NoVertTubes+1)+2,...
(NoHorizTubes-2)*(NoVertTubes+1)+LayerNo+[2,1]];
end
for VertNo = 1:NoHorizTubes
for LayerNo = 1:NoVertTubes
XS_Mask{LayerNo,VertNo} = MaskTubes(VertNo,1);
end
end
% ks values in the nodes (only the lower layer can be interpreted)
XS_ks_Node = [XS_ks(TubeEdgePos_L); XS_ks(TubeEdgePos_L(end))];
%% Plot the XS (testing purposes only)
% % Plot vs Y
% plot(XS_Y,-[XS_Depth(1);(XS_Depth(1:end-1)+XS_Depth(2:end))/2;XS_Depth(end)],'r.-')
% hold on
% plot(XS_Nodes(:,2),XS_Nodes(:,3),'go')
% for TubeNo = 1:NoHorizTubes*NoVertTubes
% plot(XS_Nodes(XS_Tubes{TubeNo},2),XS_Nodes(XS_Tubes{TubeNo},3),'b-')
% end
% hold off
%
% % Plot vs X
% plot(XS_X,-[XS_Depth(1);(XS_Depth(1:end-1)+XS_Depth(2:end))/2;XS_Depth(end)],'r.-')
% hold on
% plot(XS_Nodes(:,1),XS_Nodes(:,3),'go')
% for TubeNo = 1:NoHorizTubes*NoVertTubes
% plot(XS_Nodes(XS_Tubes{TubeNo},1),XS_Nodes(XS_Tubes{TubeNo},3),'b-')
% end
% hold off
end