潮汐特征值计算及其预报

现磨咖啡不放糖

已于 2022-05-11 21:01:42 修改

阅读量2k

点赞数 12

分类专栏：海洋学潮汐 matlab 文章标签： matlab 个人开发交互线性代数矩阵

于 2022-05-10 10:15:00 首次发布

本文链接：https://blog.csdn.net/qq_67298810/article/details/124565573

版权

海洋学同时被 3 个专栏收录

1 篇文章 2 订阅

订阅专栏

潮汐

1 篇文章 1 订阅

订阅专栏

matlab

1 篇文章 0 订阅

订阅专栏

本人为海洋大学大二在读学生，由于作业需求。自行编写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:7
    if basic_value(i) > 2*pi
        while basic_value(i) > 2*pi
            basic_value(i) = basic_value(i)-2*pi;
        end
    elseif basic_value(i) < (-2)*pi
        while basic_value(i) < (-2)*pi
            basic_value(i) = basic_value(i)+2*pi;
        end
    end
    if 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:8
    if all_v0(i) > 2*pi
        while all_v0(i) > 2*pi
            all_v0(i) = all_v0(i)-2*pi;
        end
    elseif all_v0(i) < (-2)*pi
        while all_v0(i) < (-2)*pi
            all_v0(i) = all_v0(i)+2*pi;
        end
    end
    if 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:9
    for j = 2:9
        A(j,j) = (1/2) * (N + (sin(N*Omega(j-1)*time_interval)) / (sin(Omega(j-1)*time_interval)) );
        if i ~= j
            A(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:9
    A(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:8
    for j = 1:8
        B(j,j) = 1/2*(N - (sin(N*Omega(j)*time_interval))/(sin(Omega(j)*time_interval)));
        if i ~= j
            B(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:9
    for 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:8
    if  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 == 1
    n = D - 1;
elseif M == 2
    n = D +30;
elseif M == 3
    n = D + 58;
elseif M == 4
    n = D + 89;
elseif M == 5
    n = D + 119;
elseif M == 6
    n = D + 150;
elseif M == 7
    n = D + 180;
elseif M == 8
    n = D + 211;
elseif M == 9
    n = D + 242;
elseif M == 10
    n = D + 272;
elseif M == 11
    n = D + 303;
elseif M == 12
    n = D + 333;
end

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]';    %返回基本天文元素矩阵

%此函数用来计算分潮的焦点因子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);
end

f = 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

朝朝夕夕，日夜不停。海洋学子在探索海洋的道路上也不会停歇，前有王品先院士，文圣常院士等开创我国海洋学的奠基。未来在维护国家海洋权利，建设海洋强国的道路上，仍然需要大量海洋及其他各个专业学生前赴后继的沉淀。再次希望引起大家的思考和讨论，逐渐完善出可以供大家学习的、标准的算法。