目前網(wǎng)上似乎沒有公開的代碼計(jì)算潮流橢圓參數(shù)，給潮汐潮流分析帶來了不便。這里，我介紹一下我本人開發(fā)的新型潮汐潮流調(diào)和分析MATLAB工具包S_TIDE（?https://www.researchgate.net/project/A-non-stationary-tidal-analysis-toolbox-S-TIDE）里的潮流橢圓的有關(guān)函數(shù)，希望能幫助到大家。

s_plot_tidal_ellipse.m是繪制潮流橢圓的函數(shù)，有5個(gè)輸入?yún)?shù)uam,uph,vam,vph,tides，分別代表u方向分潮流流速振幅，u方向分潮流流速遲角，v方向分潮流流速振幅，v方向分潮流流速遲角，以及對(duì)應(yīng)的分潮名字。s_estimate_tidal_ellipse.m是計(jì)算潮流橢圓參數(shù)的函數(shù)，輸入?yún)?shù)同上，有6個(gè)輸出參數(shù)maxc,minc,maxt,mint,maxp,minp,分別代表最大分潮流流速（即橢圓長軸），最小分潮流流速（即潮流橢圓短軸），最大分潮流流速時(shí)刻（從0開始計(jì)時(shí)），最小分潮流流速時(shí)刻（從0開始計(jì)時(shí)），最大分潮流流速方向（即橢圓長軸與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') ;

%下面是我提供的一個(gè)動(dòng)態(tài)實(shí)例(可以直接運(yùn)行程序，然后你會(huì)看到潮流橢圓的生成過程）來幫助大家理解潮流橢圓要素的具體含義

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時(shí)刻到36小時(shí)

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)