clc;
clear all;
close all;
datas=textread('2.txt');
data=datas(:,1:2);
%%時間計算程序： 第一個數對應的時間為2003年3月3日0時 見程序time
%%%%%%Ans： 361= 2003/3/18/0

%%%%%%%%%%%%
%先用3月時間序列求f u vo 再結合R sita ?獲取定值H g
%再用五月的時間序列求新的f u vo 代入預報公式
rad=pi/180;%弧度制轉換
Y=2003;
M=3;
D=18;?
t=0;
k=D+58;%n為1.1日起的日期數

%計算六個天文元素（PPT第9課時第7、8頁 /上機-第3頁）
i=floor((Y-1900)/4);%i為1900年至Y年的閏年數
nit=(k+i+t/24);
YY=Y-1900;
s=277.02+129.3848*YY+13.1764*nit;
h_pie=280.19-0.2387*YY+0.9857*nit;
p=334.39+40.6625*YY+0.1114*nit;
N_pie=100.84+19.3282*YY+0.0530*nit;
p_pie=281.22+0.0172*YY+0.00005*nit;
tao=15*t-s+h_pie;
Astro=[tao,s,h_pie,p,N_pie,p_pie,pi/2]*rad;

%%%%%% 角速度驗證//選取8分潮并通過6個天文元素計算其角速度
%6個天文元素隨時間的變化速率如下：#書本p34
tao_dt=14.49205211*rad;
s_dt=0.54901653*rad;
h_pie_dt=0.04106864*rad;
p_dt=0.00464183*rad;
N_pie_dt=0.00220641*rad;
p_pie_dt=0.00000196*rad;

%%%M2 由u1至u6-u0 下同
u_M2=[2,0,0,0,0,0,0];
sig_M2=tao_dt*u_M2(1)+u_M2(2)*s_dt+u_M2(3)*h_pie_dt+u_M2(4)*p_dt+u_M2(5)*N_pie_dt+u_M2(6)*p_pie_dt;
? ? ?%-得sig_M2=0.5059 ?查表得M2分潮角速度為：28.9841042*pi/180=0.5059驗算成立
? ? ? %-則接下來為了節約時間，分潮角速度直接查表得出?
? ? ? %2.0記錄 太坑了Q1和表給的不一樣需要重新算
v0_M2=sum(Astro.*u_M2); %M2分潮的初始相位（上機PPT第3頁）
for i=1:10000
?if( v0_M2 <= 0 )
? ? v0_M2=v0_M2+2*pi;
?elseif(v0_M2>= 2*pi)
? ? v0_M2=v0_M2-2*pi;
% ?elseif( v0_M2>=0 && v0_M2<= 2*pi)
% ? ? ?return
?end
end
%計算M2分潮的交點因子f &交點訂正角uu（上機PPT第4頁）：
rou_M2=[0.0005,-0.0373,1,0.0006,0.0002];%書本92
du4_M2=[0,0,0,2,2];
du5_M2=[-2,-1,0,0,1];
f_M2_cosu=sum(rou_M2.*cos(du4_M2*p + du5_M2*N_pie));%點乘求和
f_M2_sinu=sum(rou_M2.*sin(du4_M2*p + du5_M2*N_pie));
f_M2=((f_M2_cosu^2)+(f_M2_sinu^2))^(1/2);
uu_M2=atan(f_M2_sinu/f_M2_cosu);
for i=1:10000
?if( uu_M2 <= 0 )
? ? uu_M2=uu_M2+2*pi;
?elseif(uu_M2>= 2*pi)
? ? uu_M2=uu_M2-2*pi;
% ?elseif( uu_M2>=0 && uu_M2<= 2*pi)
% ? ? ?return
?end
end

%%%S2?
u_S2=[2,2,-2,0,0,0,0];
sig_S2=30*rad;
v0_S2=sum(Astro.*u_S2);
for i=1:10000
?if( v0_S2 <= 0 )
? ? v0_S2=v0_S2+2*pi;
?elseif(v0_M2>= 2*pi)
? ? v0_S2=v0_S2-2*pi;
% ?elseif( v0_S2>=0 && v0_S2<= 2*pi)
% ? ? ?return
?end
end
%交點因子1 訂正角0
f_S2=1.0;
uu_S2=0;

%%%N2?
u_N2=[2,-1,0,1,0,0,0];
sig_N2=28.43973*rad;
v0_N2=sum(Astro.*u_N2);
for i=1:10000
?if( v0_N2 <= 0 )
? ? v0_N2=v0_N2+2*pi;
?elseif(v0_M2>= 2*pi)
? ? v0_N2=v0_N2-2*pi;
% ?elseif( v0_N2>=0 && v0_N2<= 2*pi)
% ? ? ?return
?end
end
%M2計算結果賦值給其:
f_N2=f_M2;
uu_N2=uu_M2;

