function icrs_usd;
%Demonstration of photolithography and
%diffusion process steps on the example
%of making an integrated-circuit resisistor
%(Section 1.3.3 in Understanding Semiconductor Devices)

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


global con_exit_usd2; global con_go_usd2; 
global con_sl_brght_usd2; 
global con_pr_list_usd2;
con_go_usd2=0;

%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.51 0.02 0.25]; % 3 maroon     (Positive resist)
dopcolmap(4,:)=[0.51 0.25 0.4]; % 4 changed maroon (Exposed 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)
dopcolmap(9,:)=[0.55 0.22 0.79];%  9 Violet
dopcolmap(10,:)=[0 0 0]; %        10 Black
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','Section 1.3.3: Making an IC Resistor',...
   '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              
introtext ={'UNDERSTANDING SEMICONDUCTOR DEVICES ';...
      'Section 1.3.3: Making an IC Resistor' ;...
      ' ';...
      'DEMONSTRATION of photolithography and diffusion processing steps,';...
      'by example of making an integrated-circuit resistor. '};
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_go_usd2=1;');
con_go_usd2=0;
while con_go_usd2<1; 
      if ishandle(Hf)==0;
          return; 
      end
   pause(0.2);   
end
con_go_usd2=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);
Xox=1;
ymask1=2; ymask2=6; ymask3=8.25; ymask4=9; 
ymask5=1.75; ymask6=2.5; ymask7=5.5; ymask8=6.25; ymask9=8.25; ymask10=9; %Contacts
ymask11=1.5; ymask12=2.75; ymask13=5.25; ymask14=6.5; ymask15=8; ymask16=9.25; %Metal
N_PR_end=18;
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;
bright=14.4; %Brightness of doping colors - doping below log10(bright) is black
         
%THE PROCESS LIST
prlist(1,1)={' SUBSTRATE'};         
prlist(1,2)={'    N-type'};
prlist(1,3)={'    Nsub=1e+16 cm-3'};
prlist(2,1)={' INITIAL OXIDATION'};
prlist(3,1)={'      Photoresist application'};
prlist(3,2)={'         Spinning'};
prlist(3,3)={'         Soft-Baking'};
prlist(4,1)={'      Exposure'};
prlist(5,1)={'      Developing'};
prlist(5,2)={'         followed by hard-baking'};
prlist(6,1)={'      OXIDE ETCHING'};
prlist(7,1)={' FIRST PHOTOLITHOGRAPHY COMPLETE'};
prlist(7,2)={'      after resist stripping'};
prlist(8,1)={' P-TYPE DIFFUSION (Boron deposition)'};
prlist(8,2)={'      Temperature=900oC'};
prlist(8,3)={'      Time=2 min'};
prlist(9,1)={' DRIVE-IN '};
prlist(9,2)={'      Temperature=1100oC'};
prlist(9,3)={'      Time=120 min'};
prlist(10,1)={' ETCH MASKING OXIDE'};
prlist(11,1)={' DEPOSIT NEW MASKING OXIDE'};
prlist(12,1)={' SECOND PHOTOLITHOGRAPHY'};
prlist(12,2)={'      to open window for'};
prlist(12,3)={'      N-type diffusion'};
prlist(13,1)={' N-TYPE DIFFUSION (Phosporus deposition)'};
prlist(13,2)={'      Tepmperature=900oC'};
prlist(13,3)={'      Time=5 min'};
prlist(14,1)={' DRIVE-IN '};
prlist(14,2)={'      Temperature=1100oC'};
prlist(14,3)={'      Time=5 min'};
prlist(15,1)={' ETCH MASKING OXIDE'};
prlist(16,1)={' DEPOSIT ISOLATING OXIDE'};
prlist(17,1)={' THIRD PHOTOLITHOGRAPHY'};
prlist(17,2)={'      to open contact windows'};
prlist(18,1)={' METALLIZATION'};
prlist(19,1)={' FOURTH PHOTOLITHOGRAPHY'};
prlist(19,2)={'      Metal patterning '};
         
%SETTING THE CONTROL PANEL              
Hframe_out=uicontrol(Hf,'Units','normalized','Position',[0.025 0.025 0.95 0.325],...
   'Style','frame');
