使用S_TIDE绘制潮流椭圆并计算椭圆参数

目前网上似乎没有公开的代码计算潮流椭圆参数，给潮汐潮流分析带来了不便。这里，我介绍一下我本人开发的新型潮汐潮流调和分析MATLAB工具包S_TIDE（ https://www.researchgate.net/project/A-non-stationary-tidal-analysis-toolbox-S-TIDE）里的潮流椭圆的有关函数，希望能帮助到大家。

s_plot_tidal_ellipse.m是绘制潮流椭圆的函数，有5个输入参数uam,uph,vam,vph,tides，分别代表u方向分潮流流速振幅，u方向分潮流流速迟角，v方向分潮流流速振幅，v方向分潮流流速迟角，以及对应的分潮名字。s_estimate_tidal_ellipse.m是计算潮流椭圆参数的函数，输入参数同上，有6个输出参数maxc,minc,maxt,mint,maxp,minp,分别代表最大分潮流流速（即椭圆长轴），最小分潮流流速（即潮流椭圆短轴），最大分潮流流速时刻（从0开始计时），最小分潮流流速时刻（从0开始计时），最大分潮流流速方向（即椭圆长轴与x轴正方向的夹角），最小分潮流流速方向（即椭圆短轴与x轴正方向的夹角）

例子：

uam=38.6;uph=247;

vam=25.0;vph=222;

s_plot_tidal_ellipse(uam,uph,vam,vph,{'M2'})

[maxc,minc,maxt,mint,maxp,minp]=s_estimate_tidal_ellipse(uam,uph,vam,vph,'M2') ;

%下面是我提供的一个动态实例(可以直接运行程序，然后你会看到潮流椭圆的生成过程）来帮助大家理解潮流椭圆要素的具体含义

uam=38.6;uph=247;

vam=25.0;vph=222;

tides={'M2'};

load t_constituents

ju=strmatch(upper(tides),const.name);

fu=const.freq(ju);

t=[0:1/12:36]; %from 0 hour to 36 hours 从0时刻到36小时

u=uam*cos(fu*t*2*pi-uph/180*pi);

v=vam*cos(fu*t*2*pi-vph/180*pi);

L=sqrt(u.^2+v.^2);

for i=1:160

subplot(2,1,1);plot(u(i),v(i),'k.');

title([ 'time=',num2str(t(i)),'hours'])

hold on;xlim([-40 40]);ylim([-30 30]);

xlabel('u(cm/s)','fontsize',12);ylabel('v(cm/s)','fontsize',12);

pause(0.04)

subplot(2,1,2);plot(t(i),L(i),'k.');hold on

xlim([t(1) t(160)]);ylim([-10 60])

xlabel('time(hours)','fontsize',12);ylabel('the speed of tidal current vector(cm/s)','fontsize',12)

pause(0.04)

end

[a1,a2]=findpeaks(L);maxt=t(a2);

[a3,a4]=findpeaks(-1*L);mint=t(a4);

maxc=max(L);minc=min(L);

line([t(a2(1)) t(a2(1))],[-10 maxc],'color','r');

%line([0 t(a2(1))],[maxc maxc],'color','r');

text(0.5,20,'maxt = 2.083 hour');

text(0.5,25,'maxc = 45.09')

line([t(a2(2)) t(a2(2))],[-10 maxc],'color','r');

%line([0 t(a2(2))],[maxc maxc],'color','r');

text(8.4,22,'maxt = 8.29 hour');

text(8.4,27,'maxc = 45.09')

line([t(a4(1)) t(a4(1))],[-10 minc],'color','r');

%line([0 t(a4(1))],[minc minc],'color','r');

text(4.5,22,'mint = 5.17 hour');

text(4.5,27,'minc = 9.05')

line([t(a4(2)) t(a4(2))],[-10 minc],'color','r');

%line([0 t(a4(2))],[minc minc],'color','r');

text(11,21,'mint = 11.38 hour');

text(11,26,'minc = 9.05')

subplot(2,1,1);line([-40 40],[0 0],'color','r');

line([0 0],[-30 30],'color','r');

line([u(a2(1)) u(a2(2))],[v(a2(1)) v(a2(2))]);

text(10,3,'31.76°')

text(5,-13,'There are two angles for maximum tidal currents','Fontsize',12)

text(7,-18,'One is 31.76°, the other is 31.76°+ 180°','Fontsize',12)