function vdr_usd;
%Animation of electron trajectory in a semiconductor crystal,
%illustrating the concept of drift velocity.
%(Section 1.4.1 in Understanding Semiconductor Devices)

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


global con_m_usd4; global con_t_usd4; global con_d_usd4; global con_start_usd4; 
global con_stop_usd4; global con_close_usd4; global con_e_usd4; global con_go_usd4;
con_go_usd4=0;

%OPEN PARENT FIGURE
Hf=figure('NumberTitle','off','Name','Section 1.4.1: Drift Velocity',...
   'Color','k');
set(gcf,'InvertHardCopy','off');

Ha0=axes('Xlim',[0 10],'Ylim',[0 10],'Visible','off','Position',[0 0 1 1]);
Hp_bckgr=patch([0 10 10 0],[0 0 10 10],[0.8 0.8 0.8]); 
%The Figure color is black to produce a black trace (trajectory). 
%It is covered by  Hp_bckgr to convert it into normal color otherwise

%READ THE INTRO TEXT              
introtext ={'UNDERSTANDING SEMICONDUCTOR DEVICES ';...
      'Section 1.4.1: Drift Velocity' ;...
      ' ';...
      'This animation shows the electron trajectory in a semiconductor crystal.';...
      'Set the electric field to different values to see how it influences ';...
      'the electron trajectory, resulting in effective drift of the electron. '};
Ht_intro=text(0.1,0.6,char(introtext),'Units','normalized');
Hp_cont=uicontrol(Hf,'Units','normalized','Position',[0.1 0.2 0.7 0.1],...
   'Style','pushbutton','String','CLICK HERE TO CONTINUE!','CallBack','con_go_usd4=1;');
con_go_usd4=0;
while con_go_usd4<1; 
      if ishandle(Hf)==0;
         return; 
      end
   pause(0.2);
end
con_go_usd4=0;
      if ishandle(Hf)==0;
         return; 
      end
set(Ht_intro,'Visible','off');
set(Hp_cont,'Visible','off','Enable','off');
drawnow;

con_m_usd4=0; con_t_usd4=0; con_d_usd4=0; con_start_usd4=0; 
con_stop_usd4=0; con_close_usd4=0; con_e_usd4=0;
%SETTING THE CONTROL PANEL              
Hframe=uicontrol(Hf,'Units','normalized','Position',[0.77 0.01 0.22 0.98],...
   'Style','frame');

%Electric field
Ht_efu=uicontrol(Hf,'Units','normalized','Position',[0.78 0.90 0.20 0.05],...
   'Style','text','String','ELECTRIC FIELD','BackgroundColor',[.7 0 0],...
   'ForegroundColor','w');
Hsl_ef=uicontrol(Hf,'Units','normalized','Position',[0.78 0.85 0.20 0.05],...
   'Style','slider','min',0,'max',4,'Value',0,'CallBack','con_e_usd4=1;');
Ht_efd=uicontrol(Hf,'Units','normalized','Position',[0.78 0.80 0.20 0.05],...
   'Style','text','String',[num2str(round(10*(get(Hsl_ef,'Value')))/10) ' V/um'],...
   'BackgroundColor',[0.7 0 0],'ForegroundColor','w');
Ef=round(10*(get(Hsl_ef,'Value')))/10;

%m_eff slider
Ht_mu=uicontrol(Hf,'Units','normalized','Position',[0.78 0.75 0.20 0.03],...
   'Style','text','String','MATERIAL (EFFECTIVE MASS)');
Hsl_m=uicontrol(Hf,'Units','normalized','Position',[0.78 0.70 0.20 0.05],...
   'Style','slider','min',0.5,'max',2.0,'Value',0.5,'CallBack','con_m_usd4=1;');
Ht_md=uicontrol(Hf,'Units','normalized','Position',[0.78 0.65 0.20 0.05],...
   'Style','text','String',['meff= ' num2str(round(10*get(Hsl_m,'Value'))/10) '*m0']);
meff=round(10*get(Hsl_m,'Value'))/10;

%Pushbuttons
Hb_close=uicontrol(Hf,'Units','normalized','Position',[0.78 0.02 0.20 0.07],...
   'Style','pushbutton','String','CLOSE','CallBack','con_close_usd4=1;');
Hb_stop=uicontrol(Hf,'Units','normalized','Position',[0.78 0.11 0.20 0.07],...
   'Style','pushbutton','String','STOP','CallBack','con_stop_usd4=1;','Enable','off');
Hb_start=uicontrol(Hf,'Units','normalized','Position',[0.78 0.20 0.20 0.07],...
   'Style','pushbutton','String','START','CallBack','con_start_usd4=1;');