Hb_last=uicontrol(Hf,'Units','normalized','Position',[0.05 0.05 0.15 0.05],...
   'Style','pushbutton','String','LAST','CallBack','con_go_usd2=4;');
Hb_first=uicontrol(Hf,'Units','normalized','Position',[0.05 0.125 0.15 0.05],...
   'Style','pushbutton','String','FIRST','CallBack','con_go_usd2=3;','Enable','off');
Hb_previous=uicontrol(Hf,'Units','normalized','Position',[0.05 0.2 0.15 0.05],...
   'Style','pushbutton','String','PREVIOUS','CallBack','con_go_usd2=2;','Enable','off');
Hb_next=uicontrol(Hf,'Units','normalized','Position',[0.05 0.275 0.15 0.05],...
   'Style','pushbutton','String','NEXT','CallBack','con_go_usd2=1;');
Hb_exit=uicontrol(Hf,'Units','normalized','Position',[0.85 0.05 0.1 0.275],...
   'Style','pushbutton','String','CLOSE','CallBack','con_exit_usd2=1;');
H_list=uicontrol(Hf,'Style','listbox','Background',[0.8 0.8 0.8],'String',prlist(:,1),...
  'Units','normalized','Position',[0.25 0.05 0.55 0.25],'CallBack','con_pr_list_usd2=1;');
Ht_prlist=uicontrol(Hf,'Units','normalized','Position',[0.25 0.3 0.55 0.03],...
   'Style','text','String',['Click to select a process step'],...
                         'BackgroundColor',[0.8 0.8 0.8]);

Ht_msg=uicontrol(Hf,'Units','normalized','Position',[0.05 0.375 0.5 0.1],...
   'Style','text','String',prlist(1,:),'HorizontalAlignment','left');

Hsl_brght=uicontrol(Hf,'Units','normalized','Position',[0.7 0.4 0.25 0.04],...
   'Style','slider','min',6,'max',16,'Value',14.4,'CallBack','con_sl_brght_usd2=1;');
Ht_brght=uicontrol(Hf,'Units','normalized','Position',[0.7 0.44 0.25 0.03],...
   'Style','text','String','BRIGHTNESS');

%SET THE 2D GRAPH
%%The Grid
nxsteps=80; nysteps=150; %X is inside, Y is lateral
[gridx,gridy]=meshgrid(linspace(-2,4,nxsteps),linspace(0,8,nysteps));
for ii=1:16;
   gridx(:,ii)=-3+3/16*(ii-1);
end
for ii=17:nxsteps;
   gridx(:,ii)=(ii-17)*4/(nxsteps-17);
end
ymask1_i=round(ymask1/9.5*nysteps);
ymask2_i=round(ymask2/9.5*nysteps);
ymask3_i=round(ymask3/9.5*nysteps);
ymask4_i=round(ymask4/9.5*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):ymask3_i;
   gridy(ii,:)=(ii-ymask2_i)*(ymask3-ymask2)/(ymask3_i-ymask2_i)+ymask2;
end
for ii=(ymask3_i+1):ymask4_i;
   gridy(ii,:)=(ii-ymask3_i)*(ymask4-ymask3)/(ymask4_i-ymask3_i)+ymask3;
end
for ii=(ymask4_i+1):nysteps;
   gridy(ii,:)=(ii-ymask4_i)*(9.5-ymask4)/(nysteps-ymask4_i)+ymask4;
end

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

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

%%2D Graph
       if ishandle(Hf)==0;
          return; 
       end
Ha2=axes('Position',[0.025 0.50 0.95 0.475]);
%%axes(Ha2);
Hs=surface(gridx,gridy,z2,'EdgeColor','none');
...%,'FaceColor','interp'
set(Ha2,'Xlim',[-3 4],'Ylim',[0 9.5],'Zlim',[-z2_max z2_max],...
   'Color',[0.8 0.8 0.8],'Visible','off','CLim',[1 64],...
   'View',[90 90]);
