function mosc_usd;
%2D space and 1D energy-band illustrations/animations
%accompany C-V curve of a MOS capacitor.
%(Sections 2.5.3-2.5.5 in Understanding Semiconductor Devices)

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


global uf_usd11; global con_vg_usd11; 
global con_close_usd11; global con_doping_usd11;
global con_ox_usd11; global con_gate_usd11;
global con_chrg_usd11; global con_lfhf_usd11;
global con_stp_usd11; global con_strt_usd11;
global con_br_usd11;


%SET THE DOPING COLOR MAP
a=([1:64].')*0; dopcolmap=[a a a];
dopcolmap(1,:)=[0.8 0.8 0.8]; %    1  light grey  (Air)
dopcolmap(2,:)=[0.4 0.4 0.4]; %    2  grey (Aluminum)
dopcolmap(3,:)=[0.8 0.8 0.8]; % 3 grey not to be seen
dopcolmap(4,:)=[0.6 0.6 0.6]; % 4 grey for the E-B diagram
dopcolmap(5,:)=[0.51 0.02 0.25]; % 5 purple     (Nitride)
dopcolmap(6,:)=[0.8 0.8 1.0]; %    6 light 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.7*(ii-11)/16+0.3;
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.7*(64-ii)/16+0.3;
end   

%OPEN PARENT FIGURE
Hf=figure('NumberTitle','off','Name','Sections 2.5.3-2.5.5: MOS Capacitor',...
   'Color',[0.8 0.8 0.8]);
set(gcf,'InvertHardCopy','off');

set(gca,'Visible','off','Xlim',[0 15],'Ylim',[0 7.5],'Position',[0.01 0.15 .75 .825]);
set(Hf,'Colormap',dopcolmap);

%READ THE INTRO TEXT              
introtext ={'UNDERSTANDING SEMICONDUCTOR DEVICES ';...
      'Sections 2.5.3-2.5.5: MOS Capacitor' ;...
      ' ';...
      '2D-space and energy-band diagrams accompany the C-V curves of';...
      'a MOS capacitor to illustrate: (1) different modes of operation and ';...
      '(2) the effects of different technological parameters (substrate ';
      'doping, oxide thickness, type of gate, and oxide charge). ';...
      'The red and green colors indicate P-type and N-type areas, respectively ';...
      '(darker the color, higher the concentration). The white area indicates  ';...
      'the depletion layer. NOTE: Depletion/accumulation layers in the gate ';...
      'are too narrow, and therefore, not visible. ';...
      '          ';...
      'Slow-motion view of the effects of LF- or HF-AC signal, applied to the ';...
      'capacitor, is available. Note that the HF-AC signal, as opposed to the ';...
      'LF-AC signal, cannot change the charge in the inversion layer. ';...
      'Accordingly, the depletion-layer width and charge have to change as a ';....
      'response to the HF-AC signal.'}; 

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



con_vg_usd11=0; con_close_usd11=0; con_doping_usd11=0; con_ox_usd11=0; con_gate_usd11=0;
con_chrg_usd11=0; con_hflf_usd11=0; con_stp_usd11=0; con_strt_usd11=0;
%SETTING THE CONTROL PANEL              
Hframe=uicontrol(Hf,'Units','normalized','Position',[0.77 0.01 0.22 0.98],...
   'Style','frame');
%Vgq
HframeVg=uicontrol(Hf,'Units','normalized','Position',[0.775 0.875 0.21 0.11],...
   'Style','frame','BackgroundColor','y');
Hsl_vg=uicontrol(Hf,'Units','normalized','Position',[0.78 0.93 0.20 0.04],...
   'Style','slider','min',-5,'max',5,'Value',0,'CallBack','con_vg_usd11=1;');
Ht_vgd=uicontrol(Hf,'Units','normalized','Position',[0.78 0.88 0.20 0.05],...
   'Style','text','String',['Vg=' num2str(round(get(Hsl_vg,'Value')*10)/10) ' V'],...
   'BackgroundColor','y','ForegroundColor','k','FontSize',10,'FontWeight','bold');
Vgq=round(get(Hsl_vg,'Value')*10)/10;

%Doping
Ht_dop=uicontrol(Hf,'Units','normalized','Position',[0.78 0.82 0.20 0.03],...
   'Style','text','String','SUBSTRATE DOPING');
Hr_dop(1)=uicontrol(Hf,'Units','normalized','Position',[0.78 0.78 0.10 0.03],...
   'Style','radio','String','P type','CallBack','con_doping_usd11=1;','Value',1);
Hr_dop(2)=uicontrol(Hf,'Units','normalized','Position',[0.88 0.78 0.10 0.03],...
   'Style','radio','String','N type','CallBack','con_doping_usd11=2;','Value',0);
Hframe_dop=uicontrol(Hf,'Units','normalized','Position',[0.775 0.69 0.21 0.09],...
   'Style','frame','BackgroundColor','r');
Hsl_dop=uicontrol(Hf,'Units','normalized','Position',[0.78 0.73 0.20 0.04],...
   'Style','slider','min',16,'max',19,'Value',17,'CallBack','con_doping_usd11=3;');
Ht_dop=uicontrol(Hf,'Units','normalized','Position',[0.78 0.70 0.20 0.03],...
   'Style','text','String',['Na=' num2str(round(10^(get(Hsl_dop,'Value')))) ' cm-3'],...
   'BackgroundColor','r','ForegroundColor','w');
Na=round(10^(get(Hsl_dop,'Value')));
Nd=0;

%Oxide thickness
Ht_ox=uicontrol(Hf,'Units','normalized','Position',[0.78 0.64 0.20 0.03],...
   'Style','text','String','OXIDE THICKNESS');
Hsl_ox=uicontrol(Hf,'Units','normalized','Position',[0.78 0.60 0.20 0.04],...
   'Style','slider','min',3,'max',50,'Value',10,'CallBack','con_ox_usd11=1;');
Ht_ox=uicontrol(Hf,'Units','normalized','Position',[0.78 0.57 0.20 0.03],...
   'Style','text','String',['Tox=' num2str(round(get(Hsl_ox,'Value'))) ' nm']);
tox=round(get(Hsl_ox,'Value'))*1e-9;

%Gate material
Ht_gate=uicontrol(Hf,'Units','normalized','Position',[0.78 0.52 0.20 0.03],...
   'Style','text','String','GATE MATERIAL');
Hr_gate(1)=uicontrol(Hf,'Units','normalized','Position',[0.78 0.48 0.10 0.03],...
   'Style','radio','String','N+ poly','CallBack','con_gate_usd11=1;','Value',1);
Hr_gate(2)=uicontrol(Hf,'Units','normalized','Position',[0.88 0.48 0.10 0.03],...
   'Style','radio','String','P+ poly','CallBack','con_gate_usd11=2;','Value',0);
Hr_gate(3)=uicontrol(Hf,'Units','normalized','Position',[0.78 0.45 0.10 0.03],...
   'Style','radio','String','Metal','CallBack','con_gate_usd11=3;','Value',0);
Hsl_gate=uicontrol(Hf,'Units','normalized','Position',[0.78 0.41 0.20 0.04],...
   'Style','slider','min',4,'max',6,'Value',4.1,'CallBack','con_gate_usd11=4;',...
   'Enable','off');
Ht_gate=uicontrol(Hf,'Units','normalized','Position',[0.78 0.38 0.20 0.03],...
   'Style','text','String',['Work function=4.05 eV']);

%Oxide charge
Ht_chrgu=uicontrol(Hf,'Units','normalized','Position',[0.78 0.33 0.20 0.03],...
   'Style','text','String','OXIDE CHARGE');
Hsl_chrg=uicontrol(Hf,'Units','normalized','Position',[0.78 0.28 0.20 0.04],...
   'Style','slider','min',0,'max',90,'Value',0,'CallBack','con_chrg_usd11=1;');
Ht_chrgd=uicontrol(Hf,'Units','normalized','Position',[0.78 0.24 0.20 0.03],...
   'Style','text','String',['Nox=0']);

%LF-HF radiobuttons
Hr_lfhf(1)=uicontrol(Hf,'Units','normalized','Position',[0.54 0.31 0.10 0.03],...
   'Style','radio','String','LF','CallBack','con_lfhf_usd11=1;','Value',1,...
   'BackgroundColor',[0.8 0.8 0.8]);
Hr_lfhf(2)=uicontrol(Hf,'Units','normalized','Position',[0.64 0.31 0.10 0.03],...
   'Style','radio','String','HF','CallBack','con_lfhf_usd11=2;','Value',0,...
   'BackgroundColor',[0.8 0.8 0.8]);

%CLOSE
Hb_close=uicontrol(Hf,'Units','normalized','Position',[0.78 0.02 0.20 0.10],...
   'Style','pushbutton','String','CLOSE','CallBack','con_close_usd11=1;');

%AC signal
Hframe_ac=uicontrol(Hf,'Units','normalized','Position',[0.53 0.01 0.22 0.23],...
   'Style','frame');
Hb_stop=uicontrol(Hf,'Units','normalized','Position',[0.54 0.02 0.20 0.07],...
   'Style','pushbutton','String','STOP','CallBack','con_stp_usd11=1;',...
   'Enable','off');
Hb_start=uicontrol(Hf,'Units','normalized','Position',[0.54 0.11 0.20 0.07],...
   'Style','pushbutton','String','START','CallBack','con_strt_usd11=1;');
Ht_ac=uicontrol(Hf,'Units','normalized','Position',[0.54 0.19 0.20 0.04],...
   'Style','text','String','LF AC SIGNAL');

%Brightness
Ht_bru=uicontrol(Hf,'Units','normalized','Position',[0.15 0.37 0.20 0.03],...
   'Style','text','String','brightness','BackgroundColor',[0.8 0.8 0.8]);
Hsl_br=uicontrol(Hf,'Units','normalized','Position',[0.15 0.34 0.20 0.03],...
   'Style','slider','min',13,'max',15.5,'Value',14.1,'CallBack','con_br_usd11=1;');
bright=get(Hsl_br,'Value');


%CONSTANTS
kt=0.02586;
ni=1.02e10;
epsox=3.9*8.85e-12;
epssi=11.8*8.85e-12;
q=1.6e-19;
lfhf=1; %LF active; lfhf=0 - HF active
Cox=epsox/tox;
Nd=0; %Has to be 0 if Na>0; to be adjusted by a slider
Qox=0;
Npoly=-20; %log of concentration; Npoly=0 for metals
phim=4.05;
z2_max=log10(2e21); z2_min=log10(1e14);
if Na>0; 
   uf_usd11=log(Na/ni);
else
   uf_usd11=-log(Nd/ni);
end   
   Vfb=phim-(4.05+1.12/2+uf_usd11*kt)-Qox/Cox;
ld=sqrt(epssi*kt/(2*q*ni*1e6));
%the grid
Xmax=250;
Ymax=400; ygu=300; ygl=Ymax-ygu;  
scal=0.2;
bright=get(Hsl_br,'Value'); %While color for doping <10^bright
Vac=0.2;
Nsteps=8; %Number of steps in the AC signal
num_x=round((Xmax-tox*1e9-20)*scal); %Number of lateral points - 3...
    %(three lines outside silicon will be added)
num_yhalf=round((Ymax-ygu)*scal); %Number of vertical points above the gate
[xs,ys]=meshgrid(([0:num_x].^2/num_x),num_yhalf:-1:1);
%Set the grid inside the silicon
xs=xs/scal; ys=ys/scal+ygu; %The grid above ygu
ysl=ys-ygu;
ys=[ys' [ys(1,:)*0+Ymax/2]' ysl']';
xs=[xs' [xs(1,:)]' xs']';
%Set the grid points outside the silicon
gatethick=40;
xa(:,1)=xs(:,1)*0;
xa(:,2)=xa(:,1)+gatethick; %The metal/oxide interface
xs=xs+tox*1e9+xa(1,2); %%Shift the silicon grid inside the silicon...
  % to the oxide/Si interface
x=[xa xs]; %Merge the ambient and silicon x-grid points
y=[ys(:,1) ys(:,1) ys]; %Expand the y-grid points to include the ambient

%Establish the ends of the gate at grid lines
ygu=ys((num_yhalf),1);
ygl=ys((num_yhalf+2),1);

%Set material matrix 
material=[(xa(:,1)*0+1) (xa(:,2)*0+6 ) (xs*0+20)]; %air=1, oxide=6, substrate=20
material=(material==1).*(y>=ygl).*(y<=ygu).*11+... %Poly-gate=11
   (material~=1).*material;

%Set "distance from the gate interface" matrix;
rg=((x-x(1,3)).*(y>=ygl).*(y<=ygu)+...
      (y<ygl).*((x-x(1,3)).^2+(ygl-y).^2).^0.5+...
      (y>ygu).*((x-x(1,3)).^2+(y-ygu).^2).^0.5).*(material>11);


%SET THE GRAPHS

%C-V graph
   if ishandle(Hf)==0;
      return; 
   end;
Ha_cv=axes('Xlim',[-5 5],'Ylim',[0 4],'Position',[0.45 0.48 .3 .5],...
   'Color',[0.8 0.8 0.8],'YAxisLocation','left','FontAngle','italic');
hold on;
xlabel('GATE VOLTAGE (V)'); ylabel('CAPACITANCE (mF/mm^{2})');
Vgq=0;
Clf=1; Chf=1;
Hp_cvhf=plot(Vgq,Chf,'k-','LineWidth',2,'Visible','off');
set(Hp_cvhf,'Color',[0.55 0.55 0.55]);
Hp_cvlf=plot(Vgq,Clf,'k-','LineWidth',2,'Visible','off');
Hp_cvq=plot(Vgq,Clf,'y.','MarkerSize',20,'Visible','off');
hold off;

%CALCULATE C-V CURVE  %%%%%NEED TO ADD Vfb
%Calculate the Clf-vs-us
us=[uf_usd11-25:0.5:uf_usd11+25];
fus=sqrt(exp(uf_usd11)*(exp(-us)+us-1)+exp(-uf_usd11)*(exp(us)-us-1));
if Na>0;
   negl=-0.1; posl=0.1;
   usacc=(us<negl).*us+(us>=negl);
   usinv=(us>posl).*us+(us<=posl);
   weffacc=sign(us).*ld.*(2*fus./(exp(uf_usd11)*(1-exp(-usacc))+exp(-uf_usd11)*(exp(usacc)-1)));
   weff0=sqrt(2)*ld/sqrt(exp(uf_usd11)+exp(-uf_usd11));
   weffinv=sign(us).*ld.*(2*fus./(exp(uf_usd11)*(1-exp(-usinv))+exp(-uf_usd11)*(exp(usinv)-1)));
   weff=weffacc.*(us<negl)+weff0.*(us>=negl).*(us<=posl)+weffinv.*(us>posl);
else
   negl=-0.1; posl=0.1;
   usacc=(us>posl).*us+(us<=posl);
   usinv=(us<negl).*us+(us>=negl);
   weffacc=sign(us).*ld.*(2*fus./(exp(uf_usd11)*(1-exp(-usacc))+exp(-uf_usd11)*(exp(usacc)-1)));
   weff0=sqrt(2)*ld/sqrt(exp(uf_usd11)+exp(-uf_usd11));
   weffinv=sign(us).*ld.*(2*fus./(exp(uf_usd11)*(1-exp(-usinv))+exp(-uf_usd11)*(exp(usinv)-1)));
   weff=weffacc.*(us>posl)+weff0.*(us>=negl).*(us<=posl)+weffinv.*(us<negl);
end
Clf=Cox./(1+(epsox*weff./(epssi*tox)));

%Calculate Vg-vs-us
Vg=kt*(us+sign(us).*(epssi*tox/(epsox*ld)).*fus)+Vfb;

%Calculate Chf
[Cmin,imin]=min(Clf);
weffmax=weff(imin);
Vgmin=Vg(imin);
nv=length(us);
weff_hf=weff;
if Na>0;
   weff_hfmax=ld*quad('deplintg',sign(us(nv)),us(nv));
   mdepl=log(weff_hfmax/weffmax)/log(us(nv)/us(imin));
   adepl=weffmax/(us(imin)^mdepl);
   for i=(imin+1):nv;
   weff_hf(i)=adepl*us(i)^mdepl;
   end   
   Chf=Cox./(1+(epsox*weff_hf./(epssi*tox)));
else
   weff_hfmax=-ld*quad('deplintg',sign(us(1)),us(1));
   mdepl=log(weff_hfmax/weffmax)/log(us(1)/us(imin));
   adepl=weffmax/((-us(imin))^mdepl);
   for i=(imin-1):-1:1;
   weff_hf(i)=adepl*(-us(i))^mdepl;
   end   
   Chf=Cox./(1+(epsox*weff_hf./(epssi*tox)));
end;

%plot C-vs-Vg
   if ishandle(Hf)==0;
      return; 
   end;
set(Hp_cvlf,'XData',Vg,'YData',Clf*1000,'Visible','on');
set(Hp_cvhf,'XData',Vg,'YData',Chf*1000,'Visible','on');

%Find the quiescent point
Vgq=0; %To be adjusted by a slider
[Vgqcl,iminq]=min(abs(Vg-Vgq));
usq=us(iminq);
if lfhf>0;
   Cq=Clf(iminq);
else
   Cq=Chf(iminq);
end   
set(Hp_cvq,'XData',Vgq,'YData',Cq*1000,'Visible','on');

%2D Graph
if Na>0;
   Nsub=Na;
else
   Nsub=-Nd;
end
if Npoly==0; 
   material=(material==11).*2+(material~=11).*material; %Change gate to metal
   z2=(material>11).*sign(Nsub).*log10(abs(Nsub))+...
      (material==2).*2;
else
   z2=(material>11).*sign(Nsub).*log10(abs(Nsub))+...
      (material==11).*Npoly;
end   
   if ishandle(Hf)==0;
      return; 
   end;
Ha2=axes('Position',[0 (0.98-0.35*Ymax/Xmax) 0.35 0.35*Ymax/Xmax]);
Hs=surface(x,y,z2,'EdgeColor','none','Visible','off');
set(Ha2,'Xlim',[0 Xmax],'Ylim',[0 Ymax],'Zlim',[-z2_max z2_max],...
   'Color',[0.8 0.8 0.8],'FontAngle','italic','CLim',[1 64],...
   'View',[0 90],'Visible','off');
ylabel('Lateral Distance (nm)'); xlabel('Depth (nm)');
set(Hs,'CData',(material<=10).*material+...
   (material>10).*((abs(z2)>bright).*(z2-sign(z2).*bright)*26.5/(z2_max-bright)+37.5));
   
%1D Graph (dummy numbers)
[nysteps nxsteps]=size(x);
x1=x(fix(nxsteps/2),:); y1=z2(fix(nxsteps/2),:); %Initial data
z1=0*x1; 
Ha1=axes('Position',[0 (0.6-0.35*Ymax/Xmax) 0.35 0.20*Ymax/Xmax],...
   'Visible','off');
set(Ha1,'Xlim',[0 Xmax],'YLim',[-1.7 1.7],'ZLim',[z2_min z2_max],...
   'Color',[0.8 0.8 0.8],'CLim',[1 64]);
%Valence band
Hp_vb=surface([x1;x1],[abs(y1);abs(y1)],[z1;z1],'LineWidth',3,...
      'FaceColor','none','EdgeColor','flat','Visible','off');
   hold on;
%Conduction band   
Hp_cb=surface([x1;x1],[abs(y1);abs(y1)],[z1;z1],'LineWidth',3,...
   'FaceColor','none','EdgeColor','flat','Visible','off');
%Oxide
Hph_oxeb=patch([x1(2) x1(3) x1(3) x1(2)],...
   [-1.7 -1.7 1.7 1.7],'r','Visible','off'); %Oxide on EB diagram
%Metal
Hph_met=patch([x1(1) x1(2) x1(2) x1(1)],...
   [-1.7 -1.7 -uf_usd11*kt -uf_usd11*kt],'r','Visible','off'); %Oxide on EB diagram
set(Hph_met,'FaceColor',[0.4 0.4 0.4],'EdgeColor',[0.4 0.4 0.4],'Visible','on');
%Ef in substrate
[xsiz1 xsiz2]=size(x1);
Hp_efsub=plot([x1(3) x1(xsiz2)],[-uf_usd11*kt -uf_usd11*kt],'k-.','Visible','off');
%Ef in gate
Hp_efgate=plot([x1(1) x1(2)],[-uf_usd11*kt -uf_usd11*kt],'w-.','Visible','off'); %Dummy numbers

%The text
Ht_Ec=text(Xmax,1.12/2,'E_{C}','Color','r','FontSize',12,...
   'FontAngle','italic','FontWeight','bold','Visible','off');
Ht_Ev=text(Xmax,-1.12/2,'E_{V}','Color','r','FontSize',12,...
   'FontAngle','italic','FontWeight','bold','Visible','off');
Ht_Ef=text(Xmax,0,'E_{F}','Color','k','FontSize',12,...
   'FontAngle','italic','FontWeight','bold','Visible','off');
hold off;

%CALCULATE THE CONCENTRATION SURFACE FOR 2D GRAPH
%Calculate wdepl/(wacc=-wdepl);
%ld*1e6
%usq=-35; %%%THIS A TEST LINE ONLY
if abs(usq)<0.1;
   wdepl=0;
else   
if Na>0;
   wdepl=ld*quad('deplintg',sign(usq),usq);
else
   wdepl=-ld*quad('deplintg',sign(usq),usq);
end   
end

%Set the potential matrix assuming a parabolic change inside the depl. layer
if wdepl==0;
   u=rg*0;
else   
u=usq*((1-1e-9*rg./abs(wdepl)).^2).*(1e-9*rg<abs(wdepl));
end

%Calculate the concentration
Ncarr=ni*sign(uf_usd11-u).*exp(abs(uf_usd11-u));
Ncarr=(abs(Ncarr)<2e10).*1e10+(abs(Ncarr)>=2e10).*Ncarr;
%   if abs(Ncarr)<1e10; Ncarr=1e10; end;


%Set the 2D graph
if Npoly==0; 
   material=(material==11).*2+(material~=11).*material; %Change gate to metal
   z2=(material>11).*sign(Ncarr).*log10(abs(Ncarr))+...
      (material==2).*2;
else
   z2=(material>11).*sign(Ncarr).*log10(abs(Ncarr))+...
      (material==11).*Npoly;
end 
   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),...
   'Visible','on');

%Reset 1D graph
set(Hph_oxeb,'FaceColor',[0.8 0.8 1],'EdgeColor',[0.8 0.8 1],'Visible','on');
yen1=-u(fix(nxsteps/2),:)*kt+1.12/2; %Energy band numbers
yen1(1)=-uf_usd11*kt-Vgq;
yen1(2)=yen1(1);
y1=z2(fix(nxsteps/2),:); %Log concentration numbers for coloring
y1(1)=Npoly; y1(2)=Npoly;
y1c=(abs(y1)>bright).*(y1-sign(y1).*bright)*26.5/(z2_max-bright)+37.5;  
y1c(2)=6; y1c(3)=6; %same color as ox for gate<x<Si surface
set(Hp_cb,'YData',[yen1;yen1],'Visible','on','CData',[y1c;y1c]);   
set(Hp_vb,'YData',[(yen1-1.12);(yen1-1.12)],'Visible','on','CData',[y1c;y1c]);
set(Hph_met,'Visible','off');
set(Hp_efsub,'Visible','on');
set(Hp_efgate,'YData',[yen1(1) yen1(2)],'Visible','on');
set(Ht_Ec,'Position',[Xmax 1.12/2+0.35],'Visible','on');
set(Ht_Ev,'Position',[Xmax -1.12/2-0.25],'Visible','on');
set(Ht_Ef,'Position',[Xmax -uf_usd11*kt],'Visible','on');


%%%%%%%%%%%%%%%%%%%%%%
%%%THE MAIN LOOP
%%%%%%%%%%%%%%%%%%%%%%%
while con_close_usd11<1;
   
%BRIGHTNESS
%%%%%%%%%%%%%%%%%%%%%%%%
if con_br_usd11>0;
   bright=get(Hsl_br,'Value');
   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 Npoly==0;
      y1c=(abs(y1)>bright).*(y1-sign(y1).*bright)*26.5/(z2_max-bright)+37.5;  
      y1c(1)=3;
      y1c(2)=6; y1c(3)=6; %same color as ox for gate<x<Si surface
   else   
      y1c=(abs(y1)>bright).*(y1-sign(y1).*bright)*26.5/(z2_max-bright)+37.5;  
      y1c(2)=6; y1c(3)=6; %same color as ox for gate<x<Si surface
   end
   set(Hp_cb,'CData',[y1c;y1c]);   
   set(Hp_vb,'CData',[y1c;y1c]);
end; %if con_br_usd11>0;
con_br_usd11=0;
   
%LF - HF RADIOBUTTONS   
%%%%%%%%%%%%%%%%%%%%%%%%
if con_lfhf_usd11>0;
   if ishandle(Hf)==0;
      return; 
   end;
   if con_lfhf_usd11==1;
      lfhf=1;
      Cq=Clf(iminq);
      set(Hr_lfhf(1),'Value',1); 
      set(Hr_lfhf(2),'Value',0); 
      set(Ht_ac,'String','LF AC SIGNAL');
   else
      lfhf=0;
      Cq=Chf(iminq);
      set(Hr_lfhf(1),'Value',0); 
      set(Hr_lfhf(2),'Value',1); 
      set(Ht_ac,'String','HF AC SIGNAL');
   end
   %Set the quiescent point
   set(Hp_cvq,'XData',Vgq,'YData',Cq*1000,'Visible','on');
end; %if con_lfhf_usd11>0;   
con_lfhf_usd11=0;

%GATE VOLTAGE IS CHANGED   
%%%%%%%%%%%%%%%%%%%%%%%%   
if con_vg_usd11>0;   
   if ishandle(Hf)==0;
      return; 
   end;
   Vgq=round(get(Hsl_vg,'Value')*10)/10;
   set(Ht_vgd,'String',['Vg=' num2str(Vgq) ' V']);
   
   %Find the quiescent point
   [Vgqcl,iminq]=min(abs(Vg-Vgq));
   usq=us(iminq);
   if lfhf>0;
      Cq=Clf(iminq);
   else
      Cq=Chf(iminq);
   end   
   set(Hp_cvq,'XData',Vgq,'YData',Cq*1000);

   %Calculate the concentration surface for the 2D graph
   %Calculate wdepl/(wacc=-wdepl);
   if abs(usq)<0.1;
      wdepl=0;
   else   
      if Na>0;
         wdepl=ld*quad('deplintg',sign(usq),usq);
      else
         wdepl=-ld*quad('deplintg',sign(usq),usq);
      end   
   end

   %Set the potential matrix assuming a parabolic change inside the depl. layer
   if wdepl==0;
      u=rg*0;
   else   
      u=usq*((1-1e-9*rg./abs(wdepl)).^2).*(1e-9*rg<abs(wdepl));
   end

   %Calculate the concentration
   Ncarr=ni*sign(uf_usd11-u).*exp(abs(uf_usd11-u));
   Ncarr=(abs(Ncarr)<2e10).*1e10+(abs(Ncarr)>=2e10).*Ncarr;

   %Set the 2D graph
   if Npoly==0; 
      material=(material==11).*2+(material~=11).*material; %Change gate to metal
      z2=(material>11).*sign(Ncarr).*log10(abs(Ncarr))+...
         (material==2).*2;
   else
      z2=(material>11).*sign(Ncarr).*log10(abs(Ncarr))+...
         (material==11).*Npoly;
   end  
   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));
    
   %Reset 1D graph
   yen1=-u(fix(nxsteps/2),:)*kt+1.12/2; %Energy band numbers
   yen1(1)=-uf_usd11*kt-Vgq+(Npoly>0)*1.12;
   yen1(2)=yen1(1);
   y1=z2(fix(nxsteps/2),:); %Log concentration numbers for coloring
   if Npoly==0;
      yen1(1)=-10; yen2(2)=-10; %Remove it from the graph
      y1c=(abs(y1)>bright).*(y1-sign(y1).*bright)*26.5/(z2_max-bright)+37.5;  
      y1c(1)=3;
      y1c(2)=6; y1c(3)=6; %same color as ox for gate<x<Si surface
      set(Hph_met,'YData',[-1.7 -1.7 (-uf_usd11*kt-Vgq) (-uf_usd11*kt-Vgq)],'Visible','on');
   else   
      y1(1)=Npoly; y1(2)=Npoly;
      y1c=(abs(y1)>bright).*(y1-sign(y1).*bright)*26.5/(z2_max-bright)+37.5;  
      y1c(2)=6; y1c(3)=6; %same color as ox for gate<x<Si surface
      set(Hph_met,'Visible','off');
   end
   set(Hp_efsub,'YData',[-uf_usd11*kt -uf_usd11*kt]);
   set(Hp_cb,'YData',[yen1;yen1],'CData',[y1c;y1c]);   
   set(Hp_vb,'YData',[(yen1-1.12);(yen1-1.12)],'CData',[y1c;y1c]);
   set(Hp_efsub,'Visible','on');
   set(Hp_efgate,'YData',[(-uf_usd11*kt-Vgq) (-uf_usd11*kt-Vgq)],'Visible','on');

