
clear
clc
close all

N=100;                                                                        % N nodi su x carico (da 1 a 1000 g)
M=100;                                                                        % M nodi su z densità (da 70 a 160 J/mm3)
max_carico=1000;

x=linspace(1,max_carico,N);                                                 %vettore carico in grammi

y1_af = -26.82*log(x) + 745.44;                                             %HV provino C1
y2_af = -24.37*log(x) + 721.12;                                             %HV provino C2
y3_af = -17.96*log(x) + 625.37;                                             %HV provino C3

y4_af = -7.738*log(x) + 559.56;                                             %HV provino C4
y5_af = -14.66*log(x) + 661.72;                                             %HV provino C5
y6_af = -25.77*log(x) + 679.74;                                             %HV provino C6
y7_af = -18.67*log(x) + 636.34;                                             %HV provino C7

y1_tt = -15.73*log(x) + 421.18;                                              %HV provino C1_tt
y2_tt = -19.37*log(x) + 456.98;                                              %HV provino C2_tt
y3_tt = -20.47*log(x) +542.55;                                               %HV provino C3_tt

y4_tt = -13.2*log(x) + 525.57;                                               %HV provino C4_tt
y5_tt = -12.59*log(x) + 401.69;                                              %HV provino C5_tt
y6_tt = -17.02*log(x) + 464.88;                                              %HV provino C6_tt
y7_tt = -23.06*log(x) + 479.74;                                              %HV provino C7_tt

z1=400;  % 152.381 J/mm3; 192W 400mm/s                                     %densità energetica (ParProcAM) per provinoC1
z2=600;  % 101.587 J/mm3; 192W 600mm/s                                     %densità energetica (ParProcAM) per provinoC2
z3=800;  % 76.1905 J/mm3; 192W 800mm/s                                     %densità energetica (ParProcAM) per provinoC3

z4=500;  % 146.032 J/mm3; 230W 500mm/s                                     %densità energetica (ParProcAM) per provinoC4
z5=600;  % 121.693 J/mm3; 230W 600mm/s                                     %densità energetica (ParProcAM) per provinoC5
z6=700;  % 104.308 J/mm3; 230W 700mm/s                                     %densità energetica (ParProcAM) per provinoC6
z7=900;  % 81.1287 J/mm3; 230W 900mm/s

% scelta As Fabricated o Trattam Termico
provini = input('if As_Fabricated press " 1 "; if Heat_Treatment press "0" ');
 
if (provini == 1)
    y1=y1_af;
    y2=y2_af;
    y3=y3_af;
    y4=y4_af;
    y5=y5_af;
    y6=y6_af;
    y7=y7_af;

SIGMA_SN= 530;
rp1=  -0.00008.*x.^2 + 0.0946.*x + 521.1; 
rp2= SIGMA_SN*ones(1,N);
rp3= -0.00001.*x.^2 + 0.1455.*x + 671.36;  
rp4= -0.00007.*x.^2 + 0.1379.*x + 687.8; 
rp5= -0.00003.*x.^2 + 0.0575.*x + 515.07;  
rp6= -0.000009.*x.^2 + 0.0098.*x + 585.82;  
rp7= -0.00007.*x.^2 + 0.0467.*x + 566.41; 
    
end
if (provini == 0)
    y1=y1_tt;
    y2=y2_tt;
    y3=y3_tt;
    y4=y4_tt;
    y5=y5_tt;
    y6=y6_tt;
    y7=y7_tt;

SIGMA_SN= 900;
rp1=  -0.00004.*x.^2 + 0.0093.*x + 948.22; 
rp2= SIGMA_SN*ones(1,N);
rp3= -0.00007.*x.^2 + 0.0519.*x + 819.24;  
rp4= -0.0002.*x.^2 + 0.2574.*x + 770.95; 
rp5= -0.0001.*x.^2 + 0.1306.*x + 888.42;  
rp6= -0.00005.*x.^2 + 0.00972.*x + 863.4;  
rp7= -0.00004.*x.^2 + 0.0186.*x + 835.12; 
    
end
 

