function bjt_usd;
%The concept of diffusion 
%(Sections 6.1 and 6.3 in "Understanding Semiconductor Devices")

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

global convbe_usd12; global convcb_usd12; global conr_usd12; 
global con_start_usd12; global con_close_usd12; global con_n_usd12; 
global con_m_usd12; global con_q_usd12; global con_go_usd12;


%OPEN PARENT FIGURE
Hf=figure('NumberTitle','off','Name','Sections 6.1 and 6.3: BJT Principles and Modeling',...
   'Color',[0.8 0.8 0.8]);
Ha0=axes('Xlim',[0 15],'Ylim',[0 7.5],'Zlim',[-7.5 1.5],'Visible','off',...
  'Position',[0.01 0.0 .75 .975]);
   % 'Position',[0.01 0.15 .75 .825]);
set(gcf,'InvertHardCopy','off');
   
%READ THE INTRO TEXT              
introtext ={'UNDERSTANDING SEMICONDUCTOR DEVICES ';...
      'Sections 6.1 and 6.3: BJT Principles and Modeling' ;...
      ' ';...
      'Illustrations of the internal currents and the energy-band diagram of an NPN BJT';...
      'are presented along with the input and output characteristics. ';...
      'Base-to-emitter and collector-to-base voltages can be set at different values, ';
      'and the main technology parameters (emitter doping, and base doping and width) can ';...
      'be changed, to enable a study of the BJT principles. The BJT model (equations ';...
      'used to perform the calculations) can be viewed, and some questions to test/trigger ';...
      'understanding of the BJT effects are provideed. ';};

