function mfet_usd;
%3D and related 2D illustrations of the conduction energy band
%of a MOSFET. Allowing gate and drain voltages to be set at different values,
%it illustrates the depletion layers and current flow.
%(Sections 5.1.2 and 5.1.3 in Understanding Semiconductor Devices)

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


global con_vg_usd10; global con_vds_usd10; 
global con_r_usd10; global con_start_usd10; 
global con_close_usd10; global con_br_usd10;
global con_go_usd10;


%SET THE DOPING COLOR MAP
a=([1:64].')*0; mosfetcol=[a a a];
%the green part
for ii=1:25; %the dark green
   mosfetcol(ii,2)=0.5*(ii-1)^2.5/24^2.5+0.5;
end
for ii=26:32; %the light green
   mosfetcol(ii,1)=(ii-25)^0.7/7^0.7;;
   mosfetcol(ii,2)=1;
   mosfetcol(ii,3)=(ii-25)^0.7/7^0.7;
end   

%the red part
for ii=33:44;
   mosfetcol(ii,2)=(44-ii)^0.5/12^0.5;
   mosfetcol(ii,3)=(44-ii)^0.5/12^0.5;
   mosfetcol(ii,1)=1;
end 
for ii=45:64;
   mosfetcol(ii,1)=(64-ii)^0.5/19^0.5;
end


%CONSTANTS
Xmax=15; Ymax=7.5;
Xsource=1; Xdrain=14;
Ysource=2.0; Ydrain=2.0;
scal=10;
Eg=1.12;
ld=0.4; %Depletion layer weighting due to Vbi of S/D
lvd=3; %Depletion layer weighting due to Vd
lsurf=0.5; %Depletion-layer weighting due to surface potential 
phif_n=0.42; %Equilibrium N-type Fermi potential
phif_p=0.4; %Equilibrium P-type Fermi potential
ksl=45; %A larger ksl whitens the depletion layer in the gradient area
kT=0.6; %A smaller kT sharpens the colors/color map less uniform
phi2f=0.73; %Threshold surface potential
Vfb=0.5; %Flat-band voltage
spcoeff=1.0; %Coefficient in phis-vs-Vg equation
bright=0.2; %bright_max=1.5
Num_el=50; %Number of electrons in the source+that many in the drain
kTel=0.1; %Energy of electrons
num_x=Xmax*scal; %Number of lateral points
xGreedStep=Xmax/num_x;
LL=6; %Steps of every trace
num_y=Ymax*scal; %Number of vertical points
[x,y]=meshgrid(0:num_x,((num_y-[0:num_y]).^2/num_y));
x=x/scal; y=y/scal;

%Find the indeces for the random generation of the channel electrons
[amin,imin]=min(abs(Xsource+Ysource-x'));
ix_source=imin(1);
[amin,imin]=min(abs(Xdrain-Ydrain-x'));
ix_drain=imin(1);
[amin,imin]=min(abs(2*lsurf-y)); %MJax channel thickness is 2*lsurf
jy_ch=imin(1);
%The number of channel electrons, assuming same density as S/D
%Nch=round(2*lsurf*(Xdrain-Xsource-2*Ysource)*Num_el/(Ysource*(Xsource+Ysource)));
max_tr=16; %Maximum number of electron traces

%OPEN PARENT FIGURE
Hf=figure('NumberTitle','off','Name','Sections 5.1.2/5.1.3: MOSFET Principles',...
   'Color',[0.8 0.8 0.8]);
set(gcf,'InvertHardCopy','off');

Ha0=axes('Xlim',[0 15],'Ylim',[0 7.5],'Zlim',[-7.5 1.5],'Visible','off',...
   'Position',[0.01 0.15 .75 .825],...
   'CLim',[1 64],'View',[25 25]);
set(Hf,'Colormap',mosfetcol);


%READ THE INTRO TEXT              
introtext ={'UNDERSTANDING SEMICONDUCTOR DEVICES ';...
      'Sections 5.1.2/5.1.3: MOSFET Principles' ;...
      ' ';...
      '3D and related 2D presentations of the CONDUCTION ENERGY BAND (Ec) ';...
      'of a MOSFET are used to visualize the MOSFET operation. ';...
      '(NOTE: The valence band is omitted for clarity). Gate and drain';
      'voltages can be set at different values to see the energy and current effects.';...
      ' ';...
      'The red surface shows Ec in P-type areas, the green surfaces show Ec in';...
      'N-type areas, whereas the white areas indicate the depletion layers. ';...
      'Note that the gate voltage can create green N-type channel between ';...
      'the source and the drain.';...
      '          ';...
      'Changing the viewing angle can help to visualize and grasp the 3D shape ';...
      'of the Ec surface. It is very useful to relate the 3D view to ';... 
      'the commonly used 2D "x-y diagram". There is a demo rotation to enable you ';...
      'to make a direct link between these views. It is also very useful ';...
      'to view the 2D E-y and E-x diagrams.'};   

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

con_vg_usd10=0; con_vds_usd10=0; con_r_usd10=0; con_start_usd10=0; 
con_close_usd10=0; con_br_usd10=0;
%SETTING THE CONTROL PANEL              
Hframe=uicontrol(Hf,'Units','normalized','Position',[0.77 0.01 0.22 0.98],...
   'Style','frame');
%Vgs
Vg=0.5;
HframeVg=uicontrol(Hf,'Units','normalized','Position',[0.775 0.775 0.21 0.21],...
   'Style','frame','BackgroundColor','y');
Hrbg(1)=uicontrol(Hf,'Units','normalized','Position',[0.78 0.91 0.025 0.03],...
   'Style','radio','CallBack','con_vg_usd10=1;','Value',0,'BackgroundColor','y',...
   'HorizontalAlignment','center');
Hrbg(2)=uicontrol(Hf,'Units','normalized','Position',[0.805 0.91 0.025 0.03],...
   'Style','radio','CallBack','con_vg_usd10=2;','Value',1,'BackgroundColor','y',...
   'HorizontalAlignment','center');
Hrbg(3)=uicontrol(Hf,'Units','normalized','Position',[0.83 0.91 0.03 0.03],...
   'Style','radio','CallBack','con_vg_usd10=3;','Value',0,'BackgroundColor','y',...
   'HorizontalAlignment','center');
Hrbg(4)=uicontrol(Hf,'Units','normalized','Position',[0.86 0.91 0.03 0.03],...
   'Style','radio','CallBack','con_vg_usd10=4;','Value',0,'BackgroundColor','y',...
   'HorizontalAlignment','center');
Hrbg(5)=uicontrol(Hf,'Units','normalized','Position',[0.89 0.91 0.03 0.03],...
   'Style','radio','CallBack','con_vg_usd10=5;','Value',0,'BackgroundColor','y',...
   'HorizontalAlignment','center');
Hrbg(6)=uicontrol(Hf,'Units','normalized','Position',[0.92 0.91 0.03 0.03],...
   'Style','radio','CallBack','con_vg_usd10=6;','Value',0,'BackgroundColor','y',...
   'HorizontalAlignment','center');
Hrbg(7)=uicontrol(Hf,'Units','normalized','Position',[0.95 0.91 0.03 0.03],...
   'Style','radio','CallBack','con_vg_usd10=7;','Value',0,'BackgroundColor','y',...
   'HorizontalAlignment','center');
Htg(1)=uicontrol(Hf,'Units','normalized','Position',[0.79 0.885 0.055 0.03],...
   'Style','text','BackgroundColor','y','ForegroundColor','k','String','Vfb',...
   'HorizontalAlignment','center');
Htg(2)=uicontrol(Hf,'Units','normalized','Position',[0.835 0.885 0.02 0.03],...
   'Style','text','BackgroundColor','y','ForegroundColor','k','String','1',...
   'HorizontalAlignment','center');
Htg(3)=uicontrol(Hf,'Units','normalized','Position',[0.855 0.885 0.03 0.03],...
   'Style','text','BackgroundColor','y','ForegroundColor','k','String','2',...
   'HorizontalAlignment','center');
Htg(4)=uicontrol(Hf,'Units','normalized','Position',[0.885 0.885 0.03 0.03],...
   'Style','text','BackgroundColor','y','ForegroundColor','k','String','3',...
   'HorizontalAlignment','center');
Htg(5)=uicontrol(Hf,'Units','normalized','Position',[0.915 0.885 0.03 0.03],...
   'Style','text','BackgroundColor','y','ForegroundColor','k','String','4',...
   'HorizontalAlignment','center');
Htg(6)=uicontrol(Hf,'Units','normalized','Position',[0.945 0.885 0.03 0.03],...
   'Style','text','BackgroundColor','y','ForegroundColor','k','String','5',...
   'HorizontalAlignment','center');
Htg(7)=uicontrol(Hf,'Units','normalized','Position',[0.78 0.94 0.20 0.04],...
   'Style','text','String','Vg (V)','BackgroundColor','y','ForegroundColor','k');

%get(Hrbg(1))

%Vds
%Ht_vdu=uicontrol(Hf,'Units','normalized','Position',[0.78 0.85 0.20 0.01],...
%   'Style','text','BackgroundColor','y');
Hsl_vd=uicontrol(Hf,'Units','normalized','Position',[0.78 0.81 0.20 0.04],...
   'Style','slider','min',0,'max',7,'Value',0,'CallBack','con_vds_usd10=1;');
Ht_vdd=uicontrol(Hf,'Units','normalized','Position',[0.78 0.78 0.20 0.03],...
   'Style','text','String',['Vds=' num2str(round(get(Hsl_vd,'Value')*10)/10) ' V'],...
   'BackgroundColor','y','ForegroundColor','k');
Vds=round(get(Hsl_vd,'Value')*10)/10;
%3D
%The radiobutton
Hr(1)=uicontrol(Hf,'Units','normalized','Position',[0.78 0.71 0.20 0.03],...
   'Style','radio','String','3D (reset)','CallBack','con_r_usd10=1;','Value',1,...
   'ForegroundColor','b','FontSize',10,'FontWeight','bold');
%Azimut
Hsl_az=uicontrol(Hf,'Units','normalized','Position',[0.78 0.66 0.20 0.04],...
   'Style','slider','min',0,'max',90,'Value',25,'CallBack','con_r_usd10=5;');
Ht_azd=uicontrol(Hf,'Units','normalized','Position',[0.78 0.63 0.20 0.03],...
   'Style','text','String',['Azimuth=' num2str(round(get(Hsl_az,'Value'))) ' deg.']);
Az=round(get(Hsl_az,'Value'));
%Elevation
Hsl_el=uicontrol(Hf,'Units','normalized','Position',[0.78 0.58 0.20 0.04],...
   'Style','slider','min',-90,'max',45,'Value',25,'CallBack','con_r_usd10=6;');
Ht_eld=uicontrol(Hf,'Units','normalized','Position',[0.78 0.55 0.20 0.03],...
   'Style','text','String',['Elevation=' num2str(round(get(Hsl_el,'Value'))) ' deg.']);
El=round(get(Hsl_el,'Value'));
%Demo rotations
Ht_start=uicontrol(Hf,'Units','normalized','Position',[0.78 0.315 0.20 0.03],...
   'Style','text','String','DEMO  ROTATIONS');
Hb_start1=uicontrol(Hf,'Units','normalized','Position',[0.78 0.26 0.20 0.05],...
   'Style','pushbutton','String','3D->x-y','CallBack','con_start_usd10=1;');
Hb_start2=uicontrol(Hf,'Units','normalized','Position',[0.78 0.20 0.20 0.05],...
   'Style','pushbutton','String','3D->E-y (lateral)','CallBack','con_start_usd10=2;');
Hb_start3=uicontrol(Hf,'Units','normalized','Position',[0.78 0.14 0.20 0.05],...
   'Style','pushbutton','String','3D->E-x (vertical)','CallBack','con_start_usd10=3;');
%2D
Ht_2d=uicontrol(Hf,'Units','normalized','Position',[0.78 0.495 0.20 0.03],...
   'Style','text','String','2D');
%E-y (lat)
Hr(2)=uicontrol(Hf,'Units','normalized','Position',[0.78 0.46 0.20 0.03],...
   'Style','radio','String','x - y diagram','CallBack','con_r_usd10=2;','Value',0,...
   'ForegroundColor','b','FontSize',8,'FontWeight','bold');
Hr(3)=uicontrol(Hf,'Units','normalized','Position',[0.78 0.415 0.20 0.03],...
   'Style','radio','String','E - y (lateral)','CallBack','con_r_usd10=3;','Value',0);
Hr(4)=uicontrol(Hf,'Units','normalized','Position',[0.78 0.37 0.20 0.03],...
   'Style','radio','String','E - x (vertical)','CallBack','con_r_usd10=4;','Value',0);
%CLOSE
Hb_close=uicontrol(Hf,'Units','normalized','Position',[0.78 0.02 0.20 0.10],...
   'Style','pushbutton','String','CLOSE','CallBack','con_close_usd10=1;');
%Brightness
Ht_bru=uicontrol(Hf,'Units','normalized','Position',[0.02 0.465 0.20 0.03],...
   'Style','text','String','brightness','BackgroundColor',[0.8 0.8 0.8]);
Hsl_br=uicontrol(Hf,'Units','normalized','Position',[0.02 0.43 0.20 0.03],...
   'Style','slider','min',0.1,'max',0.7,'Value',0.3,'CallBack','con_br_usd10=1;');
bright=get(Hsl_br,'Value');



%SET Id-Vds GRAPH
%%%%%%%%%%%%%%%%%%%
Ha_idvd=axes('Xlim',[0 7],'Ylim',[0 5],'Position',[0.06 0.08 0.32 0.3],...
   'Color',[0.8 0.8 0.8],'FontAngle','italic');
hold on;
xlabel('DRAIN VOLTAGE (V)'); ylabel('DRAIN CURRENT (mA)');

%THE Id-Vds MODEL
Vdso=[0:0.1:7];
beta0=1.23;
theta=0.1;
sigma=0.1;
Fb=0.7;
Vt0=1.3;
for Vgo=0:5;
Vt=Vt0-sigma*Vdso;
beta=beta0./(1+theta*(Vgo-Vt));
Vdsat=(Vgo-Vt0)/(1+Fb);
Ido=0*Vdso.*(Vgo<Vt)+... %sub-threshold
    beta.*(Vgo-Vt-0.5*(1+Fb)*Vdso).*Vdso.*(Vdsat>Vdso).*(Vgo>Vt)+... %Triode
    0.5*beta/(1+Fb).*(Vgo-Vt).^2.*(Vdsat<=Vdso).*(Vgo>Vt); %Saturation
    
%Set the Id-Vd characteristics   
   if ishandle(Hf)==0;
      return; 
   end;
Hp_o(Vgo+1)=plot(Vdso,Ido,'-k','LineWidth',1);
set(Hp_o(Vgo+1),'Color',[0.5 0.5 0.5]);
end
Vgo=round(Vg); %to select a line for the current Vg
set(Hp_o(Vgo+1),'Color','y','LineWidth',2);

%Set the operating point
Vtq=Vt0-sigma*Vds;
betaq=beta0/(1+theta*(Vgo-Vtq));
Vdsat=(Vgo-Vt0)/(1+Fb);
if Vgo<Vtq;
   Idsoq=0;
else if Vds<Vdsat;
      Idsoq=betaq*(Vgo-Vtq-0.5*(1+Fb)*Vds)*Vds;
      else
         0.5*betaq/(1+Fb)*(Vgo-Vtq)^2;
      end   
end
   if ishandle(Hf)==0;
      return; 
   end;
Hp_oq=plot(Vds,Idsoq,'r.','MarkerSize',20);
set(Hp_oq,'Color','g');
hold off;

%SET Id-Vgs GRAPH
%%%%%%%%%%%%%%%%%%%
Ha_idvg=axes('Xlim',[0 5],'Ylim',[0 5],'Position',[0.5 0.08 0.2 0.3],...
   'Color',[0.8 0.8 0.8],'YAxisLocation','right','FontAngle','italic');
hold on;
xlabel('GATE VOLTAGE (V)'); ylabel('DRAIN CURRENT (mA)');

%THE Id-Vgs MODEL
Vgst=[0:0.1:5];
betat=beta0./(1+theta*(Vgst-Vtq));
Vdsatt=(Vgst-Vt0)/(1+Fb);
Idt=betat.*(Vgst-Vtq-0.5*(1+Fb)*Vds).*Vds.*(Vdsatt>Vds).*(Vgst>Vtq)+... %Triode
    0.5*betat/(1+Fb).*(Vgst-Vtq).^2.*(Vdsatt<=Vds).*(Vgst>Vtq); %Saturation
   if ishandle(Hf)==0;
      return; 
   end;
Hp_t=plot(Vgst,Idt,'-k','LineWidth',1);
set(Hp_t,'Color',[0.5 0.5 0.5]);

Hp_tq=plot(Vgo,Idsoq,'r.','MarkerSize',20);
set(Hp_tq,'Color','g');
hold off;


%SETTING THE STRUCTURE
%set r_s: Distance from the source, 0 inside
r_s0=sqrt((x-Xsource).^2+y.^2)-Ysource;
r_s=(x<=Xsource).*(y>Ysource).*(y-Ysource)...
   +(x>Xsource).*(r_s0>=0).*r_s0;  
Sour_xy=(r_s==0); %Defines the source region

%set r_d: Distance from the drain, 0 inside
r_d0=sqrt((x-Xdrain).^2+y.^2)-Ydrain;
r_d=(x>=Xdrain).*(y>Ydrain).*(y-Ydrain)...
   +(x<Xdrain).*(r_d0>=0).*r_d0;  
Drain_xy=(r_d==0); %Defines the drain region

Bulk_xy=1-Sour_xy-Drain_xy; %Defines the bulk

%Define Ec and Ef; 
Vg=0.5;
phis0=phi2f*(1-exp(-spcoeff*(Vg-Vfb)));
phis=phis0.*exp(-y/lsurf).*(x>=(Xsource+Ysource)).*(x<=(Xdrain+Ydrain));

Ef=x*0;

Ef=-Vds.*Drain_xy-Vds.*Bulk_xy.*(exp(-(Vds^0)*r_d/lvd)-...
      exp(-(Vds^0)*(Xdrain-Xsource-2.2*Ydrain)/lvd) ); %This term sets Ef at the source at 0 (a correction)
[gefx,gefy]=gradient(Ef);
gef=sqrt(gefx.^2+gefy.^2);
dphi_p=exp(-ksl*gef).*Bulk_xy+(1-Bulk_xy);

phin0=Eg/2-phif_n; 
Ec=Ef+Sour_xy.*phin0+Drain_xy.*phin0+...
   0.5.*(Eg/2+phif_p-phis).*Bulk_xy.*((1-exp(-(r_s-0.00)/ld))+...
   (1-exp(-(r_d-0.00)/ld)));

%Ec_f=tan(bright*(2*(dphi_p.*(Ec-Ef-Eg/2)+Eg/2)/Eg-1))/tan(bright)*32+32;
Ec_f=tan(bright*dphi_p.*(-exp((Ef-Ec)/kT)+exp((Ec-Eg-Ef)/kT) ))/tan(bright)*32+32;

axes(Ha0);
Hs_ec=surface(x,y,Ec,'EdgeColor','none','CData',Ec_f);

hold on;
%SET THE ELECTRONS
%The channel
xel_ch=0; yel_ch=0; Eel_ch=0; Nch=0;  
   if ishandle(Hf)==0;
      return; 
   end;
Hp3_ch=plot3(xel_ch,yel_ch,Eel_ch,'.k');
set(Hp3_ch,'EraseMode','xor','MarkerSize',10,'Visible','off');

%Source
    i=0;
    while i<Num_el;
       i=i+1;
       xel_s(i)=rand(1,1)*(Xsource+Ysource);
       yel_s(i)=rand(1,1)*Ysource;
       if xel_s(i)>Xsource;
       if sqrt((xel_s(i)-Xsource).^2+yel_s(i).^2)>=Ysource;
          i=i-1; %This is to elimiate the electrons outside the source
       end
       end
    end
    Eel_s=-(kTel*log(rand(1,Num_el)))+phin0;
%Drain    
    i=0;
    while i<Num_el;
       i=i+1;
       xel_d(i)=Xmax-rand(1,1)*(Xmax-Xdrain+Ydrain);
       yel_d(i)=rand(1,1)*Ydrain;
       if xel_d(i)<Xdrain;
       if sqrt((Xdrain-xel_d(i)).^2+yel_d(i).^2)>=Ydrain;
          i=i-1; %This is to elimiate the electrons outside the source
       end
       end
    end
    Efd=-Vds+...
    exp(-(Vds^0)*(Xdrain-Xsource-2.2*Ydrain)/lvd); %This term sets Ef at the source at 0 (a correction)
 Eel_d=Efd-(kTel*log(rand(1,Num_el)))+phin0;
    if ishandle(Hf)==0;
      return; 
    end;
Hp3_sd=plot3([xel_s xel_d],[yel_s yel_d],[Eel_s Eel_d],'.k');
set(Hp3_sd,'EraseMode','xor','MarkerSize',10);

%SET THE PLOTS FOR THE ELECTRON TRACES
xtr=[1:max_tr]*0;
ytr=[1:max_tr]*0;
Etr=[1:max_tr]*0;
for k=1:max_tr;
   Hpl3tr(k)=plot3(xtr(k),ytr(k),Etr(k),'-b','EraseMode','xor','Visible','off');
end   
hold off;

%THE TEXT
Ht_source=text('Position',[0.1 0.3 -0.7],'String','SOURCE','EraseMode','none');
Ht_drain=text('Position',[Xmax-0.1 0.3 -Vds-0.7],'String','DRAIN',...
   'HorizontalAlignment','right','EraseMode','none');
Ht_bulk1=text('Position',[Xmax/2 0.6*Ymax sign(El)*0.8+0.4],'String','BULK',...
   'EraseMode','none');


%***********************************
while con_close_usd10<1 %MAIN LOOP
   
%SET THE SLIDERS AND THE RADIO-BUTTONS 
   
%BRIGHTNESS
if con_br_usd10>0;
   bright=get(Hsl_br,'Value');
   Ec_f=tan(bright*dphi_p.*(-exp((Ef-Ec)/kT)+exp((Ec-Eg-Ef)/kT) ))/tan(bright)*32+32;
   if ishandle(Hf)==0;
      return; 
   end;
   set(Hs_ec,'ZData',Ec,'CData',Ec_f);
   con_br_usd10=0;
end   



%Vgs %%%%%%%%%%%%%%%%
if con_vg_usd10>0; 
   if ishandle(Hf)==0;
      return; 
   end;
   set(Hrbg,'Value',0);
   switch con_vg_usd10;
      case 1; Vg=0; set(Hrbg(1),'Value',1);
      case 2; Vg=0.5; set(Hrbg(2),'Value',1);
      case 3; Vg=1; set(Hrbg(3),'Value',1);
      case 4; Vg=2; set(Hrbg(4),'Value',1);
      case 5; Vg=3; set(Hrbg(5),'Value',1);
      case 6; Vg=4; set(Hrbg(6),'Value',1);
      case 7; Vg=5; set(Hrbg(7),'Value',1);
   end   
   drawnow;   
phis0=phi2f*(1-exp(-spcoeff*(Vg-Vfb)));

%Do the recalculations
phis=phis0.*exp(-y/lsurf).*(x>=(Xsource+Ysource)).*(x<=(Xdrain+Ydrain));

Ec=Ef+Sour_xy.*phin0+Drain_xy.*phin0+...
   0.5.*(Eg/2+phif_p-phis).*Bulk_xy.*((1-exp(-(r_s-0.00)/ld))+...
   (1-exp(-(r_d-0.00)/ld)));
Ec_f=tan(bright*dphi_p.*(-exp((Ef-Ec)/kT)+exp((Ec-Eg-Ef)/kT) ))/tan(bright)*32+32;
   if ishandle(Hf)==0;
      return; 
   end;
set(Hs_ec,'ZData',Ec,'CData',Ec_f);

%Regenerate S/D electrons
%Source electrons
    i=0;
    while i<Num_el;
       i=i+1;
       xel_s(i)=rand(1,1)*(Xsource+Ysource);
       yel_s(i)=rand(1,1)*Ysource;
       if xel_s(i)>Xsource;
       if sqrt((xel_s(i)-Xsource).^2+yel_s(i).^2)>=Ysource;
          i=i-1; %This is to elimiate the electrons outside the source
       end
       end
    end
    Eel_s(i)=-(kTel*log(rand(1,1)))+phin0;

%    Drain/source electrons
    i=0;
    while i<Num_el;
       i=i+1;
       xel_d(i)=Xmax-rand(1,1)*(Xmax-Xdrain+Ydrain);
       yel_d(i)=rand(1,1)*Ydrain;
       if xel_d(i)<Xdrain;
       if sqrt((Xdrain-xel_d(i)).^2+yel_d(i).^2)>=Ydrain;
          i=i-1; %This is to elimiate the electrons outside the source
       end
       end
    end
    Efd=-Vds+...
    exp(-(Vds^0)*(Xdrain-Xsource-2.2*Ydrain)/lvd); %This term sets Ef at the source at 0 (a correction)
    Eel_d=Efd-(kTel*log(rand(1,Num_el)))+phin0;
   if ishandle(Hf)==0;
      return; 
   end;
   
      if (Az>70)*(El>-30)*(El<10);

      set(Hp3_sd,'XData',xel_d,'YData',yel_d,'ZData',Eel_d);
      else 
      set(Hp3_sd,'XData',[xel_s xel_d],'YData',[yel_s yel_d],'ZData',[Eel_s Eel_d]);
      end   
    
    %The channel electrons
Nch=0;
switch round(Vg);
case 2; Nch=2;
case 3; Nch=5;
case 4; Nch=9;
case 5; Nch=16;   
end
ix=0; jy=0; Eel_chi=0; Eel_ch=0; spd=0; ii=0; xel_ch=0; yel_ch=0; spd=0;
   if ishandle(Hf)==0;
      return; 
   end;
set(Hpl3tr,'Visible','off');
if Nch>0;
  delta_ix=(ix_drain-ix_source)/Nch; 
  while ii<=(Nch-1);
   ii=ii+1;
   ix(ii)=round(delta_ix*(ii-0.5+0.8*(rand(1,1)-0.5))+ix_source); 
   %0.8*(rand(1,1)-0.5)) is a weak randomization of x pos.
   jy(ii)=round(rand(1,1)*(num_y+1-jy_ch)+jy_ch);
      
  Eel_chi(ii)=-(kTel*log(rand(1,1)))+Ef(num_y+1,ix(ii))+phin0;
  if Eel_chi(ii)<Ec(jy(ii),ix(ii));
      ii=ii-1; %Replace this electron as it is below Ec
   else
      if dphi_p(num_y+1,ix(ii))<0.2;
         ii=ii-1; Nch=Nch-1; %Randomly remove electrons from the slopes
      end
   end   
  end     

    %Convert indeces into x,y
       for i=1:Nch;
          xel_ch(i)=x(1,ix(i));
          yel_ch(i)=y(jy(i),1);
          Eel_ch(i)=Eel_chi(i);
          spd(i)=1-dphi_p(num_y+1,ix(i)); %Electron speed for the trace 
             if ishandle(Hf)==0;
                return; 
             end;
          set(Hpl3tr(i),'XData',xel_ch(i),'YData',yel_ch(i),'ZData',Eel_ch(i),'Visible','on');
       end   
    set(Hp3_ch,'XData',xel_ch,'YData',yel_ch,'ZData',Eel_ch,'Visible','on');
 else
    if ishandle(Hf)==0;
      return; 
    end;
   set(Hp3_ch,'Visible','off');
end


   %Set the initial positions of all Nch traces
xtr=0; ytr=0; Etr=0;   
for k=1:Nch;
   xtr(k)=xel_ch(k);
   ytr(k)=yel_ch(k);
   Etr(k)=Eel_chi(k);
end   

%Set the Id-Vd characteristic   
   if ishandle(Hf)==0;
      return; 
   end;
set([Hp_o(1) Hp_o(2) Hp_o(3) Hp_o(4) Hp_o(5) Hp_o(6)],'Color',[0.5 0.5 0.5],'LineWidth',1);
Vgo=round(Vg); %to select a line for the current Vg
set(Hp_o(Vgo+1),'Color','y','LineWidth',2);

%Set the output operating point
Vtq=Vt0-sigma*Vds;
betaq=beta0/(1+theta*(Vgo-Vtq));
Vdsat=(Vgo-Vt0)/(1+Fb);
if Vgo<Vtq;
   Idsoq=0;
else if Vds<Vdsat;
      Idsoq=betaq*(Vgo-Vtq-0.5*(1+Fb)*Vds)*Vds;
      else
         Idsoq=0.5*betaq/(1+Fb)*(Vgo-Vtq)^2;
      end   
end
   if ishandle(Hf)==0;
      return; 
   end;
set(Hp_oq,'Xdata',Vds,'YData',Idsoq);

% Set the Id-Vgs point

set(Hp_tq,'XData',Vgo,'YData',Idsoq);

con_vg_usd10=0;
end

%Vds %%%%%%%%%%%%%%%%%%
if con_vds_usd10>0; 
   if ishandle(Hf)==0;
      return; 
   end;
set(Ht_vdd,'String',['Vds=' num2str(round(get(Hsl_vd,'Value')*10)/10) ' V']);
Vds=round(get(Hsl_vd,'Value')*10)/10;

%Do the recalculations
Ef=-Vds.*Drain_xy-Vds.*Bulk_xy.*(exp(-(Vds^0)*r_d/lvd)-...
   exp(-(Vds^0)*(Xdrain-Xsource-2.2*Ydrain)/lvd) ); %This term sets Ef at the source at 0 (a correction)
[gefx,gefy]=gradient(Ef);
gef=sqrt(gefx.^2+gefy.^2);
dphi_p=exp(-ksl*gef).*Bulk_xy+(1-Bulk_xy);

Ec=Ef+Sour_xy.*phin0+Drain_xy.*phin0+...
   0.5.*(Eg/2+phif_p-phis).*Bulk_xy.*((1-exp(-(r_s-0.00)/ld))+...
   (1-exp(-(r_d-0.00)/ld)));
Ec_f=tan(bright*dphi_p.*(-exp((Ef-Ec)/kT)+exp((Ec-Eg-Ef)/kT) ))/tan(bright)*32+32;

   if ishandle(Hf)==0;
      return; 
   end;
set(Hs_ec,'ZData',Ec,'CData',Ec_f);
set(Ht_drain,'Position',[Xmax-0.1 0.3 -Vds-0.7]);

%Regenerate S/D electrons
%Source electrons
    i=0;
    while i<Num_el;
       i=i+1;
       xel_s(i)=rand(1,1)*(Xsource+Ysource);
       yel_s(i)=rand(1,1)*Ysource;
       if xel_s(i)>Xsource;
       if sqrt((xel_s(i)-Xsource).^2+yel_s(i).^2)>=Ysource;
          i=i-1; %This is to elimiate the electrons outside the source
       end
       end
    end
    Eel_s(i)=-(kTel*log(rand(1,1)))+phin0;

%    Drain/source electrons
    i=0;
    while i<Num_el;
       i=i+1;
       xel_d(i)=Xmax-rand(1,1)*(Xmax-Xdrain+Ydrain);
       yel_d(i)=rand(1,1)*Ydrain;
       if xel_d(i)<Xdrain;
       if sqrt((Xdrain-xel_d(i)).^2+yel_d(i).^2)>=Ydrain;
          i=i-1; %This is to elimiate the electrons outside the source
       end
       end
    end
    Efd=-Vds+...
    exp(-(Vds^0)*(Xdrain-Xsource-2.2*Ydrain)/lvd); %This term sets Ef at the source at 0 (a correction)
    Eel_d=Efd-(kTel*log(rand(1,Num_el)))+phin0;
   if ishandle(Hf)==0;
      return; 
   end;
   if (Az>70)*(El>-30)*(El<10);
      set(Hp3_sd,'XData',xel_d,'YData',yel_d,'ZData',Eel_d);
   else 
      set(Hp3_sd,'XData',[xel_s xel_d],'YData',[yel_s yel_d],'ZData',[Eel_s Eel_d]);
   end   
    
    %%The channel electrons
Nch=0;
switch round(Vg);
case 2; Nch=2;
case 3; Nch=5;
case 4; Nch=9;
case 5; Nch=16;   
end
ix=0; jy=0; Eel_chi=0; Eel_ch=0; spd=0; ii=0; xel_ch=0; yel_ch=0; spd=0;
   if ishandle(Hf)==0;
      return; 
   end;
set(Hpl3tr,'Visible','off');
if Nch>0;
  delta_ix=(ix_drain-ix_source)/Nch; 
  while ii<=(Nch-1);
   ii=ii+1;
   ix(ii)=round(delta_ix*(ii-0.5+0.8*(rand(1,1)-0.5))+ix_source); 
   %0.8*(rand(1,1)-0.5)) is a weak randomization of x pos.
   jy(ii)=round(rand(1,1)*(num_y+1-jy_ch)+jy_ch);
      
  Eel_chi(ii)=-(kTel*log(rand(1,1)))+Ef(num_y+1,ix(ii))+phin0;
  if Eel_chi(ii)<Ec(jy(ii),ix(ii));
      ii=ii-1; %Replace this electron as it is below Ec
   else
      if dphi_p(num_y+1,ix(ii))<0.2;
         ii=ii-1; Nch=Nch-1; %Randomly remove electrons from the slopes
      end
   end   
  end     

    %Convert indeces into x,y
       for i=1:Nch;
          xel_ch(i)=x(1,ix(i));
          yel_ch(i)=y(jy(i),1);
          Eel_ch(i)=Eel_chi(i);
          spd(i)=1-dphi_p(num_y+1,ix(i)); %Electron speed for the trace 
             if ishandle(Hf)==0;
                return; 
             end;
          set(Hpl3tr(i),'XData',xel_ch(i),'YData',yel_ch(i),'ZData',Eel_ch(i),'Visible','on');
       end   
    set(Hp3_ch,'XData',xel_ch,'YData',yel_ch,'ZData',Eel_ch,'Visible','on');
 else
       if ishandle(Hf)==0;
          return; 
       end;
   set(Hp3_ch,'Visible','off');
end
    
%Set the operating point
Vtq=Vt0-sigma*Vds;
betaq=beta0/(1+theta*(Vgo-Vtq));
Vdsat=(Vgo-Vt0)/(1+Fb);
if Vgo<Vtq;
   Idsoq=0;
else if Vds<Vdsat;
      Idsoq=betaq*(Vgo-Vtq-0.5*(1+Fb)*Vds)*Vds;
      else
         Idsoq=0.5*betaq/(1+Fb)*(Vgo-Vtq)^2;
      end   
end
   if ishandle(Hf)==0;
      return; 
   end;
set(Hp_oq,'Xdata',Vds,'YData',Idsoq);

%Set the Id-Vgs characteristic
betat=beta0./(1+theta*(Vgst-Vtq));
Vdsatt=(Vgst-Vt0)/(1+Fb);
Idt=betat.*(Vgst-Vtq-0.5*(1+Fb)*Vds).*Vds.*(Vdsatt>Vds).*(Vgst>Vtq)+... %Triode
    0.5*betat/(1+Fb).*(Vgst-Vtq).^2.*(Vdsatt<=Vds).*(Vgst>Vtq); %Saturation
   if ishandle(Hf)==0;
      return; 
   end;
set(Hp_t,'YData',Idt);

set(Hp_tq,'YData',Idsoq);



con_vds_usd10=0;
end

%AZIMUTH ANGLE
if con_r_usd10==5; 
   if ishandle(Hf)==0;
      return; 
   end;
   set(Ha0,'Position',[0.01 0.15 .75 .825])
Az=get(Hsl_az,'Value');
set(Ha0,'View',[Az El]);
set(Ht_azd,'String',['Azimuth=' num2str(round(get(Hsl_az,'Value'))) ' deg.']);
   for i=1:4;
      set(Hr(i),'Value',0);
   end   
   if (Az>72)*(El<3);
      set([Ht_source Ht_drain],'Visible','off');
   else set([Ht_source Ht_drain],'Visible','on');
   end
   if El<-20;
   set(Ha0,'Position',[0.01 0.15-((El+20)/70)*0.35 .75 .825+((El+20)/70)*0.35]);   
   end
   if (Az>70)*(El>-30)*(El<10);
      set(Hp3_sd,'XData',xel_d,'YData',yel_d,'ZData',Eel_d);
   else 
      set(Hp3_sd,'XData',[xel_s xel_d],'YData',[yel_s yel_d],'ZData',[Eel_s Eel_d]);
   end   
con_r_usd10=0;
end

%ELEVATION ANGLE
if con_r_usd10==6; 
   if ishandle(Hf)==0;
      return; 
   end;
   set(Ha0,'Position',[0.01 0.15 .75 .825])
El=get(Hsl_el,'Value');
set(Ha0,'View',[Az El]);
set(Ht_eld,'String',['Elevation=' num2str(round(get(Hsl_el,'Value'))) ' deg.']);
   for i=1:4;
   set(Hr(i),'Value',0);
   end   
   if (Az>72)*(El<3);
      set([Ht_source Ht_drain],'Visible','off');
   else set([Ht_source Ht_drain],'Visible','on');
   end
   if El<-20;
   set(Ha0,'Position',[0.01 0.15-((El+20)/70)*0.35 .75 .825+((El+20)/70)*0.35]);   
   end   
   if (El>-30)*(El<10)*(Az>70);
      set(Hp3_sd,'XData',xel_d,'YData',yel_d,'ZData',Eel_d);
   else 
      set(Hp3_sd,'XData',[xel_s xel_d],'YData',[yel_s yel_d],'ZData',[Eel_s Eel_d]);
   end   

con_r_usd10=0;
end

%3D (RESET)
if con_r_usd10==1;
   if ishandle(Hf)==0;
      return; 
   end;
  for i=1:4;
  if i==con_r_usd10;
     set(Hr(i),'Value',1);
  else set(Hr(i),'Value',0);
  end   
  end   
Az=25; El=25;
set([Ht_source Ht_drain],'Visible','on');  
set(Ha0,'Position',[0.01 0.15 .75 .825]);
set(Ha0,'View',[Az El]);
set(Hsl_az,'Value',Az);   
set(Hsl_el,'Value',El);   
set(Ht_azd,'String',['Azimuth=' num2str(round(get(Hsl_az,'Value'))) ' deg.']);
set(Ht_eld,'String',['Elevation=' num2str(round(get(Hsl_el,'Value'))) ' deg.']);
set(Hp3_sd,'XData',[xel_s xel_d],'YData',[yel_s yel_d],'ZData',[Eel_s Eel_d]);
con_r_usd10=0;
end
   
%2D X-Y
if con_r_usd10==2;
   if ishandle(Hf)==0;
      return; 
   end;
   for i=1:4;
      if i==con_r_usd10;   
         set(Hr(i),'Value',1);
      else set(Hr(i),'Value',0);
      end   
   end
Az=0; El=-90;   
set([Ht_source Ht_drain],'Visible','on');  
set(Ha0,'Position',[0.01 0.50 .75 .475]);  
set(Hsl_az,'Value',Az);   
set(Hsl_el,'Value',El);   
set(Ha0,'View',[Az El]);
set(Ht_azd,'String',['Azimuth=' num2str(round(get(Hsl_az,'Value'))) ' deg.']);
set(Ht_eld,'String',['Elevation=' num2str(round(get(Hsl_el,'Value'))) ' deg.']);
set(Hp3_sd,'XData',[xel_s xel_d],'YData',[yel_s yel_d],'ZData',[Eel_s Eel_d]);
con_r_usd10=0;
end   

%2D E-Y (LATERAL)
if con_r_usd10==3;
   if ishandle(Hf)==0;
      return; 
   end;
   for i=1:4;
      if i==con_r_usd10;   
         set(Hr(i),'Value',1);
      else set(Hr(i),'Value',0);
      end   
   end
Az=0; El=1;   
set([Ht_source Ht_drain],'Visible','on');  
set(Ha0,'Position',[0.01 0.15 .75 .825])   
set(Hsl_az,'Value',Az);   
set(Hsl_el,'Value',El);   
set(Ha0,'View',[Az El]);
set(Ht_azd,'String',['Azimuth=' num2str(round(get(Hsl_az,'Value'))) ' deg.']);
set(Ht_eld,'String',['Elevation=' num2str(round(get(Hsl_el,'Value'))) ' deg.']);
set(Hp3_sd,'XData',[xel_s xel_d],'YData',[yel_s yel_d],'ZData',[Eel_s Eel_d]);
con_r_usd10=0;
end   

%2D E-X (VERTICAL)
if con_r_usd10==4;
   if ishandle(Hf)==0;
      return; 
   end;
   for i=1:4;
      if i==con_r_usd10;   
         set(Hr(i),'Value',1);
      else set(Hr(i),'Value',0);
      end   
   end
Az=90; El=1;   
set([Ht_source Ht_drain],'Visible','off');
set(Ha0,'Position',[0.01 0.15 .75 .825]);   
set(Hsl_az,'Value',Az);   
set(Hsl_el,'Value',El);   
set(Ha0,'View',[Az El]);
set(Ht_azd,'String',['Azimuth=' num2str(round(get(Hsl_az,'Value'))) ' deg.']);
set(Ht_eld,'String',['Elevation=' num2str(round(get(Hsl_el,'Value'))) ' deg.']);
set(Hp3_sd,'XData',xel_d,'YData',yel_d,'ZData',Eel_d);
con_r_usd10=0;
end   

%3D-> X-Y ROTATION
if con_start_usd10==1;
   if ishandle(Hf)==0;
      return; 
   end;
   set([Ht_source Ht_drain],'Visible','on');
   set(Ha0,'Position',[0.01 0.15 .75 .825]);
   set(Hsl_vd,'Enable','off'); 
   set(Hrbg,'Enable','off');
   set(Htg,'Enable','off');
   set(Hr,'Value',0,'Enable','off');
   set([Hsl_az Hsl_el],'Enable','off');
   set([Hb_start1 Hb_start2 Hb_start3 Hb_close],'Enable','off');
   set(Hsl_br,'Enable','off');
   set([Ht_azd Ht_eld],'Visible','off')
   set(Hp3_sd,'XData',[xel_s xel_d],'YData',[yel_s yel_d],'ZData',[Eel_s Eel_d]);

Az=25; El=25;   
for i=25:-2:-19;
El=i; 
   if ishandle(Hf)==0;
      return; 
   end;
set(Ht_bulk1,'Position',[Xmax/2 0.6*Ymax sign(El)*0.8+0.4]);
set(Ha0,'View',[Az El]);
tic;
while toc==0.2;
   pause(0.1);
end   
   if ishandle(Hf)==0;
      return; 
   end;
drawnow;
end 
   if ishandle(Hf)==0;
      return; 
   end;
   pause(1);
   if ishandle(Hf)==0;
      return; 
   end;
for i=24:-2:0;
Az=i;
   if ishandle(Hf)==0;
      return; 
   end;
set(Ha0,'View',[Az El]);
tic;
while toc==0.2;
   pause(0.1);
end  
   if ishandle(Hf)==0;
      return; 
   end;
drawnow;
end 
   if ishandle(Hf)==0;
      return; 
   end;
   pause(1);
   if ishandle(Hf)==0;
      return; 
   end;
for i=-20:-5:-90;
   El=i;
   if ishandle(Hf)==0;
      return; 
   end;
set(Ha0,'View',[Az El],...
   'Position',[0.01 0.15-((i+20)/70)*0.35 .75 .825+((i+20)/70)*0.35]);
tic;
while toc==0.2;
   pause(0.1);
end   
   if ishandle(Hf)==0;
      return; 
   end;
drawnow;
end 
 
con_start_usd10=0;
   if ishandle(Hf)==0;
      return; 
   end;
   set(Hsl_vd,'Enable','on'); 
   set(Hr,'Enable','on');
   set(Hrbg,'Enable','on');
   set(Htg,'Enable','on');
   set([Hsl_az Hsl_el],'Enable','on');
   set([Hb_start1 Hb_start2 Hb_start3 Hb_close],'Enable','on');
   set(Hsl_br,'Enable','on');
   Az=0; El=-90;
   set(Hsl_az,'Value',Az);   
   set(Hsl_el,'Value',El);   
   set(Ht_azd,'Visible','on','String',...
      ['Azimuth=' num2str(round(get(Hsl_az,'Value'))) ' deg.']);
   set(Ht_eld,'Visible','on','String',...
      ['Elevation=' num2str(round(get(Hsl_el,'Value'))) ' deg.']);
end

%3D-> E-Y (LATERAL) ROTATION
if con_start_usd10==2;
   if ishandle(Hf)==0;
      return; 
   end;
   set([Ht_source Ht_drain],'Visible','on')
   set(Ha0,'Position',[0.01 0.15 .75 .825])
   set(Hsl_vd,'Enable','off'); 
   set(Hrbg,'Enable','off');
   set(Htg,'Enable','off');
   set([Hr(1) Hr(2) Hr(3) Hr(4)],'Value',0,'Enable','off');
   set([Hsl_az Hsl_el],'Enable','off');
   set([Hb_start1 Hb_start2 Hb_start3 Hb_close],'Enable','off');
   set([Ht_azd Ht_eld],'Visible','off')
   set(Hsl_br,'Enable','off');
   set(Hp3_sd,'XData',[xel_s xel_d],'YData',[yel_s yel_d],'ZData',[Eel_s Eel_d]);

Az=25; El=25;   
for i=25:-2:1;
   El=i; Az=i; 
   if ishandle(Hf)==0;
      return; 
   end;
set(Ha0,'View',[Az El]);

tic;
while toc==0.2;
   pause(0.1);
end 
   if ishandle(Hf)==0;
      return; 
   end;
drawnow;
end 
 
con_start_usd10=0;
   if ishandle(Hf)==0;
      return; 
   end;
   set(Hsl_vd,'Enable','on');
   set(Hrbg,'Enable','on');
   set(Htg,'Enable','on');
   set([Hr(1) Hr(2) Hr(3) Hr(4)],'Enable','on');
   set([Hsl_az Hsl_el],'Enable','on');
   set([Hb_start1 Hb_start2 Hb_start3 Hb_close],'Enable','on');
   set(Hsl_br,'Enable','on');
   Az=0; El=1;
   set(Ha0,'View',[Az El]);
   set(Hsl_az,'Value',Az);   
   set(Hsl_el,'Value',El);   
   set(Ht_azd,'Visible','on','String',...
      ['Azimuth=' num2str(round(get(Hsl_az,'Value'))) ' deg.']);
   set(Ht_eld,'Visible','on','String',...
      ['Elevation=' num2str(round(get(Hsl_el,'Value'))) ' deg.']);
end

%3D-> E-X (VERTICAL) ROTATION
if con_start_usd10==3;
   if ishandle(Hf)==0;
      return; 
   end;
   set([Ht_source Ht_drain],'Visible','on')
   set(Ha0,'Position',[0.01 0.15 .75 .825])
   set(Hsl_vd,'Enable','off'); 
   set(Hrbg,'Enable','off');
   set(Htg,'Enable','off');
   set([Hr(1) Hr(2) Hr(3) Hr(4)],'Value',0,'Enable','off');
   set([Hsl_az Hsl_el],'Enable','off');
   set([Hb_start1 Hb_start2 Hb_start3 Hb_close],'Enable','off');
   set([Ht_azd Ht_eld],'Visible','off')
   set(Hsl_br,'Enable','off');
   set(Hp3_sd,'XData',xel_d,'YData',yel_d,'ZData',Eel_d);

Az=25; El=25;   
for i=25:-2:1;
El=i; 
Az=(90-25)*(26-i)/25+25;   
   if ishandle(Hf)==0;
      return; 
   end;
set(Ha0,'View',[Az El]);
   if (Az>72)*(El<3);
      set([Ht_source Ht_drain],'Visible','off');
   else set([Ht_source Ht_drain],'Visible','on');
   end

tic;
while toc==0.2;
   pause(0.1);
end   
   if ishandle(Hf)==0;
      return; 
   end;
drawnow;
end 
 
con_start_usd10=0;
   if ishandle(Hf)==0;
      return; 
   end;
   set(Hsl_vd,'Enable','on'); 
   set(Hrbg,'Enable','on');
   set(Htg,'Enable','on');
   set([Hr(1) Hr(2) Hr(3) Hr(4)],'Enable','on');
   set([Hsl_az Hsl_el],'Enable','on');
   set([Hb_start1 Hb_start2 Hb_start3 Hb_close],'Enable','on');
   set(Hsl_br,'Enable','on');
   Az=90; El=1;
   set(Ha0,'View',[Az El]);
   set(Hsl_az,'Value',Az);   
   set(Hsl_el,'Value',El);   
   set(Ht_azd,'Visible','on','String',...
      ['Azimuth=' num2str(round(get(Hsl_az,'Value'))) ' deg.']);
   set(Ht_eld,'Visible','on','String',...
      ['Elevation=' num2str(round(get(Hsl_el,'Value'))) ' deg.']);
end


%THE OUTER ANIMATION LOOP
   if ishandle(Hf)==0;
      return; 
   end;
set(Ht_bulk1,'Position',[Xmax/2 0.6*Ymax sign(El)*0.8+0.4]);     
while (con_vg_usd10+con_vds_usd10+con_r_usd10+con_start_usd10+con_close_usd10+con_br_usd10)<1;
       
if Nch>0;
%Set the initial positions of all Nch traces
xtr=0; ytr=0; Etr=0;   
for k=1:Nch;
   xtr(k)=xel_ch(k);
   ytr(k)=yel_ch(k);
   Etr(k)=Eel_ch(k);
end   

%THE INNER LOOP (to create the traces)
l=1; 
while (((con_vg_usd10+con_vds_usd10+con_r_usd10+con_start_usd10+con_close_usd10+con_br_usd10)<1)+...
      l<=LL);
   l=l+1;
   
   for k=1:Nch;
      if Vds>0;
         deltax=spd(k)*1;
      else
         deltax=0.2*(rand(1,1)-0.5);
      end   
      xtr(l,k)=xtr(l-1,k)+deltax;
      ytr(l,k)=ytr(1,k);
      ixstt(k)=ix(k)+round((xtr(l,k)-xtr(1,k))/xGreedStep);
      if xtr(l,k)>=0.97*(Xdrain-Ydrain);
            ix(k)=round(delta_ix*(0.2+0.2*(rand(1,1)-0.5))+ix_source); 
            spd(k)=1-dphi_p(num_y+1,ix(k)); %Electron speed for the trace 
            %0.2*(rand(1,1)-0.5)) is a weak randomization of x pos.
            ixstt(k)=ix(k);
            xel_ch(k)=x(1,ix(k));

            jy(k)=round(rand(1,1)*(num_y+1-jy_ch)+jy_ch);
            Eel_ch(k)=-(kTel*log(rand(1,1)))+Ef(num_y+1,ix(k))+phin0;
            while Eel_ch(k)<Ec(jy(k),ix(k));
               %Find an electron that is above Ec
            jy(k)=round(rand(1,1)*(num_y+1-jy_ch)+jy_ch);
            Eel_ch(k)=-(kTel*log(rand(1,1)))+Ef(num_y+1,ix(k))+phin0;
            end   
            %spd(k)=1-dphi_p(num_y+1,ix(k));
            for li=1:l;
               %set all the previous trace points to the starting one
               xtr(li,k)=xel_ch(k);
               ytr(li,k)=yel_ch(k);
               Etr(li,k)=Eel_ch(k);
            end   
      else   
           Etr(l,k)=Etr(1,k)-Ef(num_y+1,ix(k))+Ef(num_y+1,ixstt(k));
           Eel_ch(k)=Etr(l,k);
           xel_ch(k)=xtr(l,k);
           spd(k)=1-dphi_p(num_y+1,ixstt(k)); %Electron speed for the trace 
      end
         if ishandle(Hf)==0;
            return; 
         end;
        set(Hpl3tr(k),'XData',xtr(:,k),'YData',ytr(:,k),'ZData',Etr(:,k));
    end
    if ishandle(Hf)==0;
      return; 
    end;
    set(Hp3_ch,'XData',xel_ch,'YData',yel_ch,'ZData',Eel_ch);
    
    tic;
    while toc<0.2;
       pause(0.1);
    end   
end; %while of the inner (trace) animation loop

ix=ixstt;


else 
   if ishandle(Hf)==0;
      return; 
   end;
   pause(1); %To delay for the time the tracing would take
   if ishandle(Hf)==0;
      return; 
   end;

end; %if Nch>0
end; %while of outer animation loop
   
end; %while
close