function[]=dpng_usd(introtext);
%Semiconductor Doping Simulator 
%(Sections 1.3.4/5.2.1 in Understanding Semiconductor Devices)

%   Sima Dimitrijev, September 27, 1999.
%   Copyright (c) 2000, Oxford University Press.


global con_exit_usd3; global con_stop_usd3; global con_start_usd3; 
global con_sl_ener_usd3; global con_r_pr_usd3;
global con_sl_temp_usd3; 
global con_sl_sub_usd3; global con_r_el_usd3;
global con_sl_brght_usd3; global con_sl_time_usd3;
global con_pr_list_usd3;

%SET THE DOPING COLOR MAP
a=([1:64].')*0; dopcolmap=[a a a];
dopcolmap(1,:)=[0.95 0.95 0.95]; %    1  light grey  (Air)
dopcolmap(2,:)=[0.4 0.4 0.4]; %    2  grey (Aluminum)
dopcolmap(3,:)=[0.31 0.57 0.35]; % 3 dark green (Negative resist)
dopcolmap(4,:)=[0.72 0.68 0.35]; % 4 dark khaki (Positive resist)
dopcolmap(5,:)=[0.51 0.02 0.25]; % 5 purple     (Nitride)
dopcolmap(6,:)=[0.8 0.8 1.0]; %    6 blue       (Oxide)
dopcolmap(7,:)=[0.7 1 0.7]; %      7 greenish   (i/i Donor source)
dopcolmap(8,:)=[1 0.7 0.7]; %      8 redish     (i/i Acceptor source)
for ii=11:27; %the dark green part
   dopcolmap(ii,2)=0.4*(ii-11)/16+0.6;
end   
for ii=28:35; %the light green part  
   dopcolmap(ii,2)=1;
   dopcolmap(ii,1)=(ii-27)/8;
   dopcolmap(ii,3)=(ii-27)/8;
end
for ii=36:39; %The white part in the middle
dopcolmap(ii,1)=1; dopcolmap(ii,2)=1; dopcolmap(ii,3)=1;   
end

for ii=40:47; %the light red part
   dopcolmap(ii,2)=1-(ii-40)/8;
   dopcolmap(ii,3)=1-(ii-40)/8;
   dopcolmap(ii,1)=1;
end   
for ii=48:64; %the dark red part
   dopcolmap(ii,1)=0.4*(64-ii)/16+0.6;
end   

%OPEN PARENT FIGURE
Hf=figure('NumberTitle','off','Name','Sections 1.3.4/5.2.1: Semiconductor Doping Simulator',...
   'Color',[0.8 0.8 0.8]);
set(gcf,'InvertHardCopy','off');

Ha0=axes('Xlim',[0 10],'Ylim',[0 10],'Visible','off','Position',[0 0 1 1]);
set(Hf,'Colormap',dopcolmap);


%READ THE INTRO TEXT              
Ht_intro=text(0.1,0.6,char(introtext),'Units','normalized');
Hp_cont=uicontrol(Hf,'Units','normalized','Position',[0.1 0.2 0.7 0.1],...
   'Style','pushbutton','String','CLICK HERE TO CONTINUE!','CallBack','con_start_usd3=1;');
con_start_usd3=0;
while con_start_usd3<1; 
   if ishandle(Hf)==0;
          return; 
   end
   pause(0.2);   
end
con_start_usd3=0;
if ishandle(Hf)==0;
   return; 
end
set(Ht_intro,'Visible','off');
set(Hp_cont,'Visible','off','Enable','off');
drawnow;

%CONSTANTS
z2_max=log10(2e21); z2_min=log10(1e14);
Delta_t=60; time_max=180*60; %in [seconds]
Delta_Phi=1e11; Phi_max_lcurr=1e13; 
Delta_Phim=1e13; Phi_max_mcurr=1e15;
acc_fac=30; %Time acceleration factor
ymask1=2; ymask2=6; Xox=1;
Sb=1; %Needs to be established
N0=[4e20, 1e21, 1.8e21, Sb]; %Surface concentration of [B, P, As, Sb]
EA=[3.46, 3.75, 3.90, Sb ;% in Si
   3.46, 3.66, 3.53, Sb; % in poly-Si
   3.53, 4.44, 4.89, Sb; % in SiO2
   3.53, 4.44, 4.89, Sb]; % in Si3N4
D0=[0.76, 19.7, 6.0, Sb; %in Si
   210, 3000, 55, Sb; % in poly-Si
   3.16e-4, 7.23, 1.75, Sb; % in SiO2
   3.16e-6, 0.0723, 0.0175, Sb]; % in Si3N4
kb=8.625e-5;
%Rp_ii_coeff=a0+a1*E^2+a2*E^2
Rp_ii_coef=[27.5934, 2.6199, -3.23e-4; %a0, a1 and a2 for Boron
            3.0793, 1.24786, 6.0233e-4; % Phosphorus
            6.4916, 0.60656, -1.8822e-4; %Arsenic
            Sb, Sb, Sb ]; % Anthimony
         
st_ii_coef=[13.46, 0.8025, -2.162e-3; %a0, a1 and a2 for Boron
            3.0032, 0.5322, -5.5981e-4; % Phosphorus
            2.3082, 0.2245, -1.8085e-4; %Arsenic
            Sb, Sb, Sb ]; % Anthimony
bright=13.5; %Brightness of doping colors - doping below log10(bright) is white  

%SETTING THE CONTROL PANEL              
Hb_exit=uicontrol(Hf,'Units','normalized','Position',[0.82 0.02 0.16 0.10],...
   'Style','pushbutton','String','CLOSE','CallBack','con_exit_usd3=1;');
Hb_stop=uicontrol(Hf,'Units','normalized','Position',[0.91 0.18 0.07 0.05],...
   'Style','pushbutton','String','STOP','CallBack','con_stop_usd3=1;','Enable','off');
Hb_start=uicontrol(Hf,'Units','normalized','Position',[0.82 0.18 0.07 0.05],...
   'Style','pushbutton','String','START','CallBack','con_start_usd3=1;');
Hframe_out=uicontrol(Hf,'Units','normalized','Position',[0.815 0.24 0.17 0.55],...
   'Style','frame');
Hframe_in=uicontrol(Hf,'Units','normalized','Position',[0.815 0.66 0.17 0.13],...
   'Style','frame','BackgroundColor','k');
Hsl_ener=uicontrol(Hf,'Units','normalized','Position',[0.82 0.28 0.16 0.03],...
   'Style','slider','min',20,'max',200,'Value',50,'CallBack','con_sl_ener_usd3=1;');
Ht_ener=uicontrol(Hf,'Units','normalized','Position',[0.82 0.305 0.16 0.03],...
   'Style','text','String',[num2str(round(get(Hsl_ener,'Value'))) ' keV']);
Hr_pr(4)=uicontrol(Hf,'Units','normalized','Position',[0.82 0.34 0.16 0.03],...
   'Style','radio','String','Medium Current','CallBack','con_r_pr_usd3=4;','Value',0);
Hr_pr(3)=uicontrol(Hf,'Units','normalized','Position',[0.82 0.37 0.16 0.03],...
   'Style','radio','String','Low Current','CallBack','con_r_pr_usd3=3;','Value',0);
Ht_ion=uicontrol(Hf,'Units','normalized','Position',[0.82 0.40 0.16 0.03],...
   'Style','text','String','ION IMPLANTATION');
Hsl_temp=uicontrol(Hf,'Units','normalized','Position',[0.82 0.49 0.16 0.03],...
   'Style','slider','min',900,'max',1200,'Value',1000,'CallBack',...
   'con_sl_temp_usd3=1;');
Ht_temp=uicontrol(Hf,'Units','normalized','Position',[0.82 0.515 0.16 0.03],...
   'Style','text','String',[num2str(round(get(Hsl_temp,'Value'))) ' oC']);
Hr_pr(2)=uicontrol(Hf,'Units','normalized','Position',[0.82 0.55 0.16 0.03],...
   'Style','radio','String','Drive-In','CallBack','con_r_pr_usd3=2;','Value',0);
Hr_pr(1)=uicontrol(Hf,'Units','normalized','Position',[0.82 0.58 0.16 0.03],...
   'Style','radio','String','Deposition','CallBack','con_r_pr_usd3=1;','Value',1);
Ht_dif=uicontrol(Hf,'Units','normalized','Position',[0.82 0.61 0.16 0.03],...
   'Style','text','String','DIFFUSION');
Hsl_sub=uicontrol(Hf,'Units','normalized','Position',[0.82 0.68 0.16 0.03],...
   'Style','slider','min',15,'max',20,'Value',15,...
   'CallBack','con_sl_sub_usd3=1;','BackgroundColor',[0.83 0.63 0.09]);
Ht_sub2=uicontrol(Hf,'Units','normalized','Position',[0.82 0.71 0.16 0.03],...
   'Style','text','String',[num2str(round(10^(get(Hsl_sub,'Value')))) ' cm-3'],...
   'BackgroundColor','k','ForegroundColor','y');
Ht_sub1=uicontrol(Hf,'Units','normalized','Position',[0.82 0.75 0.16 0.03],...
   'Style','text','String','RESET/ SUBSTRATE','BackgroundColor','k',...
   'ForegroundColor','y');
Ht_proc=uicontrol(Hf,'Units','normalized','Position',[0.815 0.79 0.17 0.03],...
   'Style','text','String','PROCESSING');
Hr_el(4)=uicontrol(Hf,'Units','normalized','Position',[0.90 0.84 0.08 0.03],...
   'Style','radio','String','Sb','Enable','off','CallBack','con_r_el_usd3=4;',...
   'Value',0);
Hr_el(3)=uicontrol(Hf,'Units','normalized','Position',[0.82 0.84 0.08 0.03],...
   'Style','radio','String','As','CallBack','con_r_el_usd3=3;','Value',0);
Hr_el(2)=uicontrol(Hf,'Units','normalized','Position',[0.90 0.87 0.08 0.03],...
   'Style','radio','String','P','CallBack','con_r_el_usd3=2;','Value',0);
Hr_el(1)=uicontrol(Hf,'Units','normalized','Position',[0.82 0.87 0.08 0.03],...
   'Style','radio','String','B','CallBack','con_r_el_usd3=1;','Value',1);
Ht_dop=uicontrol(Hf,'Units','normalized','Position',[0.82 0.90 0.16 0.03],...
   'Style','text','String','DOPING ELEMENT');
Hsl_brght=uicontrol(Hf,'Units','normalized','Position',[0.55 0.35 0.21 0.03],...
   'Style','slider','min',6,'max',16,'Value',13.5,'CallBack','con_sl_brght_usd3=1;');
Ht_brght=uicontrol(Hf,'Units','normalized','Position',[0.55 0.385 0.21 0.03],...
   'Style','text','String','BRIGHTNESS');
Hsl_time=uicontrol(Hf,'Units','normalized','Position',[0.5 0.25 0.26 0.05],...
   'Style','slider','min',0,'max',time_max/60,'Value',0,'CallBack','con_sl_time_usd3=1;',...
   'Visible','off','Enable','off');
Ht_time=uicontrol(Hf,'Units','normalized','Position',[0.5 0.3 0.26 0.05],...
   'Style','text','String',['TIME=' num2str(round(get(Hsl_time,'Value'))) ' min'],...
   'Visible','off','FontSize',12,'FontWeight','bold');

C_Panel=... %Set of buttons/sliders/etc. making the Control Panel
[Hsl_ener Ht_ener Hr_pr Ht_ion Hsl_temp Ht_temp Ht_dif Hsl_sub Ht_sub2 Ht_sub1...
   Ht_proc Hr_el(1) Hr_el(2) Hr_el(3) Ht_dop]; 

%SET THE GRAPHS AND THE PROCESS LIST
%Structural level A - THE GRID
nxsteps=80; nysteps=120; %X is inside, Y is lateral
[gridx,gridy]=meshgrid(linspace(-2,4,nxsteps),linspace(0,8,nysteps));
for ii=1:8;
   gridx(:,ii)=-2+2/8*(ii-1);
end
for ii=9:nxsteps;
   gridx(:,ii)=(ii-9)*4/(nxsteps-9);
end
ymask1_i=round(ymask1/8*nysteps);
ymask2_i=round(ymask2/8*nysteps);
for ii=1:ymask1_i; 
   gridy(ii,:)=(ii-1)*ymask1/(ymask1_i-1);
end   
for ii=(ymask1_i+1):ymask2_i;
   gridy(ii,:)=(ii-ymask1_i)*(ymask2-ymask1)/(ymask2_i-ymask1_i)+ymask1;
end
for ii=(ymask2_i+1):nysteps;
   gridy(ii,:)=(ii-ymask2_i)*(8-ymask2)/(nysteps-ymask2_i)+ymask2;
end

%Structural level B - Default material matrix
material(1:nysteps,1:nxsteps)=1; %Create the matrix
material=material.*(1+(gridx>=0).*36); %Silicon (index=37) for x>0

%Structural level C - Default doping matrix
Nsub(1:nysteps,1:nxsteps)=1; %Create the matrix
Nsub=(material>10)*1e15+Nsub;
N1=Nsub*0; N2=Nsub*0; N3=Nsub*0; 
ntot=Nsub; %Total Concentration
z2=(material>10).*sign(ntot).*log10(abs(ntot));

%1D Graph
if ishandle(Hf)==0;
   return; 
end
x1=gridx(fix(nxsteps/2),:); y1=z2(fix(nxsteps/2),:); z1=0*x1; %Initial data
Ha1=axes('Position',[0.1 0.1 0.4 0.3]);
%axes(Ha1);
Hp(1)=surface([x1;x1],[abs(y1);abs(y1)],[z1;z1],'LineWidth',3,...
      'FaceColor','none','EdgeColor','flat');
set(Ha1,'Xlim',[0 4],'YLim',[z2_min z2_max],'ZLim',[z2_min z2_max],...
      'Color',[0.8 0.8 0.8],'FontAngle','italic','Box','off',...
      'CLim',[1 64]);
xlabel('Depth (um)'); ylabel('log (Concentration in cm^{-3})');
y1c=(abs(y1)>bright).*(y1-sign(y1).*bright)*26.5/(z2_max-bright)+37.5;  
set(Hp(1),'CData',[y1c;y1c]);   
%colorbar('vert');

%2D Graph
if ishandle(Hf)==0;
   return; 
end
Ha2=axes('Position',[0.1 0.50 0.66 0.396]);
%axes(Ha2);
Hs=surface(gridx,gridy,z2,'EdgeColor','none');
...%,'FaceColor','interp'
set(Ha2,'Xlim',[-2 4],'Ylim',[0 8],'Zlim',[-z2_max z2_max],...
   'Color',[0.8 0.8 0.8],'FontAngle','italic','CLim',[1 64],...
   'View',[90 90]);
ylabel('Lateral Distance (um)'); xlabel('Depth (um)');
set(Hs,'CData',(material<=10).*material+...
   (material>10).*((abs(z2)>bright).*(z2-sign(z2).*bright)*26.5/(z2_max-bright)+37.5));

%Process list
Ht_prlist=uicontrol(Hf,'Units','normalized','Position',[0.55 0.235 0.21 0.09],...
   'Style','text','String',['PROCESS LIST         ';...
                            'Click on process step';...
                            'to see the parameters'],...
                         'BackgroundColor',[0.8 0.8 0.8]);
prlist_num=0; prlist(prlist_num+1,:)={[num2str(prlist_num),' Substrate   ']};                       
%prlist=['0. Substrate         '];
H_list=uicontrol(Hf,'Style','listbox','Background',[0.8 0.8 0.8],'String',prlist,...
  'Units','normalized','Position',[0.55 0.02 0.21 0.21],'CallBack','con_pr_list_usd3=1;');

%H_msg=msgbox('0','Process Parameters');
%close(H_msg);

%***********************************
%Define the oxide mask
ar=((gridy<ymask1)+(gridy>ymask2)).*(gridx<0).*(gridx>-Xox);
material=(1-ar).*material+ar*6; 
if ishandle(Hf)==0;
   return; 
end
set(Hs,'Cdata',(1-ar).*get(Hs,'Cdata')+ar*6);
%***********************************

%***********************************
con_exit_usd3=0; con_stop_usd3=0; con_start_usd3=0; 
con_sl_ener_usd3=0; con_r_pr_usd3=0; con_sl_temp_usd3=0; 
con_sl_sub_usd3=0; con_sl_brght_usd3=0; con_sl_time_usd3=0;
con_pr_list_usd3=0;
element=1; NofL=0; process=1; 
pr_step={'0 Substrate';'    BORON';...
      ['    Nsub=' num2str(round(10^(get(Hsl_sub,'Value')))) ' cm-3']};
pr_seq{1}=pr_step; time=0;
while con_exit_usd3<1 %MAIN LOOP

%SET THE SLIDER AND THE RADIOBUTTON VALUES   
   while (con_start_usd3+con_exit_usd3)<1;       
      if con_r_el_usd3>0;
         if ishandle(Hf)==0;
           return;  
         end
         for ii=1:4
         set(Hr_el(ii),'Value',ii==con_r_el_usd3);
         end   
         element=con_r_el_usd3; 
         con_r_el_usd3=0;
      end   
      if con_r_pr_usd3>0; 
         if ishandle(Hf)==0;
            return; 
         end
         for ii=1:4
         set(Hr_pr(ii),'Value',ii==con_r_pr_usd3);
         end   
         process=con_r_pr_usd3; con_r_pr_usd3=0;
      end
      if con_sl_sub_usd3>0;
         if ishandle(Hf)==0;
            return; 
         end
         sub_conc=round(10^(get(Hsl_sub,'Value')));
         set(Ht_sub2,'String',[num2str(sub_conc) ' cm-3']);
         con_sl_sub_usd3=0;
         %Reset the graphs
         prlist_num=0;
         clear prlist;
         prlist={['0 Substrate   ']};
         set(H_list,'String',prlist,'Value',1);
         clear pr_step;
         switch element;
            case 1; elmn='    BORON';
            case 2; elmn='    PPHOSPHORUS';
            case 3; elmn='    ARSENIC';
            case 4; elmn='    ANTHIMONY';   
            end
         pr_step={'0 Substrate'; elmn;...
         ['    Nsub=' num2str(round(10^(get(Hsl_sub,'Value')))) ' cm-3']};
         pr_seq{1}=pr_step;
         NofL=0;
         pm=2*(element<2)-1;
         Nsub=pm*sub_conc*(material>10)+(material<=10).*Nsub;
         N1=Nsub*0; N2=Nsub*0; N3=Nsub*0; 
         ntot=Nsub; %Total Concentration
         z2=(material>10).*sign(ntot).*log10(abs(ntot));
         if ishandle(Hf)==0;
            return; 
         end
         set(Hs,'CData',(material<=10).*material+...
            (material>10).*((abs(z2)>bright).*(z2-sign(z2).*bright)*26.5/(z2_max-bright)+37.5));;
         y1=z2(fix(nxsteps/2),:);
         y1c=(abs(y1)>bright).*(y1-sign(y1).*bright)*26.5/(z2_max-bright)+37.5;
         set(Hp(1),'YData',[abs(y1);abs(y1)],'CData',[y1c;y1c]);
         drawnow;
      end
      if con_sl_temp_usd3>0; 
         if ishandle(Hf)==0;
            return; 
         end
         temperature=round(get(Hsl_temp,'Value'));
         set(Ht_temp,'String',[num2str(temperature) ' oC']);
         con_sl_temp_usd3=0;
      end
      if con_sl_ener_usd3>0;
         if ishandle(Hf)==0;
            return; 
         end
         energy=round(get(Hsl_ener,'Value'));
         set(Ht_ener,'String',[num2str(energy) ' keV']);
         con_sl_ener_usd3=0;
      end
      if con_sl_brght_usd3>0;
         if ishandle(Hf)==0;
            return; 
         end
         bright=(get(Hsl_brght,'Value'));
         y1c=(abs(y1)>bright).*(y1-sign(y1).*bright)*26.5/(z2_max-bright)+37.5;   
         set(Hp(1),'CData',[y1c;y1c]);
         set(Hs,'CData',(material<=10).*material+...
            (material>10).*((abs(z2)>bright).*(z2-sign(z2).*bright)*26.5/(z2_max-bright)+37.5));
         con_sl_brght_usd3=0;
      end   
      if con_pr_list_usd3>0;
         con_pr_list_usd3=0;
         i_pr_seq=get(H_list,'Value');
%         msgbox(pr_seq{i_pr_seq},'Process Parameters');
         H_msg=msgbox(pr_seq{i_pr_seq},'Process Parameters');
      end   
           if ishandle(Hf)==0;
            return; 
           end
      pause(0.1);
   end % while - to set the slider and radiobutton values or view process prms.
   
   con_start_usd3=0;
         if ishandle(Hf)==0;
            return; 
         end
   set(Hb_start,'Enable','off'); set(Hb_stop,'Enable','on');
   set([Hsl_brght H_list],'Enable','off','Visible','off');
   set([Ht_prlist Ht_brght],'Visible','off');
   drawnow;
   
   
if con_exit_usd3<1;
%START IS PRESSED    
         if ishandle(Hf)==0;
            return; 
         end
   set(Hb_exit,'Enable','off');
   set(C_Panel,'Enable','off');   
      
  %%SETTING THE STRUCTURE/PARAMETERS ACCORDING TO THE SELECTED BUTTONS/SLIDERS
  prlist_num=prlist_num+1; %Display ERROR
  switch process;
  case 1; %Pre-deposition diff. - Add a diffusion layer
  pr_item=' Depos. diff.';
         NofL=NofL+1; 
         if NofL>3; %DISPLAY ERROR
            errordlg({'Cannot have more than 3';'deposition diffusion and/or';...
                      'ion implant process steps.'},'STOP');   
         end          
      switch NofL %To set the new layer
      case 1;
      %Structural level C - create X_Rp1 and Y_Rp1
      X_Rp1=gridx.*(material>10)+(material==6)*10;          
             
      Y_Rp1=gridy*0;
      for j=1:nysteps;
          if gridy(j,1)<=ymask1;
             Y_Rp1(j,:)=ymask1-gridy(j,1);
          else if gridy(j,1)<=ymask2;
             Y_Rp1(j,:)=0; %the 1D section
               else Y_Rp1(j,:)=gridy(j,1)-ymask2; 
               end
          end 
      end
      Y_Rp1=Y_Rp1.*(material>10)+(material==6)*10; 
      
         %Structural level C 
         elm(1)=element;
         kT=kb*(round(get(Hsl_temp,'Value'))+273.15); %the thermal energy         
         D_Si(1)=D0(1,element)*exp(-EA(1,element)/kT);
         D_SiO2(1)=D0(3,element)*exp(-EA(1,element)/kT);
         time=0; %in seconds
         sigma_Si_0(1)=0.001; 
         
         Nmax(1)=N0(element);
         pm(1)=2*(element<2)-1;
                          
      case 2; 
      %Structural level C - create X_Rp1 and Y_Rp1
      X_Rp2=gridx.*(material>10)+(material==6)*10;          
             
      Y_Rp2=gridy*0;
      for j=1:nysteps;
          if gridy(j,1)<=ymask1;
             Y_Rp2(j,:)=ymask1-gridy(j,1);
          else if gridy(j,1)<=ymask2;
             Y_Rp2(j,:)=0; %the 1D section
               else Y_Rp2(j,:)=gridy(j,1)-ymask2; 
               end
          end 
      end
      Y_Rp2=Y_Rp2.*(material>10)+(material==6)*10; 
      
         %Structural level C 
         elm(2)=element;
         kT=kb*(round(get(Hsl_temp,'Value'))+273.15); %the thermal energy         
         D_Si(2)=D0(1,element)*exp(-EA(1,element)/kT);
         time=0; %in seconds
         sigma_Si_0(2)=0.001; 
         
         Nmax(2)=N0(element);
         pm(2)=2*(element<2)-1;
      case 3; 
      %Structural level C - create X_Rp1 and Y_Rp1
      X_Rp3=gridx.*(material>10)+(material==6)*10;          
             
      Y_Rp3=gridy*0;
      for j=1:nysteps;
          if gridy(j,1)<=ymask1;
             Y_Rp3(j,:)=ymask1-gridy(j,1);
          else if gridy(j,1)<=ymask2;
             Y_Rp3(j,:)=0; %the 1D section
               else Y_Rp3(j,:)=gridy(j,1)-ymask2; 
               end
          end 
      end
      Y_Rp3=Y_Rp3.*(material>10)+(material==6)*10; 
      
         %Structural level C 
         elm(3)=element;
         kT=kb*(round(get(Hsl_temp,'Value'))+273.15); %the thermal energy         
         D_Si(3)=D0(1,element)*exp(-EA(1,element)/kT);
         time=0; %in seconds
         sigma_Si_0(3)=0.001; 
         
         Nmax(3)=N0(element);
         pm(3)=2*(element<2)-1;
      end  %Switch NofL 
    
   case 2; %Drive-in diff.
   pr_item=' Drive-in    ';   
      %Change the temperature and the assoc. papameters)
      if NofL>0;  
      for ii=1:NofL;   
         kT=kb*(round(get(Hsl_temp,'Value'))+273.15); %the thermal energy         
         D_Si(ii)=D0(1,elm(ii))*exp(-EA(1,elm(ii))/kT);
         time=0; %in seconds
      end   
      end %if NofL>0
         
   case 3; %low current I/I - Add a doping layer
   pr_item=' Ion implant ';  
         NofL=NofL+1; 
         if NofL>3; %DISPLAY ERROR
         errordlg({'Cannot have more than 3';'deposition diffusion and/or';...
                      'ion implant process steps.'},'STOP');   
         end 
         
      switch NofL %To set the new layer
      case 1;
      %Structural level C - create X_Rp1 and Y_Rp1
      energy=get(Hsl_ener,'Value'); %the i/i energy
      Rp=0.001*(Rp_ii_coef(element,1)+Rp_ii_coef(element,2)*energy...
         +Rp_ii_coef(element,3)*energy^2);
      
      X_Rp1=(gridx-Rp).*(material>10)+(material==6)*10;          
            
      Y_Rp1=gridy*0;
      for j=1:nysteps;
          if gridy(j,1)<=ymask1;
             Y_Rp1(j,:)=ymask1-gridy(j,1);
          else if gridy(j,1)<=ymask2;
             Y_Rp1(j,:)=0; %the 1D section
               else Y_Rp1(j,:)=gridy(j,1)-ymask2; 
               end
          end 
      end
      Y_Rp1=Y_Rp1.*(material>10)+(material==6)*10; 
      
         %Structural level C 
         elm(1)=element;
         sigma_Si(1)=0.001*(st_ii_coef(element,1)+st_ii_coef(element,2)*energy...
                       +st_ii_coef(element,3)*energy^2);
         Phi(1)=0; %in seconds
         pm(1)=2*(element<2)-1;
         
         Xnor1=(material>10).*X_Rp1/sigma_Si(1);
         Ynor1=(material>10).*Y_Rp1/(0.8*sigma_Si(1));

         
      case 2; 
         
         %Structural level C - create X_Rp2 and Y_Rp2
      energy=get(Hsl_ener,'Value'); %the i/i energy
      Rp=0.001*(Rp_ii_coef(element,1)+Rp_ii_coef(element,2)*energy...
         +Rp_ii_coef(element,3)*energy^2);
         
      X_Rp2=(gridx-Rp).*(material>10)+(material==6)*10;          
             
      Y_Rp2=gridy*0;
      for j=1:nysteps;
          if gridy(j,1)<=ymask1;
             Y_Rp2(j,:)=ymask1-gridy(j,1);
          else if gridy(j,1)<=ymask2;
             Y_Rp2(j,:)=0; %the 1D section
               else Y_Rp2(j,:)=gridy(j,1)-ymask2; 
               end
          end 
      end
      Y_Rp2=Y_Rp2.*(material>10)+(material==6)*10; 
      
         %Structural level C 
         elm(2)=element;
         sigma_Si(2)=0.001*(st_ii_coef(element,1)+st_ii_coef(element,2)*energy...
                       +st_ii_coef(element,3)*energy^2);
         Phi(2)=0; %in seconds
         pm(2)=2*(element<2)-1;
         
         Xnor2=(material>10).*X_Rp2/sigma_Si(2);
         Ynor2=(material>10).*Y_Rp2/(0.8*sigma_Si(2));

      case 3; 
      %Structural level C - create X_Rp3 and Y_Rp3
      energy=get(Hsl_ener,'Value'); %the i/i energy
      Rp=0.001*(Rp_ii_coef(element,1)+Rp_ii_coef(element,2)*energy...
         +Rp_ii_coef(element,3)*energy^2);
         
      X_Rp3=(gridx-Rp).*(material>10)+(material==6)*10;          
             
      Y_Rp3=gridy*0;
      for j=1:nysteps;
          if gridy(j,1)<=ymask1;
             Y_Rp3(j,:)=ymask1-gridy(j,1);
          else if gridy(j,1)<=ymask2;
             Y_Rp3(j,:)=0; %the 1D section
               else Y_Rp3(j,:)=gridy(j,1)-ymask2; 
               end
          end 
      end
      Y_Rp3=Y_Rp3.*(material>10)+(material==6)*10; 
      
         %Structural level C 
         elm(3)=element;
         sigma_Si(3)=0.001*(st_ii_coef(element,1)+st_ii_coef(element,2)*energy...
                       +st_ii_coef(element,3)*energy^2);
         Phi(3)=0; %in seconds
         pm(3)=2*(element<2)-1;
         
         Xnor3=(material>10).*X_Rp3/sigma_Si(3);
         Ynor3=(material>10).*Y_Rp3/(0.8*sigma_Si(3));
         
      end %switch NofL   
      
      case 4; %medium current I/I - Add a doping layer
   pr_item=' Ion implant ';  
         NofL=NofL+1; 
         if NofL>3; %DISPLAY ERROR
         errordlg({'Cannot have more than 3';'deposition diffusion and/or';...
                      'ion implant process steps.'},'STOP');   
         end 
         
      switch NofL %To set the new layer
      case 1;
      %Structural level C - create X_Rp1 and Y_Rp1
      energy=get(Hsl_ener,'Value'); %the i/i energy
      Rp=0.001*(Rp_ii_coef(element,1)+Rp_ii_coef(element,2)*energy...
         +Rp_ii_coef(element,3)*energy^2);
      
      X_Rp1=(gridx-Rp).*(material>10)+(material==6)*10;          
            
      Y_Rp1=gridy*0;
      for j=1:nysteps;
          if gridy(j,1)<=ymask1;
             Y_Rp1(j,:)=ymask1-gridy(j,1);
          else if gridy(j,1)<=ymask2;
             Y_Rp1(j,:)=0; %the 1D section
               else Y_Rp1(j,:)=gridy(j,1)-ymask2; 
               end
          end 
      end
      Y_Rp1=Y_Rp1.*(material>10)+(material==6)*10; 
      
         %Structural level C 
         elm(1)=element;
         sigma_Si(1)=0.001*(st_ii_coef(element,1)+st_ii_coef(element,2)*energy...
                       +st_ii_coef(element,3)*energy^2);
         Phi(1)=0; %in seconds
         pm(1)=2*(element<2)-1;
         
         Xnor1=(material>10).*X_Rp1/sigma_Si(1);
         Ynor1=(material>10).*Y_Rp1/(0.8*sigma_Si(1));

         
      case 2; 
         
         %Structural level C - create X_Rp2 and Y_Rp2
      energy=get(Hsl_ener,'Value'); %the i/i energy
      Rp=0.001*(Rp_ii_coef(element,1)+Rp_ii_coef(element,2)*energy...
         +Rp_ii_coef(element,3)*energy^2);
         
      X_Rp2=(gridx-Rp).*(material>10)+(material==6)*10;          
             
      Y_Rp2=gridy*0;
      for j=1:nysteps;
          if gridy(j,1)<=ymask1;
             Y_Rp2(j,:)=ymask1-gridy(j,1);
          else if gridy(j,1)<=ymask2;
             Y_Rp2(j,:)=0; %the 1D section
               else Y_Rp2(j,:)=gridy(j,1)-ymask2; 
               end
          end 
      end
      Y_Rp2=Y_Rp2.*(material>10)+(material==6)*10; 
      
         %Structural level C 
         elm(2)=element;
         sigma_Si(2)=0.001*(st_ii_coef(element,1)+st_ii_coef(element,2)*energy...
                       +st_ii_coef(element,3)*energy^2);
         Phi(2)=0; %in seconds
         pm(2)=2*(element<2)-1;
         
         Xnor2=(material>10).*X_Rp2/sigma_Si(2);
         Ynor2=(material>10).*Y_Rp2/(0.8*sigma_Si(2));

      case 3; 
      %Structural level C - create X_Rp3 and Y_Rp3
      energy=get(Hsl_ener,'Value'); %the i/i energy
      Rp=0.001*(Rp_ii_coef(element,1)+Rp_ii_coef(element,2)*energy...
         +Rp_ii_coef(element,3)*energy^2);
         
      X_Rp3=(gridx-Rp).*(material>10)+(material==6)*10;          
             
      Y_Rp3=gridy*0;
      for j=1:nysteps;
          if gridy(j,1)<=ymask1;
             Y_Rp3(j,:)=ymask1-gridy(j,1);
          else if gridy(j,1)<=ymask2;
             Y_Rp3(j,:)=0; %the 1D section
               else Y_Rp3(j,:)=gridy(j,1)-ymask2; 
               end
          end 
      end
      Y_Rp3=Y_Rp3.*(material>10)+(material==6)*10; 
      
         %Structural level C 
         elm(3)=element;
         sigma_Si(3)=0.001*(st_ii_coef(element,1)+st_ii_coef(element,2)*energy...
                       +st_ii_coef(element,3)*energy^2);
         Phi(3)=0; %in seconds
         pm(3)=2*(element<2)-1;
         
         Xnor3=(material>10).*X_Rp3/sigma_Si(3);
         Ynor3=(material>10).*Y_Rp3/(0.8*sigma_Si(3));
         
      end %switch NofL
      
   end %switch process
   
   %%PERFORMING THE ANIMATION/4 switch branches 
   if NofL<4;
   switch process;
   case 1; %Deposition diffusion 
      if ishandle(Hf)==0;
         return; 
      end
      set(Hsl_time,'Visible','on','Enable','on',...
         'min',0,'max',time_max/60,'Value',0); %Time indicator
      set(Ht_time,'Visible','on',...
         'String',['TIME=' num2str(round(get(Hsl_time,'Value'))) ' min']);
   
   if element>1 dop_sour_ind=11; %Index for the dopins source color
   else dop_sour_ind=64;
   end  
   material=(material==1).*dop_sour_ind+(material>1).*material; %Converts air to doping source
      
   while ((con_stop_usd3<1)*(time<time_max)); %CASE/PROCESS 1 animation loop itself
   tic; 
    
      time=time+Delta_t;       
           
      %Set/show the time
      if con_sl_time_usd3>0; 
      time=round(get(Hsl_time,'Value'))*60;
      con_sl_time_usd3=0;
      end
      if ishandle(Hf)==0;
         return; 
      end
      set(Hsl_time,'Value',time/60); 
      set(Ht_time,'String',['TIME=' num2str(round(get(Hsl_time,'Value'))) ' min']);
      
      %Perform the calculations as the time is increased by Delta_t
      for ii=1:NofL
         sigma_Si(ii)=10000*sqrt(2*time*D_Si(ii))+...
            sigma_Si_0(ii); %10000 converts cm into um;
      switch ii;
      case 1;   
         Xnor1=(material>10).*X_Rp1/sigma_Si(1);
         Ynor1=(material>10).*Y_Rp1/(0.8*sigma_Si(1));
         if NofL>1; %This step appears as drive-in for this layer
         Nmax(1)=sqrt(2/3.1415)*Phi(1)/sigma_Si(1); 
         end
         N1=pm(1)*Nmax(1)*exp(-0.5*(Xnor1.*Xnor1+Ynor1.*Ynor1));
      case 2;
         Xnor2=(material>10).*X_Rp2/sigma_Si(2);
         Ynor2=(material>10).*Y_Rp2/(0.8*sigma_Si(2));      
         if NofL>2; %This step appears as drive-in for this layer
         Nmax(2)=sqrt(2/3.1415)*Phi(2)/sigma_Si(2); 
         end
         N2=pm(2)*Nmax(2)*exp(-0.5*(Xnor2.*Xnor2+Ynor2.*Ynor2));
      case 3;
         Xnor3=(material>10).*X_Rp3/sigma_Si(3);
         Ynor3=(material>10).*Y_Rp3/(0.8*sigma_Si(3));      
         N3=pm(3)*Nmax(3)*exp(-0.5*(Xnor3.*Xnor3+Ynor3.*Ynor3));
      end %switch ii
      end %for NofL
      %Ploting
      
      ntot=N1+N2+N3+Nsub; %Total Concentration
      z2=(material>10).*sign(ntot).*log10(abs(ntot));
      z2=(material>10).*sign(ntot).*log10(abs(ntot));
      if ishandle(Hf)==0;
         return; 
      end
      set(Hs,'CData',(material<=10).*material+...
        (material>10).*((abs(z2)>bright).*(z2-sign(z2).*bright)*26.5/(z2_max-bright)+37.5));;
      y1=z2(fix(nxsteps/2),:);
      y1c=(abs(y1)>bright).*(y1-sign(y1).*bright)*26.5/(z2_max-bright)+37.5;
      set(Hp(1),'YData',[abs(y1);abs(y1)],'CData',[y1c;y1c]);
      drawnow;
   %Time control
   while toc<Delta_t/acc_fac;
      pause(0.1);
   end   
   if ishandle(Hf)==0;
     return; 
   end
   end %while of animation loop
   material=(material==dop_sour_ind)+(material~=dop_sour_ind).*material;
   %Converts doping source to air
   if ishandle(Hf)==0;
      return; 
   end
   set(Hs,'CData',(material<=10).*material+...
        (material>10).*((abs(z2)>bright).*(z2-sign(z2).*bright)*26.5/(z2_max-bright)+37.5));;
   drawnow; 
   % Remember sigmas, to add them up at the forthcoming high-T steps
   % Also, calculate the dose
   for ii=1:NofL; 
   sigma_Si_0(ii)=sigma_Si(ii); 
   Phi(ii)=sqrt(2/3.1415)*Nmax(ii)*sigma_Si_0(ii);
   end
   %Create pr_seq item
            switch element;
            case 1; elmn='    BORON';
            case 2; elmn='    PHOSPHORUS';
            case 3; elmn='    ARSENIC';
            case 4; elmn='    ANTHIMONY';   
            end
            clear pr_step;
            pr_step{2}=elmn;
            pr_step{3}=['    TEMP=' num2str(round(get(Hsl_temp,'Value'))) ' oC'];
            pr_step{4}=['    TIME=' num2str(time/60) ' min'];  
            pr_step{5}=['        Dose=' num2str(Phi(NofL)/10000) 'cm-2'];
            pr_step{6}=['        Nmax=' num2str(Nmax(NofL)) 'cm-3'];
               
      case 2; %Drive-in diffusion
         if ishandle(Hf)==0;
         return; 
         end
      set(Hsl_time,'Visible','on','Enable','on',...
         'min',0,'max',time_max/60,'Value',0); %Time indicator
      set(Ht_time,'Visible','on',...
         'String',['TIME=' num2str(round(get(Hsl_time,'Value'))) ' min']);
   while ((con_stop_usd3<1)*(time<time_max)); %CASE/PROCESS 2 animation loop itself
   tic; 
   
      time=time+Delta_t;       
           
      %Set/show the time
      if con_sl_time_usd3>0; 
      time=round(get(Hsl_time,'Value'))*60;
      con_sl_time_usd3=0;
      end
      if ishandle(Hf)==0;
         return; 
      end
      set(Hsl_time,'Value',time/60); 
      set(Ht_time,'String',['TIME=' num2str(round(get(Hsl_time,'Value'))) ' min']);
      
      %Perform the calculations as the time is increased by Delta_t
      for ii=1:NofL
         sigma_Si(ii)=10000*sqrt(2*time*D_Si(ii))+...
            sigma_Si_0(ii); %10000 converts cm into um;
      
      switch ii;
      case 1;   
         Xnor1=(material>10).*X_Rp1/sigma_Si(1);
         Ynor1=(material>10).*Y_Rp1/(0.8*sigma_Si(1));
         Nmax(1)=sqrt(2/3.1415)*Phi(1)/sigma_Si(1);       
         N1=pm(1)*Nmax(1)*exp(-0.5*(Xnor1.*Xnor1+Ynor1.*Ynor1));
      case 2;
         Xnor2=(material>10).*X_Rp2/sigma_Si(2);
         Ynor2=(material>10).*Y_Rp2/(0.8*sigma_Si(2));
         Nmax(2)=sqrt(2/3.1415)*Phi(2)/sigma_Si(2); 
         N2=pm(2)*Nmax(2)*exp(-0.5*(Xnor2.*Xnor2+Ynor2.*Ynor2));
      case 3;
         Xnor3=(material>10).*X_Rp3/sigma_Si(3);
         Ynor3=(material>10).*Y_Rp3/(0.8*sigma_Si(3));
         Nmax(3)=sqrt(2/3.1415)*Phi(3)/sigma_Si(3); 
         N3=pm(3)*Nmax(3)*exp(-0.5*(Xnor3.*Xnor3+Ynor3.*Ynor3));
      end %switch ii
      end %for NofL
      %Ploting
      
      ntot=N1+N2+N3+Nsub; %Total Concentration
      z2=(material>10).*sign(ntot).*log10(abs(ntot));
      z2=(material>10).*sign(ntot).*log10(abs(ntot));
      if ishandle(Hf)==0;
         return; 
      end
      set(Hs,'CData',(material<=10).*material+...
        (material>10).*((abs(z2)>bright).*(z2-sign(z2).*bright)*26.5/(z2_max-bright)+37.5));;
      y1=z2(fix(nxsteps/2),:);
      y1c=(abs(y1)>bright).*(y1-sign(y1).*bright)*26.5/(z2_max-bright)+37.5;
      set(Hp(1),'YData',[abs(y1);abs(y1)],'CData',[y1c;y1c]);
      drawnow;
   %Time control
   while toc<Delta_t/acc_fac;
      pause(0.1);
   end   
     if ishandle(Hf)==0;
         return; 
     end
   end %while of animation loop
      
   % Remember st's as st_p's, to add them up at the forthcoming high-T steps
   % Also, calculate the dose
   for ii=1:NofL; 
   sigma_Si_0(ii)=sigma_Si(ii); 
   end
   
   %Create pr_seq item
   clear pr_step;
   pr_step{2}=['    TEMP=' num2str(round(get(Hsl_temp,'Value'))) ' oC'];
   pr_step{3}=['    TIME=' num2str(time/60) ' min'];  
      
   case 3; %Low-current implant
      if ishandle(Hf)==0;
         return; 
      end
      set(Hsl_time,'Visible','on','Enable','on',...
         'min',0,'max',Phi_max_lcurr,'Value',0); %Dose indicator
      set(Ht_time,'Visible','on',...
         'String',['DOSE=' num2str(round(get(Hsl_time,'Value')/1e11)/10) 'e+12 cm-2']);
   
   if element>1 dop_sour_ind=7; %Index for the dopins source color
   else dop_sour_ind=8;
   end  
   material=(material==1).*dop_sour_ind+(material>1).*material; %Converts air to doping source
      
   while ((con_stop_usd3<1)*(Phi(NofL)<(Phi_max_lcurr-Delta_Phi))); %CASE/PROCESS 3 animation loop itself
   tic; 
    
      Phi(NofL)=Phi(NofL)+Delta_Phi;       
      
      
      %Set/show the time
      if con_sl_time_usd3>0; 
      Phi(NofL)=get(Hsl_time,'Value');
      con_sl_time_usd3=0;
      end 
      if ishandle(Hf)==0;
         return; 
      end
      set(Hsl_time,'Value',Phi(NofL)); 
      set(Ht_time,'String',['DOSE=' num2str(round(get(Hsl_time,'Value')/1e11)/10) 'e+12 cm-2']);
      
      %Perform the calculations as the dose is increased by Delta_Phi
      switch NofL;
      case 1;   
         Nmax(1)=sqrt(2/3.1415)*Phi(1)*1e4/sigma_Si(1); 
         N1=pm(1)*Nmax(1)*exp(-0.5*(Xnor1.*Xnor1+Ynor1.*Ynor1));
      case 2;
         Nmax(2)=sqrt(2/3.1415)*Phi(2)*1e4/sigma_Si(2); 
         N2=pm(2)*Nmax(2)*exp(-0.5*(Xnor2.*Xnor2+Ynor2.*Ynor2));
      case 3;
         Nmax(3)=sqrt(2/3.1415)*Phi(3)*1e4/sigma_Si(3);      
         N3=pm(3)*Nmax(3)*exp(-0.5*(Xnor3.*Xnor3+Ynor3.*Ynor3));
      end %switch NofL
      
      %Ploting
      
      ntot=N1+N2+N3+Nsub; %Total Concentration
      z2=(material>10).*sign(ntot).*log10(abs(ntot));
      z2=(material>10).*sign(ntot).*log10(abs(ntot));
      if ishandle(Hf)==0;
         return; 
      end
      set(Hs,'CData',(material<=10).*material+...
        (material>10).*((abs(z2)>bright).*(z2-sign(z2).*bright)*26.5/(z2_max-bright)+37.5));;
      y1=z2(fix(nxsteps/2),:);
      y1c=(abs(y1)>bright).*(y1-sign(y1).*bright)*26.5/(z2_max-bright)+37.5;
      set(Hp(1),'YData',[abs(y1);abs(y1)],'CData',[y1c;y1c]);
      drawnow;
   %Time control
   while toc<1;
      pause(0.1);
   end   
   if ishandle(Hf)==0;
      return; 
   end
   end %while of animation loop
   material=(material==dop_sour_ind)+(material~=dop_sour_ind).*material;
   %Converts doping source to air
   if ishandle(Hf)==0;
      return; 
   end
   set(Hs,'CData',(material<=10).*material+...
        (material>10).*((abs(z2)>bright).*(z2-sign(z2).*bright)*26.5/(z2_max-bright)+37.5));;
   drawnow; 
   % Remember sigmas, to add them up at the forthcoming high-T steps
   % and normalize the dose
   sigma_Si_0(NofL)=sigma_Si(NofL); 
   Phi(NofL)=1e4*Phi(NofL);   
   
   %Create pr_seq item
            switch element;
            case 1; elmn='    BORON';
            case 2; elmn='    PHOSPHORUS';
            case 3; elmn='    ARSENIC';
            case 4; elmn='    ANTHIMONY';   
            end
            clear pr_step;
            pr_step{2}=elmn;
            pr_step{3}=['    ENERGY=' num2str(energy) ' keV'];  
            pr_step{4}=['    DOSE=' num2str(Phi(NofL)/10000) 'cm-2'];
            pr_step{5}=['        X_av=' num2str(Rp) 'um'];
            pr_step{6}=['        sigma_X=' num2str(sigma_Si(NofL)) 'um'];

   
    case 4; %Medium current implant
      if ishandle(Hf)==0;
         return; 
      end
      set(Hsl_time,'Visible','on','Enable','on',...
         'min',0,'max',Phi_max_mcurr,'Value',0); %Dose indicator
      set(Ht_time,'Visible','on',...
         'String',['DOSE=' num2str(round(get(Hsl_time,'Value')/1e13)/10) 'e+14 cm-2']);
   
   if element>1 dop_sour_ind=7; %Index for the dopins source color
   else dop_sour_ind=8;
   end  
   material=(material==1).*dop_sour_ind+(material>1).*material; %Converts air to doping source
      
   while ((con_stop_usd3<1)*(Phi(NofL)<(Phi_max_mcurr-Delta_Phim))); %CASE/PROCESS 4 animation loop itself
   tic; 
    
      Phi(NofL)=Phi(NofL)+Delta_Phim;       
      
      
      %Set/show the time
      if con_sl_time_usd3>0; 
      Phi(NofL)=get(Hsl_time,'Value');
      con_sl_time_usd3=0;
      end
      if ishandle(Hf)==0;
         return; 
      end
      set(Hsl_time,'Value',Phi(NofL)); 
      set(Ht_time,'String',['DOSE=' num2str(round(get(Hsl_time,'Value')/1e13)/10) 'e+14 cm-2']);
      
      %Perform the calculations as the dose is increased by Delta_Phi
      switch NofL;
      case 1;   
         Nmax(1)=sqrt(2/3.1415)*Phi(1)*1e4/sigma_Si(1); 
         N1=pm(1)*Nmax(1)*exp(-0.5*(Xnor1.*Xnor1+Ynor1.*Ynor1));
      case 2;
         Nmax(2)=sqrt(2/3.1415)*Phi(2)*1e4/sigma_Si(2); 
         N2=pm(2)*Nmax(2)*exp(-0.5*(Xnor2.*Xnor2+Ynor2.*Ynor2));
      case 3;
         Nmax(3)=sqrt(2/3.1415)*Phi(3)*1e4/sigma_Si(3);      
         N3=pm(3)*Nmax(3)*exp(-0.5*(Xnor3.*Xnor3+Ynor3.*Ynor3));
      end %switch NofL
      
      %Ploting
      
      ntot=N1+N2+N3+Nsub; %Total Concentration
      z2=(material>10).*sign(ntot).*log10(abs(ntot));
      z2=(material>10).*sign(ntot).*log10(abs(ntot));
      if ishandle(Hf)==0;
         return; 
      end
      set(Hs,'CData',(material<=10).*material+...
        (material>10).*((abs(z2)>bright).*(z2-sign(z2).*bright)*26.5/(z2_max-bright)+37.5));;
      y1=z2(fix(nxsteps/2),:);
      y1c=(abs(y1)>bright).*(y1-sign(y1).*bright)*26.5/(z2_max-bright)+37.5;
      set(Hp(1),'YData',[abs(y1);abs(y1)],'CData',[y1c;y1c]);
      drawnow;
   %Time control
   while toc<1;
      pause(0.1);
   end   
   if ishandle(Hf)==0;
      return; 
   end
   end %while of animation loop
   material=(material==dop_sour_ind)+(material~=dop_sour_ind).*material;
   %Converts doping source to air
   if ishandle(Hf)==0;
      return; 
   end
   set(Hs,'CData',(material<=10).*material+...
        (material>10).*((abs(z2)>bright).*(z2-sign(z2).*bright)*26.5/(z2_max-bright)+37.5));;
   drawnow; 
   % Remember sigmas, to add them up at the forthcoming high-T steps
   % and normalize the dose
   sigma_Si_0(NofL)=sigma_Si(NofL); 
   Phi(NofL)=1e4*Phi(NofL);   
   
   %Create pr_seq item
            switch element;
            case 1; elmn='    BORON';
            case 2; elmn='    PHOSPHORUS';
            case 3; elmn='    ARSENIC';
            case 4; elmn='    ANTHIMONY';   
            end
            clear pr_step;
            pr_step{2}=elmn;
            pr_step{3}=['    ENERGY=' num2str(energy) ' keV'];  
            pr_step{4}=['    DOSE=' num2str(Phi(NofL)/10000) 'cm-2'];
            pr_step{5}=['        X_av=' num2str(Rp) 'um'];
            pr_step{6}=['        sigma_X=' num2str(sigma_Si(NofL)) 'um'];

   end %Switch  
end; %If NofL<4 
   con_stop_usd3=0;
      if ishandle(Hf)==0;
         return; 
      end
   set(Hb_stop,'Enable','off');
   set(Hb_exit,'Enable','on'); set(Hb_start,'Enable','on'); 
   set(Hsl_time,'Visible','off','Enable','off'); 
   set(Ht_time,'Visible','off');
   set(C_Panel,'Enable','on');
   set(Hsl_brght,'Enable','on','Visible','on');
   set([Ht_prlist Ht_brght],'Visible','on');
   if NofL<4;
      prlist(prlist_num+1,:)={[num2str(prlist_num),pr_item]};
      pr_step{1}=[num2str(prlist_num),pr_item]; %pr_step{2}, etc. defined earlier
      pr_seq{prlist_num+1}=pr_step;
   end   
      set(H_list,'String',prlist,'Enable','on','Visible','on');
   if NofL>3; NofL=3; prlist_num=prlist_num-1; end

end; end % while and if of main loop (con_exit_usd3=1)
close