Ht_intro=text(0.1,0.6,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_go_usd12=1;');
con_go_usd12=0;
while con_go_usd12<1; 
   pause(0.2);  
   if ishandle(Hf)==0;
      return; 
   end;
end
con_go_usd12=0;
set(Ht_intro,'Visible','off');
set(Hp_cont,'Visible','off','Enable','off');
drawnow;

convbe_usd12=0; convcb_usd12=0; conr_usd12=0; con_start_usd12=0; 
con_close_usd12=0; con_n_usd12=0; con_m_usd12=0; con_q_usd12=0;

%SETTING THE CONTROL PANEL              
Hframe=uicontrol(Hf,'Units','normalized','Position',[0.77 0.01 0.22 0.98],...
   'Style','frame');
%Vbe
HframeVbe=uicontrol(Hf,'Units','normalized','Position',[0.775 0.775 0.21 0.21],...
   'Style','frame','BackgroundColor','y');
Hrbb(1)=uicontrol(Hf,'Units','normalized','Position',[0.78 0.91 0.025 0.03],...
   'Style','radio','CallBack','convbe_usd12=1;','Value',0,'BackgroundColor','y');
Hrbb(2)=uicontrol(Hf,'Units','normalized','Position',[0.81 0.91 0.03 0.03],...
   'Style','radio','CallBack','convbe_usd12=2;','Value',0,'BackgroundColor','y');
Hrbb(3)=uicontrol(Hf,'Units','normalized','Position',[0.86 0.91 0.03 0.03],...
   'Style','radio','CallBack','convbe_usd12=3;','Value',0,'BackgroundColor','y');
Hrbb(4)=uicontrol(Hf,'Units','normalized','Position',[0.895 0.91 0.03 0.03],...
   'Style','radio','CallBack','convbe_usd12=4;','Value',1,'BackgroundColor','y');
Hrbb(5)=uicontrol(Hf,'Units','normalized','Position',[0.925 0.91 0.03 0.03],...
   'Style','radio','CallBack','convbe_usd12=5;','Value',0,'BackgroundColor','y');
Hrbb(6)=uicontrol(Hf,'Units','normalized','Position',[0.95 0.91 0.03 0.03],...
   'Style','radio','CallBack','convbe_usd12=6;','Value',0,'BackgroundColor','y');
Htb(1)=uicontrol(Hf,'Units','normalized','Position',[0.78 0.885 0.025 0.027],...
   'Style','text','BackgroundColor','y','ForegroundColor','k','String','0',...
   'HorizontalAlignment','center');
Htb(2)=uicontrol(Hf,'Units','normalized','Position',[0.92 0.885 0.03 0.027],...
   'Style','text','BackgroundColor','y','ForegroundColor','k','String','0.7',...
      'HorizontalAlignment','center');
Vbe=[0 0.3 0.65 0.68 0.7 0.715];
iVbeq=4; Vbeq=Vbe(iVbeq);
Htb(3)=uicontrol(Hf,'Units','normalized','Position',[0.78 0.94 0.20 0.04],...
   'Style','text','String',['Vbe=' num2str(Vbeq) ' V'],'BackgroundColor','y','ForegroundColor','k');

%Vcb
Ht_vcbu=uicontrol(Hf,'Units','normalized','Position',[0.78 0.85 0.20 0.01],...
   'Style','text','BackgroundColor','y');
Hsl_vcb=uicontrol(Hf,'Units','normalized','Position',[0.78 0.81 0.20 0.04],...
   'Style','slider','min',-4*Vbeq,'max',4.35,'Value',2,'CallBack','convcb_usd12=1;');
Vcbq=(get(Hsl_vcb,'Value')*1000)/1000;
if Vcbq<0; Vcbq=Vcbq/4; end
Ht_vcbd=uicontrol(Hf,'Units','normalized','Position',[0.78 0.78 0.20 0.03],...
   'Style','text','String',['Vcb=' num2str(Vcbq) ' V'],...
   'BackgroundColor','y','ForegroundColor','k');

%Ne
Ht_nat=uicontrol(Hf,'Units','normalized','Position',[0.78 0.71 0.20 0.03],...
   'Style','text','String',['EMITTER  DOPING']);
Hsl_ne=uicontrol(Hf,'Units','normalized','Position',[0.78 0.67 0.20 0.04],...
   'Style','slider','min',16.7,'max',20,'Value',20,'CallBack','con_n_usd12=1;');
Ht_nad=uicontrol(Hf,'Units','normalized','Position',[0.78 0.64 0.20 0.03],...
   'Style','text','String',[num2str(round(10^(get(Hsl_ne,'Value')))) ' cm-3']);
Ne=round(10^(get(Hsl_ne,'Value')))*1e6;

%Nb
Ht_nbt=uicontrol(Hf,'Units','normalized','Position',[0.78 0.59 0.20 0.03],...
   'Style','text','String',['BASE  DOPING']);
Hsl_nb=uicontrol(Hf,'Units','normalized','Position',[0.78 0.55 0.20 0.04],...
   'Style','slider','min',16.7,'max',19,'Value',17,'CallBack','con_n_usd12=2;');
Ht_nbd=uicontrol(Hf,'Units','normalized','Position',[0.78 0.52 0.20 0.03],...
   'Style','text','String',[num2str(round(10^(get(Hsl_nb,'Value')))) ' cm-3']);
Nb=round(10^(get(Hsl_nb,'Value')))*1e6;

%Wb
Ht_wbt=uicontrol(Hf,'Units','normalized','Position',[0.78 0.47 0.20 0.03],...
   'Style','text','String',['BASE  WIDTH']);
Hsl_wb=uicontrol(Hf,'Units','normalized','Position',[0.78 0.43 0.20 0.04],...
   'Style','slider','min',0.5,'max',1.3,'Value',1,'CallBack','con_n_usd12=3;');
Ht_wbd=uicontrol(Hf,'Units','normalized','Position',[0.78 0.40 0.20 0.03],...
   'Style','text','String',[num2str(round(get(Hsl_wb,'Value')*100)/100) ' um']);
Wb=round(get(Hsl_wb,'Value')*100)/100*1e-6;

%MODEL
Hb_model=uicontrol(Hf,'Units','normalized','Position',[0.78 0.30 0.20 0.08],...
   'Style','pushbutton','String','VIEW BJT MODEL','CallBack','con_m_usd12=1;');

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

%QUESTIONS
Ht_qu=uicontrol(Hf,'Units','normalized','Position',[0.83 0.25 0.15 0.03],...
   'Style','text','String','QUESTIONS');
Hrbq(1)=uicontrol(Hf,'Units','normalized','Position',[0.78 0.22 0.025 0.03],...
   'Style','radio','CallBack','con_q_usd12=1;','Value',1);
Hrbq(2)=uicontrol(Hf,'Units','normalized','Position',[0.83 0.22 0.03 0.03],...
   'Style','radio','CallBack','con_q_usd12=2;','Value',0);
Hrbq(3)=uicontrol(Hf,'Units','normalized','Position',[0.86 0.22 0.03 0.03],...
   'Style','radio','CallBack','con_q_usd12=3;','Value',0);
Hrbq(4)=uicontrol(Hf,'Units','normalized','Position',[0.89 0.22 0.03 0.03],...
   'Style','radio','CallBack','con_q_usd12=4;','Value',0);
Hrbq(5)=uicontrol(Hf,'Units','normalized','Position',[0.92 0.22 0.03 0.03],...
   'Style','radio','CallBack','con_q_usd12=5;','Value',0);
Hrbq(6)=uicontrol(Hf,'Units','normalized','Position',[0.95 0.22 0.03 0.03],...
   'Style','radio','CallBack','con_q_usd12=6;','Value',0);
Htbq(1)=uicontrol(Hf,'Units','normalized','Position',[0.83 0.19 0.03 0.027],...
   'Style','text','ForegroundColor','k','String','1',...
      'HorizontalAlignment','center');
Htbq(2)=uicontrol(Hf,'Units','normalized','Position',[0.86 0.19 0.03 0.027],...
   'Style','text','ForegroundColor','k','String','2',...
      'HorizontalAlignment','center');
Htbq(3)=uicontrol(Hf,'Units','normalized','Position',[0.89 0.19 0.03 0.027],...
   'Style','text','ForegroundColor','k','String','3',...
      'HorizontalAlignment','center');
Htbq(4)=uicontrol(Hf,'Units','normalized','Position',[0.92 0.19 0.03 0.027],...
   'Style','text','ForegroundColor','k','String','4',...
      'HorizontalAlignment','center');
Htbq(5)=uicontrol(Hf,'Units','normalized','Position',[0.95 0.19 0.03 0.027],...
   'Style','text','ForegroundColor','k','String','5',...
      'HorizontalAlignment','center');
   
  
%THE COLOR MAP FOR THE BANDS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Nbands=32;
a=([1:Nbands].')*0; cbcol=[a a a]; vbcol=[a a a];
%CB - green
for i=1:Nbands/2; %the dark green part
   cbcol(i,1)=0; cbcol(i,3)=0;
   cbcol(i,2)=0.8*log(i)/log(Nbands/2);
end
for i=(Nbands/2+1):Nbands; %the light green part
   cbcol(i,1)=0.8*log(i-Nbands/2)/log(Nbands/2); 
   cbcol(i,3)=cbcol(i,1);
   cbcol(i,2)=0.8;
end
%VB - red
cbreak=round(Nbands/3);
for i=1:cbreak; %the dark red part
   vbcol(i,2)=0; cbcol(i,3)=0;
   vbcol(i,1)=0.8*log(i)/log(cbreak);
end
for i=(cbreak+1):Nbands; %the light green part
   vbcol(i,2)=0.8*log(i-cbreak)/log(Nbands-cbreak); 
   vbcol(i,3)=cbcol(i,1);
   vbcol(i,1)=0.8;
end

%CONSTANTS
%%%%%%%%%%%%%%%%%%
Vt=0.02585;
epsnor=sqrt(2*11.8*8.85e-12/1.6e-19);
Ae=1e-4; Ac=1e-4;
q=1.6e-19;
ni=1.02e16;
Dn=2e-3;
Dp=8e-4; 
%Lp=100e-6; 
kt=0.026;
Eg=1.12;

Xmax=5; Ymax=7.0;
Icmax=60; %To be reset later

alphat=0.993;
We=1e-6; 
Nc=1e21; %To be set by slider
Vbi_eb=kt*log(Ne*Nb/ni^2);
Vbi_cb=kt*log(Nc*Nb/ni^2);
Wmin=0.1e-6;
Webe=epsnor*sqrt((Vbi_eb+Vbeq)/(Ne*(1+Ne/Nb)));
Webb=epsnor*sqrt((Vbi_eb+Vbeq)/(Nb*(1+Nb/Ne)))+Wmin;
Wcbc=epsnor*sqrt((Vbi_cb+Vcbq)/(Nc*(1+Nc/Nb)));
Wcbb=epsnor*sqrt((Vbi_cb+Vcbq)/(Nb*(1+Nb/Nc)))+Wmin;
phife=kt*log(Ne/ni);
phifb=kt*log(Nb/ni);
phifc=kt*log(Nc/ni);

%SET THE ENERGY-BAND DIAGRAMS
Ha_eb=axes('XLim',[0 Xmax],'Ylim',[0 Ymax],'Position',[0.05 0.0 0.35 0.6],...
   'Color',[0.8 0.8 0.8],'Visible','off');
deltay=1.25*Eg/Nbands;
hold on;
%Base holes
Eve=5.5-Eg/2-phife; Evc=(5.5-Eg/2-phifc-(Vcbq+Vbeq));
ybho=[1:(Nbands+1)].*0; yheq=[1:(Nbands+1)].*0;
ybho(1)=5.5-Vbeq; yheq(1)=5.5;
ycho=[1:(Nbands+1)].*0; ycho(1)=5.5-Vcbq-Vbeq;
for i=1:Nbands;
   ybho(i+1)=ybho(i)-deltay;
   yheq(i+1)=yheq(i)-deltay;
   ycho(i+1)=ycho(i)-deltay;
   if Vcbq>0;
   deltax=((5.5-yheq(i+1))/(Vcbq+phifb))*(Wcbb+Wcbc);
   xhb=[(We+Wb-Wcbb+deltax)*1e6 (We-Webe)*1e6 0 0 (We-Webe)*1e6 (We+Wb-Wcbb+deltax)*1e6];
         if yheq(i+1)<Eve;
            Hp_vbb(i)=patch(xhb,[ybho(i+1) ybho(i+1) yheq(i+1) yheq(i) ybho(i) ybho(i)],'r');
         else   
            Hp_vbb(i)=patch(xhb,[ybho(i+1) ybho(i+1) (Eve-deltay) Eve ybho(i) ybho(i)],'r');
         end
   else
   xhb=[Xmax (We+Wb-Wcbb)*1e6 (We-Webe)*1e6 0 0 (We-Webe)*1e6 (We+Wb-Wcbb)*1e6 Xmax];
         if yheq(i+1)<Eve;
            ylft=yheq(i+1);
         else   
            ylft=Eve-deltay;
         end
         if ycho(i+1)<Evc;
            yrght=ycho(i+1);
         else
            yrght=Evc-deltay;
         end
         Hp_vbb(i)=patch(xhb,[yrght ybho(i+1) ybho(i+1) ylft (ylft+deltay) ybho(i) ybho(i)...
            (yrght+deltay)],'r');
   end      
   set(Hp_vbb(i),'FaceColor',vbcol(i,:),'EdgeColor',vbcol(i,:));
   % set(Hp_vbb(i)) - After Coll. holes;
end   
%Collector holes 
Evc=(5.5-Eg/2-phifc-(Vcbq+Vbeq));
xhc=[(We+Wb+Wcbc)*1e6 Xmax Xmax (We+Wb+Wcbc)*1e6];
ycho=[1:(Nbands+1)].*0; 
ycho(1)=5.5-Vcbq-Vbeq;
for i=1:Nbands;
   ycho(i+1)=ycho(i)-deltay;
   if ybho(i+1)<Evc;
      Hp_vbc(i)=patch(xhc,[ybho(i) ycho(i) ycho(i+1) ybho(i+1)],'r');
   else
      Hp_vbc(i)=patch(xhc,[Evc ycho(i) ycho(i+1) (Evc-deltay)],'r');
   end
   set(Hp_vbc(i),'FaceColor',vbcol(i,:),'EdgeColor',vbcol(i,:));
end   
if Vcbq<=0; 
   set(Hp_vbc,'Visible','off'); 
else
   set(Hp_vbc,'Visible','on'); 
end

%Collector electrons
ycel=[1:(Nbands+1)].*0;
ycel(1)=5.5-Vcbq-Vbeq; %Ef in E is set at 5.5;
for i=1:Nbands;
   ycel(i+1)=ycel(i)+deltay;
   if Vcbq>=0;
   deltax=((ycel(i+1)-5.5+Vcbq+Vbeq)/(Vcbq+phifc))*(Wcbb+Wcbc);
   xleft=(We+Wb+Wcbc-deltax);
   if xleft<(We+Wb-Wcbb); xleft=We+Wb-Wcbb; end
   xec=[Xmax xleft*1e6 xleft*1e6 Xmax];
   else
   xec=[Xmax (We+Wb-Wcbb) (We+Wb-Wcbb) Xmax];
   end
   Hp_cbc(i)=patch(xec,[ycel(i) ycel(i) ycel(i+1) ycel(i+1)],'g');
   set(Hp_cbc(i),'FaceColor',cbcol(i,:),'EdgeColor',cbcol(i,:));
end   

%Emitter electrons
xeb=[0 (We+Webb)*1e6 (We+Wb-Wcbb)*1e6 (We+Wb-Wcbb)*1e6 (We+Webb)*1e6 0];
Ecb=5.5+Eg/2+phifb-Vbeq;
yeel=[1:(Nbands+1)].*0; 
yeel(1)=5.5; %Ef in E is set at 5.5;
for i=1:Nbands;
   yeel(i+1)=yeel(i)+deltay;
   if ycel(i+1)>Ecb;
      Hp_cbe(i)=patch(xeb,[yeel(i) yeel(i) ycel(i) ycel(i+1) yeel(i+1) yeel(i+1)],'g');
   else
      Hp_cbe(i)=patch(xeb,[yeel(i) yeel(i) Ecb (Ecb+deltay) yeel(i+1) yeel(i+1)],'g');
   end   
   set(Hp_cbe(i),'FaceColor',cbcol(i,:),'EdgeColor',cbcol(i,:));
end   

%The energy gap
xeb=[0 (We-Webe) (We+Webb) (We+Wb-Wcbb) (We+Wb+Wcbc) Xmax/1e6...
    Xmax/1e6 (We+Wb+Wcbc) (We+Wb-Wcbb) (We+Webb) (We-Webe) 0];
%Ef in E is set at 5.5;
yeb1=5.5+Eg/2+[-phife -phife (phifb-Vbeq) (phifb-Vbeq) (-phifc-(Vcbq+Vbeq)) (-phifc-(Vcbq+Vbeq))];
yeb2=5.5-Eg/2+[(-phifc-(Vcbq+Vbeq)) (-phifc-(Vcbq+Vbeq)) (phifb-Vbeq) (phifb-Vbeq) -phife -phife];
yeb=[yeb1 yeb2];
Hp_ebm=patch(xeb*1e6,yeb,'w');
set(Hp_ebm,'EdgeColor','w');

%The Fermi levels
H_cb(1)=plot([0 (We-Webe)*1e6],[5.5 5.5],'-.y');
H_cb(2)=plot([(We+Webb)*1e6 (We+Wb-Wcbb)*1e6],[(5.5-Vbeq) (5.5-Vbeq)],'-.y');
H_cb(3)=plot([(We+Wb+Wcbc)*1e6 Xmax],[(5.5-Vcbq-Vbeq) (5.5-Vcbq-Vbeq)],'-.y');
set(H_cb,'Color',[0.8 0.8 0]);

Ht_eb(1)=text((Xmax+0.05),(5.5+Eg/2-phifc-(Vcbq+Vbeq)),'E_{C}','Color','g',...
   'FontSize',11,'FontAngle','italic','FontWeight','bold',...
   'HorizontalAlignment','left','VerticalAlignment','bottom');
Ht_eb(2)=text((Xmax+0.05),(5.5-(Vcbq+Vbeq)),'E_{F}','Color','y',...
   'FontSize',11,'FontAngle','italic','FontWeight','bold',...
   'HorizontalAlignment','left','VerticalAlignment','middle');
Ht_eb(3)=text((Xmax+0.05),(5.5-Eg/2-phifc-(Vcbq+Vbeq)),'E_{V}','Color','r',...
   'FontSize',11,'FontAngle','italic','FontWeight','bold',...
   'HorizontalAlignment','left','VerticalAlignment','top');

%SET THE ELECTRONS/HOLE AND THEIR TRACES
max_tr=10;
xtr=[1:max_tr]*0;
ytr=[1:max_tr]*0;
xel=0; yel=0; Nel=0;
for k=1:max_tr;
   Hp_tr(k)=plot(xtr(k),ytr(k),'-g','EraseMode','xor','Visible','off');
   Hp_el(k)=plot(xel,yel,'.g');
end   
set(Hp_el,'EraseMode','xor','MarkerSize',10,'Visible','off','Color',[0 0.7 0]);
set(Hp_tr,'Color',[0 0.7 0]);
xtrh=0; ytrh=0; xho=0; yho=0;
Hp_trh=plot(xtrh,ytrh,'-r','EraseMode','xor','Visible','off');
Hp_ho=plot(xho,yho,'.r');
set(Hp_ho,'EraseMode','xor','MarkerSize',10,'Visible','off');
set(Hp_ho,'Color',[1 0.2 0.2]);
set(Hp_trh,'Color',[1 0.2 0.2]);
hold off;

%SET Ic-Vce graph
%%%%%%%%%%%%%%%%%%%
Ha_icvce=axes('Xlim',[0 5],'Ylim',[0 Icmax],'Position',[0.51 0.12 0.25 0.40],...
   'Color',[0.8 0.8 0.8],'FontAngle','italic');
hold on;
xlabel('V_{CE} (V)'); ylabel('I_{C} (A/cm^{2})');

%THE Ic-Vce MODEL
Vce=[0:0.1:5];
for i=1:6;
   Vbc=Vbe(i)-Vce;
   Webea=epsnor*sqrt((Vbi_eb+Vbe(i))./(Ne*(1+Ne/Nb)));
   Wcbba=sign(Vbi_cb-Vbc).*epsnor.*sqrt(abs(Vbi_cb-Vbc)./(Nb*(1+Nb/Nc)))+Wmin;
   Wcbca=sign(Vbi_cb-Vbc).*epsnor.*sqrt(abs(Vbi_cb-Vbc)./(Nc*(1+Nc/Nb)));

   
   Wbb=(Wb-Webb-Wcbba);
   Wee=(We-Webea);
   Wcc=(Xmax/1e6-We-Wb-Wcbca);
   
   Ines=Ae*q*ni^2*Dn./(Wbb*Nb);
   Ipes=Ae*q*ni^2*Dp./(Wee*Ne);
   Incs=Ac*q*ni^2*Dn./(Wbb*Nb);
   Ipcs=Ac*q*ni^2*Dp./(Wcc*Nc);
   Ic=alphat*Ines.*(exp(Vbe(i)./Vt)-1)-(Incs+Ipcs).*(exp(Vbc./Vt)-1);
    
%Set the Ic-Vce characteristics   
   if ishandle(Hf)==0;
      return; 
   end;
Hp_o(i)=plot(Vce,Ic,'-k','LineWidth',1);
set(Hp_o(i),'Color',[0.5 0.5 0.5]);
end

%Adjust the Ic axis automatically
Icmax=max(Ic);
Ndig=10^floor(log10(Icmax));
Icmaxgr=Ndig*2*ceil((Icmax/Ndig)/2);
set(Ha_icvce,'Ylim',[0 Icmaxgr]);

%select a line for the current Vg
set(Hp_o(iVbeq),'Color','y','LineWidth',2);

%Set the operating point
   Vceq=Vcbq+Vbe(iVbeq);
   Wbq=Wb-Webb-Wcbb;
   Weq=We-Webe;
   Wcq=Xmax/1e6-We-Wb-Wcbc;
   Inesq=Ae*q*ni^2*Dn/(Wbq*Nb);
   Ipesq=Ae*q*ni^2*Dp/(Weq*Ne);
   Incsq=Ac*q*ni^2*Dn/(Wbq*Nb);
   Ipcsq=Ac*q*ni^2*Dp/(Wcq*Nc);
   Icq=alphat*Inesq*(exp(Vbe(iVbeq)/Vt)-1)-(Incsq+Ipcsq)*(exp(-Vcbq/Vt)-1);

   if ishandle(Hf)==0;
      return; 
   end;
hold on;
Hp_oq=plot(Vceq,Icq,'g.','MarkerSize',20);
%set(Hp_oq,'Color','g');
hold off;

%SET Ib-Vbe graph
%%%%%%%%%%%%%%%%%%%
Ha_ibvbe=axes('Xlim',[0 0.75],'Ylim',[0 Icmax/100],'Position',[0.55 0.65 0.21 0.31],...
   'Color',[0.8 0.8 0.8],'FontAngle','italic');
hold on;
xlabel('V_{BE} (V)'); ylabel('I_{B} (A/cm^{2})');

Vbe_ib=[0:0.1:0.4 0.5:0.01:0.72];
Ib=((1-alphat)*Inesq+Ipesq).*(exp(Vbe_ib./Vt)-1)+((1-alphat).*Incsq+Ipcsq).*(exp(-Vcbq/Vt)-1);
Hp_ibvbe=plot(Vbe_ib,Ib,'-k');
set(Hp_ibvbe,'Color',[0.5 0.5 0.5]);
Ibmax=max(Ib);
Ndig=10^floor(log10(Ibmax));
Ibmaxgr=Ndig*2*ceil((Ibmax/Ndig)/2);
set(Ha_ibvbe,'Ylim',[0 Ibmaxgr]);

Ibq=((1-alphat)*Inesq+Ipesq)*(exp(Vbe(iVbeq)/Vt)-1)+...
   ((1-alphat)*Ipcsq+Ipcsq)*(exp(-Vcbq/Vt)-1);
Hp_ibvbeq=plot(Vbe(iVbeq),Ibq,'y.','MarkerSize',20);
hold off;


%SET THE CROSS-SECTION DIAGRAM
Ymincs=-3; Ymaxcs=3.0; 
Ha_cs=axes('XLim',[-0.05*Xmax/0.4 (Xmax+0.05*Xmax/0.4)],...
   'Ylim',[Ymincs Ymaxcs],'Position',[0.0 0.6 0.45 0.4],...
   'Color',[0.8 0.8 0.8],'Visible','off');
hold on;
%Set the neutral regions
Hp_cs(1)=patch([0 Xmax Xmax 0],[-2 -2 2 2],'w'); %The depl-layers in the background
Hp_cs(2)=patch([0 (We-Webe) (We-Webe) 0]*1e6,[-2 -2 2 2],'g'); %The emitter
Hp_cs(3)=patch([(We+Webb) (We+Wb-Wcbb) (We+Wb-Wcbb) (We+Webb)]*1e6,...
   [-2 -2 2 2],'r');
Hp_cs(4)=patch([(We+Wb+Wcbc)*1e6 Xmax Xmax (We+Wb+Wcbc)*1e6],[-2 -2 2 2],'g');
set([Hp_cs(2) Hp_cs(4)],'FaceColor','g','EdgeColor','g');
set(Hp_cs(3),'FaceColor','r','EdgeColor','r');

Imax=max(Icmax,Ibmax);
xfrl=(-0.05*Xmax/0.4); 
xfrc=(We+Webb+0.5*(Wb-Webb-Wcbb))*1e6;
xfrr=(Xmax+0.05*Xmax/0.4);
xfr=[(xfrl+0.3) xfrl (xfrl+0.3) xfrl 0 0 Xmax Xmax xfrr Xmax...
     (Xmax+0.3) Xmax Xmax+0.3 Xmax Xmax xfrc xfrc xfrc (xfrc-0.1) xfrc (xfrc+0.1) xfrc 0 0];
yfr=[0.15 0 -0.15 0 0 2 2 0 0 0 0.15 0 -0.15 0 -2 -2 (-2-0.7) -2 -2.3 -2 -2.3 -2 -2 0];
Hp_fr=plot(xfr,yfr,'-k');

Ht_cs(1)=text(-0.25,0.6,'I_{E}','Color','k',...
   'FontSize',12,'FontAngle','italic','FontWeight','bold','HorizontalAlignment','center');
Ht_cs(2)=text((Xmax+0.45),0.6,'I_{C}','Color','k',...
   'FontSize',12,'FontAngle','italic','FontWeight','bold','HorizontalAlignment','center');
Ht_cs(3)=text((xfrc+0.3),-2.5,'I_{B}','Color','k',...
   'FontSize',12,'FontAngle','italic','FontWeight','bold','HorizontalAlignment','center');

%Set Ine current
wine=2*Inesq*(exp(Vbe(iVbeq)/Vt)-1)/Imax;
if wine>2; wine=2; end;
winc=2*Incsq*(exp(Vbc./Vt)-1)/Imax;
if winc>2; winc=2; end;
ynel=0.75-wine/2;
yneu=0.75+wine/2;
yncl=0.75-winc/2;
yncu=0.75+winc/2;
dyine=0.2*wine;
if dyine<0.15; dyine=0.15; end
dyinc=0.2*winc;
if dyinc<0.15; dyinc=0.15; end
xp1=0.5*We*1e6;
dxinc=1.3*winc/2;
if dxinc<0.3; dxinc=0.3; end
xp2=xp1+dxinc;
dxine=1.3*wine/2;
if dxine<0.3; dxine=0.3; end
%xp3=(We+Wb+1.2*Wcbc)*1e6;
xp3=Xmax-0.5*(Xmax-(We+Wb+Wcbc)*1e6)-dxine;
%if xp3>(Xmax-dxine); xp3=Xmax-dxine; end;
xp4=xp3+dxine;

xp=[xp1 xp3 xp3 xp4 xp3 xp3 xp1];
yp=[yneu yneu (yneu+dyine) 0.75 (ynel-dyine) ynel ynel];
%The following lines for saturation region
%xp=[xp3 xp2 xp2 xp1 xp1 xp3 xp3 xp4 xp3 xp3 xp1 xp1 xp2 xp2 xp3];
%yp=[yncu yncu (yncu+dyinc) 0.75 yneu yneu (yneu+dyine) 0.75...
%      (ynel-dyine) ynel ynel 0.75 (yncl-dyinc) yncl yncl];
Hp_ine=patch(xp,yp,'g');
set(Hp_ine,'FaceColor',[0 0.6 0],'EdgeColor',[0 0.6 0]);
Ht_cs(4)=text((xp1+0.1),(yneu+0.05),'I_{nE}','Color',[0 0.6 0],...
   'FontSize',12,'FontAngle','italic','FontWeight','bold',...
   'HorizontalAlignment','left','VerticalAlignment','bottom');
xalphat=[xfrc xfrc (xfrc-0.1) xfrc (xfrc+0.1) xfrc xfrc];
yalphat=[ynel (ynel-0.7) (ynel-0.7) (ynel-1) (ynel-0.7) (ynel-0.7) ynel];
Hp_alphat=patch(xalphat,yalphat,'g');
set(Hp_alphat,'FaceColor',[0 0.6 0],'EdgeColor',[0 0.6 0]);
   
%Set Ipe current
wipe=2*Ipesq*(exp(Vbe(iVbeq)/Vt)-1)/Imax;
if wipe>2; wipe=2; end;
ypel=-1.2-wipe/2;
ypeu=-1.2+wipe/2;
dyipe=0.2*wipe;
if dyipe<0.15; dyipe=0.15; end
xp1pe=xp1;
dxipe=1.3*wipe/2;
if dxipe<0.3; dxipe=0.3; end
xpp2=xp1+dxipe;
%xpp3=(We+0.5*(Wb-Wcbb))*1e6;
xpp3=xfrc-0.2;
if xpp3<(We+Webb)*1e6; xpp3=(We+Webb)*1e6; end

xpipe=[xpp3 xpp2 xpp2 xp1pe xpp2 xpp2 xpp3];
ypipe=[ypeu ypeu (ypeu+dyipe) -1.2 (ypel-dyipe) ypel ypel];
Hp_ipe=patch(xpipe,ypipe,'r');
set(Hp_ipe,'FaceColor',[0.7 0 0],'EdgeColor',[0.7 0 0]);
Ht_cs(5)=text((xpp3-0.1),(ypeu+0.05),'I_{pE}','Color',[0.7 0 0],...
   'FontSize',12,'FontAngle','italic','FontWeight','bold',...
   'HorizontalAlignment','right','VerticalAlignment','bottom');

%Set Ipc current ('Visible','off' for Vcbq<=0)
xpp4=(We+Wb-Wcbb)*1e6;
xpipc=[xp4 (xpp4+0.3) (xpp4+0.3)  xpp4 (xpp4+0.3) (xpp4+0.3) xp4];
ypipc=[-1.2 -1.2 -1.05 -1.2 -1.35 -1.2 -1.2];
Hp_ipc=patch(xpipc,ypipc,'r');
set(Hp_ipc,'FaceColor',[0.7 0 0],'EdgeColor',[0.7 0 0]);
Ht_cs(6)=text((xp3-0.1),(-1.2+0.05),'I_{pC}','Color',[0.7 0 0],...
   'FontSize',12,'FontAngle','italic','FontWeight','bold',...
   'HorizontalAlignment','right','VerticalAlignment','bottom');
hold off;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%SET THE BJT MODEL FIGURE%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Hf_model=figure('NumberTitle','off','Name','BJT Model',...
   'Units','normalized','Position',[0.2 0.2 0.55 0.4],'Color',[0.87 0.87 0.87]);
%   'MenuBar','none');
Ha_eq=axes('Xlim',[0 1],'Ylim',[0 1],'Parent',Hf_model,'Visible','off',...
   'Position',[0 0 1 1]); %Axes for equation position
%model=0;
close(Hf_model);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%SET THE QUESTION FIGURE%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Hquest=figure('NumberTitle','off','Name','Questions',...
   'Units','normalized','Position',[0.05 0.2 0.2 0.5],'Color',[0.87 0.87 0.87]);
%   'MenuBar','none');
Hqu=axes('Xlim',[0 1],'Ylim',[0 1],'Parent',Hquest,'Visible','off',...
   'Position',[0 0 1 1]); %Axes for question position
