function[]=pn_usd(introtext);
%Energy-band display of Stored-Charge Related Transient Effects
%(Section 3.2.2 in Understanding Semiconductor Devices)

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

global debcon1_usd9; global debcon3_usd9; global debcon4_usd9; global debcon5_usd9;

%OPEN PARENT FIGURE
Hf =figure('NumberTitle','off','Name','Section 3.2.2: Stored-Charge Related Transient Effects',...
   'Color',[0.8 0.8 0.8]);
set(gcf,'InvertHardCopy','off');

Ha =axes('Xlim',[0 10],'Ylim',[0 10],'Visible','off',...
                 'Position',[.02 0.0 0.7 1.0]);
Hframe =uicontrol(Hf ,'Units','normalized','Position',[0.8 0.19 0.18 0.76],'Style','frame');

%INTRO TEXT          
Ht =text(0.1,0.52,char(introtext),'Units','normalized');

%SETTING THE BUTTONS
%the slider
Hb1 =uicontrol(Hf ,'Units','normalized','Position',[0.81 0.83 0.16 0.05],...
   'Style','slider','min',-1.5,'max',0.71,'Value',0,'CallBack','debcon1_usd9=1;');
Hb1l =uicontrol(Hf ,'Units','normalized','Position',[0.81 0.90 0.05 0.03],...
   'Style','text','String','-1.5 V');
Hb1r =uicontrol(Hf ,'Units','normalized','Position',[0.92 0.90 0.05 0.03],...
   'Style','text','String','0.7 V');
Hb1c =uicontrol(Hf ,'Units','normalized','Position',[0.86 0.90 0.06 0.03],...
   'Style','text','String',num2str(round(get(Hb1 ,'Value')*10)*0.1));
%The pushbuttons
Hb3 =uicontrol(Hf ,'Units','normalized','Position',[0.81 0.52 0.16 0.1],...
   'Style','pushbutton','String','START','CallBack','debcon3_usd9=1;');
Hbs =uicontrol(Hf ,'Units','normalized','Position',[0.81 0.37 0.16 0.1],...
   'Style','pushbutton','String','STOP','CallBack','debcon5_usd9=1;','Enable','off');
Hb4 =uicontrol(Hf ,'Units','normalized','Position',[0.81 0.22 0.16 0.1],...
   'Style','pushbutton','String','CLOSE','CallBack','debcon4_usd9=1;');


%CONSTANTS
Nel_n=750; Nel_m=7; Nho_p=750; Nho_m=7;
Nel_p=Nel_m; Nho_n=Nho_m; %Initial stored-charge is zero: Nel_p-Nel_m
in_depl_el=[1:1:(Nel_n+Nel_p)]*0; in_depl_ho=[1:1:(Nho_n+Nho_p)]*0;
y_Ec=[1:1:(Nel_n+Nel_p)]*0; y_Ev=[1:1:(Nho_n+Nho_p)]*0;
Eg=1.2; Ew=1.9; %Vmax+Eg+2*Ew should be 10
N_tr_st=15; lsc=0.75; 
Rrec=0.20; %Normalized recombination rate
taus=5; Idis=-0.4;

%SETTING THE FIGURE HANDLES
%Set the energy band handles
xn=4.5; xp=5.5; %Depl. layer edges
Efn=7; Efp=7; %Sets y position of the bands; Efp is fixed
Ecn=Efn+0.2; Evn=Efn-1.8; Ecp=Efp+1.8; Evp=Efp-0.2;
hold on;
x=[0 xn xn 0]; y=[Ecn Ecn Evn Evn]; %N-type Eg (fixed)
Hph_n=patch(x,y,'c','Visible','off','EdgeColor','c'); 
x=[xn xp xp xn]; y=[Ecn Ecp Evp Evn]; %Depl. layer Ec
Hph_d=patch(x,y,'c','Visible','off','EdgeColor','c'); %Ec
x=[xp 10 10 xp]; y=[Ecp Ecp Evp Evp]; %P-type Eg(0)
Hph_p=patch(x,y,'c','Visible','off','EdgeColor','c'); 
Hp_efn=plot([0 xn],[Efn Efn],'m-.','Visible','off');
Hp_efp=plot([xp 10],[Efp Efp],'m-.','Visible','off');