%The meaning of 26.5 and 37.5 constants:
%Doping colors are between 11 and 64, 37.5 being in the middle (white), 
%11 the darkest green, and 64 the darkest red.
%26.5 is the difference between 64 (as well as 10) and the middle 37.5.
set(Hs,'CData',(material<=10).*material+...
   (material>10).*((abs(z2)>bright).*(z2-sign(z2).*bright)*26.5/(z2_max-bright)+37.5));

%***********************************
con_exit_usd2=0; con_go_usd2=0; con_sl_brght_usd2=0; con_pr_list_usd2=0;
NofL=0; N_PR_old=0; N_PR_new=0; 
while con_exit_usd2<1 %MAIN LOOP
   
%SET THE SLIDER AND THE RADIOBUTTON VALUES   
   while (con_go_usd2+con_pr_list_usd2+con_exit_usd2)<1;       
      %Can set the brightness here
      if con_sl_brght_usd2>0; 
         if ishandle(Hf)==0;
          return; 
         end
         bright=(get(Hsl_brght,'Value'));
         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_usd2=0;
      end   
        if ishandle(Hf)==0;
          return; 
        end
        pause(0.1);
   end % while - to set the slider and radiobutton values or view process prms.
      
if con_exit_usd2<1;
%NEW CROSS-SECTION (either con_go_usd2 or con_pr_list_usd2 is not zero).
   
%Analyze the request
if con_go_usd2>0; 
   switch con_go_usd2;
   case 1;
      N_PR_old=N_PR_new; N_PR_new=N_PR_new+1;
   case 2; 
      N_PR_old=N_PR_new; N_PR_new=N_PR_new-1;
   case 3;
      N_PR_old=N_PR_new; N_PR_new=0;
   case 4;
      N_PR_old=N_PR_new; N_PR_new=N_PR_end;
   end   
   
else %Must be that an item on the process list is clicked
   N_PR_old=N_PR_new;
   N_PR_new=get(H_list,'value')-1;
end   
   
%Enable/disable buttons
if N_PR_new==N_PR_end; 
       if ishandle(Hf)==0;
          return; 
       end
   set([Hb_next Hb_last], 'Enable','off');
   set([Hb_first Hb_previous], 'Enable','on');
else
       if ishandle(Hf)==0;
          return; 
       end
   set([Hb_first Hb_previous Hb_next Hb_last], 'Enable','on');
end
if N_PR_new==0; 
       if ishandle(Hf)==0;
          return; 
       end
   set([Hb_previous Hb_first], 'Enable','off');
   set([Hb_next Hb_last], 'Enable','on');
end   

%PERFORM THE CALCULATIONS

if (N_PR_new<N_PR_old)+(N_PR_new==0);
   %%material matrix
   material(1:nysteps,1:nxsteps)=1; %Reset the matrix
   material=material.*(1+(gridx>=-0.45).*36); %Silicon (index=37) for x>0

   %%doping matrix
   Nsub(1:nysteps,1:nxsteps)=1; %Reset the matrix
   Nsub=-(material>10)*1e16+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(H_list,'Value',1);
   set(Hs,'CData',(material<=10).*material+...
   (material>10).*((abs(z2)>bright).*(z2-sign(z2).*bright)*26.5/(z2_max-bright)+37.5));
   set(Ht_msg,'String',prlist(1,:));
   
   N_PR_old=0;
end

if N_PR_new>0;
       if ishandle(Hf)==0;
          return; 
       end
  set(H_list,'Value',(N_PR_new+1));
  set(Ht_msg,'String',prlist((N_PR_new+1),:));
for ii=(N_PR_old+1):N_PR_new; 
switch ii;
   
case 1; %INITIAL OXIDATION
   material(1:nysteps,1:nxsteps)=1; %Reset the matrix
   material=(gridx>=0)*37+(gridx<0); %Silicon (index=37) for x>0
   area=(gridx<0).*(gridx>=-0.8);
   material=area.*6+(1-area).*material;
        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));

case 2; %Photoresist application
   area=(gridx<-0.8).*(gridx>=-1.6);
   material=area.*3+(1-area).*material;
        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));
 
