潮汐特征值计算及其预报
2024-06-30 10:36:22
本人为海洋大学大二在读学生，由于作业需求。自行编写matlab程序以计算潮汐类型书并对潮汐进行了预报。
我国作为陆地大国，由于各种历史问题，在过去忽视了海洋及海洋产业的发展，这也导致目前我国海洋学正处于起步阶段，网上各类海洋学相关资料也是少之又少。对于海洋学子，没有外部资料无疑增加了海洋学的学习难度。而由于本人并非计算机相关专业，并未受过相关数据结构，编程风格等计算机的专业训练，因此本人代码可能会有不规范的地方。
并且，本人计算出的结果有较大误差，由于本人能力有限，无法找到错误所在。希望本篇文章抛砖引玉，能够引起大家的思考和讨论，不断完善代码，改进算法。希望大家多加注释，未后继的海洋学生提供学习的素材。
以下为程序的主函数部分：（潮汐观测文件会在后续上传）
close all;
clear;%% 导入潮汐数据
% 第一个数对应的时间为2003年3月3日0时，相邻数据的时间间隔为1个小时。
tide_arr = importdata('（本地潮汐数据）'); %此处为导入本地潮汐数据文件%% 根据已知数据绘出潮流观测曲线
x1 = 1:721; %x为观测矩阵中第一列，既数据初始编号
y1 = tide_arr(:,2); %y为观测矩阵中第二列，既观测潮位
plot(x1,y1);
xlabel('序列号');
ylabel('潮位');
title('观测潮位与自报潮位');
hold on
%% 选取8个主要分潮O1,P1,K1,M2,K2,Q1,S2,N2
% 查表得八个主要分潮的杜德森数
% 数组中分别对应u1,u2,u3,u4,u5,u6,u0
O1_du = [1,-1,0,0,0,0,-1];
P1_du = [1,1,-2,0,0,0,-1];
K1_du = [1,1,0,0,0,0,1];
M2_du = [2,0,0,0,0,0,0];
K2_du = [2,2,0,0,0,0,0];
Q1_du = [1,-2,0,1,0,0,-1];
S2_du = [2,2,-2,0,0,0,0];
N2_du = [2,-1,0,1,0,0,0];%% 计算时间原点
% 第一个数对应的时间为2003年3月3日0时，相邻数据的时间间隔为1个小时。
year=2003;
month=3;
day=3;
hour=0;
time_interval = 1;
time_calculator(year,month,day,hour,time_interval);
%计算出中间时刻对应的时间并输出%% 计算基本天文元素随时间的变化速度
omega_t =(2*pi) / (24 + (50/60)); %平太阴时角，周期一个平太阴日=24小时50分钟
omega_s = (2*pi) / (27.32*24); %月球平均经度，周期一个回归月 = 27.32个平太阳日
omega_h = (2*pi) / (365.2422*24); %太阳平均经度，周期一个回归年 = 265.2422个平太阳日
omega_p1 = (2*pi) / (8.85*365.2422*24); %月球近地点平均经度，周期为8.85年
omega_n = (2*pi) / (18.61*365.2422*24); %白道升交点在黄道上经度的负值，周期为18.61年
omega_p2 = (2*pi) / (21000*365.2422*24); %太阳近地点平均经度，周期为21000年，u6=0此项可以忽略%基本天文分潮随时间变化率矩阵
basic_value_omega = [omega_t,omega_s,omega_h,omega_p1,omega_n,omega_p2]';%% 分别计算八个分潮的角速率
omega_O1 = O1_du(1:6) * basic_value_omega;
omega_P1 = P1_du(1:6) * basic_value_omega;
omega_K1 = K1_du(1:6) * basic_value_omega;
omega_M2 = M2_du(1:6) * basic_value_omega;
omega_K2 = K2_du(1:6) * basic_value_omega;
omega_Q1 = Q1_du(1:6) * basic_value_omega;
omega_S2 = S2_du(1:6) * basic_value_omega;
omega_N2 = N2_du(1:6) * basic_value_omega;
%创建八个分潮角速度向量
%Omega = [omega_O1,omega_P1,omega_K1,omega_M2,omega_K2,omega_Q1,omega_S2,omega_N2];% fprintf('O1分潮的角速率：%frad/h\n',omega_O1);
% fprintf('P1分潮的角速率：%frad/h\n',omega_P1);
% fprintf('K1分潮的角速率：%frad/h\n',omega_K1);
% fprintf('M2分潮的角速率：%frad/h\n',omega_M2);
% fprintf('K2分潮的角速率：%frad/h\n',omega_K2);
% fprintf('Q1分潮的角速率：%frad/h\n',omega_Q1);
% fprintf('S2分潮的角速率：%frad/h\n',omega_S2);
% fprintf('N2分潮的角速率：%frad/h\n',omega_N2);% 从这往上计算均正确代码从这以上应该没有问题，以下的初相计算和标准答案有较大偏差，本人能力有限，对知识的理解不够，未能找出错误。
%% 计算天文元素（以弧度制计算）
%初项为第一个数据收集时刻的位相：2003年3月3日0时
Y = 2003; M = 3; D = 3; t = 0;
n = D + 58; i = fix((Y - 1900)/4);
%月球平均经度s
s = (pi/180)*(277.02) + (pi/180)*(129.3848)*(Y-1900) + (pi/180)*(13.1764)*(n+i+t/24);
%太阳平均经度h
h = (pi/180)*(280.19) + (pi/180)*(0.2387)*(Y-1900) + (pi/180)*(0.9857)*(n+i+t/24);
%月球近地点平均经度p1
p1 = (pi/180)*(334.39) + (pi/180)*(40.6625)*(Y-1900) + (pi/180)*(0.1114)*(n+i+t/24);
%白道升交点在黄道上经度的负值n
n = (pi/180)*(100.84) + (pi/180)*(19.3282)*(Y-1900) + (pi/180)*(0.0530)*(n+i+t/24);
%太阳近地点平均经度p2
p2 = (pi/180)*(281.22) + (pi/180)*(0.0172)*(Y-1900) + (pi/180)*(0.00005)*(n+i+t/24);
%τ=15t?s+h′ 平太阴时tau
tau=(pi/180)*15*t-s+h;
fprintf('平太阴时角τ:%f\n月球平均经度s:%f\n太阳平均经度h:%f\n月球近地点平均经度p1:%f\n白道升交点在黄道上经度的负值n:%f\n太阳近地点平均经度p2:%f\n',tau,s,h,p1,n,p2);
basic_value = [tau,s,h,p1,n,p2,pi/2]'; %基本天文元素矩阵% 调整到[0,2*pi]
%将初始位相调整到[0，2*pi]
for i = 1:7if basic_value(i) > 2*piwhile basic_value(i) > 2*pibasic_value(i) = basic_value(i)-2*pi;endelseif basic_value(i) < (-2)*piwhile basic_value(i) < (-2)*pibasic_value(i) = basic_value(i)+2*pi;endendif basic_value(i) < 0 basic_value(i) = basic_value(i) + 2*pi;end
end% s=14934.4708; s=deg2rad(s);
% h=355.1596; h=deg2rad(h);
% p1=4533.8789; p1=deg2rad(p1);
% n=2096.9976; n=deg2rad(n);
% p2=282.99665; p2=deg2rad(p2);
% tau=(pi/180)*15*0-s+h;
% basic_value = [tau,s,h,p1,n,p2,pi/2]';
%% 计算初始位相
O1_v0 = O1_du * basic_value;
P1_v0 = P1_du * basic_value;
K1_v0 = K1_du * basic_value;
M2_v0 = M2_du * basic_value;
K2_v0 = K2_du * basic_value;
Q1_v0 = Q1_du * basic_value;
S2_v0 = S2_du * basic_value;
N2_v0 = N2_du * basic_value;
% fprintf('O1分潮初始位相：%f\n',O1_v0);
% fprintf('P1分潮初始位相：%f\n',P1_v0);
% fprintf('K1分潮初始位相：%f\n',K1_v0);
% fprintf('M2分潮初始位相：%f\n',M2_v0);
% fprintf('K2分潮初始位相：%f\n',K2_v0);
% fprintf('Q1分潮初始位相：%f\n',Q1_v0);
% fprintf('S2分潮初始位相：%f\n',S2_v0);
% fprintf('N2分潮初始位相：%f\n',N2_v0);%建立初项的矩阵
all_v0 = [O1_v0,P1_v0,K1_v0,M2_v0,K2_v0,Q1_v0,S2_v0,N2_v0];%将初始位相调整到[0，2*pi]
for i = 1:8if all_v0(i) > 2*piwhile all_v0(i) > 2*piall_v0(i) = all_v0(i)-2*pi;endelseif all_v0(i) < (-2)*piwhile all_v0(i) < (-2)*piall_v0(i) = all_v0(i)+2*pi;endendif all_v0(i) < 0 all_v0(i) = all_v0(i) + 2*pi;end
end%all_v0 = [0.0240 0.0000 0.7119 6.1142 1.4863 4.8209 4.7969 5.5087];%% 计算焦点订正角
%基本方法计算方法
%K1分潮的焦点因子f_K1和焦点订正角u_K1
a = [p1,n,p2,pi/2]';
K1_fcosu = 0.0002*cos([-2,-1,0,0]*a) + 0.0001*cos([0,-2,0,0]*a)....+0.0198*cos([0,-1,0,-2]*a) + 1*cos([0,0,0,0]*a)...+0.1356*cos([0,1,0,0]*a) + 0.0029*cos([0,2,0,-2]*a);K1_fsinu = 0.0002*sin([-2,-1,0,0]*a) + 0.0001*sin([0,-2,0,0]*a)....+0.0198*sin([0,-1,0,-2]*a) + 1*sin([0,0,0,0]*a)...+0.1356*sin([0,1,0,0]*a) + 0.0029*sin([0,2,0,-2]*a);f_K1 = sqrt( (K1_fcosu)^2 + (K1_fsinu)^2 );
u_K1 = atan(K1_fsinu/K1_fcosu);
fprintf('K1分潮的焦点因子:%f\n和焦点订正角:%f\n',f_K1,u_K1);
%除K1分潮其余分潮通过自定义函数f_and_u_calculator计算%% 用自定义函数计算其余分潮的焦点因子f和焦点订正角u
%O1分潮
p_O1 = [-0.0058,0.1885,1,0.0002,-0.0064,-0.0010]; %O1分潮引潮力系数比
delta_u_01 = [0,-2,0,0;0,-1,0,0;0,0,0,0;2,-1,0,0;2,0,0,0;2,1,0,0]; %O1分潮次要分潮与主要分潮的差值Δμ
fprintf('O1分潮的');
[f_O1,u_O1] = f_and_u_calculator(p_O1,delta_u_01);%P1分潮
p_P1 = [0.0008,-0.0112,1,-0.0015,-0.0003]; %P1分潮引潮力系数比
delta_u_P1 = [0,-2,0,0;0,-1,0,0;0,0,0,0;2,0,0,0;2,1,0,0]; %P1分潮次要分潮与主要分潮的差值Δμ
fprintf('P1分潮的');
[f_P1,u_P1]= f_and_u_calculator(p_P1,delta_u_P1);%M2分潮
p_M2 = [0.0005,-0.0373,1,0.0006,0.0002]; %M2分潮引潮力系数比
delta_u_M2 = [0,-2,0,0;0,-1,0,0;0,0,0,0;2,0,0,0;2,1,0,0]; %M2分潮次要分潮与主要分潮的差值Δμ
fprintf('M2分潮的');
[f_M2,u_M2] = f_and_u_calculator(p_M2,delta_u_M2);%K2分潮
p_K2 = [-0.0128,1,0.2980,0.0324]; %K2分潮引潮力系数比
delta_u_K2 = [0,-1,0,0;0,0,0,0;0,1,0,0;0,2,0,0]; %K2分潮次要分潮与主要分潮的差值Δμ
fprintf('K2分潮的');
[f_K2,u_K2] = f_and_u_calculator(p_K2,delta_u_K2);%% Q1,S2,N2三个分潮通过查表得到
%Q1分潮
%Q1分潮的焦点因子和焦点订正角都和O1分潮相同
fprintf('Q1分潮的');
[f_Q1,u_Q1] = f_and_u_calculator(p_O1,delta_u_01); %与O1分潮的相同%S2分潮
%S2分潮的焦点因子为1，焦点订正角为0
f_S2 = 1;
u_S2 = 0;
fprintf('S2分潮的焦点因子为：1\n焦点订正角为：0\n');%N2分潮
%N2分潮的焦点因子和焦点订正角都和M2分潮相同
fprintf('N2分潮的');
[f_N2,u_N2] = f_and_u_calculator(p_M2,delta_u_M2); %与O1分潮的相同% 焦点因子矩阵
all_f = [f_O1,f_P1,f_K1,f_M2,f_K2,f_Q1,f_S2,f_N2];
%焦点订正角矩阵
all_u = [u_O1,u_P1,u_K1,u_M2,u_K2,u_Q1,u_S2,u_N2];
%% 联立矛盾方程，利用最小二乘法得到法方程并求出矛盾方程最优解
%一共选取8个分潮
%一共有721个观测记录数据，那么就有721个方程构成矛盾方程，观测水位h有721个
%每个分潮的振幅和迟角为调和常数，一单确定则不随时间变化，平均水位也未知
%因此未知量的个数有2×8+1=17个（8个分潮的Hsinθ和Hcosθ以及S0）
%直接求法方程的系数矩阵A，B，F′，F"%创建八个分潮角速度向量
N = 721 ;
Omega = [omega_O1,omega_P1,omega_K1,omega_M2,omega_K2,omega_Q1,omega_S2,omega_N2];
%此步需要注意后面分潮计算的顺序都要按照上面的顺序%% 下面的循环计算系数矩阵A
A = zeros(8,8); %初始化
%A矩阵的第一行第一个
A(1,1) = N;
%A矩阵除去第一行和第一列
for i = 2:9for j = 2:9A(j,j) = (1/2) * (N + (sin(N*Omega(j-1)*time_interval)) / (sin(Omega(j-1)*time_interval)) );if i ~= jA(i,j) = (1/2) * ((sin((N/2)*(Omega(i-1)-Omega(j-1))*time_interval)) / ...(sin((1/2)*(Omega(i-1)-Omega(j-1))*time_interval))...+ (sin((N/2)*(Omega(i-1)+Omega(j-1))*time_interval)) /...(sin((1/2)*(Omega(i-1)+Omega(j-1))*time_interval)));end
% A(j,i) = A(i,j);end
end
%计算A矩阵的第一行第二个元素到最后一个元素,以及第一列第二个元素到最后一个元素
for j = 2:9A(1,j) = (sin((N/2)*Omega(j-1)*time_interval)) / (sin((1/2)*Omega(j-1)*time_interval));A(j,1) = A(1,j);
end
%disp(A)
%% 下面的循环计算系数矩阵B
B = zeros(8,8);
for i = 1:8for j = 1:8B(j,j) = 1/2*(N - (sin(N*Omega(j)*time_interval))/(sin(Omega(j)*time_interval)));if i ~= jB(i,j) = (1/2) * ((sin((N/2)*(Omega(i)-Omega(j))*time_interval)) / ...(sin((1/2)*(Omega(i)-Omega(j))*time_interval))...- (sin((N/2)*(Omega(i)+Omega(j))*time_interval)) /...(sin((1/2)*(Omega(i)+Omega(j))*time_interval)));end
% B(j,i) = B(i,j);end
end
%disp(B)%% 下面计算F'和F"，F'在下面用向量F1表示，F"用向量F2表示
tide_arr(:,3)=[-360:0,1:360];
F1 = zeros(1,9)';
F2 = zeros(1,8)';
F1(1) = sum(tide_arr(:,2));
for i = 2:9for n = tide_arr(1,3):tide_arr(end,3)F1(i) = F1(i) + (tide_arr(n+361,2)) * cos(n*Omega(i-1)*time_interval);F2(i-1) = F2(i-1) + (tide_arr(n+361,2))*sin(n*Omega(i-1)*time_interval);end
end%计算未知量X=Hcosθ,Y=Hsinθ
X = F1'/A;
Y = F2'/B;
S0 = X(1);
%准调和分潮振幅为
R = sqrt( power(X(:,2:end),2) + power(Y,2) );%t=0时刻的位相
theta = asin(Y ./ R);
for i = 1:8if theta(:,i) < 0;theta(:,i) = pi - theta(:,i);end
end% 振幅和位相对应的八个分潮分别为
% Omega = [omega_O1,omega_P1,omega_K1,omega_M2,omega_K2,omega_Q1,omega_S2,omega_N2];
% O1,P1,K1,M2,K2,Q1,S2,N2
%% 计算调和常数
% 振幅
H = R ./ all_f;
% 迟角
g = theta + all_v0 + all_u;%% 潮汐自报
h = zeros(721,2);
for n = tide_arr(:,3)for i = 1:8 h(n+361) = h(n+361) + all_f(i)*H(i)*cos(n*Omega(i)*time_interval + all_v0(i) + all_u(i) - g(i) );end
end
h = h(:,1) + S0; %加上平均潮位
h(:,2) = 1:721; %给自报潮位编号%自报潮位作图
x2 = h(:,2);
y2 = h(:,1);
plot(x2,y2);
legend('观测','自报');%% 自报余差计算
%余差r
r = tide_arr(:,2) - h(:,1);
%自报余差的标准差delta_r_2 = (1/721) * sum(power(r,2)); %方差delta_r = sqrt(delta_r_2); %标准差text(0,190,'自报余差的方差为：185.5727 标准差为：13.6225','Color','green','FontSize',14);%% 预报潮位（5月1日0时到6月1日0时）
%已有数据第一个数对应的时间为2003年3月3日0时，相邻数据的时间间隔为1个小时。
%721个数据，从第一个观测记录数到最后一个观测记录数一共经历了720个小时
%所以一共经历了30个平太阳日
%已有的观测记录数从2003年3月3日0时到2003年4月3日0时h_forecast = zeros(721,2);% 潮汐的振幅H和迟角g为调和常量，是本身的动力学要素，不随时间改变
% S0一但确定也不会改变
% 焦点因子f和焦点订正角u在1年之内变化缓慢，近似看作不变%计算新的天文元素
%计算待预报潮位的中间时刻,一共也会有721个数据
Y_forecast_middle = 2003 ;
M_forecast_middle = 5;
t_total = 360;
D_forecast_middle = fix(t_total/24);
t_forecast_middle = rem(t_total,24);
basoc_valie_forecast = zeros(7,1);
basic_value_forecast = basic_value_calculator(Y_forecast_middle,M_forecast_middle,D_forecast_middle,t_forecast_middle);%% 计算新的初相位
O1_v0_forecast = O1_du * basic_value_forecast;
P1_v0_forecast = P1_du * basic_value_forecast;
K1_v0_forecast = K1_du * basic_value_forecast;
M2_v0_forecast = M2_du * basic_value_forecast;
K2_v0_forecast = K2_du * basic_value_forecast;
Q1_v0_forecast = Q1_du * basic_value_forecast;
S2_v0_forecast = S2_du * basic_value_forecast;
N2_v0_forecast = N2_du * basic_value_forecast;
%预报的初相矩阵
all_v0_forecast = [O1_v0_forecast,P1_v0_forecast,K1_v0_forecast,M2_v0_forecast,K2_v0_forecast,Q1_v0_forecast,S2_v0_forecast,N2_v0_forecast];for n = tide_arr(:,3)for i = 1:8 h_forecast(n+361) = h_forecast(n+361) + ...all_f(i)*H(i)*cos(n*Omega(i)*time_interval + all_v0_forecast(i) + all_u(i) - g(i) );end
end
h_forecast = h_forecast(:,1) + S0; %加上平均潮位
h_forecast(:,2)=1:721;%% 画出预报的潮位
figure(2);
x3 = h_forecast(:,2);
y3 = h_forecast(:,1);
plot(x3,y3);
xlabel('预报记录数');
ylabel('潮位');
title('预报潮位（5月1日0时到6月1日0时）');以下部分为自编函数部分，主要分为三个：
%这个函数用来计算天文元素（以弧度制计算）并返回基本天文元素矩阵
%这个函数未考虑平年和闰年的情况
%记1月1日为第0天function basic_value = basic_value_calculator(Y,M,D,t)%中间时刻的世纪时为：2003年3月18日0时
% Y = 2003; M = 3; D = 18; t = 0;
if M == 1n = D - 1;
elseif M == 2n = D +30;
elseif M == 3n = D + 58;
elseif M == 4n = D + 89;
elseif M == 5n = D + 119;
elseif M == 6n = D + 150;
elseif M == 7n = D + 180;
elseif M == 8n = D + 211;
elseif M == 9n = D + 242;
elseif M == 10n = D + 272;
elseif M == 11n = D + 303;
elseif M == 12n = D + 333;
endi = fix((Y - 1900)/4);
%月球平均经度s
s = (pi/180)*(277.02) + (pi/180)*(129.3848)*(Y-1900) + (pi/180)*(13.1764)*(n+i+t/24);
%太阳平均经度h
h = (pi/180)*(280.19) + (pi/180)*(0.2387)*(Y-1900) + (pi/180)*(0.9857)*(n+i+t/24);
%月球近地点平均经度p1
p1 = (pi/180)*(334.39) + (pi/180)*(40.6625)*(Y-1900) + (pi/180)*(0.1114)*(n+i+t/24);
%白道升交点在黄道上经度的负值n
n = (pi/180)*(100.84) + (pi/180)*(19.3282)*(Y-1900) + (pi/180)*(0.0530)*(n+i+t/24);
%太阳近地点平均经度p2
p2 = (pi/180)*(281.22) + (pi/180)*(0.0172)*(Y-1900) + (pi/180)*(0.00005)*(n+i+t/24);
%τ=15t?s+h′ 平太阴时tau
tau=(pi/180)*15*t-s+h;
%fprintf('平太阴时角τ:%f\n月球平均经度s:%f\n太阳平均经度h:%f\n月球近地点平均经度p1:%f\n白道升交点在黄道上经度的负值n:%f\n太阳近地点平均经度p2:%f\n',tau,s,h,p1,n,p2);
basic_value = [tau,s,h,p1,n,p2,pi/2]'; %返回基本天文元素矩阵%此函数用来计算分潮的焦点因子f和焦点订正角u
%函数使用时注意：引潮力系数比和杜德森数之差要一一对应
%ρ前带正负号时，既不考虑μ0时，要将μ0位置用0代替。function [a,b] = f_and_u_calculator(p,delta_u)
%设置参数：引潮力系数比为p， a = [p1,n,p2,pi/2]'为一个定值
%delta_u为一个矩阵，是同一亚群其他分潮与主要分潮u4,u5,u6,u0之差，其中u6恒为0
%写delta_u矩阵时要注意每行与引潮力系数比p的分潮要一一对应%% 由于基本天文元素不变，因此将矩阵a设置为定值
Y = 2003; D = 3; t = 0;
n = D + 58; i = fix((Y - 1900)/4);
p1 = (pi/180)*(334.39) + (pi/180)*(40.6625)*(Y-1900) + (pi/180)*(0.1114)*(n+i+t/24);
%白道升交点在黄道上经度的负值n
n = (pi/180)*(100.84) + (pi/180)*(19.3282)*(Y-1900) + (pi/180)*(0.0530)*(n+i+t/24);
%太阳近地点平均经度p2
p2 = (pi/180)*(281.22) + (pi/180)*(0.0172)*(Y-1900) + (pi/180)*(0.00005)*(n+i+t/24);
a = [p1,n,p2,pi/2]';%% 焦点因子和焦点订正角计算的主要步骤
fcosu = 0;
fsinu = 0;
for j = 1:length(p)fcosu = fcosu + p(j) * cos(delta_u(j,:)*a);fsinu = fsinu + p(j) * sin(delta_u(j,:)*a);
endf = sqrt( (fcosu^2)+(fsinu^2) );
u = atan( (fsinu)/(fcosu) );% u=acot(fcosu/fsinu);
% f=fcosu/(cos(u));
% u=u+2*pi;
% fprintf('焦点因子为：%f\n焦点订正角为：%f\n',f,u);
a=f;
b=u;
end
%这个函数用来计算分潮的中间时刻以及其中间时刻的真实时间
%中间时刻N' = 第一个观测记录数 + 1/2 *(N-1)%目前任务的第一个数对应的时间为2003年3月3日0时，相邻数据的时间间隔为1个小时。
%总共有30天的数据function time_calculator(year,month,day,hour,time_interval)
tide_arr = importdata('D:\SHOU\Matlab\bin\tide_work\21.txt');
[N,~] = size(tide_arr); %N为总的观测记录数
N_pie = 1+(1/2) * (N-1); %N_pie为中间时刻的编号,第一个观测记录数为1
fprintf('中间时刻潮汐编号为：%d\n',N_pie); hours_of_N_pie = (N_pie) * (time_interval) - 1 ; %总计经过的小时数
days_of_N_pie = fix(hours_of_N_pie/24); %总天数
real_day_of_N_pie = day + days_of_N_pie; %中间时刻时间不会到四月
real_hour_of_N_pie = hour + rem(hours_of_N_pie , 24);
fprintf('中间时刻的世纪时为：');
fprintf('%d年%d月%d日%d时\n',year,month,real_day_of_N_pie,real_hour_of_N_pie);
end
朝朝夕夕，日夜不停。海洋学子在探索海洋的道路上也不会停歇，前有王品先院士，文圣常院士等开创我国海洋学的奠基。未来在维护国家海洋权利，建设海洋强国的道路上，仍然需要大量海洋及其他各个专业学生前赴后继的沉淀。再次希望引起大家的思考和讨论，逐渐完善出可以供大家学习的、标准的算法。
潮汐特征值计算及其预报相关推荐
- 【数理知识】《数值分析》李庆扬老师-第8章-矩阵特征值计算
第7章 回到目录 第9章 第8章-矩阵特征值计算 8.1 特征值性质和估计 8.1.1 特征值问题及其性质 定理1 8.2 幂法及反幂法 8.3.1 豪斯霍尔德 (Householder) 变换 8. ...
- 潮汐观测数据调和分析及预报成图
首先,需要说明两点. 第一,本案例使用的是单月潮汐观测数据,处理方法则是基于长期观测资料的调和分析来进行处理.中期观测资料的分析需要分别求主要分潮.随从分潮,短期观测资料分析则还需要计算不同观测序列的 ...
- GPU上大规模稀疏矩阵特征值计算高效算法之三——SLEPc测试
Slepc计算矩阵特征值时间测试 注: (1)GPU集群介绍: 该集群有一个登录节点(ustcgpu)和100个计算节点(node1~node100).各计算节点配置2 颗4核的IntelE5520 ...
- hermite矩阵matlab,Hermite矩阵的特征值计算问题 - 数学 - 小木虫 - 学术 科研 互动社区...
楼主,首先要知道Hermite矩阵就是实对称矩阵的复数版本.我更愿意用:A^H=A来定义Hermite矩阵,这里A是任何一个n阶的复数矩阵,H表示共轭转置.你的这个问题有三个步骤,但是我不清楚你的意图 ...
- matlab计算潮差程序,一种基于FVCOM模型的可视化潮汐潮流预报方法与流程
本发明涉及潮汐预报技术领域,特别是涉及一种基于FVCOM模型的可视化潮汐潮流预报方法. 背景技术: 潮汐预报对一定海区在未来一定时间内的潮汐涨落情况进行的推算和预报.预报内容包括逐日的高潮和低潮高度及 ...
- C语言通过QR分解计算矩阵的特征值和特征向量
#include <stdio.h> #include <stdlib.h> #include <math.h> #include <time.h>// ...
- MATLAB特征值的计算之eig()函数存在的问题
本人这段时间一直在研究特征值计算问题,当然从理论上来说,这个问题很简单.甚至,我们自己可以通过公式,用C语言和FORTRAN语言编制相应的代码来实现,或者用PYTHON也行.结果,我选择了编程效率很高 ...
- 潮汐计算php源码,科学网—潮汐再科普与再计算 - 岳东晓的博文
我之前写过几篇潮汐的科普,并且计算了潮汐高度随角度的变化.但这个计算忽略了一项.因此,我重新计算一下.顺便再从概念上讲讲. 古代中国人就知道潮汐主要是由月球引起的.牛顿在发现其运动定律及万有引力之后, ...
- 科学计算发展简史 -- 信息与计算科学
来源 | 本文节选自"创新报国70年"大型报告文学丛书之<冯康传>引子部分,浙江教育出版社,2019年12月. 数学与计算数学 在几千年数学演变的历史长河中,东方和西方 ...
最新文章
- 跟各种诡异 Bug 打交道 13 年后的总结
- 找工作 50道编程题Java实现（32-50）
- h5精准定位_HTML5 地理定位
- java面试题十 java数组初始化
- python时间去掉t_Python的set集合详解
- 结构体定义的三钟方式
- 互联网实习笔记之30天总结
- linux内核 快速分片,linux内核学习笔记------ip报文的分片
- SQL那些事儿（八）--oracle用户、表、表空间之间的关系
- 最新 IDEA 2022.1 版本即将发布，骚操作真不少...
- app自动化之移动端测试基础知识
- 2022年武汉科技大学成人高等学历教育招生简章--学历提升、高起专、专升本
- iar stm32_基于最新5.4电机库的STM32电机控制应用实战分享
- 时域有限差分法matlab程序,时域有限差分法的Matlab仿真
- v6使用手册 天正电气t20_天正电气T20手册
- f1c100s开发笔记
- 319. 灯泡开关【我亦无他唯手熟尔】
- CSS盒子模型（内容区、边框、内外边距）
- c语言输出一行星星代码,C语言打印星星的问题
- 基于恩智浦MK60DN512Z系列单片机的智能模型车主程序与子程序集