%Set the i-v and i-t plots
orx=3.5; ory=0.5;
x=[0 0 0.25 0.25 0.5 0.5 0.75 0.75 1 1 0 0.0 1.0 1 0]*6+orx;
y=[0 1 1.00 0.00 0.0 1.0 1.00 0.00 0 1 1 0.5 0.5 0 0]*1.6+ory;
Hp_it=plot(x,y,'b-','Visible','off'); %i-t axes
x_time=[1:1]*0+orx; y_time=[1:1]*0+ory+1;
Hp_it_tr=plot(x_time,y_time,'r-','Visible','off');
x=[0.8 0.8 2.3 2.3 3 3 0.8 0.8 3.0 3 0.8];
y=[0.0 1.0 1.0 0.0 0 1 1.0 0.5 0.5 0 0.0]*1.6+ory; 
Hp_iv=plot(x,y,'b-','Visible','off'); %i-v axes
x=[0.8 2 2.3 2.35 2.4 2.45 2.5 2.55 2.6 2.65 2.7 2.75 2.8 2.85 2.9 2.95 3];
y=0.02*(exp((x-2.3)/0.2)-1)+1.3;
Hp_iv_tr=plot(x,y,'c-','Visible','off'); %The I-V line
x_q=2.3; y_q=1.5;
Hp_iv_q=plot(x_q,y_q,'c.','Markersize',20,'Visible','off'); %The "operating" point
Ht_curr=text(0.05,0.13,'Current','Units','normalized','Rotation',90,...
   'Visible','off','Color','b','HorizontalAlignment','center');
Ht_volt=text(0.19,0.03,'Voltage','Units','normalized','Visible','off',...
   'Color','b','HorizontalAlignment','center');
Ht_time=text(0.59,0.03,'Time','Units','normalized','Visible','off',...
   'Color','b','HorizontalAlignment','center');
Ht_efn=text(0.01,0.1*Efn-0.01,'E_{FN}','Units','normalized','Visible',...
   'off','Color','b','FontSize',11);
Ht_efp=text(1.01,0.1*Efp,'E_{FP}','Units','normalized','Visible','off','Color','y','FontSize',11);

%Calculate kt so that thermal equilibrium is for Efn=Efp
kt=(Ecp-Ecn)/log(Nel_n/Nel_p);

%Set the electron/hole  plots 
Hprecel=plot(-1,-1,'k.','Markersize',12,'Visible','off','EraseMode','background');
xel_n=rand(1,Nel_n)*xn; xel_p=rand(1,Nel_p)*(10-xp)+xp;
Ekin_el=-kt*log(rand(1,(Nel_n+Nel_p)));
yel_n=Ekin_el(1:Nel_n)+Ecn; yel_p=Ekin_el((Nel_n+1):(Nel_n+Nel_p))+Ecp;
Hpel =plot([xel_n xel_p],[yel_n yel_p],'k.','Markersize',12,'Visible','off');
Hpstel=plot(-1,-1,'b.','Markersize',12,'Visible','off');
Hprecho=plot(-1,-1,'w.','Markersize',12,'Visible','off','EraseMode','background');
xho_p=rand(1,Nho_p)*(10-xp)+xp; xho_n=rand(1,Nho_n)*xn; 
Ekin_ho=-kt*log(rand(1,(Nho_p+Nho_n)));
yho_p=Evp-Ekin_ho(1:Nho_p); yho_n=Evn-Ekin_ho((Nho_p+1):(Nho_p+Nho_n));
Hpho =plot([xho_n xho_p],[yho_n yho_p],'w.','Markersize',12,'Visible','off');
Hpstho=plot(-1,-1,'y.','Markersize',12,'Visible','off');