end; %if con_vg_usd11>0
con_vg_usd11=0;

% SUBSTRATE DOPING IS CHANGED
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if con_doping_usd11>0;
   if ishandle(Hf)==0;
      return; 
   end;
   switch con_doping_usd11;
   case 1; %P type
   set(Hr_dop(1),'Value',1);
   set(Hr_dop(2),'Value',0);
   set(Hframe_dop,'BackgroundColor','r');
   set(Ht_dop,'String',['Na=' num2str(round(10^(get(Hsl_dop,'Value')))) ' cm-3'],...
   'BackgroundColor','r','ForegroundColor','w');
   Na=round(10^(get(Hsl_dop,'Value')));
   Nd=0;
   case 2; %N-type
   set(Hr_dop(1),'Value',0);
   set(Hr_dop(2),'Value',1);
   set(Hframe_dop,'BackgroundColor','g');
   set(Ht_dop,'String',['Nd=' num2str(round(10^(get(Hsl_dop,'Value')))) ' cm-3'],...
   'BackgroundColor','g','ForegroundColor','k');
   Nd=round(10^(get(Hsl_dop,'Value')));
   Na=0;
   case 3;
      if Na>0; 
         Na=round(10^(get(Hsl_dop,'Value')));
         set(Ht_dop,'String',['Na=' num2str(round(10^(get(Hsl_dop,'Value')))) ' cm-3']);
      else      
         Nd=round(10^(get(Hsl_dop,'Value')));
         set(Ht_dop,'String',['Nd=' num2str(round(10^(get(Hsl_dop,'Value')))) ' cm-3']);
      end   
   end; %switch
   drawnow;
   
   if Na>0;
   Nsub=Na;
   else
   Nsub=-Nd;
   end

   %CALCULATE C-V CURVE  
   if Na>0; 
      uf_usd11=log(Na/ni);
   else
      uf_usd11=-log(Nd/ni);
   end   
   Vfb=phim-(4.05+1.12/2+uf_usd11*kt)-Qox/Cox;

   %Calculate the Clf-vs-us
   us=[uf_usd11-25:0.5:uf_usd11+25];
   fus=sqrt(exp(uf_usd11)*(exp(-us)+us-1)+exp(-uf_usd11)*(exp(us)-us-1));
   if Na>0;
      negl=-0.1; posl=0.1;
      usacc=(us<negl).*us+(us>=negl);
      usinv=(us>posl).*us+(us<=posl);
      weffacc=sign(us).*ld.*(2*fus./(exp(uf_usd11)*(1-exp(-usacc))+exp(-uf_usd11)*(exp(usacc)-1)));
      weff0=sqrt(2)*ld/sqrt(exp(uf_usd11)+exp(-uf_usd11));
      weffinv=sign(us).*ld.*(2*fus./(exp(uf_usd11)*(1-exp(-usinv))+exp(-uf_usd11)*(exp(usinv)-1)));
      weff=weffacc.*(us<negl)+weff0.*(us>=negl).*(us<=posl)+weffinv.*(us>posl);
   else
      negl=-0.1; posl=0.1;
      usacc=(us>posl).*us+(us<=posl);
      usinv=(us<negl).*us+(us>=negl);
      weffacc=sign(us).*ld.*(2*fus./(exp(uf_usd11)*(1-exp(-usacc))+exp(-uf_usd11)*(exp(usacc)-1)));
      weff0=sqrt(2)*ld/sqrt(exp(uf_usd11)+exp(-uf_usd11));
      weffinv=sign(us).*ld.*(2*fus./(exp(uf_usd11)*(1-exp(-usinv))+exp(-uf_usd11)*(exp(usinv)-1)));
      weff=weffacc.*(us>posl)+weff0.*(us>=negl).*(us<=posl)+weffinv.*(us<negl);
   end
   Clf=Cox./(1+(epsox*weff./(epssi*tox)));

   %Calculate Vg-vs-us
   Vg=kt*(us+sign(us).*(epssi*tox/(epsox*ld)).*fus)+Vfb;

   %Calculate Chf
   [Cmin,imin]=min(Clf);
   weffmax=weff(imin);
   Vgmin=Vg(imin);
   nv=length(us);
   weff_hf=weff;
   if Na>0;
      weff_hfmax=ld*quad('deplintg',sign(us(nv)),us(nv));
      mdepl=log(weff_hfmax/weffmax)/log(us(nv)/us(imin));
      adepl=weffmax/(us(imin)^mdepl);
      for i=(imin+1):nv;
         weff_hf(i)=adepl*us(i)^mdepl;
      end   
      Chf=Cox./(1+(epsox*weff_hf./(epssi*tox)));
   else
      weff_hfmax=-ld*quad('deplintg',sign(us(1)),us(1));
      mdepl=log(weff_hfmax/weffmax)/log(us(1)/us(imin));
      adepl=weffmax/((-us(imin))^mdepl);
      for i=(imin-1):-1:1;
         weff_hf(i)=adepl*(-us(i))^mdepl;
      end   
      Chf=Cox./(1+(epsox*weff_hf./(epssi*tox)));
   end;
   
   %plot C-vs-Vg
   set(Hp_cvlf,'XData',Vg,'YData',Clf*1000,'Visible','on');
   set(Hp_cvhf,'XData',Vg,'YData',Chf*1000,'Visible','on');

   %Find the quiescent point
   [Vgqcl,iminq]=min(abs(Vg-Vgq));
   usq=us(iminq);
   if lfhf>0;
      Cq=Clf(iminq);
   else
      Cq=Chf(iminq);
   end   
   set(Hp_cvq,'XData',Vgq,'YData',Cq*1000,'Visible','on');


   %CALCULATE THE CONCENTRATION SURFACE FOR 2D GRAPH
   if abs(usq)<0.1;
      wdepl=0;
   else   
      if Na>0;
         wdepl=ld*quad('deplintg',sign(usq),usq);
      else
         wdepl=-ld*quad('deplintg',sign(usq),usq);
      end   
   end

   %Set the potential matrix assuming a parabolic change inside the depl. layer
   if wdepl==0;
      u=rg*0;
   else   
      u=usq*((1-1e-9*rg./abs(wdepl)).^2).*(1e-9*rg<abs(wdepl));
   end

   %Calculate the concentration
   Ncarr=ni*sign(uf_usd11-u).*exp(abs(uf_usd11-u));
   Ncarr=(abs(Ncarr)<2e10).*1e10+(abs(Ncarr)>=2e10).*Ncarr;

   %Set the 2D graph
   if Npoly==0; 
      material=(material==11).*2+(material~=11).*material; %Change gate to metal
      z2=(material>11).*sign(Ncarr).*log10(abs(Ncarr))+...
         (material==2).*2;
   else
      material=(material==2).*11+(material~=2).*material; %Chgange metal to poly
      z2=(material>11).*sign(Ncarr).*log10(abs(Ncarr))+...
         (material==11).*Npoly;
   end   
   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),...
      'Visible','on');
   
   %Reset 1D graph
   yen1=-u(fix(nxsteps/2),:)*kt+1.12/2; %Energy band numbers
   yen1(1)=-uf_usd11*kt-Vgq+(Npoly>0)*1.12;
   yen1(2)=yen1(1);
   y1=z2(fix(nxsteps/2),:); %Log concentration numbers for coloring
   if Npoly==0;
      yen1(1)=-10; yen2(2)=-10; %Remove it from the graph
      y1c=(abs(y1)>bright).*(y1-sign(y1).*bright)*26.5/(z2_max-bright)+37.5;  
      y1c(1)=3;
      y1c(2)=6; y1c(3)=6; %same color as ox for gate<x<Si surface
      set(Hph_met,'YData',[-1.7 -1.7 (-uf_usd11*kt-Vgq) (-uf_usd11*kt-Vgq)],'Visible','on');
   else   
      y1(1)=Npoly; y1(2)=Npoly;
      y1c=(abs(y1)>bright).*(y1-sign(y1).*bright)*26.5/(z2_max-bright)+37.5;  
      y1c(2)=6; y1c(3)=6; %same color as ox for gate<x<Si surface
      set(Hph_met,'Visible','off');
   end
   set(Hp_efsub,'YData',[-uf_usd11*kt -uf_usd11*kt],'Visible','on');
   set(Hp_cb,'YData',[yen1;yen1],'CData',[y1c;y1c],'Visible','on');   
   set(Hp_vb,'YData',[(yen1-1.12);(yen1-1.12)],'CData',[y1c;y1c],'Visible','on');
   set(Hp_efsub,'Visible','on');
   set(Hp_efgate,'YData',[(-uf_usd11*kt-Vgq) (-uf_usd11*kt-Vgq)],'Visible','on');
   if Na>0;
      set([Ht_Ec Ht_Ev],'Color','r');
   else
      set([Ht_Ec Ht_Ev],'Color','g');
   end
   set(Ht_Ef,'Position',[Xmax -uf_usd11*kt]);

