% MODELLO MONODIMENSIONALE A PARAMETRI COSTANTI
% Spessore semi-infinito

clear all

P=50; % Potenza laser [W]
r=0.8; % Riflettività superficie
R0=0.2; % Raggio fascio [cm]
S=pi*R0^2; % Area sezione fascio [cm2]
F0=P*(1-r)/S; % Flusso incidente sul pezzo [W/cm2]
K=0.72; % Conducibilità termica di un acciaio a 25°C [W/cm°C]
ro=7.8; % Densità 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]
tau=12; % Tempo di riscaldamento [sec]

z=0:0.05:5; % z = profondità in cm
t=0:0.1:20; % t = tempo in sec

% Fase di riscaldamento

for i=1:101
   for j=1:tau*10+1
      if j==1
         T(i,j)=0;
      else
         D=2*(alfa*t(j))^0.5;
         ierfc=(1/pi)^0.5*exp(-(z(i)/D)^2)-(z(i)/D)*erfc(z(i)/D);
         T(i,j)=F0*D*ierfc/K;
      end   
   end
end

% Fase di raffreddamento

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

% Grafica

h=axes('Color',[.9,.9,.9],...
   'Fontname','times',...
   'FontAngle','italic',...
   'FontSize',8);
surf(t,z,T)
set(get(h,'Xlabel'),'String','Time [sec]',...
   'Fontname','times',...
   'FontAngle','italic',...
   'FontSize',10);
set(get(h,'Ylabel'),'String','Depth [cm]',...
   'Fontname','times',...
   'FontAngle','italic',...
   'FontSize',10);
set(get(h,'Zlabel'),'String','Temperature [°C]',...
   'Fontname','times',...
   'FontAngle','italic',...
   'FontSize',10);
set(get(h,'Title'),'String',['Laser power: ',num2str(P),' W; Reflection: ',num2str(r),'; Spot radius: ',num2str(R0),' cm; Heating time: ',num2str(tau),' sec.'],...
   'Fontname','arial',...
   'FontSize',10);