%quest=0;
close(Hquest);


con_close_usd12=0;
%%%%%%%%%%%%%%%%%%%%%%
%%%THE MAIN LOOP
%%%%%%%%%%%%%%%%%%%%%%%
while con_close_usd12<1;
   
%%%%%%%%%%%%%%%%%%%%%%%
%Vbe VOLTAGE IS CHANGED   
%%%%%%%%%%%%%%%%%%%%%%%
if convbe_usd12>0;
   iVbeq=convbe_usd12;
   set(Hrbb,'Value',0);
   set(Hrbb(iVbeq),'Value',1);
   Vbeq=Vbe(iVbeq);
   if (Vcbq+Vbeq)<0;
      Vcbq=-Vbeq;
      set(Hsl_vcb,'Value',4*Vcbq);
      convcb_usd12=2;
   end;   
   set(Hsl_vcb,'min',-4*Vbeq); %so that Vce-min=0;
   set(Htb(3),'String',['Vbe=' num2str(Vbeq) ' V']);

Webe=epsnor*sqrt((Vbi_eb+Vbeq)/(Ne*(1+Ne/Nb)));
Webb=epsnor*sqrt((Vbi_eb+Vbeq)/(Nb*(1+Nb/Ne)))+Wmin;
Wcbc=epsnor*sqrt((Vbi_cb+Vcbq)/(Nc*(1+Nc/Nb)));
Wcbb=epsnor*sqrt((Vbi_cb+Vcbq)/(Nb*(1+Nb/Nc)))+Wmin;
phife=kt*log(Ne/ni);
phifb=kt*log(Nb/ni);

%RESET THE ENERGY-BAND DIAGRAMS DUE TO Vbe CHANGE
%Base holes
Eve=5.5-Eg/2-phife; Evc=(5.5-Eg/2-phifc-(Vcbq+Vbeq));
ybho=[1:(Nbands+1)].*0; yheq=[1:(Nbands+1)].*0;
ybho(1)=5.5-Vbeq; yheq(1)=5.5;
ycho=[1:(Nbands+1)].*0; ycho(1)=5.5-Vcbq-Vbeq;
for i=1:Nbands;
   ybho(i+1)=ybho(i)-deltay;
   yheq(i+1)=yheq(i)-deltay;
   ycho(i+1)=ycho(i)-deltay;
   if Vcbq>0;
   deltax=((5.5-yheq(i+1))/(Vcbq+phifb))*(Wcbb+Wcbc);
   xhb=[(We+Wb-Wcbb+deltax)*1e6 (We-Webe)*1e6 0 0 (We-Webe)*1e6 (We+Wb-Wcbb+deltax)*1e6];
         if yheq(i+1)<Eve;
            set(Hp_vbb(i),'XData',xhb,...
               'YData',[ybho(i+1) ybho(i+1) yheq(i+1) yheq(i) ybho(i) ybho(i)]);
         else   
            set(Hp_vbb(i),'XData',xhb,...
               'YData',[ybho(i+1) ybho(i+1) (Eve-deltay) Eve ybho(i) ybho(i)]);
         end
   else
   xhb=[Xmax (We+Wb-Wcbb)*1e6 (We-Webe)*1e6 0 0 (We-Webe)*1e6 (We+Wb-Wcbb)*1e6 Xmax];
         if yheq(i+1)<Eve;
            ylft=yheq(i+1);
         else   
            ylft=Eve-deltay;
         end
         if ycho(i+1)<Evc;
            yrght=ycho(i+1);
         else
            yrght=Evc-deltay;
         end
         set(Hp_vbb(i),'XData',xhb,...
            'YData',[yrght ybho(i+1) ybho(i+1) ylft (ylft+deltay) ybho(i) ybho(i)...
            (yrght+deltay)]);
   end      
   %set(Hp_vbb(i),'FaceColor',vbcol(i,:),'EdgeColor',vbcol(i,:));
end   
%Collector holes 
Evc=(5.5-Eg/2-phifc-(Vcbq+Vbeq));
xhc=[(We+Wb+Wcbc)*1e6 Xmax Xmax (We+Wb+Wcbc)*1e6];
ycho=[1:(Nbands+1)].*0; 
ycho(1)=5.5-Vcbq-Vbeq;
for i=1:Nbands;
   ycho(i+1)=ycho(i)-deltay;
   if ybho(i+1)<Evc;
      set(Hp_vbc(i),'XData',xhc,'YData',[ybho(i) ycho(i) ycho(i+1) ybho(i+1)]);
   else
      set(Hp_vbc(i),'XData',xhc,'YData',[Evc ycho(i) ycho(i+1) (Evc-deltay)]);
   end
   %set(Hp_vbc(i),'FaceColor',vbcol(i,:),'EdgeColor',vbcol(i,:));
end   
if Vcbq<=0; 
   set(Hp_vbc,'Visible','off'); 
else
   set(Hp_vbc,'Visible','on'); 
end

%Collector electrons
ycel=[1:(Nbands+1)].*0;
ycel(1)=5.5-Vcbq-Vbeq; %Ef in E is set at 5.5;
for i=1:Nbands;
   ycel(i+1)=ycel(i)+deltay;
   if Vcbq>=0;
   deltax=((ycel(i+1)-5.5+Vcbq+Vbeq)/(Vcbq+phifc))*(Wcbb+Wcbc);
   xleft=(We+Wb+Wcbc-deltax);
   if xleft<(We+Wb-Wcbb); xleft=We+Wb-Wcbb; end
   xec=[Xmax xleft*1e6 xleft*1e6 Xmax];
   else
   xec=[Xmax (We+Wb-Wcbb) (We+Wb-Wcbb) Xmax];
   end
   set(Hp_cbc(i),'XData',xec,'YData',[ycel(i) ycel(i) ycel(i+1) ycel(i+1)]);
   %set(Hp_cbc(i),'FaceColor',cbcol(i,:),'EdgeColor',cbcol(i,:));
end   

%Emitter electrons
xeb=[0 (We+Webb)*1e6 (We+Wb-Wcbb)*1e6 (We+Wb-Wcbb)*1e6 (We+Webb)*1e6 0];
Ecb=5.5+Eg/2+phifb-Vbeq;
yeel=[1:(Nbands+1)].*0; 
yeel(1)=5.5; %Ef in E is set at 5.5;
for i=1:Nbands;
   yeel(i+1)=yeel(i)+deltay;
   if ycel(i+1)>Ecb;
      set(Hp_cbe(i),'XData',xeb,...
         'YData',[yeel(i) yeel(i) ycel(i) ycel(i+1) yeel(i+1) yeel(i+1)]);
   else
      set(Hp_cbe(i),'XData',xeb,...
         'YData',[yeel(i) yeel(i) Ecb (Ecb+deltay) yeel(i+1) yeel(i+1)]);
   end   
   %set(Hp_cbe(i),'FaceColor',cbcol(i,:),'EdgeColor',cbcol(i,:));