% % curve carico/HV traslate per differenti valori di densità energetica (3D non interpolato)
% figure
% plot3(x,y1,v1*ones(1,M),'LineWidth',1)
% xlabel('carico [g]')
% ylabel('durezza HV')
% zlabel('velocità laser [mm/s]')
% hold on
% plot3(x,y2,v2*ones(1,M),'LineWidth',1)
% plot3(x,y3,v3*ones(1,M),'LineWidth',1)
% title('Potenza laser 190 W')
% legend('v1 = 400 mm/s','v2 = 600 mm/s','v3 = 800 mm/s')
% figure
% plot3(x,y4,v4*ones(1,M),'LineWidth',1)
% xlabel('carico [g]')
% ylabel('durezza HV')
% zlabel('velocità laser [mm/s]')
% hold on
% plot3(x,y5,v5*ones(1,M),'LineWidth',1)
% plot3(x,y6,v6*ones(1,M),'LineWidth',1)
% plot3(x,y7,v7*ones(1,M),'LineWidth',1)
% title('Potenza laser 230 W')
% legend('v4 = 500 mm/s','v5 = 600 mm/s','v6 = 700 mm/s','v7 = 900 mm/s')

%% interpolazione curve P=192W per ottenere superficie
Y192=[y1; y2; y3]; % matrice dureza 3 x N
Z192=[z1; z2; z3]; % matrice densità energetiche 1 x 3
RP192=[rp1; rp2; rp3]; %matrice di sigma sn 3 x N

Y230=[y4; y5; y6; y7]; % matrice dureza 4 x N
Z230=[z4; z5; z6; z7]; % matrice densità energetiche 1 x 4
RP230=[rp4; rp5; rp6; rp7]; % matrice di sigma sn 4 x N

%matrice per costruire polinomio
% P 192W
mat192=zeros(3,N,3);
for i=1:3
    for j=1:N
        for k=1:3
            if(k==1)
            mat192(k,j,i)=x(j); % posiziono il carico nella matrice; es: 25g x 7prove
            end
            if(k==2)
             mat192(k,j,i)=Y192(i,j); % posiziono la durezza nella matrice x ogni carico e x ogni densità energ
            end
            if(k==3)
                mat192(k,j,i)=Z192(i); % posiziono la densità energetica che varia per ogni prova i
            end
        end
    end
end

mat1_rp192=zeros(3,N,3);
for i=1:3
    for j=1:N
        for k=1:3
            if(k==1)
            mat1_rp192(k,j,i)=x(j); % posiziono il carico nella matrice; es: 25g x 7prove
            end
            if(k==2)
             mat1_rp192(k,j,i)=RP192(i,j); % posiziono l'Rp02 nella matrice x ogni carico e x ogni densità energ
            end
            if(k==3)
                mat1_rp192(k,j,i)=Z192(i); % posiziono la densità energetica che varia per ogni prova i
            end
        end
    end
end

% P 230
mat230=zeros(3,N,4);
for i=1:4
    for j=1:N
        for k=1:3
            if(k==1)
            mat230(k,j,i)=x(j); % posiziono il carico nella matrice; es: 25g x 7prove
            end
            if(k==2)
             mat230(k,j,i)=Y230(i,j); % posiziono la durezza nella matrice x ogni carico e x ogni densità energ
            end
            if(k==3)
             mat230(k,j,i)=Z230(i); % posiziono la densità energetica che varia per ogni prova i
            end
        end
    end
end

mat1_rp230=zeros(3,N,4);
for i=1:4
    for j=1:N
        for k=1:3
            if(k==1)
            mat1_rp230(k,j,i)=x(j); % posiziono il carico nella matrice; es: 25g x 7prove
            end
            if(k==2)
             mat1_rp230(k,j,i)=RP230(i,j); % posiziono l'Rp02 nella matrice x ogni carico e x ogni densità energ
            end
            if(k==3)
             mat1_rp230(k,j,i)=Z230(i); % posiziono la densità energetica che varia per ogni prova i
            end
        end
    end
end

