forked from macvicab/MITT
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathInterpUniformChan.m
52 lines (44 loc) · 1.65 KB
/
InterpUniformChan.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
function [oneD,twoD] = InterpUniformChan(Slope,Width,Depth,Length,Sideslope,Widthgrid,Depthgrid,Lengthgrid)
% function to calculate water and bed grid for centerline (oneD) and entire
% section (twoD)
% only works for trapezoidal, rectangular, and triangular channels
% Called from OrganizeInput
% requires use of griddedInterpolant embedded function
% centerline 1D data (common to all channel types)
xmem = 0:Lengthgrid:Length;
nxtot = length(xmem);
bed = zeros(1,nxtot)-Slope*xmem;
wsurf = bed+Depth;
oneD.bedElevation = bed;
oneD.waterElevation = wsurf;
oneD.xchannel = xmem;
% 2D bed and water surface calculation (rectangular, trapezoidal, and
% triangular only)
topwidth = Width + 2*Sideslope*Depth;
ymemi = -topwidth/2:Widthgrid:topwidth/2;
nmiddle = Width/Widthgrid;
% create cross section bed elevation
zcrossside = 0:Widthgrid/Sideslope:Depth;
zcross = [fliplr(zcrossside) zeros(1,nmiddle+1) zcrossside];
% add an extra point on either end of cross sections at the channel depth
ymem = [-topwidth/2-Widthgrid ymemi topwidth/2+Widthgrid];
oneD.ymem = ymem;
% 2D bed and water surface
nytot = length(ymem);
bed2 = zeros(nxtot,nytot);
wsurf2 = zeros(nxtot,nytot);
for nx=1:nxtot
bed2(nx,:) = zcross+bed(nx);
wsurf2(nx,:) = bed(nx)+Depth;
end
[xgrid,ygrid]=ndgrid(xmem,ymem); % grid for 2d xy plane
Fbed = griddedInterpolant(xgrid,ygrid,bed2,'linear','none');
Fwater = griddedInterpolant(xgrid,ygrid,wsurf2,'linear','linear');
twoD.xchannel = xgrid;
twoD.ychannel = ygrid;
twoD.bedElevation = Fbed(xgrid,ygrid);
twoD.waterElevation = Fwater(xgrid,ygrid);
% can get the bed and water surface position of any x and y point
twoD.Fbed = Fbed;
twoD.Fwater = Fwater;
end