% skript pro vykresleni lopatky s profily rady NACA 65 s kruhovou stredni 
% krivkou
% osa rotoru se shoduje s osou x, y je tecny smer, z radialni smer
% analyticky definovana stredni krivka
% parabolicky modifikovana odtokova hrana dle Jerie/Poul 

clear all
close all

% vstupni data
% tetiva [mm]
cv = 34:-1:26; 

% relativni tloustka
trHub = 0.1;
trShroud = 0.08;
trv = trHub:((trShroud - trHub)/8):trShroud;

% polomery jednotlivych rezu [mm]
rv = 40:5:80; 

% vstupni a vystupni geometricky uhel - aproximace po vyhlazeni
alfa1dv = zeros(1,length(trv));
alfa2dv = zeros(1,length(trv));
for ii = 1:length(trv)
    rm = rv(ii);
    alfa1dv(ii) = 0.000008390876300*rm^3 - 0.005253949467800*rm^2 + 0.873154025895360*rm + 25.548458427228100;
    alfa2dv(ii) = 0.000386586498683*rm^3 - 0.091357459631276*rm^2 + 7.770914633598370*rm - 183.207846188765000;

       
end

% tloustka odtokove hrany [mm]
tOH = 0.01;

% pocet intervalu na horni/spodni hrane profilu
nbod = 2*250;

% deklarace promennych
Xlop = zeros(length(trv),2*nbod+1);
Ylop = zeros(length(trv),2*nbod+1);
Zlop = zeros(length(trv),2*nbod+1);

for rez = 1:length(trv)
% pomocne operace
alfa1 = alfa1dv(rez)*pi/180;
alfa2 = alfa2dv(rez)*pi/180;
phi = alfa1-alfa2;
c = cv(rez);
tr = trv(rez);

% stredni krivka

r = c/(2*sin(phi/2));
x0 = r*sin(alfa1);
y0 = -1*r*cos(alfa1);

theta = (pi/2-alfa1):(phi/nbod):(pi/2-alfa1+phi);

Xs = zeros(1,nbod+1);
Ys = zeros(1,nbod+1);
ts = zeros(1,nbod+1);
xr = zeros(1,nbod+1);
Xh = zeros(1,nbod+1);
Yh = zeros(1,nbod+1);
Xd = zeros(1,nbod+1);
Yd = zeros(1,nbod+1);

for ii = 1:nbod+1
    
    Xs(ii) = x0 - r*cos(theta(ii));
    Ys(ii) = y0 + r*sin(theta(ii));
    
    xrel = (ii-1)/nbod;
    xr(ii) = xrel;
    
    ts(ii) = tr*c*((1-xrel)*(1.0675*sqrt(xrel)-0.2758*(xrel) ...
        +2.4478*(xrel)^2-2.8385*(xrel)^3)/(1-0.1760*(xrel)));
        
    if xrel>0.4
        ts(ii) = ts(ii) + 0.5*tOH*((xrel-0.4)/0.6)^2;
    end
    
    Xh(ii) = x0 - (r+ts(ii))*cos(theta(ii));
    Yh(ii) = y0 + (r+ts(ii))*sin(theta(ii));
    
    Xd(ii) = x0 - (r-ts(ii))*cos(theta(ii));
    Yd(ii) = y0 + (r-ts(ii))*sin(theta(ii));
    
    
    
end

% vypocet polohy teziste
S = 0;
Mx = 0;
My = 0;

for ii = 1:nbod
    % horni trojuhelniky
    x1 = Xd(ii);
    x2 = Xh(ii);
    x3 = Xh(ii+1);
    y1 = Yd(ii);
    y2 = Yh(ii);
    y3 = Yh(ii+1);
    Si = abs(0.5*((x2-x1)*(y3-y1)-(y2-y1)*(x3-x1)));
    xti = (x1+x2+x3)/3;
    yti = (y1+y2+y3)/3;
    
    S = S + Si;
    Mx = Mx + Si*xti;
    My = My + Si*yti;
    
    % dolni trojuhelniky
    x1 = Xd(ii);
    x2 = Xd(ii+1);
    x3 = Xh(ii+1);
    y1 = Yd(ii);
    y2 = Yd(ii+1);
    y3 = Yh(ii+1);
    Si = abs(0.5*((x2-x1)*(y3-y1)-(y2-y1)*(x3-x1)));
    xti = (x1+x2+x3)/3;
    yti = (y1+y2+y3)/3;
    
    S = S + Si;
    Mx = Mx + Si*xti;
    My = My + Si*yti;
    