%Set the recomb. stars plot
Hprecst=plot(-1,-1,'y*','Markersize',12,'Visible','off','EraseMode','xor');

%Set the "trace" electron/hole plots
Hptrel=plot(-1,-1,'b.','Markersize',12,'Visible','off','EraseMode','none');
Hptrho=plot(-1,-1,'y.','Markersize',12,'Visible','off','EraseMode','none');
hold off;

% MAIN LOOP
debcon1_usd9=0; debcon3_usd9=0; debcon4_usd9=0; debcon5_usd9=0;
while debcon4_usd9==0 % Main loop
   
   while (debcon3_usd9+debcon4_usd9)<1 %Reading the Intro text and setting the starting voltage
   %Set the slider value      
   set(Hb1c ,'String',num2str(round(get(Hb1 ,'Value')*10)*0.1));
   if ishandle(Hf)==0;
   return; end;
   pause(0.2);
   if ishandle(Hf)==0;
   return; end;
   end
   volt =round(get(Hb1 ,'Value')*10)*0.1; vold=volt;
   Efn=Efp+volt*2/1.12; 
   Ecn=Efn+0.2; Evn=Efn-1.8; % P-type EBD is fixed
   dwdepl=0.5*sqrt(abs(volt)/3); xn=4.5+sign(volt)*dwdepl; xp=5.5-sign(volt)*dwdepl;
   xel_n=rand(1,Nel_n)*xn; xel_p=rand(1,Nel_p)*(10-xp)+xp; xel=[xel_n xel_p];
   xho_p=rand(1,Nho_p)*(10-xp)+xp; xho_n=rand(1,Nho_n)*xn; xho=[xho_n xho_p];
   x_q=volt+2.3; y_q=0.02*(exp(volt/0.2)-1)+1.5;
   
if debcon4_usd9==0; %CLOSE is not activated, debcon3 must be 1 (START is activated)
   
   set(Ht ,'Visible','off')
   set([Hb3 Hb4],'Enable','off');
   set(Hbs,'Enable','on');
   
