clc;clear   %先对分量进行调和分析，然后利用elli_para2程序计算椭圆要素，ap2ep程序画潮流椭圆

%%

xlsend = [3477 3897 3477 3478];

lat_point = [39+39.922/60    39+42.356/60   39+37.024/60    39+25.933/60];

tidename = ['O1  ';'K1  ';'M2  ';'S2  ';'M4  ';'MS4 '];

color={'r';'g';'b';'k';'c';'m';'r--';'g--';'b--';'k--';'c--';'m--'};

ef_l=[11 17 31 14];%有效水层

depth=[15.5      22    19.8    18.6]    %局地水深

%% 流速分析

%层循环

for i_site = 1:4            %站位循环

%     subplot(3,4,(layer-1)*4+i_site);

time = xlsread('.\QHD_ADCP.xlsx', i_site, 'B2:G2');

if i_site==2

N=xlsread('.\QHD_ADCP.xlsx',i_site,['J2:BZ',num2str(xlsend(i_site))]);

min = xlsread('.\QHD_ADCP.xlsx', i_site, ['F2:F',num2str(xlsend(i_site))]);

idx = find( min==17 | min==37 | min==57);

N = N(idx,:);

intev = 2/6;

else

N=xlsread('.\QHD_ADCP.xlsx',i_site,['J2:DJ',num2str(xlsend(i_site))]);

intev = 1/6;

end

for k=1:ef_l(i_site)

dep_v=N(1,3*k-2);

w = N(:,3*k-1)*0.1;% 换算成cm/s

dir= N(:,3*k);  u1 = w .* sind(dir);  v1 = w.* cosd(dir);%sind,cosd适用于单位为度的角度

fuc=ones(1,5)/5; u=filtfilt(fuc,1,u1); v=filtfilt(fuc,1,v1); %5点平滑

[u_tide,pout]=t_tide(u,'interval',intev, 'start',time,...

'latitude',lat_point(i_site),'output',['tideanalysis_u_',num2str(i_site),'.dat'], 'synthesis',1);

u_left = u - pout; % E--W，余流

[v_tide,pout]=t_tide(v,'interval',intev, 'start',time,...

'latitude',lat_point(i_site),'output',['tideanalysis_v_',num2str(i_site),'.dat'], 'synthesis',2);

v_left = v - pout; % N--S，余流

%%

u_ave = mean(u_left);   v_ave = mean(v_left);

umean(i_site,k)=u_ave; vmean(i_site,k)=v_ave;

[th,w_left] = cart2pol(u_ave,v_ave);

left_dir = 90-th*180/pi;

if left_dir<0

left_dir = left_dir+360;    %%转化成海流坐标系，角度制

end

save([num2str(i_site),'_left_current.mat'],'u_left','v_left','w_left','left_dir');

left(1,i_site,k) = w_left;

left(2,i_site,k) = left_dir;

%     fuc=ones(1,5)/5

%     plot(filtfilt(fuc,1,v_left));%hold on;plot(u_left);

%%

for i_tide = 1:4      %分潮循环

tname = tidename(i_tide,:);

for ii = 1:length(u_tide.name)

if strcmp(u_tide.name(ii,:),tname)==1

idx = ii;

freq = u_tide.freq(ii);

uamp = u_tide.tidecon(ii,1);

upha = u_tide.tidecon(ii,3);

end

if strcmp(u_tide.name(ii,:),tname)==1

idx = ii;

vamp = v_tide.tidecon(ii,1);

vpha = v_tide.tidecon(ii,3);

end

end

[W, w, theta, tao] = elli_para2(uamp, upha, vamp, vpha, freq);

elli(i_tide,1,k,i_site) = W;

elli(i_tide,2,k,i_site) = w/W;

elli(i_tide,3,k,i_site) = theta;

elli(i_tide,4,k,i_site) = tao;

elli(i_tide,5,k,i_site) = sind(upha-vpha);

if elli(i_tide,5,k,i_site)<0

l{i_tide}=color{i_tide};

else

l{i_tide}=color{i_tide+6};

end

[major,eccentricity,inclination,phase,w]=ap2ep(uamp,upha,vamp,vpha);

lli(i_tide,1,k,i_site) = major;

lli(i_tide,2,k,i_site) = eccentricity;

lli(i_tide,3,k,i_site) = inclination;

lli(i_tide,4,k,i_site) = phase;

aaa=squeeze(w);%画椭圆的关键数组

plot(real([aaa;aaa(1)]),imag([aaa;aaa(1)]),l{i_tide},'linewidth',1.5);hold on;%画椭圆图的函数，问题在这里

%-depth(i_site)+0.5-dep_v为对应水深

end           %分潮循环

%       if i_site==1

%            ylabel('V(cm/s)');

%       end

%       if layer==3

%            xlabel('U(cm/s)');

%       end

%       set(gca,'xlim',[-30,30],'ylim',[-30,30]); axis equal;

end           %站位循环

end              %层循环

% saveas(gca,'.\图\秦皇岛椭圆垂向分布.emf');

ellipse.png

(176.02 KB, 下载次数: 27)

2015-3-10 08:58 上传

分潮椭圆垂向结构