end   

%The energy gap
xeb=[0 (We-Webe) (We+Webb) (We+Wb-Wcbb) (We+Wb+Wcbc) Xmax/1e6...
    Xmax/1e6 (We+Wb+Wcbc) (We+Wb-Wcbb) (We+Webb) (We-Webe) 0];
%Ef in E is set at 5.5;
yeb1=5.5+Eg/2+[-phife -phife (phifb-Vbeq) (phifb-Vbeq) (-phifc-(Vcbq+Vbeq)) (-phifc-(Vcbq+Vbeq))];
yeb2=5.5-Eg/2+[(-phifc-(Vcbq+Vbeq)) (-phifc-(Vcbq+Vbeq)) (phifb-Vbeq) (phifb-Vbeq) -phife -phife];
yeb=[yeb1 yeb2];
set(Hp_ebm,'XData',xeb*1e6,'YData',yeb);

%The Fermi levels
set(H_cb(1),'XData',[0 (We-Webe)*1e6]);
set(H_cb(2),'XData',[(We+Webb)*1e6 (We+Wb-Wcbb)*1e6],'YData',[(5.5-Vbeq) (5.5-Vbeq)]);
set(H_cb(3),'Xdata',[(We+Wb+Wcbc)*1e6 Xmax],'YData',[(5.5-Vcbq-Vbeq) (5.5-Vcbq-Vbeq)]);
set(Ht_eb(1),'Position',[(Xmax+0.05) (5.5+Eg/2-phifc-(Vcbq+Vbeq))]);
set(Ht_eb(2),'Position',[(Xmax+0.05) (5.5-(Vcbq+Vbeq))]);
set(Ht_eb(3),'Position',[(Xmax+0.05) (5.5-Eg/2-phifc-(Vcbq+Vbeq))]);


%RESET Ic-Vce graph DUE TO Vbe CHANGE
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for i=1:6;
   Vbc=Vbe(i)-Vce;
   Webea=epsnor*sqrt((Vbi_eb+Vbe(i))./(Ne*(1+Ne/Nb)));
   Wcbba=sign(Vbi_cb-Vbc).*epsnor.*sqrt(abs(Vbi_cb-Vbc)./(Nb*(1+Nb/Nc)))+Wmin;
   Wcbca=sign(Vbi_cb-Vbc).*epsnor.*sqrt(abs(Vbi_cb-Vbc)./(Nc*(1+Nc/Nb)));
   
   Wbb=(Wb-Webb-Wcbba);
   Wee=(We-Webea);
   Wcc=(Xmax/1e6-We-Wb-Wcbca);
   
   Ines=Ae*q*ni^2*Dn./(Wbb*Nb);
   Ipes=Ae*q*ni^2*Dp./(Wee*Ne);
   Incs=Ac*q*ni^2*Dn./(Wbb*Nb);
   Ipcs=Ac*q*ni^2*Dp./(Wcc*Nc);
   Ic=alphat*Ines.*(exp(Vbe(i)./Vt)-1)-(Incs+Ipcs).*(exp(Vbc./Vt)-1);
    
%Set the Ic-Vce characteristics   
if ishandle(Hf)==0;
   if ishandle(Hf_model)==1; close(Hf_model); end;
   if ishandle(Hquest)==1; close(Hquest); end;
   return; 
end;
set(Hp_o(i),'YData',Ic);
end

%select a line for the current Vg
set(Hp_o,'Color',[0.5 0.5 0.5],'LineWidth',1);
set(Hp_o(iVbeq),'Color','y','LineWidth',2);
%Set the operating point
   Vceq=Vcbq+Vbe(iVbeq);
   Wbq=Wb-Webb-Wcbb;
   Weq=We-Webe;
   Wcq=Xmax/1e6-We-Wb-Wcbc;
   Inesq=Ae*q*ni^2*Dn/(Wbq*Nb);
   Ipesq=Ae*q*ni^2*Dp/(Weq*Ne);
   Incsq=Ac*q*ni^2*Dn/(Wbq*Nb);
   Ipcsq=Ac*q*ni^2*Dp/(Wcq*Nc);
   Icq=alphat*Inesq*(exp(Vbe(iVbeq)/Vt)-1)-(Incsq+Ipcsq)*(exp(-Vcbq/Vt)-1);
   set(Hp_oq,'XData',Vceq,'YData',Icq);
   
%RESET Ib-Vbe graph DUE TO Vbe CHANGE
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Ib=((1-alphat)*Inesq+Ipesq).*(exp(Vbe_ib./Vt)-1)+...
   ((1-alphat).*Incsq+Ipcsq).*(exp(-Vcbq/Vt)-1);
Ibmax=max(Ib);
Ndig=10^floor(log10(Ibmax));
Ibmaxgr=Ndig*2*ceil((Ibmax/Ndig)/2);
set(Ha_ibvbe,'Ylim',[0 Ibmaxgr]);
set(Hp_ibvbe,'YData',Ib);
Ibq=((1-alphat)*Inesq+Ipesq)*(exp(Vbe(iVbeq)/Vt)-1)+...
   ((1-alphat)*Ipcsq+Ipcsq)*(exp(-Vcbq/Vt)-1);
set(Hp_ibvbeq,'XData',Vbe(iVbeq),'YData',Ibq);

%RESET THE CROSS-SECTION DIAGRAM, DUE TO Vbe CHANGE
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Set the neutral regions
set(Hp_cs(2),'XData',[0 (We-Webe) (We-Webe) 0]*1e6); %The emitter
set(Hp_cs(3),'XData',[(We+Webb) (We+Wb-Wcbb) (We+Wb-Wcbb) (We+Webb)]*1e6); %The base
set(Hp_cs(4),'XData',[(We+Wb+Wcbc)*1e6 Xmax Xmax (We+Wb+Wcbc)*1e6]); %The collector

xfrc=(We+Webb+0.5*(Wb-Webb-Wcbb))*1e6;
xfr=[(xfrl+0.3) xfrl (xfrl+0.3) xfrl 0 0 Xmax Xmax xfrr Xmax...
     (Xmax+0.3) Xmax Xmax+0.3 Xmax Xmax xfrc xfrc xfrc (xfrc-0.1) xfrc (xfrc+0.1) xfrc 0 0];
set(Hp_fr,'XData',xfr);
set(Ht_cs(3),'Position',[(xfrc+0.3) -2.5]);

%Set Ine current
if Vbeq==0; 
   set([Hp_ine Hp_alphat],'Visible','off');
   set(Ht_cs(4),'Visible','off');
else
if (Vcbq+Vbeq)>0;
wine=2*Inesq*(exp(Vbe(iVbeq)/Vt)-1)/Imax;
if wine>2; wine=2; end;
winc=2*Incsq*(exp(-Vcbq/Vt)-1)/Imax;
if winc>2; winc=2; end;
ynel=0.75-wine/2;
yneu=0.75+wine/2;
yncl=0.75-winc/2;
yncu=0.75+winc/2;
dyine=0.2*wine;
if dyine<0.15; dyine=0.15; end
dyinc=0.2*winc;
if dyinc<0.15; dyinc=0.15; end
xp1=0.5*We*1e6;
dxinc=1.3*winc/2;
if dxinc<0.3; dxinc=0.3; end
xp2=xp1+dxinc;
dxine=1.3*wine/2;
if dxine<0.3; dxine=0.3; end
%xp3=(We+Wb+1.2*Wcbc)*1e6;
xp3=Xmax-0.5*(Xmax-(We+Wb+Wcbc)*1e6)-dxine;
%if xp3>(Xmax-dxine); xp3=Xmax-dxine; end;
xp4=xp3+dxine;

if Vcbq>=0;
xp=[xp1 xp3 xp3 xp4 xp3 xp3 xp1];
yp=[yneu yneu (yneu+dyine) 0.75 (ynel-dyine) ynel ynel];
set(Ht_cs(4),'String','I_{nE}','Position',[(xp1+0.1) (yneu+0.05)],'Visible','on');
else
%The following lines for saturation region
xp=[xp3 xp2 xp2 xp1 xp1 xp3 xp3 xp4 xp3 xp3 xp1 xp1 xp2 xp2 xp3];
yp=[yncu yncu (yncu+dyinc) 0.75 yneu yneu (yneu+dyine) 0.75...
      (ynel-dyine) ynel ynel 0.75 (yncl-dyinc) yncl yncl];
set(Ht_cs(4),'String','I_{nE}-I_{nCC}','Position',[(xp1+0.1) (yneu+0.05)],'Visible','on');
end
set(Hp_ine,'XData',xp,'YData',yp,'Visible','on');
xalphat=[xfrc xfrc (xfrc-0.1) xfrc (xfrc+0.1) xfrc xfrc];
yalphat=[ynel (ynel-0.7) (ynel-0.7) (ynel-1) (ynel-0.7) (ynel-0.7) ynel];
set(Hp_alphat,'XData',xalphat,'YData',yalphat,'Visible','on');
else
   set(Hp_ine,'Visible','off');
   set(Ht_cs(4),'Visible','off');
   set(Hp_alphat,'Visible','off');
end   
end

%Set Ipe current
if Vbeq>0;
wipe=2*Ipesq*(exp(Vbe(iVbeq)/Vt)-1)/Imax;
if wipe>2; wipe=2; end;
ypel=-1.2-wipe/2;
ypeu=-1.2+wipe/2;
dyipe=0.2*wipe;
if dyipe<0.15; dyipe=0.15; end
dxipe=1.3*wipe/2;
if dxipe<0.3; dxipe=0.3; end
xpp2=xp1+dxipe;
xpp3=xfrc-0.2;
if xpp3<(We+Webb)*1e6; xpp3=(We+Webb)*1e6; end
if xpp2<(xpp3-0.5);
   xp1pe=0.5*We*1e6;
else
   xpp2=xpp3-0.5;
   xp1pe=xpp2-dxipe;
end
xpipe=[xpp3 xpp2 xpp2 xp1pe xpp2 xpp2 xpp3];
ypipe=[ypeu ypeu (ypeu+dyipe) -1.2 (ypel-dyipe) ypel ypel];
set(Hp_ipe,'XData',xpipe,'YData',ypipe,'Visible','on');
set(Ht_cs(5),'Position',[(xpp3-0.1) (ypeu+0.05)],'Visible','on');
else
set(Hp_ipe,'Visible','off');
set(Ht_cs(5),'Visible','off');
end   

%Set Ipc current ('Visible','off' for Vcbq<=0)
xpp4=(We+Wb-Wcbb)*1e6;
if Vcbq>0;
xpipc=[xp4 (xpp4+0.3) (xpp4+0.3)  xpp4 (xpp4+0.3) (xpp4+0.3) xp4];
ypipc=[-1.2 -1.2 -1.05 -1.2 -1.35 -1.2 -1.2];
set(Hp_ipc,'XData',xpipc,'YData',ypipc,'Visible','on');
set(Ht_cs(6),'Position',[(xp3-0.1) (-1.2+0.05)],'Visible','on',...
         'HorizontalAlignment','right');
else
   if Vcbq==0;
      set(Hp_ipc,'Visible','off');
      set(Ht_cs(6),'Visible','off');
   else %Vcbq is negative, Ipc changes direction
      wipc=2*Ipcsq*(exp(-Vcbq/Vt)-1)/Imax;
      if wipc>2; wipc=2; end;
      ypcl=-1.2-wipc/2;
      ypcu=-1.2+wipc/2;
      dyipc=0.2*wipc;
      if dyipc<0.15; dyipc=0.15; end
      %xp1=0.5*We*1e6;
      dxipc=1.3*wipc/2;
      if dxipc<0.3; dxipc=0.3; end
      %xpp2=xp1+dxipe;
      xpc3=xfrc+0.2;
      if xpc3>(We+Wb-Wcbb)*1e6; xpc3=(We+Wb-Wcbb)*1e6; end
      xp33=xp4-dxipc;
      %xpipe=[xpp3 xpp2 xpp2 xp1 xpp2 xpp2 xpp3];
      %ypipe=[ypeu ypeu (ypeu+dyipe) -1.2 (ypel-dyipe) ypel ypel];
      xpipc=[xpc3 xp33 xp33 xp4 xp33 xp33 xpc3];
      ypipc=[ypcu ypcu (ypcu+dyipc) -1.2 (ypcl-dyipc) ypcl ypcl];
      set(Hp_ipc,'XData',xpipc,'YData',ypipc,'Visible','on');
      set(Ht_cs(6),'Position',[(xpc3+0.05) (ypcu+0.05)],'Visible','on',...
         'HorizontalAlignment','left');
   end   
end   

end; %if convbe_usd12>0;
convbe_usd12=0;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Vcbq IS CHANGED
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if convcb_usd12>0;
Vcbq=get(Hsl_vcb,'Value');   
if Vcbq<0; Vcbq=Vcbq/4; end   
Vcbq=round(Vcbq*1000)/1000;
if abs(Vcbq)<0.03; Vcbq=0; end;
if Vcbq>4.3; Vcbq=4.3; end;
if Vcbq<0; 
   set(Hsl_vcb,'Value',4*Vcbq);
else
   set(Hsl_vcb,'Value',Vcbq);
end   
set(Ht_vcbd,'String',['Vcb=' num2str(Vcbq) ' V']);

Webe=epsnor*sqrt((Vbi_eb+Vbeq)/(Ne*(1+Ne/Nb)));
Webb=epsnor*sqrt((Vbi_eb+Vbeq)/(Nb*(1+Nb/Ne)))+Wmin;
Wcbc=epsnor*sqrt((Vbi_cb+Vcbq)/(Nc*(1+Nc/Nb)));
Wcbb=epsnor*sqrt((Vbi_cb+Vcbq)/(Nb*(1+Nb/Nc)))+Wmin;