case 3; %Exposure
   area=(gridx<=-2).*(gridx>=-2.5).*((gridy<ymask1)+(gridy>ymask2));
   material=area.*10+(1-area).*material;  
      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));

      if con_go_usd2==1; 
         if ishandle(Hf)==0;
            return; 
         end
          pause(1); 
      end
   
   area=(gridy>=ymask1).*(gridy<=ymask2).*(gridx<-1.6);
   material=area.*9+(1-area).*material;  
   material=(gridx<-2.5).*9+(gridx>=-2.5).*material;
      area=(material==3).*(gridy>=ymask1).*(gridy<=ymask2);
   material=area.*4+(1-area).*material;
     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));
   
      if con_go_usd2==1; 
         if ishandle(Hf)==0;
          return; 
         end
         pause(1); 
      end
   
   material=(material==9)+(material~=9).*material;
      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));

case 4; %Developing
   material=(material==10)+(material~=10).*material;
   material=(material==4)+(material~=4).*material;
      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));

case 5; %Etching
   area=(material==6).*(gridy>=ymask1).*(gridy<=ymask2);
   material=area+(1-area).*material;
      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));
   
case 6; %Resist stripping
   material=(material==3)+(material~=3).*material; 
      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));

case 7; %Boron deposition
   NofL=1;
   %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; 
      
     %Calculate N1 
     time=120; 
     kT=kb*(900+273.15); %the thermal energy         
     D_Si(1)=D0(1,1)*exp(-EA(1,1)/kT);
     sigma_Si(1)=10000*sqrt(2*time*D_Si(1));
     Xnor1=(material>10).*X_Rp1/sigma_Si(1);
     Ynor1=(material>10).*Y_Rp1/(0.8*sigma_Si(1));
     N1=N0(1)*exp(-0.5*(Xnor1.*Xnor1+Ynor1.*Ynor1));
     
     %Plotting matrices
      ntot=N1+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));
     
     if con_go_usd2==1; 
        %Convert the air into doping source color 
      if ishandle(Hf)==0;
          return; 
      end
       area=get(Hs,'CData'); 
        set(Hs,'CData',(material==1).*64+(material~=1).*area);
        pause(1); 
        %Convert back to air
         if ishandle(Hf)==0;
          return; 
         end
        set(Hs,'CData',(material==1)+(material~=1).*area);
     end
     
     sigma_Si_0(1)=sigma_Si(1); 
     Phi(1)=sqrt(2/3.1415)*N0(1)*sigma_Si_0(1);
     
case 8; %Drive-in
   kT=kb*(1100+273.15); 
   time=120*60;
   D_Si(1)=D0(1,1)*exp(-EA(1,1)/kT);
   sigma_Si(1)=10000*sqrt(2*time*D_Si(1))+...
            sigma_Si_0(1); %10000 converts cm into um;
   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=Nmax(1)*exp(-0.5*(Xnor1.*Xnor1+Ynor1.*Ynor1)); 
   
   ntot=N1+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));
   
   sigma_Si_0(1)=sigma_Si(1);
   
case 9; %Etch masking oxide   
   if ishandle(Hf)==0;
          return; 
   end
   area=get(Hs,'CData');
   set(Hs,'CData',(material<=10)+(material>10).*area);   
   
case 10; %Deposit new masking oxide
   area=(gridx<0).*(gridx>=-0.8);
   material=area.*6+(1-area).*material;
      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));
   
case 11; %Second Photolithography
   if con_go_usd2==1; 
      %Apply photoresist
      area=(gridx<-0.8).*(gridx>=-1.6);
      material=area.*3+(1-area).*material;
      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));
      pause(1); 
      
      %Exposure
      area=(gridx<=-2).*(gridx>=-2.5).*((gridy<ymask3)+(gridy>ymask4));
      material=area.*10+(1-area).*material;  
      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));
      pause(1);
   
      area=(gridy>=ymask3).*(gridy<=ymask4).*(gridx<-1.6);
      material=area.*9+(1-area).*material;  
      material=(gridx<-2.5).*9+(gridx>=-2.5).*material;
      area=(material==3).*(gridy>=ymask3).*(gridy<=ymask4);
      material=area.*4+(1-area).*material;
      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));
      pause(1);
   
      material=(material==9)+(material~=9).*material;
      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));
      pause(1);
    
      %Developing
      material=(material==10)+(material~=10).*material;
      material=(material==4)+(material~=4).*material;
      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));
      pause(1);
      
      %Etching
      area=(material==6).*(gridy>=ymask3).*(gridy<=ymask4);
      material=area+(1-area).*material;
      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));
      pause(1);
     
      %Resist stripping
      material=(material==3)+(material~=3).*material;
      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));
      
   end
   area=(gridx<0).*(gridx>=-0.8).*(gridy>=ymask3).*(gridy<=ymask4);
   material=area+(1-area).*material; 
      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));
   