%punti delle 7 prove al variare del carico
% P =192W
for j=1:N
    for i=1:3
        nodiy192(i,j)=mat192(2,j,i); %nodi durezza HV
        nodiz192(i,j)=mat192(3,j,i); %nodi densità energetica
    end
end

for j=1:N
    for i = 1:3
        nodi_rp192(i,j)=mat1_rp192(2,j,i);
        %nodiz192(i,j)=mat192(3,j,i);
    end
end

% P =230W
for j=1:N
    for i=1:4
        nodiy230(i,j)=mat230(2,j,i); %nodi durezza HV
        nodiz230(i,j)=mat230(3,j,i); %nodi densità energetica
    end
end
for j=1:N
    for i=1:4
        nodi_rp230(i,j)=mat1_rp230(2,j,i); %nodi durezza HV
        nodiz230(i,j)=mat230(3,j,i); %nodi densità energetica
    end
end

%ciclo for per avere N nodi; p = polyfit(x,y,n) returns the COEFFICIENTS for a polynomial p(x) of degree n that is a best fit (in a least-squares sense) for the data in y.
% P = 192 W
for i=1:N
    p192(:,i)=polyfit(nodiz192(:,i)',nodiy192(:,i)',3); % matrice coefficienti per la durezza
end

for i=1:N
    p_rp192(:,i)=polyfit(nodiz192(:,i)',nodi_rp192(:,i)',3); % matrice coefficienti per l'Rp02
end
zz192=linspace(nodiz192(1),nodiz192(3),M);

% P = 230 W
for i=1:N
p230(:,i)=polyfit(nodiz230(:,i)',nodiy230(:,i)',3); % matrice coefficienti per la durezza
end

for i=1:N
p_rp230(:,i)=polyfit(nodiz230(:,i)',nodi_rp230(:,i)',3); % matrice coefficienti per Rp02
end
zz230=linspace(nodiz230(1),nodiz230(4),M);

%% scelta carico e velocità
xs=input('Inserire valore di carico tra 1g e 1000g: ');
errx=(x(3)-x(2))/2;
zs=input('Inserire valore di velocità tra 500 e 800 mm/s: ');
errz192=(zz192(3)-zz192(2))/2;
errz230=(zz230(3)-zz230(2))/2;

%% interpolazione e ricerca valori sulla superficie
% matrice delle durezze valuate con polyval

%P=192W
for i=1:M
    yy192(i,:)=polyval(p192(:,i),zz192); % matrice delle durezze valuate con polyval
end 
for i=1:M
    rp192(i,:)=polyval(p_rp192(:,i),zz192); % matrice degli Rp02 valuati con polyval
end 

%P=230W
for i=1:M
    yy230(i,:)=polyval(p230(:,i),zz230); % matrice delle durezze valuate con polyval
end 
for i=1:M
    rp230(i,:)=polyval(p_rp230(:,i),zz230); % matrice degli Rp02 valuati con polyval
end 
%matrice in tre dimensioni per carico durezza densità energetica

%P=192W
mat1_192=zeros(3,N,M);
for i=1:M
    for j=1:N
        for k=1:3
            if(k==1)
            mat1_192(k,j,i)=x(j);
            end
            if(k==2)
             mat1_192(k,j,i)=yy192(j,i);
            end
            if(k==3)
                mat1_192(k,j,i)=zz192(i);
            end
        end
    end
end

mat1_rp_192=zeros(3,N,M);
for i=1:M
    for j=1:N
        for k=1:3
            if(k==1)
            mat1_rp_192(k,j,i)=x(j);
            end
            if(k==2)
             mat1_rp_192(k,j,i)=rp192(j,i);
            end
            if(k==3)
                mat1_rp_192(k,j,i)=zz192(i);
            end
        end
    end
end

%P=230W
mat1_230=zeros(3,N,M);
for i=1:M
    for j=1:N
        for k=1:3
            if(k==1)
            mat1_230(k,j,i)=x(j);
            end
            if(k==2)
             mat1_230(k,j,i)=yy230(j,i);
            end
            if(k==3)
                mat1_230(k,j,i)=zz230(i);
            end
        end
    end