%RESET THE ENERGY-BAND DIAGRAMS DUE TO Vcb CHANGE
%Base holes
Eve=5.5-Eg/2-phife; Evc=(5.5-Eg/2-phifc-(Vcbq+Vbeq));
ybho=[1:(Nbands+1)].*0; yheq=[1:(Nbands+1)].*0;
ybho(1)=5.5-Vbeq; yheq(1)=5.5;
ycho=[1:(Nbands+1)].*0; ycho(1)=5.5-Vcbq-Vbeq;
for i=1:Nbands;
   ybho(i+1)=ybho(i)-deltay;
   yheq(i+1)=yheq(i)-deltay;
   ycho(i+1)=ycho(i)-deltay;
   if Vcbq>0;
   deltax=((5.5-yheq(i+1))/(Vcbq+phifb))*(Wcbb+Wcbc);
   xhb=[(We+Wb-Wcbb+deltax)*1e6 (We-Webe)*1e6 0 0 (We-Webe)*1e6 (We+Wb-Wcbb+deltax)*1e6];
         if yheq(i+1)<Eve;
            set(Hp_vbb(i),'XData',xhb,...
               'YData',[ybho(i+1) ybho(i+1) yheq(i+1) yheq(i) ybho(i) ybho(i)]);
         else   
            set(Hp_vbb(i),'XData',xhb,...
               'YData',[ybho(i+1) ybho(i+1) (Eve-deltay) Eve ybho(i) ybho(i)]);
         end
   else
   xhb=[Xmax (We+Wb-Wcbb)*1e6 (We-Webe)*1e6 0 0 (We-Webe)*1e6 (We+Wb-Wcbb)*1e6 Xmax];
         if yheq(i+1)<Eve;
            ylft=yheq(i+1);
         else   
            ylft=Eve-deltay;
         end
         if ycho(i+1)<Evc;
            yrght=ycho(i+1);
         else
            yrght=Evc-deltay;
         end
         set(Hp_vbb(i),'XData',xhb,...
            'YData',[yrght ybho(i+1) ybho(i+1) ylft (ylft+deltay) ybho(i) ybho(i)...
            (yrght+deltay)]);
   end      
   %set(Hp_vbb(i),'FaceColor',vbcol(i,:),'EdgeColor',vbcol(i,:));
end   
%Collector holes 
Evc=(5.5-Eg/2-phifc-(Vcbq+Vbeq));
xhc=[(We+Wb+Wcbc)*1e6 Xmax Xmax (We+Wb+Wcbc)*1e6];
ycho=[1:(Nbands+1)].*0; 
ycho(1)=5.5-Vcbq-Vbeq;
for i=1:Nbands;
   ycho(i+1)=ycho(i)-deltay;
   if ybho(i+1)<Evc;
      set(Hp_vbc(i),'XData',xhc,'YData',[ybho(i) ycho(i) ycho(i+1) ybho(i+1)]);
   else
      set(Hp_vbc(i),'XData',xhc,'YData',[Evc ycho(i) ycho(i+1) (Evc-deltay)]);
   end
   %set(Hp_vbc(i),'FaceColor',vbcol(i,:),'EdgeColor',vbcol(i,:));
end   
if Vcbq<=0; 
   set(Hp_vbc,'Visible','off'); 
else
   set(Hp_vbc,'Visible','on'); 
end

%Collector electrons
ycel=[1:(Nbands+1)].*0;
ycel(1)=5.5-Vcbq-Vbeq; %Ef in E is set at 5.5;
for i=1:Nbands;
   ycel(i+1)=ycel(i)+deltay;
   if Vcbq>=0;
   deltax=((ycel(i+1)-5.5+Vcbq+Vbeq)/(Vcbq+phifc))*(Wcbb+Wcbc);
   xleft=(We+Wb+Wcbc-deltax);
   if xleft<(We+Wb-Wcbb); xleft=We+Wb-Wcbb; end
   xec=[Xmax xleft*1e6 xleft*1e6 Xmax];
   else
   xec=[Xmax (We+Wb-Wcbb) (We+Wb-Wcbb) Xmax];
   end
   set(Hp_cbc(i),'XData',xec,'YData',[ycel(i) ycel(i) ycel(i+1) ycel(i+1)]);
   %set(Hp_cbc(i),'FaceColor',cbcol(i,:),'EdgeColor',cbcol(i,:));
end   

%Emitter electrons
xeb=[0 (We+Webb)*1e6 (We+Wb-Wcbb)*1e6 (We+Wb-Wcbb)*1e6 (We+Webb)*1e6 0];
Ecb=5.5+Eg/2+phifb-Vbeq;
yeel=[1:(Nbands+1)].*0; 
yeel(1)=5.5; %Ef in E is set at 5.5;
for i=1:Nbands;
   yeel(i+1)=yeel(i)+deltay;
   if ycel(i+1)>Ecb;
      set(Hp_cbe(i),'XData',xeb,...
         'YData',[yeel(i) yeel(i) ycel(i) ycel(i+1) yeel(i+1) yeel(i+1)]);
   else
      set(Hp_cbe(i),'XData',xeb,...
         'YData',[yeel(i) yeel(i) Ecb (Ecb+deltay) yeel(i+1) yeel(i+1)]);
   end   
   %set(Hp_cbe(i),'FaceColor',cbcol(i,:),'EdgeColor',cbcol(i,:));
end   

%The energy gap
xeb=[0 (We-Webe) (We+Webb) (We+Wb-Wcbb) (We+Wb+Wcbc) Xmax/1e6...
    Xmax/1e6 (We+Wb+Wcbc) (We+Wb-Wcbb) (We+Webb) (We-Webe) 0];
%Ef in E is set at 5.5;
yeb1=5.5+Eg/2+[-phife -phife (phifb-Vbeq) (phifb-Vbeq) (-phifc-(Vcbq+Vbeq)) (-phifc-(Vcbq+Vbeq))];
yeb2=5.5-Eg/2+[(-phifc-(Vcbq+Vbeq)) (-phifc-(Vcbq+Vbeq)) (phifb-Vbeq) (phifb-Vbeq) -phife -phife];
yeb=[yeb1 yeb2];
set(Hp_ebm,'XData',xeb*1e6,'YData',yeb);

%The Fermi levels
set(H_cb(1),'XData',[0 (We-Webe)*1e6]);
set(H_cb(3),'Xdata',[(We+Wb+Wcbc)*1e6 Xmax],'YData',[(5.5-Vcbq-Vbeq) (5.5-Vcbq-Vbeq)]);
set(Ht_eb(1),'Position',[(Xmax+0.05) (5.5+Eg/2-phifc-(Vcbq+Vbeq))]);
set(Ht_eb(2),'Position',[(Xmax+0.05) (5.5-(Vcbq+Vbeq))]);
set(Ht_eb(3),'Position',[(Xmax+0.05) (5.5-Eg/2-phifc-(Vcbq+Vbeq))]);

%RESET Ic-Vce graph DUE TO Vcb CHANGE
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for i=1:6;
   Vbc=Vbe(i)-Vce;
   Webea=epsnor*sqrt((Vbi_eb+Vbe(i))./(Ne*(1+Ne/Nb)));
   Wcbba=sign(Vbi_cb-Vbc).*epsnor.*sqrt(abs(Vbi_cb-Vbc)./(Nb*(1+Nb/Nc)))+Wmin;
   Wcbca=sign(Vbi_cb-Vbc).*epsnor.*sqrt(abs(Vbi_cb-Vbc)./(Nc*(1+Nc/Nb)));
   
   Wbb=(Wb-Webb-Wcbba);
   Wee=(We-Webea);
   Wcc=(Xmax/1e6-We-Wb-Wcbca);
   
   Ines=Ae*q*ni^2*Dn./(Wbb*Nb);
   Ipes=Ae*q*ni^2*Dp./(Wee*Ne);
   Incs=Ac*q*ni^2*Dn./(Wbb*Nb);
   Ipcs=Ac*q*ni^2*Dp./(Wcc*Nc);
   Ic=alphat*Ines.*(exp(Vbe(i)./Vt)-1)-(Incs+Ipcs).*(exp(Vbc./Vt)-1);
    
%Set the Ic-Vce characteristics   
if ishandle(Hf)==0;
   if ishandle(Hf_model)==1; close(Hf_model); end;
   if ishandle(Hquest)==1; close(Hquest); end;
   return; 
end;
set(Hp_o(i),'YData',Ic);
end

%select a line for the current Vg
set(Hp_o,'Color',[0.5 0.5 0.5],'LineWidth',1);
set(Hp_o(iVbeq),'Color','y','LineWidth',2);
%Set the operating point
   Vceq=Vcbq+Vbe(iVbeq);
   Wbq=Wb-Webb-Wcbb;
   Weq=We-Webe;
   Wcq=Xmax/1e6-We-Wb-Wcbc;
   Inesq=Ae*q*ni^2*Dn/(Wbq*Nb);
   Ipesq=Ae*q*ni^2*Dp/(Weq*Ne);
   Incsq=Ac*q*ni^2*Dn/(Wbq*Nb);
   Ipcsq=Ac*q*ni^2*Dp/(Wcq*Nc);
   Icq=alphat*Inesq*(exp(Vbe(iVbeq)/Vt)-1)-(Incsq+Ipcsq)*(exp(-Vcbq/Vt)-1);
   set(Hp_oq,'XData',Vceq,'YData',Icq);
   
%RESET Ib-Vbe graph DUE TO Vcb CHANGE
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Ib=((1-alphat)*Inesq+Ipesq).*(exp(Vbe_ib./Vt)-1)+...
   ((1-alphat).*Incsq+Ipcsq).*(exp(-Vcbq/Vt)-1);
Ibmax=max(Ib);
Ndig=10^floor(log10(Ibmax));
Ibmaxgr=Ndig*2*ceil((Ibmax/Ndig)/2);
set(Ha_ibvbe,'Ylim',[0 Ibmaxgr]);
set(Hp_ibvbe,'YData',Ib);
Ibq=((1-alphat)*Inesq+Ipesq)*(exp(Vbe(iVbeq)/Vt)-1)+...
   ((1-alphat)*Ipcsq+Ipcsq)*(exp(-Vcbq/Vt)-1);
set(Hp_ibvbeq,'XData',Vbe(iVbeq),'YData',Ibq);

%RESET THE CROSS-SECTION DIAGRAM, DUE TO Vcb CHANGE
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Set the neutral regions
set(Hp_cs(2),'XData',[0 (We-Webe) (We-Webe) 0]*1e6); %The emitter
set(Hp_cs(3),'XData',[(We+Webb) (We+Wb-Wcbb) (We+Wb-Wcbb) (We+Webb)]*1e6); %The base
set(Hp_cs(4),'XData',[(We+Wb+Wcbc)*1e6 Xmax Xmax (We+Wb+Wcbc)*1e6]); %The collector

xfrc=(We+Webb+0.5*(Wb-Webb-Wcbb))*1e6;
xfr=[(xfrl+0.3) xfrl (xfrl+0.3) xfrl 0 0 Xmax Xmax xfrr Xmax...
     (Xmax+0.3) Xmax Xmax+0.3 Xmax Xmax xfrc xfrc xfrc (xfrc-0.1) xfrc (xfrc+0.1) xfrc 0 0];
set(Hp_fr,'XData',xfr);
set(Ht_cs(3),'Position',[(xfrc+0.3) -2.5]);

%Set Ine current
if Vbeq==0; 
   set([Hp_ine Hp_alphat],'Visible','off');
   set(Ht_cs(4),'Visible','off');
else
if (Vcbq+Vbeq)>0;
wine=2*Inesq*(exp(Vbe(iVbeq)/Vt)-1)/Imax;
if wine>2; wine=2; end;
winc=2*Incsq*(exp(-Vcbq/Vt)-1)/Imax;
if winc>2; winc=2; end;
ynel=0.75-wine/2;
yneu=0.75+wine/2;
yncl=0.75-winc/2;
yncu=0.75+winc/2;
dyine=0.2*wine;
if dyine<0.15; dyine=0.15; end
dyinc=0.2*winc;
if dyinc<0.15; dyinc=0.15; end
xp1=0.5*We*1e6;
dxinc=1.3*winc/2;
if dxinc<0.3; dxinc=0.3; end
xp2=xp1+dxinc;
dxine=1.3*wine/2;
if dxine<0.3; dxine=0.3; end
%xp3=(We+Wb+1.2*Wcbc)*1e6;
xp3=Xmax-0.5*(Xmax-(We+Wb+Wcbc)*1e6)-dxine;
%if xp3>(Xmax-dxine); xp3=Xmax-dxine; end;
xp4=xp3+dxine;

if Vcbq>=0;
xp=[xp1 xp3 xp3 xp4 xp3 xp3 xp1];
yp=[yneu yneu (yneu+dyine) 0.75 (ynel-dyine) ynel ynel];
set(Ht_cs(4),'String','I_{nE}','Position',[(xp1+0.1) (yneu+0.05)],'Visible','on');
else
%The following lines for saturation region
xp=[xp3 xp2 xp2 xp1 xp1 xp3 xp3 xp4 xp3 xp3 xp1 xp1 xp2 xp2 xp3];
yp=[yncu yncu (yncu+dyinc) 0.75 yneu yneu (yneu+dyine) 0.75...
      (ynel-dyine) ynel ynel 0.75 (yncl-dyinc) yncl yncl];