case 12; %PHOSPORUS DEPOSITION
   NofL=2;
   %Create X_Rp2 and Y_Rp2
      X_Rp2=gridx.*(material>10)+(material==6)*10;          
      Y_Rp2=gridy*0;
      for j=1:nysteps;
          if gridy(j,1)<=ymask3;
             Y_Rp2(j,:)=ymask3-gridy(j,1);
          else if gridy(j,1)<=ymask4;
             Y_Rp2(j,:)=0; %the 1D section
               else Y_Rp2(j,:)=gridy(j,1)-ymask4; 
               end
          end 
      end
      Y_Rp2=Y_Rp2.*(material>10)+(material==6)*10; 
      
     %Calculate N2 
     time=5; 
     kT=kb*(900+273.15); %the thermal energy         
     D_Si(2)=D0(1,2)*exp(-EA(1,2)/kT);
     sigma_Si(2)=10000*sqrt(2*time*D_Si(2));
     Xnor2=(material>10).*X_Rp2/sigma_Si(2);
     Ynor2=(material>10).*Y_Rp2/(0.8*sigma_Si(2));
     N2=-N0(2)*exp(-0.5*(Xnor2.*Xnor2+Ynor2.*Ynor2));
     
     %Plotting matrices
      ntot=N1+N2+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));
     
     if con_go_usd2==1; 
        %Convert the air into doping source color
        if ishandle(Hf)==0;
          return; 
        end
        area=get(Hs,'CData'); 
        set(Hs,'CData',(material==1).*11+(material~=1).*area);
        pause(1); 
        %Convert back to air
        if ishandle(Hf)==0;
          return; 
        end
        set(Hs,'CData',(material==1)+(material~=1).*area);
     end
     
     sigma_Si_0(2)=sigma_Si(2); 
     Phi(2)=sqrt(2/3.1415)*N0(2)*sigma_Si_0(2);
     
case 13; %Drive-in
   kT=kb*(1050+273.15); 
   time=5*60;
   D_Si(1)=D0(1,1)*exp(-EA(1,1)/kT);
   sigma_Si(1)=10000*sqrt(2*time*D_Si(1))+...
            sigma_Si_0(1); %10000 converts cm into um;
   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=Nmax(1)*exp(-0.5*(Xnor1.*Xnor1+Ynor1.*Ynor1)); 
   
   D_Si(2)=D0(1,2)*exp(-EA(1,2)/kT);
   sigma_Si(2)=10000*sqrt(2*time*D_Si(2))+...
            sigma_Si_0(2); %10000 converts cm into um;
   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=-Nmax(2)*exp(-0.5*(Xnor2.*Xnor2+Ynor2.*Ynor2));
   
   ntot=N1+N2+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));
   
   sigma_Si_0(2)=sigma_Si(2);
   
case 14; %Etch masking oxide   
   if ishandle(Hf)==0;
          return; 
   end
   area=get(Hs,'CData');
   set(Hs,'CData',(material<=10)+(material>10).*area);   
   
case 15; %Deposit new isolation oxide
   area=(gridx<0).*(gridx>=-0.8);
   material=area.*6+(1-area).*material;
       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));
   