end; %if con_doping_usd11>0;
con_doping_usd11=0;

%OXIDE THICKNESS IS CHANGED
%%%%%%%%%%%%%%%%%%%%%%%%%%%
if con_ox_usd11>0;
   if ishandle(Hf)==0;
      return; 
   end;
   set(Ht_ox,'String',['Tox=' num2str(round(get(Hsl_ox,'Value'))) ' nm']);
   tox=round(get(Hsl_ox,'Value'))*1e-9;
   Cox=epsox/tox;
   Vfb=phim-(4.05+1.12/2+uf_usd11*kt)-Qox/Cox;
   
   % SET THE NEW GRID
   clear xs; clear ys; clear xa; clear z2; clear material;
   Xmax=250;
   Ymax=400; ygu=300; ygl=Ymax-ygu;  

   num_x=round((Xmax-tox*1e9-20)*scal); %Number of lateral points - 3...
   %(three lines outside silicon will be added)
   num_yhalf=round((Ymax-ygu)*scal); %Number of vertical points above the gate
   [xs,ys]=meshgrid(([0:num_x].^2/num_x),num_yhalf:-1:1);
   %Set the grid inside the silicon
   xs=xs/scal; ys=ys/scal+ygu; %The grid above ygu
   ysl=ys-ygu;
   ys=[ys' [ys(1,:)*0+Ymax/2]' ysl']';
   xs=[xs' [xs(1,:)]' xs']';
   %Set the grid points outside the silicon
   gatethick=40;
   xa(:,1)=xs(:,1)*0;
   xa(:,2)=xa(:,1)+gatethick; %The metal/oxide interface
   xs=xs+tox*1e9+xa(1,2); %%Shift the silicon grid inside the silicon...
   % to the oxide/Si interface
   x=[xa xs]; %Merge the ambient and silicon x-grid points
   y=[ys(:,1) ys(:,1) ys]; %Expand the y-grid points to include the ambient

   %Establish the ends of the gate at grid lines
   ygu=ys((num_yhalf),1);
   ygl=ys((num_yhalf+2),1);

   %Set material matrix 
   material=[(xa(:,1)*0+1) (xa(:,2)*0+6 ) (xs*0+20)]; %air=1, oxide=6, substrate=20
   material=(material==1).*(y>=ygl).*(y<=ygu).*11+... %Poly-gate=11
   (material~=1).*material;
   
   %Set "distance from the gate interface" matrix;
   rg=((x-x(1,3)).*(y>=ygl).*(y<=ygu)+...
      (y<ygl).*((x-x(1,3)).^2+(ygl-y).^2).^0.5+...
      (y>ygu).*((x-x(1,3)).^2+(y-ygu).^2).^0.5).*(material>11);

      %CALCULATE C-V CURVE  
   if Na>0; 
      uf_usd11=log(Na/ni);
   else
      uf_usd11=-log(Nd/ni);
   end   
   Vfb=phim-(4.05+1.12/2+uf_usd11*kt)-Qox/Cox;

   %Calculate the Clf-vs-us
   us=[uf_usd11-25:0.5:uf_usd11+25];
   fus=sqrt(exp(uf_usd11)*(exp(-us)+us-1)+exp(-uf_usd11)*(exp(us)-us-1));
   if Na>0;
      negl=-0.1; posl=0.1;
      usacc=(us<negl).*us+(us>=negl);
      usinv=(us>posl).*us+(us<=posl);
      weffacc=sign(us).*ld.*(2*fus./(exp(uf_usd11)*(1-exp(-usacc))+exp(-uf_usd11)*(exp(usacc)-1)));
      weff0=sqrt(2)*ld/sqrt(exp(uf_usd11)+exp(-uf_usd11));
      weffinv=sign(us).*ld.*(2*fus./(exp(uf_usd11)*(1-exp(-usinv))+exp(-uf_usd11)*(exp(usinv)-1)));
      weff=weffacc.*(us<negl)+weff0.*(us>=negl).*(us<=posl)+weffinv.*(us>posl);
   else
      negl=-0.1; posl=0.1;
      usacc=(us>posl).*us+(us<=posl);
      usinv=(us<negl).*us+(us>=negl);
      weffacc=sign(us).*ld.*(2*fus./(exp(uf_usd11)*(1-exp(-usacc))+exp(-uf_usd11)*(exp(usacc)-1)));
      weff0=sqrt(2)*ld/sqrt(exp(uf_usd11)+exp(-uf_usd11));
      weffinv=sign(us).*ld.*(2*fus./(exp(uf_usd11)*(1-exp(-usinv))+exp(-uf_usd11)*(exp(usinv)-1)));
      weff=weffacc.*(us>posl)+weff0.*(us>=negl).*(us<=posl)+weffinv.*(us<negl);
   end
   Clf=Cox./(1+(epsox*weff./(epssi*tox)));

   %Calculate Vg-vs-us
   Vg=kt*(us+sign(us).*(epssi*tox/(epsox*ld)).*fus)+Vfb;

   %Calculate Chf
   [Cmin,imin]=min(Clf);
   weffmax=weff(imin);
   Vgmin=Vg(imin);
   nv=length(us);
   weff_hf=weff;
   if Na>0;
      weff_hfmax=ld*quad('deplintg',sign(us(nv)),us(nv));
      mdepl=log(weff_hfmax/weffmax)/log(us(nv)/us(imin));
      adepl=weffmax/(us(imin)^mdepl);
      for i=(imin+1):nv;
         weff_hf(i)=adepl*us(i)^mdepl;
      end   
      Chf=Cox./(1+(epsox*weff_hf./(epssi*tox)));
   else
      weff_hfmax=-ld*quad('deplintg',sign(us(1)),us(1));
      mdepl=log(weff_hfmax/weffmax)/log(us(1)/us(imin));
      adepl=weffmax/((-us(imin))^mdepl);
      for i=(imin-1):-1:1;
         weff_hf(i)=adepl*(-us(i))^mdepl;
      end   
      Chf=Cox./(1+(epsox*weff_hf./(epssi*tox)));
   end;
   if ishandle(Hf)==0;
      return; 
   end;
   %plot C-vs-Vg
   if tox<10e-9;
   Clfmax=max(Clf*1000);
   Ndig=10^floor(log10(Clfmax));
   Clfmaxgr=Ndig*2*ceil((Clfmax/Ndig)/2);
   set(Ha_cv,'Ylim',[0 Clfmaxgr]);
   else
   set(Ha_cv,'Ylim',[0 4]);
   end;
   set(Hp_cvlf,'XData',Vg,'YData',Clf*1000);
   set(Hp_cvhf,'XData',Vg,'YData',Chf*1000);

   %Find the quiescent point
   [Vgqcl,iminq]=min(abs(Vg-Vgq));
   usq=us(iminq);
   if lfhf>0;
      Cq=Clf(iminq);
   else
      Cq=Chf(iminq);
   end   
   set(Hp_cvq,'XData',Vgq,'YData',Cq*1000,'Visible','on');


   %CALCULATE THE CONCENTRATION SURFACE FOR 2D GRAPH
   if abs(usq)<0.1;
      wdepl=0;
   else   
      if Na>0;
         wdepl=ld*quad('deplintg',sign(usq),usq);
      else
         wdepl=-ld*quad('deplintg',sign(usq),usq);
      end   
   end

   %Set the potential matrix assuming a parabolic change inside the depl. layer
   if wdepl==0;
      u=rg*0;
   else   
      u=usq*((1-1e-9*rg./abs(wdepl)).^2).*(1e-9*rg<abs(wdepl));
   end

   %Calculate the concentration
   Ncarr=ni*sign(uf_usd11-u).*exp(abs(uf_usd11-u));
   Ncarr=(abs(Ncarr)<2e10).*1e10+(abs(Ncarr)>=2e10).*Ncarr;

   %Set the 2D graph
   if Npoly==0; 
      material=(material==11).*2+(material~=11).*material; %Change gate to metal
      z2=(material>11).*sign(Ncarr).*log10(abs(Ncarr))+...
         (material==2).*2;
   else
      z2=(material>11).*sign(Ncarr).*log10(abs(Ncarr))+...
         (material==11).*Npoly;
   end   
   if ishandle(Hf)==0;
      return; 
   end;
   set(Hs,'Xdata',x,'YData',y,'ZData',z2,'CData',(material<=10).*material+...
      (material>10).*((abs(z2)>bright).*(z2-sign(z2).*bright)*26.5/(z2_max-bright)+37.5));
   
   %Reset 1D graph
   [nysteps nxsteps]=size(x);
   x1=x(fix(nxsteps/2),:); y1=z2(fix(nxsteps/2),:); %Initial data
   z1=0*x1; 
   set(Hph_oxeb,'XData',[x1(2) x1(3) x1(3) x1(2)]);
   yen1=-u(fix(nxsteps/2),:)*kt+1.12/2; %Energy band numbers
   yen1(1)=-uf_usd11*kt-Vgq+(Npoly>0)*1.12;
   yen1(2)=yen1(1);
   y1=z2(fix(nxsteps/2),:); %Log concentration numbers for coloring
   if Npoly==0;
      x1(1)=x1(3); x1(2)=x1(3); %To remove gate e.b.
      set([Hp_cb Hp_vb],'XData',[x1;x1]);
      yen1(1)=-10; yen2(2)=-10; %Remove it from the graph
      y1c=(abs(y1)>bright).*(y1-sign(y1).*bright)*26.5/(z2_max-bright)+37.5;  
      y1c(1)=3;
      y1c(2)=6; y1c(3)=6; %same color as ox for gate<x<Si surface
      set(Hph_met,'YData',[-1.7 -1.7 yen1(1) yen1(2)]);
      set(Hph_met,'YData',[-1.7 -1.7 (-uf_usd11*kt-Vgq) (-uf_usd11*kt-Vgq)],'Visible','on');
   else   
      x1=x(fix(nxsteps/2),:);
      set([Hp_cb Hp_vb],'XData',[x1;x1]);
      y1(1)=Npoly; y1(2)=Npoly;
      y1c=(abs(y1)>bright).*(y1-sign(y1).*bright)*26.5/(z2_max-bright)+37.5;  
      y1c(2)=6; y1c(3)=6; %same color as ox for gate<x<Si surface
      set(Hph_met,'Visible','off');
   end
   set(Hp_cb,'XData',[x1;x1],'YData',[yen1;yen1],'ZData',[z1;z1],'CData',[y1c;y1c]);   
   set(Hp_vb,'XData',[x1;x1],'YData',[(yen1-1.12);(yen1-1.12)],'ZData',[z1;z1],'CData',[y1c;y1c]);
   [xsiz1 xsiz2]=size(x1);
   set(Hp_efsub,'XData',[x1(3) x1(xsiz2)],'YData',[-uf_usd11*kt -uf_usd11*kt],'Visible','on');
   set(Hp_efgate,'YData',[(-uf_usd11*kt-Vgq) (-uf_usd11*kt-Vgq)],'Visible','on');
  
end; %if con_ox_usd11>0;
con_ox_usd11=0;
pause(0.1);

%GATE MATERIAL CHANGED
%%%%%%%%%%%%%%%%%%%%%%
if con_gate_usd11>0;
   if ishandle(Hf)==0;
      return; 
   end;
   switch con_gate_usd11;
   case 1;
      set([Hr_gate(2) Hr_gate(3)],'Value',0);
      set(Hr_gate(1),'Value',1);
      set(Hsl_gate,'Enable','off');
      set(Ht_gate,'String',['Work function=4.05 eV']);
      phim=4.05; Npoly=-20;
   case 2;
      set([Hr_gate(1) Hr_gate(3)],'Value',0);
      set(Hr_gate(2),'Value',1);
      set(Hsl_gate,'Enable','off');
      set(Ht_gate,'String',['Work function=5.17 eV']);
      phim=5.17; Npoly=19.5;
   case 3;
      set([Hr_gate(1) Hr_gate(2)],'Value',0);
      set(Hr_gate(3),'Value',1);
      set(Hsl_gate,'Enable','on','Value',4.10);
      set(Ht_gate,'String',['Work function=4.10 eV']);
      phim=4.1; Npoly=0;
   case 4;   
      set([Hr_gate(1) Hr_gate(2)],'Value',0);
      set(Hr_gate(3),'Value',1);
      set(Hsl_gate,'Enable','on');
      set(Ht_gate,'String',['Work function=' num2str(round(100*get(Hsl_gate,'Value'))/100) ' eV']);
      phim=round(100*get(Hsl_gate,'Value'))/100; Npoly=0;
   end; %end switch
   
      Vfb=phim-(4.05+1.12/2+uf_usd11*kt)-Qox/Cox;

   %Calculate the Clf-vs-us
   us=[uf_usd11-25:0.5:uf_usd11+25];
   fus=sqrt(exp(uf_usd11)*(exp(-us)+us-1)+exp(-uf_usd11)*(exp(us)-us-1));
   if Na>0;
      negl=-0.1; posl=0.1;
      usacc=(us<negl).*us+(us>=negl);
      usinv=(us>posl).*us+(us<=posl);
      weffacc=sign(us).*ld.*(2*fus./(exp(uf_usd11)*(1-exp(-usacc))+exp(-uf_usd11)*(exp(usacc)-1)));
      weff0=sqrt(2)*ld/sqrt(exp(uf_usd11)+exp(-uf_usd11));
      weffinv=sign(us).*ld.*(2*fus./(exp(uf_usd11)*(1-exp(-usinv))+exp(-uf_usd11)*(exp(usinv)-1)));
      weff=weffacc.*(us<negl)+weff0.*(us>=negl).*(us<=posl)+weffinv.*(us>posl);
   else
      negl=-0.1; posl=0.1;
      usacc=(us>posl).*us+(us<=posl);
      usinv=(us<negl).*us+(us>=negl);
      weffacc=sign(us).*ld.*(2*fus./(exp(uf_usd11)*(1-exp(-usacc))+exp(-uf_usd11)*(exp(usacc)-1)));
      weff0=sqrt(2)*ld/sqrt(exp(uf_usd11)+exp(-uf_usd11));
      weffinv=sign(us).*ld.*(2*fus./(exp(uf_usd11)*(1-exp(-usinv))+exp(-uf_usd11)*(exp(usinv)-1)));
      weff=weffacc.*(us>posl)+weff0.*(us>=negl).*(us<=posl)+weffinv.*(us<negl);
   end
   Clf=Cox./(1+(epsox*weff./(epssi*tox)));

   %Calculate Vg-vs-us
   Vg=kt*(us+sign(us).*(epssi*tox/(epsox*ld)).*fus)+Vfb;

   %Calculate Chf
   [Cmin,imin]=min(Clf);
   weffmax=weff(imin);
   Vgmin=Vg(imin);
   nv=length(us);
   weff_hf=weff;
   if Na>0;
      weff_hfmax=ld*quad('deplintg',sign(us(nv)),us(nv));
      mdepl=log(weff_hfmax/weffmax)/log(us(nv)/us(imin));
      adepl=weffmax/(us(imin)^mdepl);
      for i=(imin+1):nv;
         weff_hf(i)=adepl*us(i)^mdepl;
      end   
      Chf=Cox./(1+(epsox*weff_hf./(epssi*tox)));
   else
      weff_hfmax=-ld*quad('deplintg',sign(us(1)),us(1));
      mdepl=log(weff_hfmax/weffmax)/log(us(1)/us(imin));
      adepl=weffmax/((-us(imin))^mdepl);
      for i=(imin-1):-1:1;
         weff_hf(i)=adepl*(-us(i))^mdepl;
      end   
      Chf=Cox./(1+(epsox*weff_hf./(epssi*tox)));
   end;
   
   %plot C-vs-Vg
   set(Hp_cvlf,'XData',Vg,'YData',Clf*1000);
   set(Hp_cvhf,'XData',Vg,'YData',Chf*1000);

   %Find the quiescent point
   [Vgqcl,iminq]=min(abs(Vg-Vgq));
   usq=us(iminq);
   if lfhf>0;
      Cq=Clf(iminq);
   else
      Cq=Chf(iminq);
   end   
   set(Hp_cvq,'XData',Vgq,'YData',Cq*1000,'Visible','on');


   %CALCULATE THE CONCENTRATION SURFACE FOR 2D GRAPH
   if abs(usq)<0.1;
      wdepl=0;
   else   
      if Na>0;
         wdepl=ld*quad('deplintg',sign(usq),usq);
      else
         wdepl=-ld*quad('deplintg',sign(usq),usq);
      end   
   end

   %Set the potential matrix assuming a parabolic change inside the depl. layer
   if wdepl==0;
      u=rg*0;
   else   
      u=usq*((1-1e-9*rg./abs(wdepl)).^2).*(1e-9*rg<abs(wdepl));
   end

   %Calculate the concentration
   Ncarr=ni*sign(uf_usd11-u).*exp(abs(uf_usd11-u));
   Ncarr=(abs(Ncarr)<2e10).*1e10+(abs(Ncarr)>=2e10).*Ncarr;

   %Set the 2D graph
   if Npoly==0; 
      material=(material==11).*2+(material~=11).*material; %Change gate to metal
      z2=(material>11).*sign(Ncarr).*log10(abs(Ncarr))+...
         (material==2).*2;
   else
      material=(material==2).*11+(material~=2).*material; %Chgange metal to poly
      z2=(material>11).*sign(Ncarr).*log10(abs(Ncarr))+...
         (material==11).*Npoly;
   end  
   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));
   
      %Reset 1D graph
   yen1=-u(fix(nxsteps/2),:)*kt+1.12/2; %Energy band numbers
   yen1(1)=-uf_usd11*kt-Vgq+(Npoly>0)*1.12;
   yen1(2)=yen1(1);
   y1=z2(fix(nxsteps/2),:); %Log concentration numbers for coloring
   if Npoly==0;
      x1(1)=x1(3); x1(2)=x1(3); %To remove gate E.B.
      set([Hp_cb Hp_vb],'XData',[x1;x1]);
      y1c=(abs(y1)>bright).*(y1-sign(y1).*bright)*26.5/(z2_max-bright)+37.5;  
      y1c(1)=3;
      y1c(2)=6; y1c(3)=6; %same color as ox for gate<x<Si surface
      set(Hph_met,'YData',[-1.7 -1.7 (-uf_usd11*kt-Vgq) (-uf_usd11*kt-Vgq)],'Visible','on');
   else   
      x1=x(fix(nxsteps/2),:);
      set([Hp_cb Hp_vb],'XData',[x1;x1]);
      y1(1)=Npoly; y1(2)=Npoly;
      y1c=(abs(y1)>bright).*(y1-sign(y1).*bright)*26.5/(z2_max-bright)+37.5;  
      y1c(2)=6; y1c(3)=6; %same color as ox for gate<x<Si surface
      set(Hph_met,'Visible','off');
   end
   set(Hp_efsub,'YData',[-uf_usd11*kt -uf_usd11*kt]);
   set(Hp_cb,'YData',[yen1;yen1],'CData',[y1c;y1c]);   
   set(Hp_vb,'YData',[(yen1-1.12);(yen1-1.12)],'CData',[y1c;y1c]);
   set(Hp_efsub,'Visible','on');
   set(Hp_efgate,'YData',[(-uf_usd11*kt-Vgq) (-uf_usd11*kt-Vgq)],'Visible','on');