end

mat1_rp_230=zeros(3,N,M);
for i=1:M
    for j=1:N
        for k=1:3
            if(k==1)
            mat1_rp_230(k,j,i)=x(j);
            end
            if(k==2)
             mat1_rp_230(k,j,i)=rp230(j,i);
            end
            if(k==3)
                mat1_rp_230(k,j,i)=zz230(i);
            end
        end
    end
end

% in ogni casella della matrice metto le coordinate di ogni punto nello spazio: se apro questa matirce trovo le coordinate in x y z
%P=192W
for i=1:N
    for j=1:M
        x11_192(i,j)=mat1_192(1,i,j);
        y11_192(i,j)=mat1_192(2,i,j);
        z11_192(i,j)=mat1_192(3,i,j);
    end
end

for i=1:N
    for j=1:M
        x11_rp_192(i,j)=mat1_rp_192(1,i,j);
        y11_rp_192(i,j)=mat1_rp_192(2,i,j);
        z11_rp_192(i,j)=mat1_rp_192(3,i,j);
    end
end

for i=1:N
    for j=1:M
        if(x11_192(i,j)<=xs+errx)
            if(x11_192(i,j)>xs-errx)
            indice=i;
        end
    end
    end
end 

for i=1:M
    if(zz192(i)<=zs+errz192)
            if(zz192(i)>zs-errz192)
            indicez192=i;
            end
    end
end

%plot curve traslate per diversi valori di velocità del laser
figure
grid on
plot3(x,y1,z1*ones(1,M),'k','LineWidth',2)
xlabel('carico [g]')
ylabel('durezza HV [MPa]')
zlabel('velocità laser [mm/s]')
hold on
plot3(x,y2,z2*ones(1,M),'k','LineWidth',2)
plot3(x,y3,z3*ones(1,M),'k','LineWidth',2)
title('Potenza laser 190 W')
grid on
%plot superficie
surf(x11_192,y11_192,z11_192) %in serf inserisco matrici di coordinate
title('Potenza laser 192 W')
xlabel('carico [g]')
ylabel('durezza HV [MPa]')
zlabel('velocità laser [mm/s]')
hold on
%plot value
plot3(mat1_192(1,indice,indicez192),mat1_192(2,indice,indicez192),mat1_192(3,indice,indicez192),'*','LineWidth',3,'Color','g')
xlabel('carico [g]')
ylabel('durezza HV [MPa]')
zlabel('velocità laser [mm/s]')
title ('Potenza laser 192 W')
legend('v1 = 400 mm/s','v2 = 600 mm/s','v3 = 800 mm/s','surf','value')

figure
grid on
plot3(x, rp1,z1*ones(1,M),'LineWidth',2)
xlabel('carico [g]')
ylabel('Rp02 [MPa]')
zlabel('velocità laser [mm/s]')
hold on
plot3(x,rp2,z2*ones(1,M),'LineWidth',3)
plot3(x,rp3,z3*ones(1,M),'LineWidth',2)
title('Potenza laser 190 W')
grid on
%plot superficie
surf(x11_rp_192,y11_rp_192,z11_192) %in serf inserisco matrici di coordinate
title('Potenza laser 192 W')
xlabel('carico [g]')
ylabel('Rp02 [MPa]')
zlabel('velocità laser [mm/s]')
hold on
%plot value
plot3(mat1_rp_192(1,indice,indicez192),mat1_rp_192(2,indice,indicez192),mat1_rp_192(3,indice,indicez192),'*','LineWidth',3,'Color','g')
xlabel('carico [g]')
ylabel('Rp02 [MPa]')
zlabel('velocità laser [mm/s]')
title ('Potenza laser 192 W')
legend('Rp02 [MPa] c1','Rp02 [MPa] c2','Rp02 [MPa] c3','surf','value')

durezza_192W_in_MPa =mat1_192(2,indice,indicez192)
Rp02_192W =mat1_rp_192(2,indice,indicez192)

