matlab实现海洋要素计算代码

代码：

clear all;

close all;

clc

%M2,S2,N2,K2,K1,O1,P1,Q1八个分潮的杜德深常数U1、2、3、4、5、6、0

U=[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:8

sigg(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.0001

V0(k)=0;

end

if V0(k)<0

V0(k)=V0(k)+2*pi;

end

end

%% 求交点因子ff、交点订正角uu

fK1_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];