????????本人為海洋大學大二在讀學生，由于作業需求。自行編寫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

????????朝朝夕夕，日夜不停。海洋學子在探索海洋的道路上也不會停歇，前有王品先院士，文圣常院士等開創我國海洋學的奠基。未來在維護國家海洋權利，建設海洋強國的道路上，仍然需要大量海洋及其他各個專業學生前赴后繼的沉淀。再次希望引起大家的思考和討論，逐漸完善出可以供大家學習的、標準的算法。

總結

以上是生活随笔為你收集整理的潮汐特征值计算及其预报的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。