% ANIMATION LOOP - SETTINGS        
   time_step=0; Qs=0; %This is for the i-t plot
   x_time=[1:1]*0+orx; y_time=y_time*0; y_time=[1:1]*0+ory+1;
   while debcon5_usd9<1; %The animation loop
      time_step=time_step+1; 
         if time_step>100; %Reset the matrix dimensions, to start the trace from the origin
            y_fin=y_time(time_step-1);
            time_step=1; 
            x_time=[1:1]*0; y_time=[1:1]*0;
            y_time(1)=y_fin;
         end
         %Enable change of EBD during the animation
            if ishandle(Hf)==0;
            return; end;
      set(Hb1c ,'String',num2str(round(get(Hb1 ,'Value')*10)*0.1));
      %%Calculate the EB numbers 
      if debcon1_usd9>0 vold=volt; end %Remember the voltage before it is changed
      volt =round(get(Hb1 ,'Value')*10)*0.1;
      Efn=Efp+volt*2/1.12; 
      Ecn=Efn+0.2; Evn=Efn-1.8; % P-type EBD is fixed
      dwdepl=0.5*sqrt(abs(volt)/3); xn=4.5+sign(volt)*dwdepl; xp=5.5-sign(volt)*dwdepl;
      x_q=volt+2.3; y_q=0.02*(exp(volt/0.2)-1)+1.3; %Change Q point on i-v plot
      
      %% Display the EBD...
         if ishandle(Hf)==0;
         return; end;
      x=[0 xn xn 0]; y=[Ecn Ecn Evn Evn]; %N-type Eg (fixed)
      set(Hph_n,'Xdata',x,'Ydata',y,'Visible','on'); 
      x=[xn xp xp xn]; y=[Ecn Ecp Evp Evn]; %Depl. layer Ec
      set(Hph_d,'Xdata',x,'Ydata',y,'Visible','on'); 
      x=[xp 10 10 xp]; y=[Ecp Ecp Evp Evp]; %P-type Eg(0)
      set(Hph_p,'Xdata',x,'Ydata',y,'Visible','on'); 
      set(Hp_efn,'Xdata',[0 xn],'Ydata',[Efn Efn],'Visible','on'); 
      set(Hp_efp,'Xdata',[xp 10],'Visible','on'); %Efp is fixed
      % and i-t and i-v plots
      set(Hp_it,'Visible','on'); 
      set(Hp_iv,'Visible','on'); set(Hp_iv_tr,'Visible','on');
      set(Hp_iv_q,'Xdata',x_q,'Ydata',y_q,'Visible','on');
      set(Ht_curr,'Visible','on');
      set(Ht_volt,'Visible','on');
      set(Ht_time,'Visible','on');
      set(Ht_efn,'Position',[0.01 0.1*Efn-0.01 0],'Visible','on');
      set(Ht_efp,'Visible','on');

            
      %Reset the current counters
      if debcon1_usd9>0; %The voltage has been changed -> reset the counters
         debcon1_usd9=0;
      end
             
      %%Set EBD reference for y's
      for ist=1:(Nel_n+Nel_p); % set reference Ec for every electron;
            if xel(ist)<xn; y_Ec(ist)=Ecn; else
               if xel(ist)>xp; y_Ec(ist)=Ecp; else
               y_Ec(ist)=(Ecp-Ecn)*(xel(ist)-xn)/(xp-xn)+Ecn;   
               end   
            end      
      end      
      for ist=1:(Nho_n+Nho_p); % set reference Ev for every hole;
            if xho(ist)>xp; y_Ev(ist)=Evp; else
               if xho(ist)<xn; y_Ev(ist)=Evn; else
               y_Ev(ist)=(Evp-Evn)*(xho(ist)-xn)/(xp-xn)+Evn;   
               end   
            end  
      end      
          
      %Generate random scattering length,...
      lsc_el=randn(1,(Nel_n+Nel_p))*lsc;
      lsc_ho=randn(1,(Nho_n+Nho_p))*lsc;
      
      %and random kinetic energy 
      Ekin_el=-kt*log(rand(1,(Nel_n+Nel_p)));
      Ekin_ho=-kt*log(rand(1,(Nho_n+Nho_p)));
      
      %Set particles' y for the starting x
      yel=y_Ec+Ekin_el;
      yho=y_Ev-Ekin_ho;
      
  % ANIMATION LOOP - CALCULATIONS DUE TO CHANGED X POSITION    
      %Find intended final x positions
      xel_old=xel; xho_old=xho;
      xel=xel+lsc_el; xho=xho+lsc_ho;
           
      %%Find the carriers in the depl. layer (to produce "tail")
      %% and the stored-charge for recombination
      Nel_out_n=0; Nel_out_p=0; Nho_out_p=0; Nho_out_n=0; 
      Nel_rec=0; Nho_rec=0; xel_rec=-1*[1:1]; yel_rec=-1*[1:1];
      xho_rec=-1*[1:1]; yho_rec=-1*[1:1];
      Nel_st=0; xel_st=-1*[1:1]; yel_st=-1*[1:1];
      Nho_st=0; xho_st=-1*[1:1]; yho_st=-1*[1:1];
      tel_out_n=[1:1]; tel_out_p=[1:1]; tho_out_p=[1:1]; tho_out_n=[1:1];
      %%%Electrons
      for ist=1:Nel_n;
         if xel(ist)<0; xel(ist)=rand*xn; end %Keep electrons inside the figure
         if xel(ist)>xn; %A majority electron outside N-type
            if xel_old(ist)<xn; %This is a newly injected el.;
                     Nel_out_n=Nel_out_n+1; 
                     tel_out_n(1,Nel_out_n)=xel_old(ist); 
                     tel_out_n(2,Nel_out_n)=xel(ist);
                     tel_out_n(3,Nel_out_n)=yel(ist);
                  end
                  if yel(ist)>Ecp; if xel(ist)>xp; %This is a stored-charge electron
                        Nel_st=Nel_st+1;
                        xel_st(Nel_st)=xel(ist); yel_st(Nel_st)=yel(ist);
                  if rand<Rrec;  %The el. is recombined.
                     Nel_rec=Nel_rec+1;
                     xel_rec(Nel_rec)=xel(ist); yel_rec(Nel_rec)=yel(ist);
                     xel(ist)=rand*xn; yel(ist)=Ecn+Ekin_el(ist);
                  end;         
               if xel(ist)>10; %This electron closed the circuit
                  xel(ist)=rand*xn; yel(ist)=Ecn+Ekin_el(ist); 
               end;     end;  
            else 
               x_cs=(xp-xn)*(yel(ist)-Ecn)/(Ecp-Ecn)+xn; %x crossecting point
               if xel(ist)>x_cs; %This electron hit the Ec
               xel(ist)=xel(ist)-2*(xel(ist)-x_cs); end   
            end
         end
      end
      
      for ist=(Nel_n+1):(Nel_p+Nel_n);
         if xel(ist)>10; xel(ist)=rand*(10-xp)+xp; end %Keep electrons inside the figure
         if xel(ist)<xp; if xel_old(ist)>xp; %New min. elec. in the depl. layer
               Nel_out_p=Nel_out_p+1; 
               tel_out_p(1,Nel_out_p)=xel_old(ist); 
               tel_out_p(2,Nel_out_p)=xel(ist);
               tel_out_p(3,Nel_out_p)=yel(ist);
         end; end;
         if xel(ist)<xn; %This electron loses itself in N-type; place another in P-type
            xel(ist)=rand*(10-xp)+xp; yel(ist)=Ecp+Ekin_el(ist); 
         else if yel(ist)<Ecp; %A minority electron inside the depl. layer
               x_cs=(xp-xn)*(yel(ist)-Ecn)/(Ecp-Ecn)+xn; %x crossecting point
               if xel(ist)>x_cs; %This electron hit the Ec
               xel(ist)=xel(ist)-2*(xel(ist)-x_cs); end 
              end 
         end
      end
      %%%Holes
      for ist=1:Nho_n;
         if xho(ist)<0; xho(ist)=rand*xn; end %Keep holes inside the figure
         if xho(ist)>xn; if xho_old(ist)<xn; %New min. hole. in the depl. layer
               Nho_out_n=Nho_out_n+1; 
               tho_out_n(1,Nho_out_n)=xho_old(ist); 
               tho_out_n(2,Nho_out_n)=xho(ist);
               tho_out_n(3,Nho_out_n)=yho(ist);
         end; end;
         if xho(ist)>xp; %This hole loses itself inside P-type; place another in Ntype
            xho(ist)=rand*xn; yho(ist)=Evn-Ekin_ho(ist); 
         else if yho(ist)>Evn; % A minority hole inside the depl. layer
                 x_cs=xn+(xp-xn)*(yho(ist)-Evn)/(Evp-Evn); %x crossecting point
                 if xho(ist)<x_cs; %This hole hit the Ec
                 xho(ist)=xho(ist)+2*(x_cs-xho(ist)); end 
              end
         end   
      end
      for ist=(Nho_n+1):(Nho_n+Nho_p);
          if xho(ist)>10; xho(ist)=rand*(10-xp)+xp; end %Keep holes inside the figure
          if xho(ist)<xp; %A majoity hole in the depl. layer
             if xho_old(ist)>xp; %This is a newly injected hole
                Nho_out_p=Nho_out_p+1;
                tho_out_p(1,Nho_out_p)=xho_old(ist); 
                tho_out_p(2,Nho_out_p)=xho(ist);
                tho_out_p(3,Nho_out_p)=yho(ist);
             end
             if yho(ist)<Evn; if xho(ist)<xn;  %This is a stored-charge hole
                   Nho_st=Nho_st+1; 
                   xho_st(Nho_st)=xho(ist); yho_st(Nho_st)=yho(ist);
                  if rand<Rrec;   %The hole is recombined.
                     Nho_rec=Nho_rec+1;
                     xho_rec(Nho_rec)=xho(ist); yho_rec(Nho_rec)=yho(ist);
                     xho(ist)=rand*(10-xp)+xp; yho(ist)=Evp-Ekin_ho(ist);
                  end;         end;  
             else 
                x_cs=xn+(xp-xn)*(yho(ist)-Evn)/(Evp-Evn); %x crossecting point
                if xho(ist)<x_cs; %This hole hit the Ec
                xho(ist)=xho(ist)+2*(x_cs-xho(ist)); end     
             end 
          end
      end
      