set(Ht_cs(4),'String','I_{nE}-I_{nCC}','Position',[(xp1+0.1) (yneu+0.05)],'Visible','on');
end
set(Hp_ine,'XData',xp,'YData',yp,'Visible','on');
xalphat=[xfrc xfrc (xfrc-0.1) xfrc (xfrc+0.1) xfrc xfrc];
yalphat=[ynel (ynel-0.7) (ynel-0.7) (ynel-1) (ynel-0.7) (ynel-0.7) ynel];
set(Hp_alphat,'XData',xalphat,'YData',yalphat,'Visible','on');
else
   set(Hp_ine,'Visible','off');
   set(Ht_cs(4),'Visible','off');
   set(Hp_alphat,'Visible','off');
end   
end

%Set Ipe current
if Vbeq>0;
wipe=2*Ipesq*(exp(Vbe(iVbeq)/Vt)-1)/Imax;
if wipe>2; wipe=2; end;
ypel=-1.2-wipe/2;
ypeu=-1.2+wipe/2;
dyipe=0.2*wipe;
if dyipe<0.15; dyipe=0.15; end
dxipe=1.3*wipe/2;
if dxipe<0.3; dxipe=0.3; end
xpp2=xp1+dxipe;
xpp3=xfrc-0.2;
if xpp3<(We+Webb)*1e6; xpp3=(We+Webb)*1e6; end
if xpp2<(xpp3-0.5);
   xp1pe=0.5*We*1e6;
else
   xpp2=xpp3-0.5;
   xp1pe=xpp2-dxipe;
end
xpipe=[xpp3 xpp2 xpp2 xp1pe xpp2 xpp2 xpp3];
ypipe=[ypeu ypeu (ypeu+dyipe) -1.2 (ypel-dyipe) ypel ypel];
set(Hp_ipe,'XData',xpipe,'YData',ypipe,'Visible','on');
set(Ht_cs(5),'Position',[(xpp3-0.1) (ypeu+0.05)],'Visible','on');
else
set(Hp_ipe,'Visible','off');
set(Ht_cs(5),'Visible','off');
end   

%Set Ipc current ('Visible','off' for Vcbq<=0)
xpp4=(We+Wb-Wcbb)*1e6;
if Vcbq>0;
xpipc=[xp4 (xpp4+0.3) (xpp4+0.3)  xpp4 (xpp4+0.3) (xpp4+0.3) xp4];
ypipc=[-1.2 -1.2 -1.05 -1.2 -1.35 -1.2 -1.2];
set(Hp_ipc,'XData',xpipc,'YData',ypipc,'Visible','on');
set(Ht_cs(6),'Position',[(xp3-0.1) (-1.2+0.05)],'Visible','on',...
         'HorizontalAlignment','right');
else
   if Vcbq==0;
      set(Hp_ipc,'Visible','off');
      set(Ht_cs(6),'Visible','off');
   else %Vcbq is negative, Ipc changes direction
      wipc=2*Ipcsq*(exp(-Vcbq/Vt)-1)/Imax;
      if wipc>2; wipc=2; end;
      ypcl=-1.2-wipc/2;
      ypcu=-1.2+wipc/2;
      dyipc=0.2*wipc;
      if dyipc<0.15; dyipc=0.15; end
      %xp1=0.5*We*1e6;
      dxipc=1.3*wipc/2;
      if dxipc<0.3; dxipc=0.3; end
      %xpp2=xp1+dxipe;
      %xpc3=(We+0.6*(Wb-Wcbb))*1e6;
      xpc3=xfrc+0.2;
      if xpc3>(We+Wb-Wcbb)*1e6; xpc3=(We+Wb-Wcbb)*1e6; end
      xp33=xp4-dxipc;
      %xpipe=[xpp3 xpp2 xpp2 xp1 xpp2 xpp2 xpp3];
      %ypipe=[ypeu ypeu (ypeu+dyipe) -1.2 (ypel-dyipe) ypel ypel];
      xpipc=[xpc3 xp33 xp33 xp4 xp33 xp33 xpc3];
      ypipc=[ypcu ypcu (ypcu+dyipc) -1.2 (ypcl-dyipc) ypcl ypcl];
      set(Hp_ipc,'XData',xpipc,'YData',ypipc,'Visible','on');
      set(Ht_cs(6),'Position',[(xpc3+0.05) (ypcu+0.05)],'Visible','on',...
         'HorizontalAlignment','left');
   end   
end   

end; %if convcb_usd12>0;
convcb_usd12=0;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%EMITTER OR BASE DOPING, OR BASE WIDTH IS CHANGED
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if con_n_usd12>0;
   switch con_n_usd12;
   case 1;   
      set(Ht_nad,'String',[num2str(round(10^(get(Hsl_ne,'Value')))) ' cm-3']);
      Ne=round(10^(get(Hsl_ne,'Value')))*1e6;
   case 2;
      set(Ht_nbd,'String',[num2str(round(10^(get(Hsl_nb,'Value')))) ' cm-3']);
      Nb=round(10^(get(Hsl_nb,'Value')))*1e6;
   case 3;
      set(Ht_wbd,'String',[num2str(round(get(Hsl_wb,'Value')*100)/100) ' um']);
      Wb=round(get(Hsl_wb,'Value')*100)/100*1e-6;
   end;   
Vbi_eb=kt*log(Ne*Nb/ni^2);
Vbi_cb=kt*log(Nc*Nb/ni^2);
%Wmin=0.1e-6;
Webe=epsnor*sqrt((Vbi_eb+Vbeq)/(Ne*(1+Ne/Nb)));
Webb=epsnor*sqrt((Vbi_eb+Vbeq)/(Nb*(1+Nb/Ne)))+Wmin;
Wcbc=epsnor*sqrt((Vbi_cb+Vcbq)/(Nc*(1+Nc/Nb)));
Wcbb=epsnor*sqrt((Vbi_cb+Vcbq)/(Nb*(1+Nb/Nc)))+Wmin;
phife=kt*log(Ne/ni);
phifb=kt*log(Nb/ni);
%phifc=kt*log(Nc/ni);

%RESET THE ENERGY-BAND DIAGRAMS DUE TO Ne/Nb/Wb CHANGE
%Base holes
Eve=5.5-Eg/2-phife; Evc=(5.5-Eg/2-phifc-(Vcbq+Vbeq));
ybho=[1:(Nbands+1)].*0; yheq=[1:(Nbands+1)].*0;
ybho(1)=5.5-Vbeq; yheq(1)=5.5;
ycho=[1:(Nbands+1)].*0; ycho(1)=5.5-Vcbq-Vbeq;
for i=1:Nbands;
   ybho(i+1)=ybho(i)-deltay;
   yheq(i+1)=yheq(i)-deltay;
   ycho(i+1)=ycho(i)-deltay;
   if Vcbq>0;
   deltax=((5.5-yheq(i+1))/(Vcbq+phifb))*(Wcbb+Wcbc);
   xhb=[(We+Wb-Wcbb+deltax)*1e6 (We-Webe)*1e6 0 0 (We-Webe)*1e6 (We+Wb-Wcbb+deltax)*1e6];
         if yheq(i+1)<Eve;
            set(Hp_vbb(i),'XData',xhb,...
               'YData',[ybho(i+1) ybho(i+1) yheq(i+1) yheq(i) ybho(i) ybho(i)]);
         else   
            set(Hp_vbb(i),'XData',xhb,...
               'YData',[ybho(i+1) ybho(i+1) (Eve-deltay) Eve ybho(i) ybho(i)]);
         end
   else
   xhb=[Xmax (We+Wb-Wcbb)*1e6 (We-Webe)*1e6 0 0 (We-Webe)*1e6 (We+Wb-Wcbb)*1e6 Xmax];
         if yheq(i+1)<Eve;
            ylft=yheq(i+1);
         else   
            ylft=Eve-deltay;
         end
         if ycho(i+1)<Evc;
            yrght=ycho(i+1);
         else
            yrght=Evc-deltay;
         end
         set(Hp_vbb(i),'XData',xhb,...
            'YData',[yrght ybho(i+1) ybho(i+1) ylft (ylft+deltay) ybho(i) ybho(i)...
            (yrght+deltay)]);
   end      
   %set(Hp_vbb(i),'FaceColor',vbcol(i,:),'EdgeColor',vbcol(i,:));
end   
%Collector holes 
Evc=(5.5-Eg/2-phifc-(Vcbq+Vbeq));
xhc=[(We+Wb+Wcbc)*1e6 Xmax Xmax (We+Wb+Wcbc)*1e6];
ycho=[1:(Nbands+1)].*0; 
ycho(1)=5.5-Vcbq-Vbeq;
for i=1:Nbands;
   ycho(i+1)=ycho(i)-deltay;
   if ybho(i+1)<Evc;
      set(Hp_vbc(i),'XData',xhc,'YData',[ybho(i) ycho(i) ycho(i+1) ybho(i+1)]);
   else
      set(Hp_vbc(i),'XData',xhc,'YData',[Evc ycho(i) ycho(i+1) (Evc-deltay)]);
   end
   %set(Hp_vbc(i),'FaceColor',vbcol(i,:),'EdgeColor',vbcol(i,:));
end   
if Vcbq<=0; 
   set(Hp_vbc,'Visible','off'); 
else
   set(Hp_vbc,'Visible','on'); 
end

%Collector electrons
ycel=[1:(Nbands+1)].*0;
ycel(1)=5.5-Vcbq-Vbeq; %Ef in E is set at 5.5;
for i=1:Nbands;
   ycel(i+1)=ycel(i)+deltay;
   if Vcbq>=0;
   deltax=((ycel(i+1)-5.5+Vcbq+Vbeq)/(Vcbq+phifc))*(Wcbb+Wcbc);
   xleft=(We+Wb+Wcbc-deltax);
   if xleft<(We+Wb-Wcbb); xleft=We+Wb-Wcbb; end
   xec=[Xmax xleft*1e6 xleft*1e6 Xmax];
   else
   xec=[Xmax (We+Wb-Wcbb) (We+Wb-Wcbb) Xmax];
   end
   set(Hp_cbc(i),'XData',xec,'YData',[ycel(i) ycel(i) ycel(i+1) ycel(i+1)]);
   %set(Hp_cbc(i),'FaceColor',cbcol(i,:),'EdgeColor',cbcol(i,:));
end   

%Emitter electrons
xeb=[0 (We+Webb)*1e6 (We+Wb-Wcbb)*1e6 (We+Wb-Wcbb)*1e6 (We+Webb)*1e6 0];
Ecb=5.5+Eg/2+phifb-Vbeq;
yeel=[1:(Nbands+1)].*0; 
yeel(1)=5.5; %Ef in E is set at 5.5;
for i=1:Nbands;
   yeel(i+1)=yeel(i)+deltay;
   if ycel(i+1)>Ecb;
      set(Hp_cbe(i),'XData',xeb,...
         'YData',[yeel(i) yeel(i) ycel(i) ycel(i+1) yeel(i+1) yeel(i+1)]);
   else
      set(Hp_cbe(i),'XData',xeb,...
         'YData',[yeel(i) yeel(i) Ecb (Ecb+deltay) yeel(i+1) yeel(i+1)]);
   end   
   %set(Hp_cbe(i),'FaceColor',cbcol(i,:),'EdgeColor',cbcol(i,:));
end   

%The energy gap
xeb=[0 (We-Webe) (We+Webb) (We+Wb-Wcbb) (We+Wb+Wcbc) Xmax/1e6...
    Xmax/1e6 (We+Wb+Wcbc) (We+Wb-Wcbb) (We+Webb) (We-Webe) 0];
%Ef in E is set at 5.5;
yeb1=5.5+Eg/2+[-phife -phife (phifb-Vbeq) (phifb-Vbeq) (-phifc-(Vcbq+Vbeq)) (-phifc-(Vcbq+Vbeq))];
yeb2=5.5-Eg/2+[(-phifc-(Vcbq+Vbeq)) (-phifc-(Vcbq+Vbeq)) (phifb-Vbeq) (phifb-Vbeq) -phife -phife];
yeb=[yeb1 yeb2];
set(Hp_ebm,'XData',xeb*1e6,'YData',yeb);

%The Fermi levels
set(H_cb(1),'XData',[0 (We-Webe)*1e6]);
set(H_cb(2),'XData',[(We+Webb)*1e6 (We+Wb-Wcbb)*1e6],'YData',[(5.5-Vbeq) (5.5-Vbeq)]);
set(H_cb(3),'Xdata',[(We+Wb+Wcbc)*1e6 Xmax],'YData',[(5.5-Vcbq-Vbeq) (5.5-Vcbq-Vbeq)]);
set(Ht_eb(1),'Position',[(Xmax+0.05) (5.5+Eg/2-phifc-(Vcbq+Vbeq))]);
set(Ht_eb(2),'Position',[(Xmax+0.05) (5.5-(Vcbq+Vbeq))]);
set(Ht_eb(3),'Position',[(Xmax+0.05) (5.5-Eg/2-phifc-(Vcbq+Vbeq))]);

%RESET Ic-Vce graph DUE TO Ne/Nb/Wb CHANGE
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for i=1:6;
   Vbc=Vbe(i)-Vce;
   Webea=epsnor*sqrt((Vbi_eb+Vbe(i))./(Ne*(1+Ne/Nb)));
   Wcbba=sign(Vbi_cb-Vbc).*epsnor.*sqrt(abs(Vbi_cb-Vbc)./(Nb*(1+Nb/Nc)))+Wmin;
   Wcbca=sign(Vbi_cb-Vbc).*epsnor.*sqrt(abs(Vbi_cb-Vbc)./(Nc*(1+Nc/Nb)));
   
   Wbb=(Wb-Webb-Wcbba);
   Wee=(We-Webea);
   Wcc=(Xmax/1e6-We-Wb-Wcbca);
   
   Ines=Ae*q*ni^2*Dn./(Wbb*Nb);
   Ipes=Ae*q*ni^2*Dp./(Wee*Ne);
   Incs=Ac*q*ni^2*Dn./(Wbb*Nb);
   Ipcs=Ac*q*ni^2*Dp./(Wcc*Nc);
   Ic=alphat*Ines.*(exp(Vbe(i)./Vt)-1)-(Incs+Ipcs).*(exp(Vbc./Vt)-1);
    
