%%Macro MATLAB R2012a
%Diplomová práce Jaroslav Krejčí
%Analýza lopatky turbíny
%Zpracování výsledků z ANSYSU
%nebezpečné otáčky
%_______________________________________________
clear all
%close all
clc
%PARAMETRY

for d=30:1:80
    clearvars -EXCEPT d
         %DÉLKA LOPATKY pro kterou se mají zpracovat výsledky

r=78;           %počet lopatek turbíny
p=3;            %počet uzlových průměrů
%nebezpečné násobky otáček
k(1)=37;
k(2)=74;
k(3)=111;
k(4)=148;
k(5)=148+37;
k(6)=148+37+37;
t=6;
%rozmezí otáček
ot_min=8000;   %spodní hranice sledovaných otáček
ot_max=14000;   %horní hranice sledovaných otáček
%barva_minoritni=[1,0.35,0];%[1,0,1];
barva_minoritni=[1,0.5,0];
%barva_minoritni2=[1,0.5,0]; %zelená
barva_minoritni2=[0,0.9,0]; %lutá
barva_frekvence='b';
barva_provozni='k';%[0.478,0.063,0.894];
%_______________________________________________

%načtení dat
load(['vysledky' num2str(d)],'-ascii')      
data=eval(['vysledky' num2str(d)]);              
l=length(data);                 %počet sloupců v matici výsledků

otacky=data(:,1);
for i=1:l-1
    eval(['frekvence' num2str(i) '=data(:,i+1);']);
end

%interpolace
otacky2=0:max(otacky);
for i=1:l-1
    frek=eval(['frekvence' num2str(i)]);
    fre=interp1(otacky,frek,otacky2,'linear');
    eval(['fre' num2str(i) '=fre;']);
end

%Nebezpečné násobky
n1=1; n2=1; n3=1; n4=1;
m1=1; m2=1; m3=1; m4=1;
X1(1)=-2000; X2(1)=-2000; X3(1)=-2000; X4(1)=-2000;
Y1(1)=0; Y2(1)=0; Y3(1)=0; Y4(1)=0;
for j=1:length(k);
    y=k(j)*otacky2/60;
    for i=1:l-1
        frekv=eval(['fre' num2str(i)]);
        for o=1:length(otacky2)-1
            if y(o)<=frekv(o) && y(o+1)>frekv(o)
                if j==i
                    X1(n1)=otacky2(o);
                    Y1(n1)=y(o);
                    r1(m1)=Y1(n1)-0.05*Y1(n1); r1(m1+1)=Y1(n1)+0.05*Y1(n1);
                    for I=1:length(y)-1
                        if y(I)<=r1(m1) && y(I+1)>r1(m1)
                            s1(m1)=otacky2(I);
                        end
                        if y(I)<=r1(m1+1) && y(I+1)>r1(m1+1)
                            s1(m1+1)=otacky2(I);
                        end
                    end
                    if length(s1)==m1
                        s1(m1+1)=14000;
                    end                    
                    n1=n1+1;
                    m1=m1+2;
                elseif i<=6
                %elseif j==i+1 || j==i-1 || j==i+2 || j==i-2 || j==i+3 || j==i-3 || j==i+4 || j==i-4
                    X2(n2)=otacky2(o);
                    Y2(n2)=y(o);
                    r2(m2)=Y2(n2)-0.05*Y2(n2); r2(m2+1)=Y2(n2)+0.05*Y2(n2);
                    for I=1:length(y)-1
                        if y(I)<=r2(m2) && y(I+1)>r2(m2)
                            s2(m2)=otacky2(I);
                        end
                        if y(I)<=r2(m2+1) && y(I+1)>r2(m2+1)
                            s2(m2+1)=otacky2(I);
                        end
                    end
                    if length(s2)==m2
                        s2(m2+1)=14000;
                    end
                    m2=m2+2;
                    n2=n2+1;
                else
                    X4(n4)=otacky2(o);
                    Y4(n4)=y(o);
                    r4(m4)=Y4(n4)-0.05*Y4(n4); r4(m4+1)=Y4(n4)+0.05*Y4(n4);
                    for I=1:length(y)-1
                        if y(I)<=r4(m4) && y(I+1)>r4(m4)
                            s4(m4)=otacky2(I);
                        end
                        if y(I)<=r4(m4+1) && y(I+1)>r4(m4+1)
                            s4(m4+1)=otacky2(I);
                        end
                    end
                    if length(s4)==m4
                        s4(m4+1)=14000;
                    end
                    m4=m4+2;
                    n4=n4+1;   
                end
            end
        end
    end
end

%Zobrazení výsledků
fig=1;
figure(fig)
set(figure(fig),'Position',[200 100 1000 550])
set(gca,'fontsize',14)
set(gca,'fontname','Times New Roman')
set(gca,'XTickLabel',{'8 000';'9 000';'10 000';'11 000';'12 000';'13 000';'14 000'})
%set(title(['Campbellův diagram (' num2str(d) '{\fontname{Times}\fontsize{12}\it mm})']),...
%    'Position',[9120,26.1])
%set(xlabel('Otáčky [{\fontname{Times}\fontsize{12}\itot/min}]'),'Position',[10990,-1.6]);
%set(ylabel('Frekvence [{\fontname{Times}\fontsize{12}\it kHz}]'),'Position',[7762,13]);
grid on
hold on
axis([ot_min,ot_max,30,80])
title('Nebezpečné otáčky turbíny')
xlabel('Otáčky [ot/min]')
ylabel('Délka lopatky [mm]')

plot([0,14000],[d,d],'b','LineWidth',t)

if m4>1
    JJ=1;
    for J=1:length(s4)/2
        plot([s4(JJ),s4(JJ+1)],[d,d],'color',barva_minoritni2,'LineWidth',t)
        JJ=JJ+2;
    end
end
JJ=1;
for J=1:length(s2)/2
    plot([s2(JJ),s2(JJ+1)],[d,d],'color',barva_minoritni,'LineWidth',t)
    JJ=JJ+2;
end
JJ=1;
for J=1:length(s1)/2
    plot([s1(JJ),s1(JJ+1)],[d,d],'r','LineWidth',t)
    JJ=JJ+2;
end
end

%Zobrazení provozních otáček
for l=1:13
plot([l*1000,l*1000],[30,80],':','color',barva_provozni,'LineWidth',1)
end
for l=1:9
plot([0,14000],[30+5*l,30+5*l],':','color',barva_provozni,'LineWidth',1)
end
plot([11000,11000],[30,80],'color',barva_provozni,'LineWidth',4)