%CONSTANTS
r=0.2; meff=0.5; mu=1/(20*meff); vsat=0.1;
alpha_e=pi/8; 
dr0=0.2/2; % for 1/meff=2
alpha_i=(rand-0.5)*pi/3+3*pi/4; %Initial incident angle
Xmax=24; Ymax=16;

%SET THE 2-D GRAPH
Haan=axes('Xlim',[0 Xmax],'Ylim',[0 Ymax],'Visible','off',...
   'Position',[0.05 0.5 0.675 0.45]);
xe1=[-r:0.1*r:r];
ye1=sqrt(r*r-xe1.*xe1);
xe=[xe1 xe1]; %elementary phonon
ye=[ye1 -ye1]; %elementary phonon

hold on;
Hp_bckgr=patch([0 Xmax Xmax 0],[0 0 Ymax Ymax],'y','EdgeColor','y','EraseMode','none');
for i=1:25;
for j=1:17;
   Hp(i,j)=patch((i-1+xe),(j-1+ye),[1 0.75 0],'EdgeColor',[1 0.75 0],'EraseMode','none',...
      'Visible','off');
end   
end

%SET the electron
xel=6.6; yel=6.6;
Hp_el=plot(xel,yel,'r.','MarkerSize',10,'EraseMode','xor','Visible','off');

%SET the electron trace
xtr=xel; ytr=yel;
Hp_tr=plot(xtr,ytr,'k-','LineWidth',1,'EraseMode','none','Visible','off');

hold off;

%SET Vdr-E graph
efx=[0:0.1:5];
vdr=mu*efx./(1+mu*efx/vsat);
vdrq=mu*Ef./(1+mu*Ef/vsat);
Havdr=axes('Xlim',[0 4],'Ylim',[0 0.1],'Position',[0.1 0.1 0.4 0.3],...
   'FontAngle','italic');
hold on;
      if ishandle(Hf)==0;
         return; 
      end
Hp_vdr=plot(efx,vdr,'b-','LineWidth',2);
xlabel('ELECTRIC FIELD  (V/\mum)'); ylabel('DRIFT VELOCITY (\mum/ps)');
set(Havdr,'Color',[0.85 0.85 0.85]);
Hp_vdrq=plot(Ef,vdrq,'r.','MarkerSize',15);
hold off;

for i=1:25;
   for j=1:17;
      if ishandle(Hf)==0;
         return; 
      end
   set(Hp(i,j),'Visible','on');
end   
end
%***********************************
while con_close_usd4<1 %MAIN LOOP
   
%SET THE SLIDER VALUES and RENEW 
  
if con_e_usd4>0; %Electric field
      if ishandle(Hf)==0;
         return; 
      end
set(Ht_efd,'String',[num2str(round(10*(get(Hsl_ef,'Value')))/10) ' V/um']);
Ef=round(10*(get(Hsl_ef,'Value')))/10;
vdrq=mu*Ef./(1+mu*Ef/vsat);
set(Hp_vdrq,'XData',Ef,'YData',vdrq);
end
con_e_usd4=0; 

if con_m_usd4>0; %Effective mass
      if ishandle(Hf)==0;
         return; 
      end
set(Ht_md,'String',['meff= ' num2str(round(10*get(Hsl_m,'Value'))/10) '*m0']);
meff=round(10*get(Hsl_m,'Value'))/10;
mu=1/(20*meff);
vdr=mu*efx./(1+mu*efx/vsat);
vdrq=mu*Ef./(1+mu*Ef/vsat);
set(Hp_vdr,'YData',vdr);
set(Hp_vdrq,'XData',Ef,'YData',vdrq);
end;
con_m_usd4=0;
      if ishandle(Hf)==0;
         return; 
      end
pause(1);
      if ishandle(Hf)==0;
         return; 
      end