% ANIMATION LOOP - SELECT "TRACE CARRIERS"
      %Order "trace" maj. electrons by energy
      if Nel_out_n>5; 
         for isto=1:(Nel_out_n-5);
            for ist=1:(Nel_out_n-1);
               if tel_out_n(3,(ist+1))>tel_out_n(3,ist);
                  te(1)=tel_out_n(1,ist); te(2)=tel_out_n(2,ist); te(3)=tel_out_n(3,ist);
                  tel_out_n(1,ist)=tel_out_n(1,ist+1); tel_out_n(2,ist)=tel_out_n(2,ist+1); tel_out_n(3,ist)=tel_out_n(3,ist+1);
                  tel_out_n(1,ist+1)=te(1); tel_out_n(2,ist+1)=te(2); tel_out_n(3,ist+1)=te(3);
               end   
            end
         end
         Nel_out_n=5;
      end   
      %Order "trace" min. electrons by position
      if Nel_out_p>1;
         for isto=1:(Nel_out_p-1);
            for ist=1:(Nel_out_p-1);
               if tel_out_p(1,(ist+1))<tel_out_p(1,ist);
                  te(1)=tel_out_p(1,ist); te(2)=tel_out_p(2,ist); te(3)=tel_out_p(3,ist);
                  tel_out_p(1,ist)=tel_out_p(1,ist+1); tel_out_p(2,ist)=tel_out_p(2,ist+1); tel_out_p(3,ist)=tel_out_p(3,ist+1);
                  tel_out_p(1,ist+1)=te(1); tel_out_p(2,ist+1)=te(2); tel_out_p(3,ist+1)=te(3);
               end   
            end
         end
         Nel_out_p=1;
      end   
      
      %Order "trace" maj. holes by energy
      if Nho_out_p>5; 
         for isto=1:(Nho_out_p-5);
            for ist=1:(Nho_out_p-1);
               if tho_out_p(3,(ist+1))<tho_out_p(3,ist);
                  te(1)=tho_out_p(1,ist); te(2)=tho_out_p(2,ist); te(3)=tho_out_p(3,ist);
                  tho_out_p(1,ist)=tho_out_p(1,ist+1); tho_out_p(2,ist)=tho_out_p(2,ist+1); tho_out_p(3,ist)=tho_out_p(3,ist+1);
                  tho_out_p(1,ist+1)=te(1); tho_out_p(2,ist+1)=te(2); tho_out_p(3,ist+1)=te(3);
               end   
            end
         end
         Nho_out_p=5;
      end   
      %Order "trace" min. electrons by position
      if Nho_out_n>1;
         for isto=1:(Nho_out_n-1);
            for ist=1:(Nho_out_n-1);
               if tho_out_n(1,(ist+1))>tho_out_n(1,ist);
                  te(1)=tho_out_n(1,ist); te(2)=tho_out_n(2,ist); te(3)=tho_out_n(3,ist);
                  tho_out_n(1,ist)=tho_out_n(1,ist+1); tho_out_n(2,ist)=tho_out_n(2,ist+1); tho_out_n(3,ist)=tho_out_n(3,ist+1);
                  tho_out_n(1,ist+1)=te(1); tho_out_n(2,ist+1)=te(2); tho_out_n(3,ist+1)=te(3);
               end   
            end
         end
         Nho_out_n=1;
      end   