end; %con_gate_usd11>0;
con_gate_usd11=0;

%OXIDE CHARGE IS CHANGED
%%%%%%%%%%%%%%%%%%%%%%%%%
if con_chrg_usd11>0;
   if ishandle(Hf)==0;
      return; 
   end;
   if get(Hsl_chrg,'Value')==0;
      set(Ht_chrgd,'String',['Nox=0']);
      Qox=0;
   else   
      chrg=1e14*round(10*get(Hsl_chrg,'Value'))/10;
      if chrg>1e15;
         set(Ht_chrgd,'String',['Nox='num2str(round(get(Hsl_chrg,'Value'))/10) ' e+11 cm-2']);
      else
         set(Ht_chrgd,'String',['Nox='num2str(round(10*get(Hsl_chrg,'Value'))/10) ' e+10 cm-2']);
      end;   
      Qox=q*chrg;
   end
   
         Vfb=phim-(4.05+1.12/2+uf_usd11*kt)-Qox/Cox;

   %Calculate the Clf-vs-us
   us=[uf_usd11-25:0.5:uf_usd11+25];
   fus=sqrt(exp(uf_usd11)*(exp(-us)+us-1)+exp(-uf_usd11)*(exp(us)-us-1));
   if Na>0;
      negl=-0.1; posl=0.1;
      usacc=(us<negl).*us+(us>=negl);
      usinv=(us>posl).*us+(us<=posl);
      weffacc=sign(us).*ld.*(2*fus./(exp(uf_usd11)*(1-exp(-usacc))+exp(-uf_usd11)*(exp(usacc)-1)));
      weff0=sqrt(2)*ld/sqrt(exp(uf_usd11)+exp(-uf_usd11));
      weffinv=sign(us).*ld.*(2*fus./(exp(uf_usd11)*(1-exp(-usinv))+exp(-uf_usd11)*(exp(usinv)-1)));
      weff=weffacc.*(us<negl)+weff0.*(us>=negl).*(us<=posl)+weffinv.*(us>posl);
   else
      negl=-0.1; posl=0.1;
      usacc=(us>posl).*us+(us<=posl);
      usinv=(us<negl).*us+(us>=negl);
      weffacc=sign(us).*ld.*(2*fus./(exp(uf_usd11)*(1-exp(-usacc))+exp(-uf_usd11)*(exp(usacc)-1)));
      weff0=sqrt(2)*ld/sqrt(exp(uf_usd11)+exp(-uf_usd11));
      weffinv=sign(us).*ld.*(2*fus./(exp(uf_usd11)*(1-exp(-usinv))+exp(-uf_usd11)*(exp(usinv)-1)));
      weff=weffacc.*(us>posl)+weff0.*(us>=negl).*(us<=posl)+weffinv.*(us<negl);
   end
   Clf=Cox./(1+(epsox*weff./(epssi*tox)));

   %Calculate Vg-vs-us
   Vg=kt*(us+sign(us).*(epssi*tox/(epsox*ld)).*fus)+Vfb;

   %Calculate Chf
   [Cmin,imin]=min(Clf);
   weffmax=weff(imin);
   Vgmin=Vg(imin);
   nv=length(us);
   weff_hf=weff;
   if Na>0;
      weff_hfmax=ld*quad('deplintg',sign(us(nv)),us(nv));
      mdepl=log(weff_hfmax/weffmax)/log(us(nv)/us(imin));
      adepl=weffmax/(us(imin)^mdepl);
      for i=(imin+1):nv;
         weff_hf(i)=adepl*us(i)^mdepl;
      end   
      Chf=Cox./(1+(epsox*weff_hf./(epssi*tox)));
   else
      weff_hfmax=-ld*quad('deplintg',sign(us(1)),us(1));
      mdepl=log(weff_hfmax/weffmax)/log(us(1)/us(imin));
      adepl=weffmax/((-us(imin))^mdepl);
      for i=(imin-1):-1:1;
         weff_hf(i)=adepl*(-us(i))^mdepl;
      end   
      Chf=Cox./(1+(epsox*weff_hf./(epssi*tox)));
   end;
   
   %plot C-vs-Vg
   set(Hp_cvlf,'XData',Vg,'YData',Clf*1000);
   set(Hp_cvhf,'XData',Vg,'YData',Chf*1000);

   %Find the quiescent point
   [Vgqcl,iminq]=min(abs(Vg-Vgq));
   usq=us(iminq);
   if lfhf>0;
      Cq=Clf(iminq);
   else
      Cq=Chf(iminq);
   end   
   set(Hp_cvq,'XData',Vgq,'YData',Cq*1000,'Visible','on');


   %CALCULATE THE CONCENTRATION SURFACE FOR 2D GRAPH
   if abs(usq)<0.1;
      wdepl=0;
   else   
      if Na>0;
         wdepl=ld*quad('deplintg',sign(usq),usq);
      else
         wdepl=-ld*quad('deplintg',sign(usq),usq);
      end   
   end

   %Set the potential matrix assuming a parabolic change inside the depl. layer
   if wdepl==0;
      u=rg*0;
   else   
      u=usq*((1-1e-9*rg./abs(wdepl)).^2).*(1e-9*rg<abs(wdepl));
   end

   %Calculate the concentration
   Ncarr=ni*sign(uf_usd11-u).*exp(abs(uf_usd11-u));
   Ncarr=(abs(Ncarr)<2e10).*1e10+(abs(Ncarr)>=2e10).*Ncarr;


   %Set the 2D graph
   if Npoly==0; 
      material=(material==11).*2+(material~=11).*material; %Change gate to metal
      z2=(material>11).*sign(Ncarr).*log10(abs(Ncarr))+...
         (material==2).*2;
   else
      material=(material==2).*11+(material~=2).*material; %Chgange metal to poly
      z2=(material>11).*sign(Ncarr).*log10(abs(Ncarr))+...
         (material==11).*Npoly;
   end  
   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));
    %Reset 1D graph
   yen1=-u(fix(nxsteps/2),:)*kt+1.12/2; %Energy band numbers
   yen1(1)=-uf_usd11*kt-Vgq+(Npoly>0)*1.12;
   yen1(2)=yen1(1);
   y1=z2(fix(nxsteps/2),:); %Log concentration numbers for coloring
   if Npoly==0;
      yen1(1)=-10; yen2(2)=-10; %Remove it from the graph
      y1c=(abs(y1)>bright).*(y1-sign(y1).*bright)*26.5/(z2_max-bright)+37.5;  
      y1c(1)=3;
      y1c(2)=6; y1c(3)=6; %same color as ox for gate<x<Si surface
      set(Hph_met,'YData',[-1.7 -1.7 (-uf_usd11*kt-Vgq) (-uf_usd11*kt-Vgq)],'Visible','on');
   else   
      y1(1)=Npoly; y1(2)=Npoly;
      y1c=(abs(y1)>bright).*(y1-sign(y1).*bright)*26.5/(z2_max-bright)+37.5;  
      y1c(2)=6; y1c(3)=6; %same color as ox for gate<x<Si surface
      set(Hph_met,'Visible','off');
   end
   set(Hp_efsub,'YData',[-uf_usd11*kt -uf_usd11*kt]);
   set(Hp_cb,'YData',[yen1;yen1],'CData',[y1c;y1c]);   
   set(Hp_vb,'YData',[(yen1-1.12);(yen1-1.12)],'CData',[y1c;y1c]);
   set(Hp_efsub,'Visible','on');
   set(Hp_efgate,'YData',[(-uf_usd11*kt-Vgq) (-uf_usd11*kt-Vgq)],'Visible','on');
