clear
clc 

%% GESTIONE DATI
DB = importdata('db_geometry');  
x = DB(:,1);    
y = DB(:,2);    

%%
figure(1)

hold on
plot(x,y,'k')

%inlet
plot(x(1),y(1),'ko','LineWidth',1)
plot([x(1) x(1)],[y(1) 0],'r','LineWidth',1)
plot(x(1),0,'ko','LineWidth',1)

% asse ugello
plot(x(end),y(end),'ko','LineWidth',1)
plot(x(end),0,'ko','LineWidth',1)
plot([x(1) x(end)],[0 0],'r','LineWidth',1)

% tratto dello spike 0.001 [m]
plot(x(end), y(end)+0.001,'ko','LineWidth',1)
plot([x(end) x(end)],[y(end) y(end)+0.001],'r','LineWidth',1)

% tratto r_e, tratto superiore ugello
plot(x(end) - 0.02 ,y(end)+0.001,'ko','LineWidth',1)
plot([x(end) x(end)-0.02],[y(end)+0.001 y(end)+0.001],'r','LineWidth',1)

% tratto farfield
plot(x(end)-0.02, y(end)+0.001+0.2,'ko','LineWidth',1)
plot([x(end)-0.02 x(end)-0.02],[y(end)+0.001 y(end)+0.001+0.2],'r','LineWidth',1)

% tratto orizzontale superiore
plot(x(end)+0.3,y(end)+0.001+0.2,'ko','LineWidth',1)
plot([x(end)-0.02 x(end)+0.3],[y(end)+0.001+0.2 y(end)+0.001+0.2],'r','LineWidth',1)

% Tratto orizzontale inferiore
plot(x(end)+0.3,0,'ko','LineWidth',1)
plot([x(end) x(end)+0.3],[0 0],'r','LineWidth',1)

% Tratto outlet
plot([x(end)+0.3 x(end)+0.3],[0 y(end)+0.001+0.2],'r','LineWidth',1)

grid on
grid minor
xlim([-0.02 0.5])
ylim([-0.02 0.3])

%risultato = DB(2651)-0.5*(DB(2651)-DB(1274)) 
L_tot = abs(x(1))+x(end);
L_exp = x(end)-0.062;
L_primo = L_exp/L_tot;
%%
clc 
clear
%%
DB = importdata('db_geometry');  
x = DB(:,1);    
y = DB(:,2);    
%% LEGGE DELLE AREE

Ai = DB(1,2)^2;
At = 0.01^2;
Ai_At = Ai/At;

Ae = DB(end,2)^2;
Ae_At = Ae/At;

% chiamo funzione M
[M_i, M_e] = calcola_M(Ai_At,Ae_At);

disp(['Il valore di M_i è: ', num2str(M_i)]);
disp(['Il valore di M_e è: ', num2str(M_e)]);


%% DATI 

%NPR = 20;
%NPR = 25.25;
%NPR = 30;
%NPR = 35;
%NPR = 40;
%NPR = 42;
%NPR = 44;
%NPR = 45.7;
%NPR = 46;
%NPR = 47.4;
%NPR = 48;
%NPR = 50;
%NPR = 50.5;
%NPR = 52;
%NPR = 54;
%NPR = 56;
%NPR = 58;
%NPR = 60;
%NPR = 62;
%NPR = 64;
%NPR = 66;
%NPR = 68;
%NPR = 70;
%NPR = 70.5;
%NPR = 71;
%NPR = 71.5;
%NPR = 72;
%NPR = 72.5;
%NPR = 73;
%NPR = 74;
%NPR = 76;

P_a = 101325;  %[Pa]

% chiamo funzione che calcola boundary condition
[P_t_inlet, P_inlet, T_i_t, P_outlet, T_o_t] = calcola_BC(NPR, M_i, P_a);

disp('BC INLET: ');
disp(['P totale inlet [Pa]: ', num2str(P_t_inlet)]);
disp(['P statica inlet [Pa]: ', num2str(P_inlet)]);
disp(['T totale inlet [K]: ', num2str(T_i_t)]);

disp('-------------------------------------------------------------------')
disp('BC OUTLET: ');
disp(['P statica outlet [Pa]: ', num2str(P_outlet)]);
disp(['T totale outlet [K]: ', num2str(T_o_t)]);

%% GESTIONE FILE DATI PER ANSYS

%gestisco i dati in modo da avere anche
% numero linea
% numero punto
n_linea = ones(length(x),1);
n_punto = linspace(1,length(x),length(x));
z = zeros(length(x),1);

%scrivo file

fileID = fopen('DB.dat','w');

for i = 1:length(x)
    fprintf(fileID,'%d %d %f %f %d\n',n_linea(i),n_punto(i),x(i),y(i),z(i));
end

fclose(fileID);

%% Aggiungo il contenuto del file "Dominio"

dominioID = fopen("Dominio.txt",'r');
fileID = fopen('DB.dat','a');

% copio il contenuto del file dominio
while ~feof(dominioID)
    line = fgets(dominioID);
    fprintf(fileID, '%s', line);
end

fclose(dominioID);
fclose(fileID);

%% DATO CHE DA PROBLEMI PER VIA DEL NUMERO ECCESSIVO DI PUNTI DIVIDO IN SOTTOBLOCCHI

DB_part = importdata('DB.dat');

int1 = 131;
int2 = 1271;

%suddivido in dati

block1 = DB_part(1:int1,:);
block2 = DB_part(int1:int2,:);
block3 = DB_part(int2:end,:);

