潮汐预测 （长期预报的方法，数据为中期）

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:10