if con_start_usd4>0;
   con_start_usd4=0;
   if ishandle(Hf)==0;
      return; 
   end
   %Renew the background
   set([Hp_el Hp_tr],'Visible','off');
   set(Hp_bckgr,'Xdata',[0 Xmax Xmax 0]); %To renew the background
   for i=1:25; % To renew the scattering centers
      for j=1:17;
         if ishandle(Hf)==0;
            return; 
         end
     set(Hp(i,j),'Xdata',(i-1+xe));
   end   
   end
   % renew the electron position and starting angle
   dir=1; xel=6.6; yel=6.6; xtr=xel; ytr=yel;
   alpha_i=(rand-0.5)*pi/3+3*pi/4; beta=alpha_i;%Initial incident angle
   
   %set the buttons
   if ishandle(Hf)==0;
      return; 
   end
   set([Hb_start Hb_close Hsl_ef Hsl_m],'Enable','off');
   set(Hb_stop,'Enable','on');
   set(Hp_el,'Xdata',xel,'Ydata',yel','Visible','on');
   set(Hp_tr,'Xdata',[xel xtr],'Ydata',[yel ytr],'Visible','on');
   
nst=0;   
dr=dr0/meff; 
accel_y=0;
while con_stop_usd4<1; %The animation loop
   
nst=nst+1; %To accelerate the particle
   
%remember the old position, and project the new position
xel_old=xel; yel_old=yel;
dx=dir*dr*cos(alpha_i); 
dy=dir*dr*sin(alpha_i);
accel=0.2*dr*Ef;
xel=xel+dx+accel*nst; 
yel=yel+dy+accel_y*nst;  
   
   %find the step size
   stp=sqrt((xel-xel_old)^2+(yel-yel_old)^2);
   
   if stp>=r %Need to cut the step into substeps
      beta=atan((yel-yel_old)/(xel-xel_old));
      alpha_i;
      dxels=0.9*sign(xel-xel_old)*r*abs(cos(beta));
      dyels=0.9*sign(yel-yel_old)*r*abs(sin(beta));
   else
      dxels=xel-xel_old; dyels=yel-yel_old;
   end   
        
   xels=xel_old; yels=yel_old; csec=0; stps=0;
   while ((stps<stp)*(csec==0));
      xels=xels+dxels; yels=yels+dyels;   
      stps=sqrt((xels-xel_old)^2+(yels-yel_old)^2);
      if stps>stp;
         xels=xel; yels=yel;
         stps=stp;
      end   
      
   %Check if the electron is inside a circle
   xdist=xels-round(xels); ydist=yels-round(yels);
   dist=sqrt(xdist^2+ydist^2);
   if (dist<=r); %Check if the electron is inside a circle
      csec=1; %There has been a cross section
   %put it back at the cross-section with the circle
   mm=tan(beta);
   aa=1+mm*mm;
   bb=2*mm*(ydist-mm*xdist);
   cc=mm*mm*xdist*xdist-2*mm*xdist*ydist+ydist*ydist-r*r;
   sq=sqrt(bb*bb-4*aa*cc);
   x1=(-bb+sq)/(2*aa); %x1 and x2 are the two cross-secting x co-ord.
   y1=mm*(x1-xdist)+ydist;
   x2=(-bb-sq)/(2*aa);
   y2=mm*(x2-xdist)+ydist;
   x1o=x1+round(xels);
   y1o=y1+round(yels);
   x2o=x2+round(xels);
   y2o=y2+round(yels);
   if ((x1o-xel_old)^2+(y1o-yel_old)^2)<((x2o-xel_old)^2+(y2o-yel_old)^2);
      xels=x1o; yels=y1o;
      xdist=x1; ydist=y1;
   else
      xels=x2o; yels=y2o;
      xdist=x2; ydist=y2;
   end
   %ydist_new=mm*(xdist_new-xdist)+ydist;
   if csec==1;   
   accel_y=(2*rand-1)*accel*0.5;
   nst=0;
   dir=-dir; 
   alpha_i=2*atan(ydist/xdist)-alpha_i;
   beta=alpha_i;
   end
   end   
   end % while
   xel=xels;
   yel=yels;
   
if xel>Xmax; 
   if Ef==0;
   alpha_i=pi-alpha_i;
   xel=Xmax;
   else
   con_stop_usd4=1;
   end   
end   
if (xel<0);
   alpha_i=pi-alpha_i;
   xel=0;
end
if yel>Ymax;
   alpha_i=pi-alpha_i;
   dir=-dir;
   yel=Ymax;
end
if yel<0;
   alpha_i=pi-alpha_i;
   dir=-dir;
   yel=0;
end
      if ishandle(Hf)==0;
         return; 
      end
set(Hp_el,'Xdata',xel,'YData',yel);
set(Hp_tr,'Xdata',[xel xtr],'Ydata',[yel ytr]);
xtr=xel; ytr=yel;

tic;
      if ishandle(Hf)==0;
         return; 
      end
while toc<0.05;
   pause(0.01);
end   
if ishandle(Hf)==0;
   return; 
end
end; end % if/while con_stop_usd4<1

con_stop_usd4=0;
      if ishandle(Hf)==0;
         return; 
      end
set([Hb_start Hb_close Hsl_ef Hsl_m],'Enable','on');
set(Hb_stop,'Enable','off');

end; % while of main loop (con_close_usd4=1)
close