filename1 = 'DB_part1.dat';
writematrix(block1,filename1, 'delimiter',' ');

filename2 = 'DB_part2.dat';
writematrix(block2,filename2, 'delimiter',' ');

filename3 = 'DB_part3.dat';
writematrix(block3,filename3, 'delimiter',' ');

%% 
clear 
clc
%% PROVA DEL PAPER
NPR = 45.7;
gamma = 1.4;
%inlet 
P_tot_inlet = 36e05;
P_a = P_tot_inlet/NPR;
P_inlet = P_tot_inlet/(1 + 0.5*(gamma - 1)*M_i^2)^(gamma/(gamma - 1));       %[Pa]
T_i_t = 283;

% outlet
P_outlet = P_a;
T_o_t = 283;

disp('BC INLET: ');
disp(['P totale inlet [Pa]: ', num2str(P_tot_inlet)]);
disp(['P statica inlet [Pa]: ', num2str(P_inlet)]);
disp(['T totale inlet [K]: ', num2str(T_i_t)]);

disp('-------------------------------------------------------------------')
disp('BC OUTLET: ');
disp(['P statica outlet [Pa]: ', num2str(P_outlet)]);
disp(['T totale outlet [K]: ', num2str(T_o_t)]);


%% 
clear clc
%%

figure(2)
hold on
plot(x,y,'b','LineWidth',1.5)
plot([x(1) x(end)],[0 0],'k--','LineWidth',1)
grid on
grid minor
title('$\textbf{Punti della geometria}$','Interpreter','latex','FontSize',15)
xlabel('$\mathbf{X/R^*}$','Interpreter','latex','FontSize',10)
ylabel('$\mathbf{R/R^*}$','Interpreter','latex','Rotation',0,'FontSize',12)
legend ('$\textbf{Profilo}$','$\textbf{Asse di simmetria}$','Interpreter','latex','Location','northwest','FontSize',10)
ylim([-0.01 0.07])

%% 
clear 
clc
%% GESTIONE DATI
DB = importdata('db_geometry');  
x = DB(:,1);    
y = DB(:,2);    

%% 

figure(3)
hold on
plot(x,y,'b','LineWidth',1)

%inlet
%plot(x(1),y(1),'ko','LineWidth',1)
plot([x(1) x(1)],[y(1) 0],'b','LineWidth',1)
%plot(x(1),0,'ko','LineWidth',1)

% asse ugello
%plot(x(end),y(end),'ko','LineWidth',1)
%plot(x(end),0,'ko','LineWidth',1)
plot([x(1) x(end)],[0 0],'b','LineWidth',1)

% tratto dello spike 0.001 [m]
%plot(x(end), y(end)+0.001,'ko','LineWidth',1)
plot([x(end) x(end)],[y(end) y(end)+0.001],'b','LineWidth',1)

% tratto r_e, tratto superiore ugello
%plot(x(end) - 0.02 ,y(end)+0.001,'ko','LineWidth',1)
plot([x(end) x(end)-0.02],[y(end)+0.001 y(end)+0.001],'b','LineWidth',1)

% tratto farfield
%plot(x(end)-0.02, y(end)+0.001+0.2,'ko','LineWidth',1)
plot([x(end)-0.02 x(end)-0.02],[y(end)+0.001 y(end)+0.001+0.2],'b','LineWidth',1)

% tratto orizzontale superiore
%plot(x(end)+0.3,y(end)+0.001+0.2,'ko','LineWidth',1)
plot([x(end)-0.02 x(end)+0.3],[y(end)+0.001+0.2 y(end)+0.001+0.2],'b','LineWidth',1)

% Tratto orizzontale inferiore
%plot(x(end)+0.3,0,'ko','LineWidth',1)
plot([x(end) x(end)+0.3],[0 0],'b','LineWidth',1)

% Tratto outlet
plot([x(end)+0.3 x(end)+0.3],[0 y(end)+0.001+0.2],'b','LineWidth',1)

% annotazioni

annotation("doublearrow",[0.362 0.362],[y(end)+0.3 y(end)+0.765],'Color','k','LineWidth',1);
text(x(end)-0.045, y(end)+0.001+0.1, '$\mathbf{10 r_e}$','FontSize',12,'HorizontalAlignment','Center','Interpreter','latex')

annotation("doublearrow",[0.38 0.841],[y(end)+0.801 y(end)+0.801],'Color','k','LineWidth',1);
text(x(end)+0.13, y(end)+0.22+0.01, '$\mathbf{16 r_e}$','FontSize',12,'HorizontalAlignment','Center','Interpreter','latex')

annotation("doublearrow",[0.86 0.86],[y(end)+0.175 y(end)+0.765],'Color','k','LineWidth',1);
text(x(end)+0.325, y(end)+0.07, '$\mathbf{13 r_e}$','FontSize',12,'HorizontalAlignment','Center','Interpreter','latex')

annotation("doublearrow",[0.38 0.41],[y(end)+0.32 y(end)+0.32],'Color','k','LineWidth',1);
text(x(end)-0.009, y(end)+0.02, '$\mathbf{r_e}$','FontSize',12,'HorizontalAlignment','Center','Interpreter','latex')

ylim([-0.05 0.3])
xlim([-0.05 0.485911])
grid on
grid minor
title('$\textbf{Geometria con dominio esterno}$','Interpreter','latex','FontSize',15)
xlabel('$\mathbf{X/R^*}$','Interpreter','latex','FontSize',10)
ylabel('$\mathbf{R/R^*}$','Interpreter','latex','Rotation',0,'FontSize',12)