end

XT = Mx/S;
YT = My/S;

Xd = Xd - XT;
Xh = Xh - XT;
Xs = Xs - XT;
Ys = Ys - YT;
Yd = Yd - YT;
Yh = Yh - YT;

% globalni souradnice
for ii = 1:nbod+1
    theta = Yh(ii)/rv(rez);
    Xlop(rez,nbod+2-ii) = Xh(ii);
    Ylop(rez,nbod+2-ii) = Yh(ii);
    Zlop(rez,nbod+2-ii) = Yh(ii);
    
    theta = Yd(ii)/rv(rez);
    Xlop(rez,nbod+ii)   = Xd(ii);
    Ylop(rez,nbod+ii)   = Yd(ii);
    Zlop(rez,nbod+ii)   = Yd(ii);
    
    % stredni krivka
    theta = Ys(ii)/rv(rez);
    XSp(rez,ii)   = Xs(ii);
    YSp(rez,ii)   = Ys(ii);
    ZSp(rez,ii)   = Ys(ii);
    
 end   
%     
%     for ii = 1:nbod+1
%     theta = Yh(ii)/rv(rez);
%     Xlop(rez,nbod+2-ii) = Xh(ii);
%     Ylop(rez,nbod+2-ii) = rv(rez)*sin(theta);
%     Zlop(rez,nbod+2-ii) = rv(rez)*cos(theta);
%     
%     theta = Yd(ii)/rv(rez);
%     Xlop(rez,nbod+ii)   = Xd(ii);
%     Ylop(rez,nbod+ii)   = rv(rez)*sin(theta);
%     Zlop(rez,nbod+ii)   = rv(rez)*cos(theta);
%     
%     % stredni krivka
%     theta = Ys(ii)/rv(rez);
%     XSp(rez,ii)   = Xs(ii);
%     YSp(rez,ii)   = rv(rez)*sin(theta);
%     ZSp(rez,ii)   = rv(rez)*cos(theta);
%     
    
% end

% grafy
figure
plot(Xs,Ys,'-',Xh,Yh,'k',Xd,Yd,'k',[0],[0],'+')
axis equal
grid on


% overeni tloustky OH
tOHo(rez) = sqrt((Xlop(rez,1)-Xlop(rez,2*nbod+1))^2+(Ylop(rez,1)-Ylop(rez,2*nbod+1))^2);


end


% vykresleni profilu po vysce
figure
hold on
for ii = 1:length(trv)
    plot(Xlop(ii,:),Ylop(ii,:),'k')
    
end
plot([0],[0],'r+')
% set(gca,'XColor', 'none','YColor','none')
hold off

figure
surf(Xlop,Ylop,Zlop)
xlabel('x [mm]')
ylabel('y [mm]')
zlabel('z [mm]')

figure
surf(XSp,YSp,ZSp)
xlabel('x [mm]')
ylabel('y [mm]')
zlabel('z [mm]')
title('plocha strednich krivek')

figure
plot(alfa1dv,rv,alfa2dv,rv)
xlabel('\alpha [\circ]')
ylabel('radius [mm]')
legend('\alpha_1','\alpha_2','Location','northwest')
grid on

% zapis dat do souboru
% profily
for rez = 1:length(trv)
    data = zeros(2*nbod+1,3);
    for ii = 1:2*nbod+1
        data(ii,1) = Xlop(rez,ii);
        data(ii,2) = Ylop(rez,ii);
%         data(ii,3) = Zlop(rez,ii); %tohle okomentovat %%
    end
    
    if rez <10
        filename = strcat('rez_rotor_0',num2str(rez));
    else
        filename = strcat('rez_rotor_',num2str(rez));
    end
    
    save(filename,'data','-ASCII');
end

% stredni krivky
for rez = 1:length(trv)
    data = zeros(2*nbod+1,3);
    for ii = 1:nbod+1
        data(ii,1) = XSp(rez,ii);
        data(ii,2) = YSp(rez,ii);
        data(ii,3) = ZSp(rez,ii);
    end
    
    if rez <10
        filename = strcat('strKrivka_rotor_0',num2str(rez));
    else
        filename = strcat('strKrivka_rotor_',num2str(rez));
    end
    
    save(filename,'data','-ASCII');
end