%P=230W
%durezza
for i=1:N
    for j=1:M
        x11_230(i,j)=mat1_230(1,i,j);
        y11_230(i,j)=mat1_230(2,i,j);
        z11_230(i,j)=mat1_230(3,i,j);
    end
end

figure 
grid on
plot3(x,y4,z4*ones(1,M),'LineWidth',2)
xlabel('carico [g]')
ylabel('durezza HV [MPa]')
zlabel('velocità laser [mm/s]')
hold on
plot3(x,y5,z5*ones(1,M),'LineWidth',2)
plot3(x,y6,z6*ones(1,M),'LineWidth',2)
plot3(x,y7,z7*ones(1,M),'LineWidth',2)
title('Potenza laser 230 W')
grid on

surf(x11_230,y11_230,z11_230) %in serf inserisco matrici di coordinate
title('Potenza laser 230 W')
xlabel('carico [g]')
ylabel('durezza HV [MPa]')
zlabel('velocità laser [mm/s]')
hold on

for i=1:N
    for j=1:M
        if(x11_230(i,j)<=xs+errx)
            if(x11_230(i,j)>xs-errx)
            indice=i;
        end
    end
    end
end 

for i=1:M
    if(zz230(i)<=zs+errz230)
            if(zz230(i)>zs-errz230)
            indicez230=i;
            end
    end
end
durezza_230W_in_MPa=mat1_230(2,indice,indicez230)

%figure
plot3(mat1_230(1,indice,indicez230),mat1_230(2,indice,indicez230),mat1_230(3,indice,indicez230),'*','LineWidth',3,'Color','g')
xlabel('carico [g]')
ylabel('durezza HV [MPa]')
zlabel('velocità laser [mm/s]')
title ('Potenza laser 230 W')
legend('v4 = 500 mm/s','v5 = 600 mm/s','v6 = 700 mm/s','v7 = 900 mm/s','surf','value')

%Rp02
for i=1:N
    for j=1:M
        x11_c_230(i,j)=mat1_rp_230(1,i,j);
        y11_c_230(i,j)=mat1_rp_230(2,i,j);
        z11_c_230(i,j)=mat1_rp_230(3,i,j);
    end
end
% 
figure
plot3(x,rp4,z4*ones(1,M),'LineWidth',2)
xlabel('carico [g]')
ylabel('Rp02[Mpa]')
zlabel('velocità laser [mm/s]')
hold on
plot3(x,rp5,z5*ones(1,M),'LineWidth',2)
plot3(x,rp6,z6*ones(1,M),'LineWidth',2)
plot3(x,rp7,z7*ones(1,M),'LineWidth',2)
title('Potenza laser 230 W')
grid on
 
surf(x11_c_230,y11_c_230,z11_c_230) %in surf inserisco matrici di coordinate
title('Potenza laser 230 W')
xlabel('carico [g]')
ylabel('Rp02[Mpa]')
zlabel('velocità laser [mm/s]')
hold on

for i=1:N
    for j=1:M
        if(x11_c_230(i,j)<=xs+errx)
            if(x11_c_230(i,j)>xs-errx)
            indice=i;
        end
    end
    end
end 


for i=1:M
    if(zz230(i)<=zs+errz230)
            if(zz230(i)>zs-errz230)
            indicez230=i;
            end
    end
end
Rp02_230W=mat1_rp_230(2,indice,indicez230)

%figure
plot3(mat1_rp_230(1,indice,indicez230),mat1_rp_230(2,indice,indicez230),mat1_rp_230(3,indice,indicez230),'*','LineWidth',3,'Color','g')
xlabel('carico [g]')
ylabel('Rp02[MPa]')
zlabel('velocità laser [mm/s]')
title ('Potenza laser 230 W')
legend('Rp02[MPa] C4','Rp02[MPa] C5','Rp02[MPa] C6','Rp02[MPa] C7','surf','value')

%% ricerca coefficienti Rp02 = K x HV