%Set the Ic-Vce characteristics   
if ishandle(Hf)==0;
   if ishandle(Hf_model)==1; close(Hf_model); end;
   if ishandle(Hquest)==1; close(Hquest); end;
   return; 
end;
set(Hp_o(i),'YData',Ic);
end

%Adjust the Ic axis automatically
Icmax=max(Ic);
Ndig=10^floor(log10(Icmax));
Icmaxgr=Ndig*2*ceil((Icmax/Ndig)/2);
set(Ha_icvce,'Ylim',[0 Icmaxgr]);

%select a line for the current Vg
if ishandle(Hf)==0;
   if ishandle(Hf_model)==1; close(Hf_model); end;
   if ishandle(Hquest)==1; close(Hquest); end;
   return; 
end;
set(Hp_o,'Color',[0.5 0.5 0.5],'LineWidth',1);
set(Hp_o(iVbeq),'Color','y','LineWidth',2);
%Set the operating point
   Vceq=Vcbq+Vbe(iVbeq);
   Wbq=Wb-Webb-Wcbb;
   Weq=We-Webe;
   Wcq=Xmax/1e6-We-Wb-Wcbc;
   Inesq=Ae*q*ni^2*Dn/(Wbq*Nb);
   Ipesq=Ae*q*ni^2*Dp/(Weq*Ne);
   Incsq=Ac*q*ni^2*Dn/(Wbq*Nb);
   Ipcsq=Ac*q*ni^2*Dp/(Wcq*Nc);
   Icq=alphat*Inesq*(exp(Vbe(iVbeq)/Vt)-1)-(Incsq+Ipcsq)*(exp(-Vcbq/Vt)-1);
   set(Hp_oq,'XData',Vceq,'YData',Icq);
   
%RESET Ib-Vbe graph DUE TO Ne/Nb/Wb CHANGE
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Ib=((1-alphat)*Inesq+Ipesq).*(exp(Vbe_ib./Vt)-1)+...
   ((1-alphat).*Incsq+Ipcsq).*(exp(-Vcbq/Vt)-1);
Ibmax=max(Ib);
Ndig=10^floor(log10(Ibmax));
Ibmaxgr=Ndig*2*ceil((Ibmax/Ndig)/2);
set(Ha_ibvbe,'Ylim',[0 Ibmaxgr]);
set(Hp_ibvbe,'YData',Ib);
Ibq=((1-alphat)*Inesq+Ipesq)*(exp(Vbe(iVbeq)/Vt)-1)+...
   ((1-alphat)*Ipcsq+Ipcsq)*(exp(-Vcbq/Vt)-1);
set(Hp_ibvbeq,'XData',Vbe(iVbeq),'YData',Ibq);

%RESET THE CROSS-SECTION DIAGRAM, DUE TO Ne/Nb/Wb CHANGE
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Set the neutral regions
set(Hp_cs(2),'XData',[0 (We-Webe) (We-Webe) 0]*1e6); %The emitter
set(Hp_cs(3),'XData',[(We+Webb) (We+Wb-Wcbb) (We+Wb-Wcbb) (We+Webb)]*1e6); %The base
set(Hp_cs(4),'XData',[(We+Wb+Wcbc)*1e6 Xmax Xmax (We+Wb+Wcbc)*1e6]); %The collector

Imax=max(Icmax,Ibmax);

xfrc=(We+Webb+0.5*(Wb-Webb-Wcbb))*1e6;
xfr=[(xfrl+0.3) xfrl (xfrl+0.3) xfrl 0 0 Xmax Xmax xfrr Xmax...
     (Xmax+0.3) Xmax Xmax+0.3 Xmax Xmax xfrc xfrc xfrc (xfrc-0.1) xfrc (xfrc+0.1) xfrc 0 0];
set(Hp_fr,'XData',xfr);
set(Ht_cs(3),'Position',[(xfrc+0.3) -2.5]);

%Set Ine current
if (Vcbq+Vbeq)>0;
wine=2*Inesq*(exp(Vbe(iVbeq)/Vt)-1)/Imax;
if wine>2; wine=2; end;
winc=2*Incsq*(exp(-Vcbq/Vt)-1)/Imax;
if winc>2; winc=2; end;
ynel=0.75-wine/2;
yneu=0.75+wine/2;
yncl=0.75-winc/2;
yncu=0.75+winc/2;
dyine=0.2*wine;
if dyine<0.15; dyine=0.15; end
dyinc=0.2*winc;
if dyinc<0.15; dyinc=0.15; end
xp1=0.5*We*1e6;
dxinc=1.3*winc/2;
if dxinc<0.3; dxinc=0.3; end
xp2=xp1+dxinc;
dxine=1.3*wine/2;
if dxine<0.3; dxine=0.3; end
%xp3=(We+Wb+1.2*Wcbc)*1e6;
xp3=Xmax-0.5*(Xmax-(We+Wb+Wcbc)*1e6)-dxine;
%if xp3>(Xmax-dxine); xp3=Xmax-dxine; end;
xp4=xp3+dxine;

if Vcbq>=0;
xp=[xp1 xp3 xp3 xp4 xp3 xp3 xp1];
yp=[yneu yneu (yneu+dyine) 0.75 (ynel-dyine) ynel ynel];
set(Ht_cs(4),'String','I_{nE}','Position',[(xp1+0.1) (yneu+0.05)],'Visible','on');
else
%The following lines for saturation region
xp=[xp3 xp2 xp2 xp1 xp1 xp3 xp3 xp4 xp3 xp3 xp1 xp1 xp2 xp2 xp3];
yp=[yncu yncu (yncu+dyinc) 0.75 yneu yneu (yneu+dyine) 0.75...
      (ynel-dyine) ynel ynel 0.75 (yncl-dyinc) yncl yncl];
set(Ht_cs(4),'String','I_{nE}-I_{nCC}','Position',[(xp1+0.1) (yneu+0.05)],'Visible','on');
end
set(Hp_ine,'XData',xp,'YData',yp,'Visible','on');
xalphat=[xfrc xfrc (xfrc-0.1) xfrc (xfrc+0.1) xfrc xfrc];
yalphat=[ynel (ynel-0.7) (ynel-0.7) (ynel-1) (ynel-0.7) (ynel-0.7) ynel];
set(Hp_alphat,'XData',xalphat,'YData',yalphat,'Visible','on');
else
   set(Hp_ine,'Visible','off');
   set(Ht_cs(4),'Visible','off');
   set(Hp_alphat,'Visible','off');
end   

%Set Ipe current
if Vbeq>0;
wipe=2*Ipesq*(exp(Vbe(iVbeq)/Vt)-1)/Imax;
if wipe>2; wipe=2; end;
ypel=-1.2-wipe/2;
ypeu=-1.2+wipe/2;
dyipe=0.2*wipe;
if dyipe<0.15; dyipe=0.15; end
dxipe=1.3*wipe/2;
if dxipe<0.3; dxipe=0.3; end
xpp2=xp1+dxipe;
xpp3=xfrc-0.2;
if xpp3<(We+Webb)*1e6; xpp3=(We+Webb)*1e6; end
if xpp2<(xpp3-0.5);
   xp1pe=0.5*We*1e6;
else
   xpp2=xpp3-0.5;
   xp1pe=xpp2-dxipe;
end
xpipe=[xpp3 xpp2 xpp2 xp1pe xpp2 xpp2 xpp3];
ypipe=[ypeu ypeu (ypeu+dyipe) -1.2 (ypel-dyipe) ypel ypel];
set(Hp_ipe,'XData',xpipe,'YData',ypipe,'Visible','on');
set(Ht_cs(5),'Position',[(xpp3-0.1) (ypeu+0.05)],'Visible','on');
else
set(Hp_ipe,'Visible','off');
set(Ht_cs(5),'Visible','off');
end   

%Set Ipc current 
xpp4=(We+Wb-Wcbb)*1e6;
if Vcbq>0;
xpipc=[xp4 (xpp4+0.3) (xpp4+0.3)  xpp4 (xpp4+0.3) (xpp4+0.3) xp4];
ypipc=[-1.2 -1.2 -1.05 -1.2 -1.35 -1.2 -1.2];
set(Hp_ipc,'XData',xpipc,'YData',ypipc,'Visible','on');
set(Ht_cs(6),'Position',[(xp3-0.1) (-1.2+0.05)],'Visible','on',...
         'HorizontalAlignment','right');
else
   if Vcbq==0;
      set(Hp_ipc,'Visible','off');
      set(Ht_cs(6),'Visible','off');
   else %Vcbq is negative, Ipc changes direction
      wipc=2*Ipcsq*(exp(-Vcbq/Vt)-1)/Imax;
      if wipc>2; wipc=2; end;
      ypcl=-1.2-wipc/2;
      ypcu=-1.2+wipc/2;
      dyipc=0.2*wipc;
      if dyipc<0.15; dyipc=0.15; end
      %xp1=0.5*We*1e6;
      dxipc=1.3*wipc/2;
      if dxipc<0.3; dxipc=0.3; end
      %xpp2=xp1+dxipe;
      xpc3=xfrc+0.2;
      if xpc3>(We+Wb-Wcbb)*1e6; xpc3=(We+Wb-Wcbb)*1e6; end
      xp33=xp4-dxipc;
      %xpipe=[xpp3 xpp2 xpp2 xp1 xpp2 xpp2 xpp3];
      %ypipe=[ypeu ypeu (ypeu+dyipe) -1.2 (ypel-dyipe) ypel ypel];
      xpipc=[xpc3 xp33 xp33 xp4 xp33 xp33 xpc3];
      ypipc=[ypcu ypcu (ypcu+dyipc) -1.2 (ypcl-dyipc) ypcl ypcl];
      set(Hp_ipc,'XData',xpipc,'YData',ypipc,'Visible','on');
      set(Ht_cs(6),'Position',[(xpc3+0.05) (ypcu+0.05)],'Visible','on',...
         'HorizontalAlignment','left');
   end   
end   

end; % if con_n_usd12>0
con_n_usd12=0;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%THE MODEL WINDOW%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if con_m_usd12>0;
   if ishandle(Hquest)==1; 
   set(Hrbq,'Value',0);
   set(Hrbq(1),'Value',1);
   close(Hquest); 
   end
   
if ishandle(Hf_model)==0;

Hf_model=figure('NumberTitle','off','Name','BJT Model',...
   'Units','normalized','Position',[0.1 0.1 0.55 0.4],'Color',[0.8 0.8 0.8]);
%'MenuBar','none');
%axes(Ha_eq);
Ha_eq=axes('Xlim',[0 1],'Ylim',[0 1],'Visible','off',...
   'Position',[0 0 1 1]); %Axes for equation position
set(Ha_eq,'Parent',Hf_model);

hold on;
xm=[1.5 1.0 1.0 0.5 1.0 1.0 1.5]/10;
ym=[8.5 8.5 7.7 7.7 7.7 6.9 6.9]/10;
Hp_eq(1)=plot(xm,ym,'-k'); %E wires
xm=[1.50 1.50 1.50 1.80 1.80 1.50]/10;
ym=[8.25 8.75 8.50 8.75 8.25 8.50]/10;
Hp_eq(2)=plot(xm,ym,'-k'); %The E-B diode
xm=[1.90 1.30 1.45 1.30 1.45]/10;
ym=[9.00 9.00 9.07 9.00 8.93]/10;
Hp_eq(3)=plot(xm,ym,'-k'); %InE+IpE arrow
xm=[0.18 0.35]; ym=[0.85 0.85];
Hp_eq(4)=plot(xm,ym,'-k'); %E--B dioode to E--C CS wire
xm=[0.41 0.46 0.46 0.51 0.46 0.46 0.4];
ym=[0.85 0.85 0.77 0.77 0.77 0.69 0.69];
Hp_eq(5)=plot(xm,ym,'-k'); %Coll. wires
xm=[0.400 0.400 0.40 0.370 0.370 0.40];
ym=[0.715 0.665 0.69 0.715 0.665 0.69];
Hp_eq(6)=plot(xm,ym,'-k'); %C--B diode
xm=[0.37 0.28 0.28 0.28 0.28 0.21];
ym=[0.69 0.69 0.85 0.64 0.69 0.69];
Hp_eq(7)=plot(xm,ym,'-k'); %B wires
xml=[3.5 3.525 3.55 3.6 3.65 3.7 3.75 3.8 3.85 3.9 3.95 4.0 4.05 4.075 4.1]/10;
xmu=[4.1 4.075 4.05 4.0 3.95 3.9 3.85 3.8 3.75 3.7 3.65 3.6 3.55 3.525 3.5]/10;
yml=-1.3*0.03*sqrt(abs(1-((xml-0.38)/0.03).^2))+0.85;
ymu=1.3*0.03*sqrt(abs(1-((xml-0.38)/0.03).^2))+0.85;
Hp_eq(8)=plot([xml xmu],[yml ymu],'-k'); %C--B cs/circle
xm=[0.4 0.36 0.375 0.36 0.375];
ym=[0.850 0.85 0.857 0.85 0.843];
Hp_eq(9)=plot(xm,ym,'-k'); %C--B cs/arrow
xml=xml-0.38+0.18; xmu=xmu-0.38+0.18;
yml=yml-0.85+0.69; ymu=ymu-0.85+0.69;
Hp_eq(10)=plot([xml xmu],[yml ymu],'-k'); %E--B cs/circle
xm=[0.16 0.20 0.185 0.20 0.185];
ym=[0.69 0.69 0.683 0.69 0.697];
Hp_eq(11)=plot(xm,ym,'-k'); %E--B cs/arrow
xm=[0.355 0.415 0.40 0.415 0.40];
ym=[0.64 0.64 0.647 0.64 0.633];
Hp_eq(12)=plot(xm,ym,'-k'); %InCC+Ipc arrow
xm=[0.10 0.28 0.28 0.46];
ym=[0.77 0.85 0.69 0.77];
Hp_dots=plot(xm,ym,'.k'); %Dots