%%%K2
u_K2=[2,2,0,0,0,0,0];
sig_K2=30.082137*rad;
v0_K2=sum(Astro.*u_K2);
for i=1:10000
?if( v0_K2 <= 0 )
? ? v0_K2=v0_K2+2*pi;
?elseif(v0_K2>= 2*pi)
? ? v0_K2=v0_K2-2*pi;
% ?elseif( v0_K2>=0 && v0_K2<= 2*pi)
% ? ? ?return
?end
end
%計算K2分潮的交點因子f &交點訂正角uu(書本83 & 92）
rou_K2=[-0.0128,1,0.2980,0.0324];
du4_K2=[0,0,0,0];
du5_K2=[-1,0,1,2];
f_K2_cosu=sum(rou_K2.*cos(du4_K2*p + du5_K2*N_pie));
f_K2_sinu=sum(rou_K2.*sin(du4_K2*p + du5_K2*N_pie));
f_K2=((f_K2_cosu^2)+(f_K2_sinu^2))^(1/2);
uu_K2=atan(f_K2_sinu/f_K2_cosu);
for i=1:10000
?if( uu_K2 <= 0 )
? ? uu_K2=uu_K2+2*pi;
?elseif(uu_K2>= 2*pi)
? ? uu_K2=uu_K2-2*pi;
% ?elseif( uu_K2>=0 && uu_K2<= 2*pi)
% ? ? ?return
?end
end

%%%K1?
u_K1=[1,1,0,0,0,0,1];
sig_K1=15.0410686*rad;
v0_K1=sum(Astro.*u_K1);
for i=1:10000
?if( v0_K1 <= 0 )
? ? v0_K1=v0_K1+2*pi;
?elseif(v0_K1>= 2*pi)
? ? v0_K1=v0_K1-2*pi;
% ?elseif( v0_K1>=0 && v0_K1<= 2*pi)
% ? ? ?return
?end
end

%計算K1分潮的交點因子f &交點訂正角uu(書本83 & 92）
rou_K1=[0.0002,0.0001,-0.0198,1,0.1356,-0.0029];
du4_K1=[-2,0,0,0,0,0];
du5_K1=[-1,-2,-1,0,1,2];
f_K1_cosu=sum(rou_K1.*cos(du4_K1*p + du5_K1*N_pie));
f_K1_sinu=sum(rou_K1.*sin(du4_K1*p + du5_K1*N_pie));
f_K1=((f_K1_cosu^2)+(f_K1_sinu^2))^(1/2);
uu_K1=atan(f_K1_sinu/f_K1_cosu);
for i=1:10000
?if( uu_K1 <= 0 )
? ? uu_K1=uu_K1+2*pi;
?elseif(uu_K1>= 2*pi)
? ? uu_K1=uu_K1-2*pi;
% ?elseif( uu_K1>=0 && uu_K1<= 2*pi)
% ? ? ?return
?end
end

%%%O1?
u_O1=[1,-1,0,0,0,0,-1];
sig_O1=13.94303559*rad;
v0_O1=sum(Astro.*u_O1);
for i=1:10000
?if( v0_O1 <= 0 )
? ? v0_O1=v0_O1+2*pi;
?elseif(v0_O1>= 2*pi)
? ? v0_O1=v0_O1-2*pi;
% ?elseif( v0_O1>=0 && v0_O1<= 2*pi)
% ? ? ?return
?end
end
%計算O1分潮的交點因子f &交點訂正角uu(書本92）
rou_O1=[-0.0058,0.1885,1,0.0002,-0.0064,-0.0010];
du4_O1=[0,0,0,2,2,2];
du5_O1=[-2,-1,0,-1,0,1];
f_O1_cosu=sum(rou_O1.*cos(du4_O1*p + du5_O1*N_pie));
f_O1_sinu=sum(rou_O1.*sin(du4_O1*p + du5_O1*N_pie));
f_O1=((f_O1_cosu^2)+(f_O1_sinu^2))^(1/2);
uu_O1=atan(f_O1_sinu/f_O1_cosu);
for i=1:10000
?if( uu_O1 <= 0 )
? ? uu_O1=uu_O1+2*pi;
?elseif(uu_O1>= 2*pi)
? ? uu_O1=uu_O1-2*pi;
% ?elseif( uu_O1>=0 && uu_O1<= 2*pi)
% ? ? ?return
?end
end

%%%P1?
u_P1=[1,1,-2,0,0,0,-1];
sig_P1=14.958931*rad;
v0_P1=sum(Astro.*u_P1);
for i=1:10000
?if( v0_P1 <= 0 )
? ? v0_P1=v0_P1+2*pi;
?elseif(v0_P1>= 2*pi)
? ? v0_P1=v0_P1-2*pi;
% ?elseif( v0_P1>=0 && v0_P1<= 2*pi)
% ? ? ?return
?end
end
%計算P1分潮的交點因子f &交點訂正角uu(書本83 & 92）
rou_P1=[0.0008,-0.0112,1,-0.0015,-0.0003];
du4_P1=[0,0,0,2,2];
du5_P1=[-2,-1,0,0,1];
f_P1_cosu=sum(rou_P1.*cos(du4_P1*p + du5_P1*N_pie));
f_P1_sinu=sum(rou_P1.*sin(du4_P1*p + du5_P1*N_pie));
f_P1=((f_P1_cosu^2)+(f_P1_sinu^2))^(1/2);
uu_P1=atan(f_P1_sinu/f_P1_cosu);
for i=1:10000
?if( uu_P1 <= 0 )
? ? uu_P1=uu_P1+2*pi;
?elseif(uu_P1>= 2*pi)
? ? uu_P1=uu_P1-2*pi;
% ?elseif( uu_P1>=0 && uu_P1<= 2*pi)
% ? ? ?return
?end
end

%%%Q1
u_Q1=[1,-2,0,1,0,0,-1];
sig_Q1=tao_dt*u_Q1(1)+u_Q1(2)*s_dt+u_Q1(3)*h_pie_dt+u_Q1(4)*p_dt+u_Q1(5)*N_pie_dt+u_Q1(6)*p_pie_dt;
sig_Q11=13.55365176*rad;
v0_Q1=sum(Astro.*u_Q1);
for i=1:10000
?if( v0_Q1 <= 0 )
? ? v0_Q1=v0_Q1+2*pi;
?elseif(v0_Q1>= 2*pi)
? ? v0_Q1=v0_Q1-2*pi;
% ?elseif( v0_Q1>=0 && v0_Q1<= 2*pi)
% ? ? ?return
?end
end
%O1計算結果賦值給其:
f_Q1=f_O1;
uu_Q1=uu_O1;

N=721;
n=[-360:1:360];%定義新時間序列，第361個數為0
h=data(:,2);%數據2.txt的第二列潮高數據
sig=[sig_M2,sig_S2,sig_N2,sig_K2,sig_K1,sig_O1,sig_P1,sig_Q1]; % 8個分潮的角速度排序
%size(sig) %看下K=8
dt=1.0;%每隔一小時一個觀測記錄數
%%%%構建X、Y的等號右側矩陣（F’、F''）與系數矩陣A、B
? %%X、Y的等號右側矩陣（F’、F''）：%詳見書本P73 （3.42d）
F0_pie=sum(h);
F_1pie=zeros(1,8);
F_2pie=zeros(1,8);
for i=1:8
? ? Fi_pie=sum(h.*cos(n'*sig(i)*dt));
? ? F_1pie(i)=Fi_pie?
end
for i=1:8
? ? Fi_2pie=sum(h.*sin(n'*sig(i)*dt));
? ? F_2pie(i)=Fi_2pie ?%得到Y的右矩陣 F''
end
F_pie=[F0_pie,F_1pie]; %得到X的右矩陣 F'
size(F_pie)
size(F_2pie)
? %%X、Y的系數矩陣A、B：%詳見書本P72 （3.42a-b）
%A
A=zeros(9,9);
A(1,1)=N; ?%即A_00=N
for i=2:9 ?%算A_0i 和 A_i0
? ? A_0i=sin((N/2)*sig(i-1)*dt)/sin((1/2)*sig(i-1)*dt);
? ? A(1,i)=A_0i
? ? A(i,1)=A_0i
end
for i=2:9 ?%算對角線A_ii
? ? A_ii=(1/2)*(N+(sin(N*sig(i-1)*dt)/sin(sig(i-1)*dt)));
? ? A(i,i)=A_ii
end

for i=2:9 ?%算A_ij & A_ji
? ? for j=2:8
? ? ? ? if(i==j)
? ? ? ? ? ? continue
? ? ? ? elseif(i<j)
? ? ? ? ? ? continue ? ? ? ? ? ?
? ? ? ? elseif(i>j)
? ? ? ? ? ? A_ij=(1/2)*(((sin((N/2)*(sig(i-1)-sig(j-1)))*dt)/(sin((1/2)*(sig(i-1)-sig(j-1)))*dt))+...
? ? ? ? ? ? ?((sin((N/2)*(sig(i-1)+sig(j-1)))*dt)/(sin((1/2)*(sig(i-1)+sig(j-1)))*dt))); %書上要滿足i>j，待問
? ? ? ? ? ? A(i,j)=A_ij %那循環中的第一個應該是A_32 接下來是A_42、A_52....A43【即左下三角】
? ? ? ? ? ? A(j,i)=A(i,j) %A_23、A_24...
? ? ? ? end
? ? end
end

%B
B=zeros(8,8);
for i=1:8 ?%算對角線B_ii
? ? B_ii=(1/2)*(N-(sin(N*sig(i)*dt)/sin(sig(i)*dt)));
? ? B(i,i)=B_ii
end
for i=1:8 ?%算B_ij & B_ji
? ? for j=1:7
? ? ? ? if(i==j)
? ? ? ? ? ? continue
? ? ? ? elseif(i<j)
? ? ? ? ? ? continue ? ? ? ? ? ?
? ? ? ? elseif(i>j)
? ? ? ? ? ? B_ij=(1/2)*(((sin((N/2)*(sig(i)-sig(j)))*dt)/(sin((1/2)*(sig(i)-sig(j)))*dt))-...
? ? ? ? ? ? ?((sin((N/2)*(sig(i)+sig(j)))*dt)/(sin((1/2)*(sig(i)+sig(j)))*dt))); %書上要滿足i>j，待問
? ? ? ? ? ? B(i,j)=B_ij?
? ? ? ? ? ? B(j,i)=B(i,j)?
? ? ? ? end
? ? end
end
%%%%%%%%%%接下來終于到了接矩陣方程這一步了嗚嗚嗚
X=A\F_pie' ?%x=Hcosθ
? ? ? ? ? ? % X=F_pie/A' 和老師的右除答案一樣的
Y=B\F_2pie'; %y=Hsinθ
S0=X(1,:) ? %平均水位S0 ? %%%%%%%X的矩陣 ?第一個變量就是S0 ?然后才是X1 X2...
%%%%%%%%%%計算分潮的準調和振幅 R 和位相sita （詳見上機PPT第7頁）
R=(X(2:9,:).^2+Y.^2).^(1/2) %求8個分潮對應振幅H
% sita=atan(Y./X(2:9,:)) % ！是否是弧度制？
sita=zeros(8,1);
for ii=1:8
? ? if(X(ii+1)>=0)
? ? ? ? sita(ii)=asin(Y(ii)./R(ii))
? ? elseif(X(ii+1)<0)
? ? ? ? sita(ii)=pi-asin(Y(ii)./R(ii)) %求得8個分潮對應位相sita
? ? end
end

%%%%%%%%%%計算分潮的調和常數H、g （詳見上機PPT第7頁）
f=[f_M2,f_S2,f_N2,f_K2,f_K1,f_O1,f_P1,f_Q1];% 8個分潮的交點因子f排序
uu=[uu_M2,uu_S2,uu_N2,uu_K2,uu_K1,uu_O1,uu_P1,uu_Q1]; % 8個分潮的交點訂正角uu排序
v0=[v0_M2,v0_S2,v0_N2,v0_K2,v0_K1,v0_O1,v0_P1,v0_Q1];% 8個分潮的初始相位v0排序
H=R./f'
g=sita + v0' + uu'
for i=1:8
?if( g(i) <= 0 )
? ? g(i)=g(i)+2*pi;
?elseif(g(i)>= 2*pi)
? ? g(i)=g(i)-2*pi;
?elseif( g(i)>=0 && g(i)<= 2*pi)
? ? ?continue
?end
end

%%%%%潮汐的預報和誤差分析%%%%%上機PPT第8頁
%%%%潮汐預報
%%%（i）計算各個時刻的自報潮位（n：-N'...0...N'）:
%根據程序chaoxi2獲取新的f vo u? 方法同上
f_new=[0.998818888564329,1,0.998818888564329,1.025068647340836,1.017975334956810,1.031932629767377,0.999403602946857,1.031932629767377];
v0_new=[0.050178215993988,6.283185307179497,6.133926494876704,1.861075997328001,0.957953566444805,5.375409956728703,5.325231740734693,5.175972928431833];
uu_new=[0.037881982228036,0,0.037881982228036,0.309899471911240,0.153025638522877,6.093785539018169,0.009695053375056,6.093785539018169];
?
n_new=(-372:372); ?%用time程序，3.3日作為-360，推到5.1日為1081
h_self_reported=zeros(1,745);
for i=1:745
? ? h_self_reported(1,i)=S0+sum(f_new.*H'.*cos(sig*n_new(i)+v0_new+uu_new-g'));
end

figure
x=1:745;
plot(x,h_self_reported,'r+','linewidth',2)
hold on
plot(x,h_self_reported,'b-','linewidth',1)
line([1,745],[0,0],'linestyle','--')
axis( [1 745 -80 60] )?
xlabel('時間(h)','FontSize',10)
ylabel('潮位(cm)','FontSize',10)
title('2003年5月1日0時至6月1日0時潮汐預報圖','FontSize',15)

?

總結

以上是生活随笔為你收集整理的潮汐预测 （长期预报的方法，数据为中期）的全部內容，希望文章能夠幫你解決所遇到的問題。

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