end; %if con_chrg_usd11>0      
con_chrg_usd11=0;

%%THE ANIMATION LOOP
%%%%%%%%%%%%%%%%%%%%

if con_strt_usd11>0;
   con_strt_usd11=0;      
   %Determine the time step;
   wt=0; delta_wt=2*pi/Nsteps;
   Vgq0=Vgq;
   usq0=usq;
   while con_stp_usd11<1;
      %Disable all the other buttons
      if ishandle(Hf)==0;
         return; 
      end;
      set([Hr_dop Hr_gate Hr_lfhf],'Enable','off');
      set([Hsl_dop Hsl_vg Hsl_ox Hsl_gate Hsl_chrg Hsl_br],'Enable','off');
      set([Hb_close Hb_start],'Enable','off');
      set(Hb_stop,'Enable','on');
      wt=wt+delta_wt;
      if wt>2*pi; wt=wt-2*pi; end
      if ((lfhf==1)+(sign(usq-uf_usd11)*sign(Nsub)<0)); %LOW FREQUENCY CASE
         %if (lfhf==1);
         dVac=Vac*sin(wt);
         %Vgq_old=Vgq;
         Vgq=Vgq0+dVac;
         %Find the quiescent point
         [Vgqcl,iminq]=min(abs(Vg-Vgq));
         %usq_old=usq;
         usq=us(iminq);
           Cq=Clf(iminq);
         set(Hp_cvq,'XData',Vgq,'YData',Cq*1000);

         %Calculate the concentration surface for the 2D graph
         %Calculate wdepl/(wacc=-wdepl);
         if abs(usq)<0.1;
            wdepl=0;
         else   
            if Na>0;
               wdepl=ld*quad('deplintg',sign(usq),usq);
            else
               wdepl=-ld*quad('deplintg',sign(usq),usq);
            end   
         end
         
         %Set the potential matrix assuming a parabolic change inside the depl. layer
         if wdepl==0;
            u=rg*0;
         else   
            u=usq*((1-1e-9*rg./abs(wdepl)).^2).*(1e-9*rg<abs(wdepl));
         end
         
         %Calculate the concentration
         Ncarr=ni*sign(uf_usd11-u).*exp(abs(uf_usd11-u));
            Ncarr=(abs(Ncarr)<2e10).*1e10+(abs(Ncarr)>=2e10).*Ncarr;

         linv=0.3*abs(wdepl)*1e9;
         if linv==0; linv=1; end;
         plusmin=sign(usq);
         if plusmin==0; plusmin=1; end;
         %Ncarr_coeff=2*plusmin*abs(Vgq-Vgq_old-usq*kt+usq_old*kt);
         Ncarr_coeff=4*plusmin*(dVac-usq*kt+usq0*kt);
         Ncarr_c=Ncarr.*(1+Ncarr_coeff*exp(-(rg/linv)));
         
         %Set the 2D graph
         if Npoly==0; 
            material=(material==11).*2+(material~=11).*material; %Change gate to metal
            z2=(material>11).*sign(Ncarr_c).*log10(abs(Ncarr_c))+...
               (material==2).*2;
         else
            z2=(material>11).*sign(Ncarr_c).*log10(abs(Ncarr_c))+...
               (material==11).*Npoly;
         end  
         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));
  
        %Reset 1D graph
         yen1=-u(fix(nxsteps/2),:)*kt+1.12/2; %Energy band numbers
         yen1(1)=-uf_usd11*kt-Vgq+(Npoly>0)*1.12;
         yen1(2)=yen1(1);
         y1=z2(fix(nxsteps/2),:); %Log concentration numbers for coloring
         if Npoly==0;
            yen1(1)=-10; yen2(2)=-10; %Remove it from the graph
            y1c=(abs(y1)>bright).*(y1-sign(y1).*bright)*26.5/(z2_max-bright)+37.5;  
            y1c(1)=3;
            y1c(2)=6; y1c(3)=6; %same color as ox for gate<x<Si surface
            set(Hph_met,'YData',[-1.7 -1.7 (-uf_usd11*kt-Vgq) (-uf_usd11*kt-Vgq)],'Visible','on');
         else   
            y1(1)=Npoly; y1(2)=Npoly;
            y1c=(abs(y1)>bright).*(y1-sign(y1).*bright)*26.5/(z2_max-bright)+37.5;  
            y1c(2)=6; y1c(3)=6; %same color as ox for gate<x<Si surface
            set(Hph_met,'Visible','off');
         end
         set(Hp_efsub,'YData',[-uf_usd11*kt -uf_usd11*kt]);
         set(Hp_cb,'YData',[yen1;yen1],'CData',[y1c;y1c]);   
         set(Hp_vb,'YData',[(yen1-1.12);(yen1-1.12)],'CData',[y1c;y1c]);
         set(Hp_efsub,'Visible','on');
         set(Hp_efgate,'YData',[(-uf_usd11*kt-Vgq) (-uf_usd11*kt-Vgq)]);
         
      else %High-frequency case
         %Vgq_old=Vgq;
         Vgq=Vgq0+Vac*sin(wt);
         %Find the quiescent point
         [Vgqcl,iminq]=min(abs(Vg-Vgq));
         usq=us(iminq);
           Cq=Chf(iminq);
         set(Hp_cvq,'XData',Vgq,'YData',Cq*1000);

         %Calculate altered concentration surface for the 2D graph
         wdepl=weff_hf(iminq)+(Vgq-Vgq0)*Cq/(q*Nsub*1e6);
         u=usq*((1-1e-9*rg./abs(wdepl)).^2).*(1e-9*rg<abs(wdepl));

         %Calculate the altered concentration
         Ncarr_c=ni*sign(uf_usd11-u).*exp(abs(uf_usd11-u))+sign(uf_usd11-u)*10;
         arcor=(Ncarr.*Nsub<0); %1 for different sign carriers; 
         Ncarr_c=arcor.*Ncarr+(1-arcor).*Ncarr_c; %not to alter the inversion l.


         %Set the 2D graph
         if Npoly==0; 
            material=(material==11).*2+(material~=11).*material; %Change gate to metal
            z2=(material>11).*sign(Ncarr_c).*log10(abs(Ncarr_c))+...
               (material==2).*2;
         else
            z2=(material>11).*sign(Ncarr_c).*log10(abs(Ncarr_c))+...
               (material==11).*Npoly;
         end   
         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));
    
         %Reset 1D graph
         yen1=-u(fix(nxsteps/2),:)*kt+1.12/2; %Energy band numbers
         yen1(1)=-uf_usd11*kt-Vgq+(Npoly>0)*1.12;
         yen1(2)=yen1(1);
         y1=z2(fix(nxsteps/2),:); %Log concentration numbers for coloring
         if Npoly==0;
            yen1(1)=-10; yen2(2)=-10; %Remove it from the graph
            y1c=(abs(y1)>bright).*(y1-sign(y1).*bright)*26.5/(z2_max-bright)+37.5;  
            y1c(1)=3;
            y1c(2)=6; y1c(3)=6; %same color as ox for gate<x<Si surface
            set(Hph_met,'YData',[-1.7 -1.7 (-uf_usd11*kt-Vgq) (-uf_usd11*kt-Vgq)],'Visible','on');
         else   
            y1(1)=Npoly; y1(2)=Npoly;
            y1c=(abs(y1)>bright).*(y1-sign(y1).*bright)*26.5/(z2_max-bright)+37.5;  
            y1c(2)=6; y1c(3)=6; %same color as ox for gate<x<Si surface
            set(Hph_met,'Visible','off');
         end
         set(Hp_efsub,'YData',[-uf_usd11*kt -uf_usd11*kt]);
         set(Hp_cb,'YData',[yen1;yen1],'CData',[y1c;y1c]);   
         set(Hp_vb,'YData',[(yen1-1.12);(yen1-1.12)],'CData',[y1c;y1c]);
         set(Hp_efsub,'Visible','on');
         set(Hp_efgate,'YData',[(-uf_usd11*kt-Vgq) (-uf_usd11*kt-Vgq)]);
      end; %if lfhf==1; 
      if ishandle(Hf)==0;
         return; 
      end;
      drawnow;
      tic;
      while toc<0.3;
         pause(0.1);
      end;
      if ishandle(Hf)==0;
         return; 
      end;
      
   end; %while con_stp_usd11<1;
      if ishandle(Hf)==0;
         return; 
      end;
      set([Hr_dop Hr_gate Hr_lfhf],'Enable','on');
      set([Hsl_dop Hsl_vg Hsl_ox Hsl_chrg Hsl_br],'Enable','on');
      if Npoly==0;
         set(Hsl_gate,'Enable','on');
      end;   
      set([Hb_close Hb_start],'Enable','on');
      set(Hb_stop,'Enable','off');
      con_stp_usd11=0;
      con_vg_usd11=1; %To reset the graphs at the DC Vgq;
end; %if con_strt_usd11>0;   

   if ishandle(Hf)==0;
      return; 
   end;
pause(1);
   if ishandle(Hf)==0;
      return; 
   end;

end; %while of main loop
close;