Ht_eqc(1)=text(0.16,0.95,'I_{nE}+I_{pE}',...
   'HorizontalAlignment','center','FontSize',10);
Ht_eqc(2)=text(0.38,0.95,'\alpha_{T}I_{nE}=\alpha(I_{nE}+I_{pE})',...
   'HorizontalAlignment','center','FontSize',10);
Ht_eqc(3)=text(0.385,0.59,'I_{nCC}+I_{pC}',...
   'HorizontalAlignment','center','FontSize',10);
Ht_eqc(4)=text(0.18,0.60,'\alpha_{T}I_{nCC}=\alpha(I_{nCC}+I_{pC})',...
   'HorizontalAlignment','center','FontSize',10);

Ht_eq(1)=text(0.05,0.47,...
       'I_{C}=\alpha_{T} I_{nE}-I_{pC}-I_{nCC}','Color','k','FontSize',12,...
       'FontAngle','italic','FontWeight','bold');
Ht_eq(2)=text(0.05,0.33,...
       'I_{E}=I_{nE}+I_{pE}-\alpha_{T} I_{nCC}','Color','k','FontSize',12,...
       'FontAngle','italic','FontWeight','bold');
Ht_eq(3)=text(0.05,0.19,...
       'I_{B}=I_{pE}+I_{pC}','Color','k','FontSize',12,...
       'FontAngle','italic','FontWeight','bold');
Ht_eq(4)=text(0.1,0.08,...
       '+(1-\alpha_{T} )I_{nE}+(1-\alpha_{T})I_{nCC}','Color','k','FontSize',12,...
       'FontAngle','italic','FontWeight','bold');
        
Ht_eq(5)=text(0.6,0.925,...
       'I_{nE}=I_{nES} (e^{Vbe/Vt}-1)','Color','k','FontSize',12,...
       'FontAngle','italic');
Ht_eq(6)=text(0.6,0.825,...
       'I_{pE}=I_{pES} (e^{Vbe/Vt}-1)','Color','k','FontSize',12,...
       'FontAngle','italic');
Ht_eq(7)=text(0.6,0.65,...
       'I_{nCC}=I_{nCS} (e^{Vbc/Vt}-1)','Color','k','FontSize',12,...
       'FontAngle','italic');
Ht_eq(8)=text(0.6,0.55,...
       'I_{pC}  =I_{pCS} (e^{Vbc/Vt}-1)','Color','k','FontSize',12,...
       'FontAngle','italic');
Ht_eq(9)=text(0.65,0.40,...
       'I_{nES}=A_{E}qn_{i}^{2}D_{n} /(w_{Beff }N_{B} )','Color',[0.5 0.5 0.5],'FontSize',10,...
       'FontAngle','italic');
Ht_eq(10)=text(0.65,0.30,...
       'I_{pES}=A_{E}qn_{i}^{2}D_{p} /(w_{Eeff }N_{E} )','Color',[0.5 0.5 0.5],'FontSize',10,...
       'FontAngle','italic');
Ht_eq(11)=text(0.65,0.175,...
       'I_{nCS}=A_{C}qn_{i}^{2}D_{n} /(w_{Beff }N_{B} )','Color',[0.5 0.5 0.5],'FontSize',10,...
       'FontAngle','italic');
Ht_eq(12)=text(0.65,0.075,...
   'I_{pCS}=A_{C}qn_{i}^{2}D_{p} /(w_{Ceff }N_{C} )','Color',[0.5 0.5 0.5],'FontSize',10,...
   'FontAngle','italic');
hold off;    

else
   axes(Ha_eq);
end   

con_m_usd12=0;
end; %con_m_usd12 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%THE QUESTIONS %%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if con_q_usd12>1;
   if ishandle(Hf_model)==1; close(Hf_model); end
if ishandle(Hquest)==0; %Create the figure and the handles
Hquest=figure('NumberTitle','off','Name','Question',...
   'Units','normalized','Position',[0.03 0.15 0.65 0.16],'Color',[0.8 0.8 0.8],...
   'MenuBar','none');
Hqu=axes('Xlim',[0 1],'Ylim',[0 1],'Visible','off',...
   'Position',[0 0 1 1]); %Axes for question position
set(Hqu,'Parent',Hquest);

Ht_q(1)=text(0.05,0.88,...
   'Reduce emitter doping (N_{E}) and note that this does not significantly influence');
Ht_q(2)=text(0.05,0.68,...
   'I_{C} and I_{nE} currents. WHY?');
Ht_q(3)=text(0.05,0.42,...
   'HINT QUESTIONS: Is the gradient of the electrons injected into the base influenced ');
Ht_q(4)=text(0.05,0.22,...
   'by N_{E}? Is your answer consistent with the I_{nE} equation in the BJT model?');
Ht_q(5)=text(0.05,0.8,...
   'If I_{C} is independent of the emitter doping, does the emitter doping influence');
Ht_q(6)=text(0.05,0.6,...
   'any important BJT characteristic? Does any important BJT parameter depend ');
Ht_q(7)=text(0.05,0.4,...
   'on the emitter doping?');
Ht_q(8)=text(0.05,0.88,...
   'Can changes in the gradient of electron concentration in the base explain significant');
Ht_q(9)=text(0.05,0.68,...
   'influences of the base doping and width on I_{nE} and I_{C} currents? Is I_{nE} limited');
Ht_q(10)=text(0.05,0.48,...
   'by the diffusion in the base or by the drift in the depletion layer? Is the answer ');
Ht_q(11)=text(0.05,0.28,...
   'related to the fact that I_{C} does not significantly depend on V_{CB} (V_{CE})?');
Ht_q(12)=text(0.05,0.8,...
   'Reduce the base doping and width to small values.');
Ht_q(13)=text(0.05,0.6,...
   'Why does this make I_{C} significantly dependent on V_{CE} (V_{CB})?');
Ht_q(14)=text(0.05,0.8,...
   'A decrease in V_{CB} below approximately -0.5V significantly reduces');
Ht_q(15)=text(0.05,0.6,...
   'the I_{C} current (saturation region). WHY?');


set(Ht_q,'HorizontalAlignment','left','VerticalAlignment','top','Visible','off');
else %put the figure at the top
   axes(Hqu);
end; %if ishandle(Hfquest)==0; 
   set(Hquest,'Name',['Question ' num2str(con_q_usd12-1)]);
   set(Hrbq,'Value',0);
   set(Hrbq(con_q_usd12),'Value',1);
   set(Ht_q,'Visible','off');
   switch con_q_usd12;
   case 2;
      set([Ht_q(1) Ht_q(2) Ht_q(3) Ht_q(4)],'Visible','on');
   case 3;
      set([Ht_q(5) Ht_q(6) Ht_q(7)],'Visible','on');
   case 4;
      set([Ht_q(8) Ht_q(9) Ht_q(10) Ht_q(11)],'Visible','on');
   case 5;   
      set([Ht_q(12) Ht_q(13)],'Visible','on');
   case 6;   
      set([Ht_q(14) Ht_q(15)],'Visible','on');
   end   
else; %if con_q_usd12>1
   if con_q_usd12==1;
   set(Hrbq,'Value',0);
   set(Hrbq(1),'Value',1);
   if ishandle(Hquest)==1;
      close(Hquest);
   end   
   end;   
end; %if con_q_usd12>1
con_q_usd12=0;
if ishandle(Hquest)==0;
   set(Hrbq,'Value',0);
   set(Hrbq(1),'Value',1);
end;

if ishandle(Hf)==0;
   if ishandle(Hf_model)==1; close(Hf_model); end;
   if ishandle(Hquest)==1; close(Hquest); end;
   return; 
end;
drawnow;
tic;
while toc==0.4;
   pause(0.2);
end   
pause(1);
if ishandle(Hf)==0;
   if ishandle(Hf_model)==1; close(Hf_model); end; 
   if ishandle(Hquest)==1; close(Hquest); end; 
   return; 
end;

Nel=ceil(wine*5);
if ((Vcbq>0)*(Vbeq>0.3)*(Nel>0));
%%% PREPARE THE ANIMATION LOOP%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Ecc=(5.5+Eg/2-phifc-(Vcbq+Vbeq));
%Ecb=5.5+Eg/2+phifb-Vbeq;
Evb=5.5-Eg/2+phifb-Vbeq;
%Evc=(5.5-Eg/2-phifc-(Vcbq+Vbeq));
xll=(We+Wb-Wcbb)*1e6;
xrr=(We+Wb+Wcbc)*1e6;
mm=(xrr-xll)/(Ecb-Ecc);
deltayel=(Ecb-Ecc)/(Nel+1);
ddeltayel=(Ecb-Ecc)/5; %steps in a fall
ipca=-8; %coefficient to slow-down/randomize the Ipc current
xel=0; yel=0; xtr=0; ytr=0;
%Create electrons and initialize the traces
   for i=1:Nel;
      yel(i)=Ecb-(i-0.5)*deltayel*(0.9+rand(1,1)*0.2);
      xel(i)=xll+mm*(Ecb-yel(i))+rand(1,1)*mm*0.5;
      xtr(1,i)=xel(i);
      ytr(1,i)=yel(i);
      set(Hp_tr(i),'Visible','on');
      set(Hp_el(i),'Visible','on');
   end   
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  %%%THE  ANIMATION LOOP
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  while (convbe_usd12+convcb_usd12+conr_usd12+con_close_usd12+con_n_usd12+...
        con_m_usd12+con_q_usd12)<1; 
   
   for ist=2:3; %5 steps in a trace
      for i=1:Nel;
      xel(i)=xel(i)+mm*ddeltayel;
      yel(i)=yel(i)-ddeltayel;
      xtr(ist,i)=xel(i);
      ytr(ist,i)=yel(i);
      if ytr(ist,i)<(Ecc+0.3); %El. into the coll. - put back
         yel(i)=Ecb+rand(1,1)*0.5;         
         xel(i)=xll-mm*0.5*(yel(i)-Ecb);
         for ii=1:ist;
            xtr(ii,i)=xel(i);
            ytr(ii,i)=yel(i);
         end   
      end
      set(Hp_tr(i),'XData',xtr(:,i),'YData',ytr(:,i));
      set(Hp_el(i),'XData',xel(i),'YData',yel(i));
      end
      
      if ipca>0; %Ipc hole to be shown;
         xho=xho-mm*ddeltayel;
         yho=yho+ddeltayel;
         xtrh(ist)=xho;
         ytrh(ist)=yho;
         set(Hp_trh,'XData',xtrh,'YData',ytrh,'Visible','on');
         set(Hp_ho,'XData',xho,'YData',yho,'Visible','on');
         if yho>(Evb-0.2);
            ipca=0;
            set(Hp_trh,'Visible','off');
            set(Hp_ho,'Visible','off');
         end   
      end
              if (ishandle(Hquest)==0)*(get(Hrbq(1),'Value')==0);
              set(Hrbq,'Value',0);
              set(Hrbq(1),'Value',1);
           end;
       if ishandle(Hf)==0;
          if ishandle(Hf_model)==1; close(Hf_model); end;
          if ishandle(Hquest)==1; close(Hquest); end;
          return; 
       end;
          drawnow;
      tic;
      while toc<0.2;
         pause(0.1);
      end
      if ishandle(Hf)==0;
         if ishandle(Hf_model)==1; close(Hf_model); end;
         if ishandle(Hquest)==1; close(Hquest); end;
         return; 
      end;

   end
   xtr=0; ytr=0;
   for i=1:Nel;
   xtr(1,i)=xel(i);
   ytr(1,i)=yel(i);
   end   
   
   if (ipca<1);
      ipca=ipca-1;
      if ipca<-10; %Ipc hole to be shown - initialize it here
      ipca=1;
      xtrh=0; ytrh=0; xho=0; yho=0;
      yho=Evc-rand(1,1)*0.3-0.1;         
      xho=xrr+mm*0.5*(Evc-yho);
      xtrh(1)=xho;
      ytrh(1)=yho;
      end
   else
      if ipca>0;
      ytrh=0; xtrh=0;
      xtrh(1)=xho; ytrh(1)=yho;
      end
   end
   
   if ishandle(Hf)==0;
   if ishandle(Hf_model)==1; close(Hf_model); end;
   if ishandle(Hquest)==1; close(Hquest); end;
   return; 
   end;
   tic;
  while toc<0.2;
     pause(0.1);
  end   
   if ishandle(Hf)==0;
   if ishandle(Hf_model)==1; close(Hf_model); end;
   if ishandle(Hquest)==1; close(Hquest); end;
   return; 
end;

end; %while of inner animation loop
if ishandle(Hf)==0;
   if ishandle(Hf_model)==1; close(Hf_model); end;
   if ishandle(Hquest)==1; close(Hquest); end;
   return; 
end;

set(Hp_tr,'Visible','off');
set(Hp_el,'Visible','off');
set(Hp_trh,'Visible','off');
set(Hp_ho,'Visible','off');

drawnow;

end; %if Vcbq>0;


end; %while con_close_us12<1; MAIN LOOP
if ishandle(Hf_model)==1;
   close(Hf_model);
end; 
close;