% ANIMATION LOOP - CALCULATIONS FOR i-v AND i-t PLOTS
              
      %Current calculation fot the i-t plot
      I=0.02*(exp(volt/0.2)-1);
      if volt<vold; %There may be discharging current
         Qs=Qs+Idis;
         if Qs>0; I=I+Idis; 
         else vold=volt; end
      else Qs=taus*I;   
      end   
          
      x_time(time_step)=0.06*time_step+3.5;
      y_time(time_step)=I+1.3; 
      
% ANIMATION LOOP - SETTING THE PLOTS AND PLOTTING "NON-TRACE" ELEMENTS     
      
%Set the new  positions on the plots
      if ishandle(Hf)==0;
      return; end;
      set(Hprecel, 'Xdata',xel_rec,'Ydata',yel_rec,'Visible','on');
      set(Hprecho, 'Xdata',xho_rec,'Ydata',yho_rec,'Visible','on');
      set(Hprecst, 'Xdata',[xel_rec xho_rec],'Ydata',[yel_rec yho_rec],'Visible','on');
      set(Hpel,'Xdata',xel,'Ydata',yel,'Visible','on');
      set(Hpho,'Xdata',xho,'Ydata',yho,'Visible','on');
      set(Hpstel, 'Xdata',xel_st,'Ydata',yel_st,'Visible','on');
      set(Hpstho, 'Xdata',xho_st,'Ydata',yho_st,'Visible','on');
      set(Hp_it_tr,'Xdata',x_time,'Ydata',y_time,'Visible','on');
      set(Hptrel,'Visible','off'); set(Hptrho,'Visible','off');
         if ishandle(Hf)==0;
         return; end;
      drawnow;
      
