K=0.72; % Conducibilità termica di un acciaio a 25°C [W/cm°C]
ro=7.8; % Peso specifico di un acciaio [g/cm3]
Cp=0.46; % Calore specifico di un acciaio a 25 °C [J/g°C]
alfa=K/(ro*Cp); % Diffusività termica del materiale [cm2/sec]
Tmax=1100; % Temperatura massima sopportabile dal materiale [°C]
Taust=750; % Temperatura di austenitizzazione del materiale [°C]
zprof=0.5; % Profondità di tempra [cm]

% Calcolo profondità normalizzata

delta0=100;
for zn=0:0.01:2
   ierfc=(1/pi)^0.5*exp(-zn^2)-zn*erfc(zn);
   Ttau=Tmax*pi^0.5*ierfc;
   deltaT=abs(Ttau-Taust);
   if deltaT<delta0
      delta0=deltaT;
      zn0=zn;
   end
end

% Calcolo tempo di interazione

tau=(1/(4*alfa))*(zprof/zn0)^2;

% Calcolo del flusso termico necessario

F0=Tmax*(pi)^0.5*K/(2*(alfa*tau)^0.5);

% Analisi curve di raffreddamento

t=tau:0.05:tau+5;

z=0;

D=2*(alfa*tau)^0.5;
ierfc=(1/pi)^0.5*exp(-(z/D)^2)-(z/D)*erfc(z/D);
T(1)=Tmax*pi^0.5*ierfc;

for i=2:101
   D=2*(alfa*t(i))^0.5;
   Dprimo=2*(alfa*(t(i)-tau))^0.5;
   ierfc=(1/pi)^0.5*exp(-(z/D)^2)-(z/D)*erfc(z/D);
   ierfcprimo=(1/pi)^0.5*exp(-(z/Dprimo)^2)-(z/Dprimo)*erfc(z/Dprimo);
   T(i)=F0*(D*ierfc-Dprimo*ierfcprimo)/K;
end

plot(t,T)
text(t(2),T(2),[num2str(z),' cm'])


z=zprof/2;

D=2*(alfa*tau)^0.5;
ierfc=(1/pi)^0.5*exp(-(z/D)^2)-(z/D)*erfc(z/D);
T(1)=Tmax*pi^0.5*ierfc;

for i=2:101
   D=2*(alfa*t(i))^0.5;
   Dprimo=2*(alfa*(t(i)-tau))^0.5;
   ierfc=(1/pi)^0.5*exp(-(z/D)^2)-(z/D)*erfc(z/D);
   ierfcprimo=(1/pi)^0.5*exp(-(z/Dprimo)^2)-(z/Dprimo)*erfc(z/Dprimo);
   T(i)=F0*(D*ierfc-Dprimo*ierfcprimo)/K;
end

xlabel('Time [sec]');
ylabel('Temperature [°C]');
hold on
plot(t,T)
text(t(2),T(2),[num2str(z),' cm'])

z=zprof;

D=2*(alfa*tau)^0.5;
ierfc=(1/pi)^0.5*exp(-(z/D)^2)-(z/D)*erfc(z/D);
T(1)=Tmax*pi^0.5*ierfc;

for i=2:101
   D=2*(alfa*t(i))^0.5;
   Dprimo=2*(alfa*(t(i)-tau))^0.5;
   ierfc=(1/pi)^0.5*exp(-(z/D)^2)-(z/D)*erfc(z/D);
   ierfcprimo=(1/pi)^0.5*exp(-(z/Dprimo)^2)-(z/Dprimo)*erfc(z/Dprimo);
   T(i)=F0*(D*ierfc-Dprimo*ierfcprimo)/K;
end

plot(t,T)
text(t(2),T(2),[num2str(z),' cm'])

% Curva CCT

ttt = [	6		5		4		3		2.5	2.4	2.3	2.2	2.1	2.05	2.1	2.2	2.3	2.4	2.5	3	4	5	6	];
TTT = [	390	410	430	450	460	465	470	480	490	510	530	540	545	550	555	575	610	630	650];
ttt=ttt+tau;

plot(ttt,TTT,'r--')
hold off