case 16; %Third Photolithography
   if con_go_usd2==1; 
      %Apply photoresist
      area=(gridx<-0.8).*(gridx>=-1.6);
      material=area.*3+(1-area).*material;
      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));
      pause(1); 
      
      %Exposure
      area=(gridx<=-2).*(gridx>=-2.5).*((gridy<ymask5)+(gridy>ymask10)+...
         (gridy>ymask6).*(gridy<ymask7)+(gridy>ymask8).*(gridy<ymask9));
      material=area.*10+(1-area).*material;  
      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));
      pause(1);
   
      area=((gridy>=ymask5).*(gridy<=ymask6)+(gridy>=ymask7).*(gridy<=ymask8)+...
         (gridy>=ymask9).*(gridy<=ymask10)).*(gridx<-1.6);
      material=area.*9+(1-area).*material;  
      material=(gridx<-2.5).*9+(gridx>=-2.5).*material;
      area=(material==3).*((gridy>=ymask5).*(gridy<=ymask6)+(gridy>=ymask7).*(gridy<=ymask8)+...
         (gridy>=ymask9).*(gridy<=ymask10));
      material=area.*4+(1-area).*material;
      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));
      pause(1);
   
      material=(material==9)+(material~=9).*material;
      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));
      pause(1);
    
      %Developing
      material=(material==10)+(material~=10).*material;
      material=(material==4)+(material~=4).*material;
      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));
      pause(1);
      
      %Etching
      area=(material==6).*((gridy>=ymask5).*(gridy<=ymask6)+(gridy>=ymask7).*(gridy<=ymask8)+...
         (gridy>=ymask9).*(gridy<=ymask10));
      material=area+(1-area).*material;
      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));
      pause(1);
     
      %Resist stripping
      material=(material==3)+(material~=3).*material;
      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));
      
   end
   area=(gridx<0).*(gridx>=-0.8).*((gridy>=ymask5).*(gridy<=ymask6)+(gridy>=ymask7).*(gridy<=ymask8)+...
         (gridy>=ymask9).*(gridy<=ymask10));
      material=area+(1-area).*material;
      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));
   
case 17; %Metallization   
   area=(material==1).*(gridx>=-1.4);
      material=area.*2+(1-area).*material; 
      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));
   
case 18; %Fourth photolithography
   if con_go_usd2==1; 
      %Apply photoresist
      area=(gridx<-1.4).*(gridx>=-2.2);
      material=area.*3+(1-area).*material; 
      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));
      pause(1); 
      
      %Exposure
      area=(gridx<=-2.4).*(gridx>=-2.8).*((gridy>ymask11).*(gridy<ymask12)+...
         (gridy>ymask13).*(gridy<ymask14)+(gridy>ymask15).*(gridy<ymask16));
      material=area.*10+(1-area).*material;  
      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));
      pause(1);
   
      area=((gridy<ymask11)+(gridy>=ymask12).*(gridy<=ymask13)+(gridy>=ymask14).*...
         (gridy<=ymask15)+(gridy>ymask16)).*(gridx<-2.2);
      material=area.*9+(1-area).*material;  
      material=(gridx<-2.8).*9+(gridx>=-2.8).*material;
      area=(material==3).*((gridy<ymask11)+(gridy>=ymask12).*(gridy<=ymask13)+(gridy>=ymask14).*...
         (gridy<=ymask15)+(gridy>ymask16));
      material=area.*4+(1-area).*material; 
      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));
      pause(1);
      
      material=(material==9)+(material~=9).*material; 
      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));
      pause(1);
    
      %Developing
      material=(material==10)+(material~=10).*material;
      material=(material==4)+(material~=4).*material; 
      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));
      pause(1);
      
      %Etching
      area=(material==2).*((gridy<ymask11)+(gridy>ymask12).*(gridy<ymask13)+(gridy>ymask14).*...
         (gridy<ymask15)+(gridy>ymask16));
      material=area+(1-area).*material; 
      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));
      pause(1);
     
      %Resist stripping
      material=(material==3)+(material~=3).*material;
      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));
      
   end
   area=(gridx<-0.8).*(gridx>=-1.4).*((gridy<ymask11)+(gridy>ymask12).*(gridy<ymask13)+(gridy>ymask14).*...
         (gridy<ymask15)+(gridy>ymask16));
      material=area+(1-area).*material; 
      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));
   
  
end; %switch  
end; %for      
end; %if N_PR_new>0   
con_go_usd2=0; con_pr_list_usd2=0;
end; end % while and if of main loop (con_exit_usd2=1)
close