% ANIMATION LOOP - CALCULATING, SETTING AND PLOTTING "TRACES"     
%Set the constants and the initial positions  - Electrons      
ist=0; xtr=[1:1]*0; ytr=[1:1]*0;
if Nel_out_n>0;
for ist=1:Nel_out_n; 
xtr(ist)=tel_out_n(1,ist);
ytr(ist)=tel_out_n(3,ist);
  if ytr(ist)>Ecp; 
   dx(ist)=(tel_out_n(2,ist)-tel_out_n(1,ist))/N_tr_st;
   x_cs(ist)=10;
  else 
   x_cs(ist)=(xp-xn)*(ytr(ist)-Ecn)/(Ecp-Ecn)+xn;
   dx(ist)=(x_cs(ist)-tel_out_n(1,ist)+x_cs(ist)-tel_out_n(2,ist))/N_tr_st;
  end
end
end
if Nel_out_p>0; 
   xtr(ist+1)=tel_out_p(1,1);
   ytr(ist+1)=tel_out_p(3,1);
   dx(ist+1)=(tel_out_p(1,1)-tel_out_p(2,1))/N_tr_st; %The min. el.
end   

%Set the constants and the initial positions - Holes      
ist=0; xtr_ho=[1:1]*0; ytr_ho=[1:1]*0;
if Nho_out_p>0;
for ist=1:Nho_out_p; 
xtr_ho(ist)=tho_out_p(1,ist);
ytr_ho(ist)=tho_out_p(3,ist);
  if ytr_ho(ist)<Evn; 
   dx_ho(ist)=(tho_out_p(1,ist)-tho_out_p(2,ist))/N_tr_st;
   x_cs_ho(ist)=0;
  else 
     x_cs_ho(ist)=xn+(xp-xn)*(ytr_ho(ist)-Evn)/(Evp-Evn);
     dx_ho(ist)=(-2*x_cs_ho(ist)+tho_out_p(1,ist)+tho_out_p(2,ist))/N_tr_st;
  end
