function Southsaint
% 洪水模型（圣维南方程的差分求解方法）；
% All Right Reserved 城市公共安全应急系统中一维洪水演进及其可视化研究 ;
% Modified by lilongduzhi@yahoo.com.cn hehe:)
% RiverLength = 5089.2
clear all ;
zd = [11.018 ,10.915 ,10.810 ,10.710 ,10.610 , ...
10.510 ,10.410 ,10.300 ,10.200 ,10.100 ,10.000 ];
Zo = [19.518 ,19.515 ,19.510 ,19.510 ,19.510 , ...
19.510 ,19.510 ,19.500 ,19.500 ,19.500 ,19.500 ];
Qo = [30 ,30 ,30 ,30 ,30 ,30 ,30 ,30 ,30 ,30 ,30 ] ;
b = 5; % 底宽 ；
m = 3; % 边坡 ；
I = 0.0002; % 底坡 ；
n = 0.013 ; % 河床糙率 ；
g = 9.8 ;
% 初始赋值 ；
j = 1 ;
z(j ,:) = Zo ;
Q(j ,:) = Qo ;
% 时间步长 ；
t = 0 ;
Times = 30 ;
numOfs = 11 ;
delta_s = 500;
delta_t = 60 ;
% 测试：
Wm =[];
Np = [];
zpp = [];
zll =[];
while j < Times
t = t + delta_t;
% 计算各断面在t时刻的水面宽和过水面积(河道断面为梯形) ；
% 计算断面面积及水面宽 ；
h = z(j ,:) - zd ; % 水深 ；
A = (b + m.*h).*h ; % 过水断面面积 ；
p = b +(1+m^2.*h).^0.5; % 湿周 ；
B = b + 2*m.*h ; % 水面宽 ；
j = j +1 ; % 时间数 ；
for i=1:numOfs
% 计算M点的值（应当是P点的前一时间段的值） ；
AM = A(i);
BM = B(i);
QM = Q(j-1 ,i) ;
zM = z(j-1 ,i) ;
RM = A(i)/p(i) ;
vM = QM/AM ;
M = [vM ,AM ,BM ,QM ,RM ,n ,g];
wminus = vM - (g*AM/BM)^0.5;
wplus = vM + (g*AM/BM)^0.5;
N = BM*(QM/AM)^2*I-g*n^2*QM^2/(AM*RM^(4/3));
% 计算P点的左右节点的值 ；
if i == 1 % 上游 ；
zP = 19.518 ; % 上游边界条件，定值 ；
QE = Q(j-1 ,i+1); zE = z(j-1 ,i+1) ;
zR = (delta_t/delta_s)*wminus*(zM-zE)+zM;
QR = (delta_t/delta_s)*wminus*(QM-QE)+QM ;
QP = QR + BM*wplus*(zP-zR) + N*delta_t ;
elseif i == numOfs % 下游 ；
QD = Q(j-1 ,i-1); zD = z(j-1 ,i-1);
zL = (delta_t/delta_s)*wplus*(zD-zM)+zM ;
QL = (delta_t/delta_s)*wplus*(QD-QM)+QM ;
QP = 30+ t/60*(630-30)/20 ;
QP = min(QP ,630) ;
% 水位计算公式：Bw-(M)*(zP - zL)-QP+QL = N*delta_t ;
zP = (N*delta_t -QL +QP)/(BM*wminus) + zL ;
else % 中游节点
QD = Q(j-1 ,i-1); zD = z(j-1 ,i-1);
QE = Q(j-1 ,i+1); zE = z(j-1 ,i+1);
zL = (delta_t/delta_s)*wplus*(zD-zM)+zM ;
QL = (delta_t/delta_s)*wplus*(QD-QM)+QM ;
zR = (delta_t/delta_s)*wminus*(zM-zE)+zM;
QR = (delta_t/delta_s)*wminus*(QM-QE)+QM;
% 计算P点的值 ；
zP = (QR - QL + BM*wminus*zL-BM*wplus*zR)...
/ (BM *wminus - BM *wplus);
QP = QR + BM*wplus*(zP-zR) + N*delta_t;
end
% 记录数据：
Q(j ,i) = QP ;
z(j ,i) = zP ;
end
end
i = 7 ; % 取第七个切面做观测 ，可任选：
subplot(1,2,1)
plot(z(: ,i)); % 上游流量变化 ；
title(['第',num2str(i),'断面水位变化']);
xlabel('时间');
ylabel('水位');
subplot(1,2,2)
plot(Q(: ,i)); % 上游流量变化 ；
title(['第',num2str(i),'断面流量变化']);
xlabel('时间');
ylabel('流量');
% ------------------------------------------- %
%PS： 文献所述的有限，当初编这个程序有很多的疑问，后解决了，相信对大家了解圣维南方程的求解会有帮助
%研究洪水模型源于参加的研究生数模，因文献的算例有一定的特殊性（给出上下游水位或流量过程线），没有用到提交的论文上