% if provini == 1
%     SIGMA_SN
%     verifica_coeff_angol_192W = 730/durezza_192W_in_MPa
%     verifica_coeff_angol_230W = 730/durezza_230W_in_MPa
%     
% end
% if provini == 0
%     SIGMA_SN
%     verifica_coeff_angol_192W = 1303/durezza_192W_in_MPa
%     verifica_coeff_angol_230W = 1303/durezza_230W_in_MPa
% end
% 

% SIGMA_SN
%     verifica_coeff_angol_192W = SIGMA_SN/durezza_192W_in_MPa
%     verifica_coeff_angol_230W = SIGMA_SN/durezza_230W_in_MPa

%% ricerca valori di durezza


% xs=input('Inserire valore di carico tra 1g e 1kg: ');
% errx=(x(3)-x(2))/2;
% zs=input('Inserire valore di velocità: ');
% errz192=(zz192(3)-zz192(2))/2;
% errz230=(zz230(3)-zz230(2))/2;

%P=192W
% for i=1:N
%     for j=1:M
%         if(x11_192(i,j)<=xs+errx)
%             if(x11_192(i,j)>xs-errx)
%             indice=i;
%         end
%     end
%     end
% end 
% 
% 
% for i=1:M
%     if(zz192(i)<=zs+errz192)
%             if(zz192(i)>zs-errz192)
%             indicez192=i;
%             end
%     end
% end
% 
% ys192=mat1_192(2,indice,indicez192)
% 
% %figure
% plot3(mat1_192(1,indice,indicez192),mat1_192(2,indice,indicez192),mat1_192(3,indice,indicez192),'*','LineWidth',3,'Color','g')
% xlabel('carico [g]')
% ylabel('durezza HV [MPa]')
% zlabel('velocità laser [mm/s]')
% title ('Potenza laser 192 W')

%P=230W
% for i=1:N
%     for j=1:M
%         if(x11_230(i,j)<=xs+errx)
%             if(x11_230(i,j)>xs-errx)
%             indice=i;
%         end
%     end
%     end
% end 
% 
% 
% for i=1:M
%     if(zz230(i)<=zs+errz230)
%             if(zz230(i)>zs-errz230)
%             indicez230=i;
%             end
%     end
% end
% 
% ys230=mat1_230(2,indice,indicez230)
% 
% %figure
% plot3(mat1_230(1,indice,indicez230),mat1_230(2,indice,indicez230),mat1_230(3,indice,indicez230),'*','LineWidth',3,'Color','g')
% xlabel('carico [g]')
% ylabel('durezza HV [MPa]')
% zlabel('velocità laser [mm/s]')
% title ('Potenza laser 230 W')

%% prova su un'interpolazione di un valore di carico : vediamo come ha interpolato per valutare il grado del polinomio

% yy192_1 = yy192(80,:);
% figure
% plot(nodiy192(:,80), nodiz192(:,80),'*')  %punti notevoli delle 7 prove che già ho
% xlabel('[mm/s]')
% ylabel('durezza HV [MPa]')
% hold on
% pause
% plot(yy192_1,zz192)
% title ('Potenza laser 192 W')
% 
% yy230_1 = yy230(80,:);
% figure
% plot(nodiy230(:,80), nodiz230(:,80),'*')  %punti notevoli delle 7 prove che già ho
% xlabel('[mm/s]')
% ylabel('durezza HV [MPa]')
% hold on
% pause
% plot(yy230_1,zz230)
% title ('Potenza laser 230 W')
% 
% yy192_1 = yy192(80,:);
% figure
% plot(nodiy192(:,80), nodiz192(:,80),'*')  %punti notevoli delle 7 prove che già ho
% xlabel('[mm/s]')
% ylabel('durezza HV [MPa]')
% hold on
% pause
% plot(yy192_1,zz192)
% title ('Potenza laser 192 W')
% 
% yy230_1 = yy230(80,:);
% figure
% plot(nodiy230(:,80), nodiz230(:,80),'*')  %punti notevoli delle 7 prove che già ho
% xlabel('[mm/s]')
% ylabel('durezza HV [MPa]')
% hold on
% pause
% plot(yy230_1,zz230)
% title ('Potenza laser 230 W')
