代码：clear all;close all;clc%M2,S2,N2,K2,K1,O1,P1,Q1八个分潮的杜德深常数U1、2、3、4、5、6、0U=[2,0,0,0,0,0,0;2,2,-2,0,0,0,0;2,-1,0,1,0,0,0;2,2,0,0,0,0,0;1,1,0,0,0,0,1;1,-1,0,0,0,0,-1;1,1,-2,0,0,0,-1;1,-2,0,1,0,0,-1;];%% 天文元素的变化速度twys_[]twys_=[14.4921,0.549,0.0411,0.0046,0.0022,0.0000000196].*(pi/180); %% 求天文要素twys[tao,s,h_,p,N_,p_]year=2003;mon=3;day=3+15;time=0;i=floor((year-1900)/4);ddd=76;dy=year-1900;di=ddd+i+time/24;s=277.02+129.3848*dy+13.1764*di;h_=280.19-0.2387*dy+0.9857*di;p=334.39+40.6625*dy+0.1114*di;N_=100.84+19.2382*dy+0.053*di;p_=281.22+0.0172*dy+0.00005*di;tao=15*time-s+h_;twys=[tao,s,h_,p,N_,p_].*pi/180;%% 求[M2,S2,N2,K2,K1,O1,P1,Q1]八个分潮角速度sig、初始相位V sigg=zeros(1,8);V0=zeros(1,8);for k =1:8sigg(k)=sum(U(k,1:6).*twys_);V0(k)=sum(U(k,1:7).*[twys,pi/2]);V0(k)=rem(V0(k),(2*pi));if abs(V0(k))<0.0001V0(k)=0;endif V0(k)<0V0(k)=V0(k)+2*pi;endend%% 求交点因子ff、交点订正角uufK1_cosuK1=1+0.002*cos(-2*twys(4)-twys(5))+0.0001*cos(-2*twys(5)) -0.0198*cos(-twys(5))+...0.1356*cos(twys(5))-0.0029*cos(2*twys(5));fK1_sinuK1=0.002*sin(-2*twys(4)-twys(5))+0.0001*sin(-2*twys(5))-0.0 198*sin(-twys(5))+...0.1356*sin(twys(5))-0.0029*sin(2*twys(5));uK1=atan(fK1_sinuK1/fK1_cosuK1)+2*pi;fK1=abs(sqrt(fK1_cosuK1^2+fK1_sinuK1^2));%%%计算M2分潮的交点因子和交点订正角fM2_cosuM2=1+0.0005*cos(-2*twys(5))-0.0373*cos(-twys(5))+0.0006* cos(2*twys(4))+...0.0002*cos(2*twys(4)+twys(5));fM2_sinuM2=0.0005*sin(-2*twys(5))-0.0373*sin(-twys(5))+0.0006*sin (2*twys(4))+...0.0002*sin(2*twys(4)+twys(5));uM2=atan(fM2_sinuM2/fM2_cosuM2)+2*pi;fM2=abs(sqrt(fM2_cosuM2^2+fM2_sinuM2^2));%计算K2分潮的交点因子和交点订正角fK2_cosuK2=1-0.0128*cos(-twys(5))+0.2890*cos(twys(5))+0.0324*cos (2*twys(5));fK2_sinuK2=-0.0128*sin(-twys(5))+0.2890*sin(twys(5))+0.0324*sin(2*t wys(5));uK2=atan(fK2_sinuK2/fK2_cosuK2)+2*pi;fK2=abs(sqrt(fK2_cosuK2^2+fK2_sinuK2^2));%%%计算P1分潮的交点因子和交点订正角fP1_cosuP1=1+0.0008*cos(-2*twys(5))-0.0112*cos(-twys(5))-0.0015*co s(2*twys(4))...-0.0003*cos(2*twys(4)+twys(5));fP1_sinuP1=0.0008*cos(-2*twys(5))-0.0112*cos(-twys(5))-0.0015*cos(2 *twys(4))...-0.0003*cos(2*twys(4)+twys(5));uP1=atan(fP1_sinuP1/fP1_cosuP1)+2*pi;fP1=abs(sqrt(fP1_cosuP1^2+fP1_sinuP1^2));%%%计算O1分潮的交点因子和交点订正角fO1_cosuO1=1-0.0058*cos(-2*twys(5))+0.1885*cos(-twys(5))+0.0002*c os(2*twys(4)-twys(5))...-0.0064*cos(2*twys(4))-0.0010*cos(2*twys(4)+twys(5));fO1_sinuO1=-0.0058*sin(-2*twys(5))+0.1885*sin(-twys(5))+0.0002*sin (2*twys(4)-twys(5))...-0.0064*sin(2*twys(4))-0.0010*sin(2*twys(4)+twys(5));uO1=atan(fO1_sinuO1/fO1_cosuO1);fO1=abs(sqrt(fO1_cosuO1^2+fO1_sinuO1^2));fS2=1;uS2=0;fQ1=fO1;uQ1=uO1;fN2=fM2;uN2=uM2;ff=[fM2,fS2,fN2,fK2,fK1,fO1,fP1,fQ1];ff=[0.9837,1,0.9837,1.511,1.0632,0.0968,0.9936,1.0968];uu=[uM2,uS2,uN2,uK2,uK1,uO1,uP1,uQ1];uu=[6.2504,0,6.2504,6.0167,6.1549,0.144,6.2724,0.144];%% 最小二乘法计算法方程date=load('E:\magicalcity\2020.3\海洋要素计算\上机实验\3.txt'); %构造系数矩阵AN=721;N_=360;A=zeros(9,9);for f=2:9aa=(sin(sigg(f-1)*N/2));bb=(sin(sigg(f-1)/2));A(1,f)=aa/bb;A(f,1)=A(1,f);endA(1,1)=N;for k=2:9for j=2:9if k==jA(k,j)=(N+sin(sigg(k-1)*N)/sin(sigg(k-1)))/2;endif k~=ja=(sin((sigg(k-1)-sigg(j-1))*N/2));b=(sin((sigg(k-1)-sigg(j-1))/2));c=(sin((sigg(k-1)+sigg(j-1))*N/2));d=(sin((sigg(k-1)+sigg(j-1))/2));A(k,j)=(a/b+c/d)/2;endendend%构造系数矩阵BB=zeros(8,8);for k=1:8for j=1:8if k==jB(k,j)=(N-sin(sigg(k)*N)/sin(sigg(k)))/2;endif k~=ja=(sin((sigg(k)-sigg(j))*N/2));b=(sin((sigg(k)-sigg(j))/2));c=(sin((sigg(k)+sigg(j))*N/2));d=(sin((sigg(k)+sigg(j))/2));B(k,j)=(a/b-c/d)/2;endendend%构造F1、F2F1=zeros(9,1);F2=zeros(8,1);F1(1,1)=sum(date(:,2));for k=2:9for j=1:NF1(k,1)=F1(k,1)+date(j,2)*cos(sigg(k-1)*(j-361));F2(k-1,1)=F2(k-1,1)+date(j,2)*sin(sigg(k-1)*(j-361));endend%法方程的解x，yx=A\F1;y=B\F2;%% 计算准调和振幅R和位相sitafor k=2:9R(k-1)=(x(k)*x(k)+y(k-1)*y(k-1))^(1/2);if x(k)>=0sita(k-1)=asin(y(k-1)/R(k-1));endif x(k)<0sita(k-1)=pi-asin(y(k-1)/R(k-1));endend%% 计算调和常数H、gfor k=1:8H(k)=R(k)./ff(k);g(k)=sita(k)+V0(k)+uu(k);g(k)=rem(g(k),2*pi);if g(k)<0g(k)=g(k)+2*pi;endg(k)=g(k)*180/pi;end%% 潮汐自报hhdelt=0;for k=1:Nhh(k)=0;for j=1:8hh(k)=hh(k)+ff(j)*H(j)*cos(sigg(j)*(k-361)+V0(j)+uu(j)-g(j));endhh(k)=hh(k)+x(1);delt=delt+(hh(k)-date(k,2))^2;enddelt=(delt/N)^(1/2);figure(1)plot(date(:,2),'r')hold on;plot(hh(1,:),'b')title('实测（红色）、自报（蓝色）潮汐曲线','fontsize',12)%% 5.1-6.1的潮汐预报hhpfor k=1:745hhp(k)=0;for j=1:8hhp(k)=hhp(k)+ff(j)*H(j)*cos(sigg(j)*(k+1056)+V0(j)+uu(j)-g (j));endendfigure(2)plot(hhp(1,:),'b')title('预报5.1-6.1潮汐','fontsize',12)。