end
end
if Nho_out_n>0; 
   xtr_ho(ist+1)=tho_out_n(1,1);
   ytr_ho(ist+1)=tho_out_n(3,1);
   dx_ho(ist+1)=(tho_out_n(2,1)-tho_out_n(1,1))/N_tr_st; %The min. ho.
end  

%The loop for the trace animations
dir_f=[1:Nel_out_n]*0+1; dir_f_ho=[1:Nho_out_p]*0+1;
for iist=1:N_tr_st; 
   
ist=0; %If Nel_out_n=0, ist+1 index for el_out_p will be 1 
if Nel_out_n>0;
   for ist=1:Nel_out_n;
      xtr(ist)=xtr(ist)+dir_f(ist)*dx(ist);
      if ytr(ist)<Ecp; if xtr(ist)>x_cs(ist); %The electron hit the Ec
         dir_f(ist)=-1; xtr(ist)=x_cs(ist);
      end; end 
   end %for ist
end %if
if Nel_out_p>0;  
   xtr(ist+1)=xtr(ist+1)-dx(ist+1); 
end   
   if ishandle(Hf)==0;
   return; end;
set(Hptrel,'Xdata',xtr,'Ydata',ytr,'Visible','on');

ist=0; %If Nho_out_p=0, ist+1 index for ho_out_n will be 1 
if Nho_out_p>0;
   for ist=1:Nho_out_p;
      xtr_ho(ist)=xtr_ho(ist)-dir_f_ho(ist)*dx_ho(ist);
      if ytr_ho(ist)>Evn; if xtr_ho(ist)<x_cs_ho(ist); %The hole hit the Ev
         dir_f_ho(ist)=-1; xtr_ho(ist)=x_cs_ho(ist);
      end; end 
   end %for ist
end %if
if Nho_out_n>0;  
   xtr_ho(ist+1)=xtr_ho(ist+1)+dx_ho(ist+1); 
end 
   if ishandle(Hf)==0;
   return; end;
set(Hptrho,'Xdata',xtr_ho,'Ydata',ytr_ho,'Visible','on');

   if ishandle(Hf)==0;
   return; end;
tic; % Delay for less <1s
while toc<0.001;
pause(0.001);
end   
   if ishandle(Hf)==0;
   return; end;
end %iist
          
end % while/ END ANIMATION LOOP

   debcon5_usd9=0; debcon3_usd9=0; 
   if ishandle(Hf)==0;
   return; end;
   set(Hph_n,'Visible','off'); set(Hph_d,'Visible','off'); set(Hph_p,'Visible','off');
   set(Hpel ,'Visible','off'); 
   set(Hpho ,'Visible','off'); 
   set(Hpstel,'Visible','off');
   set(Hpstho,'Visible','off');
   set(Hp_efn,'Visible','off'); set(Hp_efp,'Visible','off');
   set(Hp_it,'Visible','off');  
   set(Hp_it_tr,'Visible','off');
   set(Hp_iv,'Visible','off'); set(Hp_iv_tr,'Visible','off'); 
   set(Hp_iv_q,'Visible','off'); 
   set(Ht_curr,'Visible','off');
   set(Ht_volt,'Visible','off');
   set(Ht_time,'Visible','off');
   set(Hptrel,'Visible','off'); set(Hptrho,'Visible','off');
   set(Hprecel,'Visible','off'); set(Hprecho,'Visible','off');
   set(Hprecst,'Visible','off');
   set(Ht_efn,'Visible','off');
   set(Ht_efp,'Visible','off');
   set(Ht ,'Visible','on')
   set([Hb3 Hb4],'Enable','on');
   set(Hbs,'Enable','off');
  
end %if debcon4_usd9==0;
end % END MAIN LOOP
close
