⭐ 欢迎来到虫虫下载站！ | 📦 资源下载 📁 资源专辑 ℹ️ 关于我们
⭐ 虫虫下载站
📤 上传资源

📄 tidepredict_main.m

📁 MATLAB语言下的(有源码):Mike Zero系列软件的前后处理 GE中岸线与MIKE岸线相互转换 Mapinfo与MIKE中数据相互转换.(对学习MIKE与学习MATLAB编程人皆有帮助,且有说
💻 M
字号:
🔍 
function [Result]=TidePredict_Main(DataType,StartTime,Interval,TotalNum,TideConst,DataAverage)
%[Result]=TidePredict(DataType,StartTime,Interval,TotalNum,TideConst,DataAverage)
%功能基于O1、K1、M2、S2、M4、MS4、P1、Q1、N2、K2十个分潮的调和常数进行潮位/潮流预报
%DataType，取值为2时为潮位预报，取值为1时为潮流预报
%StartTime，开始计算时间，一维数组，数组单元数为4，分别为年、月、日和时，其中时可以为小数
%Interval,数据时间间隔，单位为小时，可以为小数
%TotalNum,总预报数据个数
%TideConst，调和常数，DataType=2时，二维数组，第一列为振幅，第二列为迟角（单位°，正东为90，正北为0）
%TideConst，调和常数，DataType=1时，二维数组，第一列为东分量振幅，第二列为东分量迟角（单位°，正东为90，正北为0）
%第三列为北分量振幅，第三列为北分量迟角（单位°，正东为90，正北为0）
%TideConst，各行顺序依次为O1、K1、M2、S2、M4、MS4、P1、Q1、N2、K2
%DataAverage，当DataType=2时为余水位；当DataType=1时为余流，第一个记录为余流东分量，第二个记录为余流北分量
%Result，DataType=2时，二维数组，第一列为时间（datestr(Result(:,1))既可转换为时间字串），第二列为预报潮位
%Result，DataType=1时，二维数组，第一列为时间（datestr(Result(:,1))既可转换为时间字串），第二列为预报流速
%，第三列为预报流向
%%
TideConst(:,2)=TideConst(:,2).*pi./180;
%计算分潮角速率
D_tshpp=[14.49205211 0.54901653 0.04106864 0.00464183 0.00000196]; %式10.5中的系数
D_tshpp=D_tshpp.*pi./180; %系数转换为弧度
%
Eta_O1 =(1)*D_tshpp(1)+(-1)*D_tshpp(2)+ (0)*D_tshpp(3)+(0)*D_tshpp(4)+(0)*D_tshpp(5);
Eta_K1 =(1)*D_tshpp(1)+ (1)*D_tshpp(2)+ (0)*D_tshpp(3)+(0)*D_tshpp(4)+(0)*D_tshpp(5);
Eta_M2 =(2)*D_tshpp(1)+ (0)*D_tshpp(2)+ (0)*D_tshpp(3)+(0)*D_tshpp(4)+(0)*D_tshpp(5);
Eta_S2 =(2)*D_tshpp(1)+ (2)*D_tshpp(2)+(-2)*D_tshpp(3)+(0)*D_tshpp(4)+(0)*D_tshpp(5);
Eta_M4 =(4)*D_tshpp(1)+ (0)*D_tshpp(2)+ (0)*D_tshpp(3)+(0)*D_tshpp(4)+(0)*D_tshpp(5);
Eta_MS4=(4)*D_tshpp(1)+ (2)*D_tshpp(2)+(-2)*D_tshpp(3)+(0)*D_tshpp(4)+(0)*D_tshpp(5);
Eta_P1 =(1)*D_tshpp(1)+ (1)*D_tshpp(2)+(-2)*D_tshpp(3)+(0)*D_tshpp(4)+(0)*D_tshpp(5);
Eta_Q1 =(1)*D_tshpp(1)+(-2)*D_tshpp(2)+ (0)*D_tshpp(3)+(1)*D_tshpp(4)+(0)*D_tshpp(5);
Eta_N2 =(2)*D_tshpp(1)+(-1)*D_tshpp(2)+ (0)*D_tshpp(3)+(1)*D_tshpp(4)+(0)*D_tshpp(5);
Eta_K2 =(2)*D_tshpp(1)+ (2)*D_tshpp(2)+ (0)*D_tshpp(3)+(0)*D_tshpp(4)+(0)*D_tshpp(5);
%
%计算天文要素,tao,s,h,p,N,ps
[s,h,p,N,ps]=shpnps_predict(StartTime(1),StartTime(2),StartTime(3),0);%当日零时的天文要素
Tao=pi/12.*0-s+h;
%计算分潮天文相角
v0_O1 =(-1)*pi/2 + (1)*Tao + (-1)*s + (0)*h + (0)*p + (0)*N + (0)*ps;
v0_K1 = (1)*pi/2 + (1)*Tao +  (1)*s + (0)*h + (0)*p + (0)*N + (0)*ps;
v0_M2 = (0)*pi/2 + (2)*Tao +  (0)*s + (0)*h + (0)*p + (0)*N + (0)*ps;
v0_S2 = (0)*pi/2 + (2)*Tao +  (2)*s +(-2)*h + (0)*p + (0)*N + (0)*ps;
v0_M4 = (0)*pi/2 + (4)*Tao +  (0)*s + (0)*h + (0)*p + (0)*N + (0)*ps;
v0_MS4= (0)*pi/2 + (4)*Tao +  (2)*s +(-2)*h + (0)*p + (0)*N + (0)*ps;
v0_P1 =(-1)*pi/2 + (1)*Tao +  (1)*s +(-2)*h + (0)*p + (0)*N + (0)*ps;
v0_Q1 =(-1)*pi/2 + (1)*Tao + (-2)*s + (0)*h + (1)*p + (0)*N + (0)*ps;
v0_N2 = (0)*pi/2 + (2)*Tao + (-1)*s + (0)*h + (1)*p + (0)*N + (0)*ps;
v0_K2 = (0)*pi/2 + (2)*Tao +  (2)*s + (0)*h + (0)*p + (0)*N + (0)*ps;
%
for c1=1:TotalNum
    CurrentTime=StartTime(4)+Interval*(c1-1);
    %%
    [s,h,p,N,ps]=shpnps_predict(StartTime(1),StartTime(2),StartTime(3),CurrentTime);
    %%
    %计算分潮f u
    FCos_O1=-0.0058*cos( 0*p-2*N)+0.1885*cos(0*p-1*N)+1*cos(0*p-0*N)+0.0002*cos(2*p-1*N)-0.0064*cos(2*p+0*N)-0.001*cos(2*p+1*N);
    FSin_O1=-0.0058*sin( 0*p-2*N)+0.1885*sin(0*p-1*N)+1*sin(0*p-0*N)+0.0002*sin(2*p-1*N)-0.0064*sin(2*p+0*N)-0.001*sin(2*p+1*N);
    FCos_K1= 0.0002*cos(-2*p-1*N)+0.0001*cos(0*p-2*N)+1*cos(0*p-0*N)-0.0198*cos(0*p-1*N)+0.1356*cos(0*p+1*N)-0.0029*cos(0*p+2*N);
    FSin_K1= 0.0002*sin(-2*p-1*N)+0.0001*sin(0*p-2*N)+1*sin(0*p-0*N)-0.0198*sin(0*p-1*N)+0.1356*sin(0*p+1*N)-0.0029*sin(0*p+2*N);
    FCos_M2= 0.0005*cos(0*p-2*N) -0.0373*cos(0*p-1*N)+1*cos(0*p-0*N)+0.0006*cos(2*p+0*N)+0.002*cos(2*p+1*N);
    FSin_M2= 0.0005*sin(0*p-2*N) -0.0373*sin(0*p-1*N)+1*sin(0*p-0*N)+0.0006*sin(2*p+0*N)+0.002*sin(2*p+1*N);
    FCos_P1=0.0008*cos( 0*p-2*N)-0.0112*cos(0*p-1*N)+1*cos(0*p-0*N)-0.0015*cos(2*p+0*N)-0.0003*cos(2*p+1*N);
    FSin_P1=0.0008*sin( 0*p-2*N)-0.0112*sin(0*p-1*N)+1*sin(0*p-0*N)-0.0015*sin(2*p+0*N)-0.0003*sin(2*p+1*N);
    FCos_K2=1*cos(0*p-0*N)-0.0128*cos(0*p-1*N)+0.2980*cos(0*p+1*N)+0.0324*cos(0*p+2*N);
    FSin_K2=1*sin(0*p-0*N)-0.0128*sin(0*p-1*N)+0.2980*sin(0*p+1*N)+0.0324*sin(0*p+2*N);
    %
    f_O1=sqrt(FCos_O1.^2+FSin_O1.^2);
    f_K1=sqrt(FCos_K1.^2+FSin_K1.^2);
    f_M2=sqrt(FCos_M2.^2+FSin_M2.^2);
    f_P1=sqrt(FCos_P1.^2+FSin_P1.^2);
    f_K2=sqrt(FCos_K2.^2+FSin_K2.^2);
    u_O1=asin(FSin_O1./f_O1);
    u_K1=asin(FSin_K1./f_K1);
    u_M2=asin(FSin_M2./f_M2);
    u_P1=asin(FSin_P1./f_P1);
    u_K2=asin(FSin_K2./f_K2);
    f_S2=1;
    u_S2=0;
    f_M4=f_M2*f_M2;
    u_M4=2*u_M2;
    f_MS4=f_M2;
    u_MS4=u_M2;
    f_Q1=f_O1;
    u_Q1=u_O1;
    f_N2=f_M2;
    u_N2=u_M2;
    %%%%%%%%%%%%%%%
    if DataType==2 %潮位预报
        E_Resudial=DataAverage(1);
        E_H_O1=TideConst(1,1);
        E_g_O1=TideConst(1,2);
        E_H_K1=TideConst(2,1);
        E_g_K1=TideConst(2,2);
        E_H_M2=TideConst(3,1);
        E_g_M2=TideConst(3,2);
        E_H_S2=TideConst(4,1);
        E_g_S2=TideConst(4,2);
        E_H_M4=TideConst(5,1);
        E_g_M4=TideConst(5,2);
        E_H_MS4=TideConst(6,1);
        E_g_MS4=TideConst(6,2);
        E_H_P1=TideConst(7,1);
        E_g_P1=TideConst(7,2);
        E_H_Q1=TideConst(8,1);
        E_g_Q1=TideConst(8,2);
        E_H_N2=TideConst(9,1);
        E_g_N2=TideConst(9,2);
        E_H_K2=TideConst(10,1);
        E_g_K2=TideConst(10,2);
        %
        E_Part_O1=f_O1.*E_H_O1.*cos(Eta_O1.*CurrentTime+v0_O1+u_O1-E_g_O1);
        E_Part_K1=f_K1.*E_H_K1.*cos(Eta_K1.*CurrentTime+v0_K1+u_K1-E_g_K1);
        E_Part_M2=f_M2.*E_H_M2.*cos(Eta_M2.*CurrentTime+v0_M2+u_M2-E_g_M2);
        E_Part_S2=f_S2.*E_H_S2.*cos(Eta_S2.*CurrentTime+v0_S2+u_S2-E_g_S2);
        E_Part_M4=f_M4.*E_H_M4.*cos(Eta_M4.*CurrentTime+v0_M4+u_M4-E_g_M4);
        E_Part_MS4=f_MS4.*E_H_MS4.*cos(Eta_MS4.*CurrentTime+v0_MS4+u_MS4-E_g_MS4);
        E_Part_P1=f_P1.*E_H_P1.*cos(Eta_P1.*CurrentTime+v0_P1+u_P1-E_g_P1);
        E_Part_Q1=f_Q1.*E_H_Q1.*cos(Eta_Q1.*CurrentTime+v0_Q1+u_Q1-E_g_Q1);
        E_Part_N2=f_N2.*E_H_N2.*cos(Eta_N2.*CurrentTime+v0_N2+u_N2-E_g_N2);
        E_Part_K2=f_K2.*E_H_K2.*cos(Eta_K2.*CurrentTime+v0_K2+u_K2-E_g_K2);
        Result(c1,1)=datenum([StartTime(1:3),CurrentTime,0,0]);
        Result(c1,2)=E_Resudial+E_Part_O1+E_Part_K1+E_Part_M2+E_Part_S2+E_Part_M4+E_Part_MS4+E_Part_P1+E_Part_Q1+E_Part_N2+E_Part_K2;
    else    %潮流预报
        %%东分量预报
        E_Resudial=DataAverage(1);
        E_H_O1=TideConst(1,1);
        E_g_O1=TideConst(1,2);
        E_H_K1=TideConst(2,1);
        E_g_K1=TideConst(2,2);
        E_H_M2=TideConst(3,1);
        E_g_M2=TideConst(3,2);
        E_H_S2=TideConst(4,1);
        E_g_S2=TideConst(4,2);
        E_H_M4=TideConst(5,1);
        E_g_M4=TideConst(5,2);
        E_H_MS4=TideConst(6,1);
        E_g_MS4=TideConst(6,2);
        E_H_P1=TideConst(7,1);
        E_g_P1=TideConst(7,2);
        E_H_Q1=TideConst(8,1);
        E_g_Q1=TideConst(8,2);
        E_H_N2=TideConst(9,1);
        E_g_N2=TideConst(9,2);
        E_H_K2=TideConst(10,1);
        E_g_K2=TideConst(10,2);
        %
        E_Part_O1=f_O1.*E_H_O1.*cos(Eta_O1.*CurrentTime+v0_O1+u_O1-E_g_O1);
        E_Part_K1=f_K1.*E_H_K1.*cos(Eta_K1.*CurrentTime+v0_K1+u_K1-E_g_K1);
        E_Part_M2=f_M2.*E_H_M2.*cos(Eta_M2.*CurrentTime+v0_M2+u_M2-E_g_M2);
        E_Part_S2=f_S2.*E_H_S2.*cos(Eta_S2.*CurrentTime+v0_S2+u_S2-E_g_S2);
        E_Part_M4=f_M4.*E_H_M4.*cos(Eta_M4.*CurrentTime+v0_M4+u_M4-E_g_M4);
        E_Part_MS4=f_MS4.*E_H_MS4.*cos(Eta_MS4.*CurrentTime+v0_MS4+u_MS4-E_g_MS4);
        E_Part_P1=f_P1.*E_H_P1.*cos(Eta_P1.*CurrentTime+v0_P1+u_P1-E_g_P1);
        E_Part_Q1=f_Q1.*E_H_Q1.*cos(Eta_Q1.*CurrentTime+v0_Q1+u_Q1-E_g_Q1);
        E_Part_N2=f_N2.*E_H_N2.*cos(Eta_N2.*CurrentTime+v0_N2+u_N2-E_g_N2);
        E_Part_K2=f_K2.*E_H_K2.*cos(Eta_K2.*CurrentTime+v0_K2+u_K2-E_g_K2);
        E_Result=E_Resudial+E_Part_O1+E_Part_K1+E_Part_M2+E_Part_S2+E_Part_M4+E_Part_MS4+E_Part_P1+E_Part_Q1+E_Part_N2+E_Part_K2;
        %%
        %北分量计算
        N_Resudial=DataAverage(2);
        N_H_O1=TideConst(1,1);
        N_g_O1=TideConst(1,2);
        N_H_K1=TideConst(2,1);
        N_g_K1=TideConst(2,2);
        N_H_M2=TideConst(3,1);
        N_g_M2=TideConst(3,2);
        N_H_S2=TideConst(4,1);
        N_g_S2=TideConst(4,2);
        N_H_M4=TideConst(5,1);
        N_g_M4=TideConst(5,2);
        N_H_MS4=TideConst(6,1);
        N_g_MS4=TideConst(6,2);
        N_H_P1=TideConst(7,1);
        N_g_P1=TideConst(7,2);
        N_H_Q1=TideConst(8,1);
        N_g_Q1=TideConst(8,2);
        N_H_N2=TideConst(9,1);
        N_g_N2=TideConst(9,2);
        N_H_K2=TideConst(10,1);
        N_g_K2=TideConst(10,2);
        %
        N_Part_O1=f_O1.*N_H_O1.*cos(Eta_O1.*CurrentTime+v0_O1+u_O1-N_g_O1);
        N_Part_K1=f_K1.*N_H_K1.*cos(Eta_K1.*CurrentTime+v0_K1+u_K1-N_g_K1);
        N_Part_M2=f_M2.*N_H_M2.*cos(Eta_M2.*CurrentTime+v0_M2+u_M2-N_g_M2);
        N_Part_S2=f_S2.*N_H_S2.*cos(Eta_S2.*CurrentTime+v0_S2+u_S2-N_g_S2);
        N_Part_M4=f_M4.*N_H_M4.*cos(Eta_M4.*CurrentTime+v0_M4+u_M4-N_g_M4);
        N_Part_MS4=f_MS4.*N_H_MS4.*cos(Eta_MS4.*CurrentTime+v0_MS4+u_MS4-N_g_MS4);
        N_Part_P1=f_P1.*N_H_P1.*cos(Eta_P1.*CurrentTime+v0_P1+u_P1-N_g_P1);
        N_Part_Q1=f_Q1.*N_H_Q1.*cos(Eta_Q1.*CurrentTime+v0_Q1+u_Q1-N_g_Q1);
        N_Part_N2=f_N2.*N_H_N2.*cos(Eta_N2.*CurrentTime+v0_N2+u_N2-N_g_N2);
        N_Part_K2=f_K2.*N_H_K2.*cos(Eta_K2.*CurrentTime+v0_K2+u_K2-N_g_K2);
        N_Result=N_Resudial+N_Part_O1+N_Part_K1+N_Part_M2+N_Part_S2+N_Part_M4+N_Part_MS4+N_Part_P1+N_Part_Q1+N_Part_N2+N_Part_K2;
        %%
        %潮流合成
        [CurrentSpeed,CurrentDirection]=EN2VD(E_Result,N_Result);
        Result(c1,1)=datenum([StartTime(1:3),CurrentTime,0,0]);
        Result(c1,2)=CurrentSpeed;
        Result(c1,3)=CurrentDirection;
    end
end
%%
function [s,h,p,N,ps]=shpnps_predict(year,month,day,time)
% 计算天文变量
% [s,h,p,N,ps]=shpnps_predict(year,month,day,time)
ruinian=fix(0.25.*(year-1901));
if mod(year,4)==0
    switch month
        case 1,n=day-1;
        case 2,n=day+30;
        case 3,n=day+59;
        case 4,n=day+90;
        case 5,n=day+120;
        case 6,n=day+151;
        case 7,n=day+181;
        case 8,n=day+212;
        case 9,n=day+243;
        case 10,n=day+273;
        case 11,n=day+304;
        case 12,n=day+334;
        otherwise,error('Error Date');
    end
else    
    switch month
        case 1,n=day-1;
        case 2,n=day+30;
        case 3,n=day+58;
        case 4,n=day+89;
        case 5,n=day+119;
        case 6,n=day+150;
        case 7,n=day+180;
        case 8,n=day+211;
        case 9,n=day+242;
        case 10,n=day+272;
        case 11,n=day+303;
        case 12,n=day+333;
        otherwise,error('Error Date');
    end
end
% s、h'、p、n、ps
s= 277.02+129.3848.*(year-1900)+13.1764.*(n+ruinian+time./24);
h= 280.19-0.238700.*(year-1900)+0.98570.*(n+ruinian+time./24);
p= 334.39+40.66250.*(year-1900)+0.11140.*(n+ruinian+time./24);
N= 100.84+19.32820.*(year-1900)+0.05300.*(n+ruinian+time./24);
ps=281.22+0.017200.*(year-1900)+0.00005.*(n+ruinian+time./24);
%转换为弧度
s=s.*pi./180;
h=h.*pi./180;
p=p.*pi./180;
N=N.*pi./180;
ps=ps.*pi./180;
%%
function [V,D]=EN2VD(V_E,V_N)
%由流速的东/北分量计算流速流向
V=sqrt(V_E.^2+V_N.^2);
c2=find((V_E>=0)&(V_N>=0));  %第一象限内流向
D(c2)=asin(abs(V_E(c2))./V(c2)).*180./pi;
clear c2;
c2=find((V_E>=0)&(V_N<0));  %第四象限内流向
D(c2)=180-asin(abs(V_E(c2))./V(c2)).*180./pi;
clear c2;
c2=find((V_E<0)&(V_N>=0));  %第二象限内流向
D(c2)=360-asin(abs(V_E(c2))./V(c2)).*180./pi;
clear c2;
c2=find((V_E<0)&(V_N<0));  %第三象限内流向
D(c2)=180+asin(abs(V_E(c2))./V(c2)).*180./pi;
clear c2;
D=D';

⌨️ 快捷键说明

复制代码 Ctrl + C
搜索代码 Ctrl + F
全屏模式 F11
切换主题 Ctrl + Shift + D
显示快捷键 ?
增大字号 Ctrl + =
减小字号 Ctrl + -