%%Macro MATLAB R2012a
%Diplomová práce Jaroslav Krejčí
%Analýza lopatky turbíny
%Zpracování výsledků z ANSYSU
%tvorba Campbellova diagramu
%_______________________________________________
clear all
%close all
clc
%PARAMETRY
d=60;           %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)=185;
%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=[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
% if m4>1
%     for II=1:m4-2%length(s4)    
%         if s4(II)==0
%             s4(II)=14000;
%         end
%     end
% end

%Zobrazení výsledků
graf=2;
figure(graf)
set(figure(graf),'Position',[200 100 700 400])
%set(figure(graf),'Position',[200 100 700 400])
plot(X1,Y1/1000,'r.','MarkerSize',30)
hold on
plot(X2,Y2/1000,'.','color',barva_minoritni,'MarkerSize',23)
%legend('Kritické otáčky','Provozní otáčky')
%plot(X4,Y4/1000,'g.','MarkerSize',30)
set(gca,'fontsize',12)
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
axis([ot_min,ot_max,0,26])
plot(otacky2,fre1/1000,'color',barva_frekvence,'LineWidth',2)
y=k(1)*otacky2/60;
plot(otacky2,y/1000,'k','LineWidth',2)
% legend('majoritní rezonance','minoritní rezonance',...'minoritní rezonance',...
%             'Vlastní frekvence','Náběhová přímka')


for J=1:length(r1)
    plot([s1(J),s1(J)],[1000,r1(J)]/1000,'r')
end
% for J=1:length(r2)
%     plot([s2(J),s2(J)],[3000,r2(J)]/1000,'color',barva_minoritni)
% end


% for J=1:length(r3)
%     plot([s3(J),s3(J)],[3000,r3(J)]/1000,'m')
% end
%plot(s1,r1/1000,'r*')
%plot(s2,r2/1000,'*','color',barva_minoritni)
%plot(s4,r4/1000,'m*')
if m4>1
    JJ=1;
    for J=1:length(s4)/2
        plot([s4(JJ),s4(JJ+1)],[0,0],'color',barva_minoritni2,'LineWidth',16)
        JJ=JJ+2;
    end
end
JJ=1;
for J=1:length(s2)/2
    plot([s2(JJ),s2(JJ+1)],[0,0],'color',barva_minoritni,'LineWidth',16)
    JJ=JJ+2;
end
JJ=1;
for J=1:length(s1)/2
    plot([s1(JJ),s1(JJ+1)],[0,0],'r','LineWidth',16)
    JJ=JJ+2;
end

% Zobrazení vlastních frekvencí


%Zobrazení náběhových přímek
for j=2:length(k);
    y=k(j)*otacky2/60;
    plot(otacky2,y/1000,'k','LineWidth',2)
end

%Zobrazení provozních otáček
plot([11000,11000],[0,1.2],'color',barva_provozni,'LineWidth',4)

if d==30
    text(13400,9.3,'\fontname{Times}\fontsize{12}\bf\itk=1x37')
    text(13400,18.3,'\fontname{Times}\fontsize{12}\bf\itk=2x37')
elseif d==40
    text(13400,10.6,'\fontname{Times}\fontsize{12}\bf\itk=1x37')
    text(13400,18,'\fontname{Times}\fontsize{12}\bf\itk=2x37')
elseif d==50
    text(13400,9.3,'\fontname{Times}\fontsize{12}\bf\itk=1x37')
    text(13400,18,'\fontname{Times}\fontsize{12}\bf\itk=2x37')
elseif d==60
    text(13400,9.3,'\fontname{Times}\fontsize{12}\bf\itk=1x37')
    text(13400,18,'\fontname{Times}\fontsize{12}\bf\itk=2x37')
elseif d==70
    text(13400,9.3,'\fontname{Times}\fontsize{12}\bf\itk=1x37')
    text(13400,18,'\fontname{Times}\fontsize{12}\bf\itk=2x37')
else
    text(13400,9.3,'\fontname{Times}\fontsize{12}\bf\itk=1x37')
    text(13400,18,'\fontname{Times}\fontsize{12}\bf\itk=2x37')
end
text(13400,27,'\fontname{Times}\fontsize{12}\bf\itk=3x37')
text(10640,27,'\fontname{Times}\fontsize{12}\bf\itk=4x37')

for i=1:l-1
    plot(otacky2,eval(['fre' num2str(i) '/1000']),'color',barva_frekvence,'LineWidth',2)
end
plot(X4,Y4/1000,'.','color',barva_minoritni2,'MarkerSize',23)
plot(X2,Y2/1000,'.','color',barva_minoritni,'MarkerSize',23)
plot(X1,Y1/1000,'r.','MarkerSize',30)

