function ef_usd;
%Fermi-Dirac distribution and Fermi level:
%Carrier concentration dependence on 
%material properties, temperature and doping level.
%(Section 1.5.2 in Understanding Semiconductor Devices)

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


global con_sl_usd7; global con_b_usd7; global con_go_usd7; global con_b_ex_usd7;

%OPEN PARENT FIGURE
Hf_ef=figure('NumberTitle','off','Name','Section 1.5.2: Fermi-Dirac Distribution and Fermi level',...
   'Color',[0.8 0.8 0.8]);
set(gcf,'InvertHardCopy','off');

Ha_ex=axes('Xlim',[0 3],'Ylim',[-5 5],'Visible','off',...
                 'Position',[0.02 0 0.3 1]); %Axes for E-x diagram
              
%READ THE INTRO TEXT    
con_go_usd7=0;              
introtext ={'UNDERSTANDING SEMICONDUCTOR DEVICES ';...
      'Section 1.3.2: Fermi-Dirac Distribution and Fermi level';...
      ' ';...
      'Set the material energy gap, temperature and doping level,';...
      'to see the influence on carrier concentrations. '};
Ht_intro=text(0.4,0.6,char(introtext),'Units','normalized');
Hp_cont=uicontrol(Hf_ef,'Units','normalized','Position',[0.1 0.2 0.7 0.1],...
   'Style','pushbutton','String','CLICK HERE TO CONTINUE!','CallBack','con_go_usd7=1;');
con_go_usd7=0;
while con_go_usd7<1; 
pause(0.2);
      if ishandle(Hf_ef)==0;
         return; 
      end;
end
con_go_usd7=0;
set(Ht_intro,'Visible','off');
set(Hp_cont,'Visible','off','Enable','off');
drawnow;

%CONSTANTS
k=8.625e-5; Ncv=2.6e+19; Vt=k*300; n=1e10; p=1e10;

%SETTING THE CONTROL PANEL
con_sl_ef_usd7=0; con_sl_eg_usd7=0; con_sl_t_usd7=0; con_b_usd7=0; con_b_ex_usd7=0;
Hframe=uicontrol(Hf_ef,'Units','normalized','Position',[0.72 0.01 0.27 0.98],...
   'Style','frame');
%Temperature slider
Ht_tu=uicontrol(Hf_ef,'Units','normalized','Position',[0.73 0.35 0.25 0.05],...
   'Style','text','String','TEMPERATURE');
Hsl_t=uicontrol(Hf_ef,'Units','normalized','Position',[0.73 0.30 0.25 0.05],...
   'Style','slider','min',100,'max',600,'Value',300,'CallBack','con_sl_usd7=4;');
Ht_td=uicontrol(Hf_ef,'Units','normalized','Position',[0.73 0.25 0.25 0.05],...
   'Style','text','String',[num2str(round(get(Hsl_t,'Value'))) ' K']);
T=(round(get(Hsl_t,'Value')));

%P-type doping
Ht_dpu=uicontrol(Hf_ef,'Units','normalized','Position',[0.73 0.55 0.25 0.05],...
   'Style','text','String','P-TYPE DOPING','BackgroundColor',[0.6 0 0],...
   'ForegroundColor','w');
Hsl_dp=uicontrol(Hf_ef,'Units','normalized','Position',[0.73 0.50 0.25 0.05],...
   'Style','slider','min',14,'max',19,'Value',14,'CallBack','con_sl_usd7=3;');
Ht_dpd=uicontrol(Hf_ef,'Units','normalized','Position',[0.73 0.45 0.25 0.05],...
   'Style','text','String',[num2str(round(10^(get(Hsl_dp,'Value')))) ' cm-3'],...
   'BackgroundColor',[0.6 0 0],'ForegroundColor','w');
Na=round(10^(get(Hsl_dp,'Value')));

%N-type doping
Ht_dnu=uicontrol(Hf_ef,'Units','normalized','Position',[0.73 0.70 0.25 0.05],...
   'Style','text','String','N-TYPE DOPING','BackgroundColor',[0 0.5 0],...
   'ForegroundColor','w');
Hsl_dn=uicontrol(Hf_ef,'Units','normalized','Position',[0.73 0.65 0.25 0.05],...
   'Style','slider','min',14,'max',19,'Value',14,'CallBack','con_sl_usd7=2;');
Ht_dnd=uicontrol(Hf_ef,'Units','normalized','Position',[0.73 0.60 0.25 0.05],...
   'Style','text','String',[num2str(round(10^(get(Hsl_dn,'Value')))) ' cm-3'],...
   'BackgroundColor',[0 .5 0],'ForegroundColor','w');
Nd=round(10^(get(Hsl_dn,'Value')));

%Eg slider
Ht_egu=uicontrol(Hf_ef,'Units','normalized','Position',[0.73 0.90 0.25 0.05],...
   'Style','text','String','MATERIAL');
Hsl_eg=uicontrol(Hf_ef,'Units','normalized','Position',[0.73 0.85 0.25 0.05],...
   'Style','slider','min',0.2,'max',3.0,'Value',1.12,'CallBack','con_sl_usd7=1;');
Ht_egd=uicontrol(Hf_ef,'Units','normalized','Position',[0.73 0.80 0.25 0.05],...
   'Style','text','String',['Eg= ' num2str(round(100*get(Hsl_eg,'Value'))/100) ' eV']);
Eg=round(100*get(Hsl_eg,'Value'))/100;

%The pushbutton
Hb_ex=uicontrol(Hf_ef,'Units','normalized','Position',[0.73 0.02 0.25 0.10],...
   'Style','pushbutton','String','CLOSE','CallBack','con_b_ex_usd7=1;');

%SET THE GRAPHS
%E-x diagram
aa=Eg/2; ww=2; Ef=0;
hold on;
Hp_Evx=patch([0 2.8 2.8 0],[-aa -aa -(aa+ww) -(aa+ww)],'r',...
   'EdgeColor','r');
Hp_Ecx=patch([0 2.8 2.8 0],[aa aa (aa+ww) (aa+ww)],'g',...
   'EdgeColor','g');
Hl_Efx=plot([0 2.8],[Ef Ef],'b-.','LineWidth',2);

xx=[0 2.8]; 
yxv=[-aa -(aa+0.08) -(aa+0.15) -(aa+0.35)];
yxc=[aa (aa+0.08) (aa+0.15) (aa+0.35)];
for ii=1:4;
   Hl_Evx(ii)=plot(xx,[yxv(ii) yxv(ii)],'g'); 
   Hl_Ecx(ii)=plot(xx,[yxc(ii) yxc(ii)],'r'); 
end
      if ishandle(Hf_ef)==0;
         return; 
      end;
set(Hl_Evx(1),'LineWidth',3'); set(Hl_Ecx(1),'LineWidth',3);

Ekin=-5*Vt*log(rand(1,3));
xx=rand(1,3)*2.8; yy=aa+Ekin;
Hp_xel=plot(xx,yy,'k.','Markersize',15);
Ekin=-5*Vt*log(rand(1,3));
xx=rand(1,3)*2.8; yy=-aa-Ekin;
Hp_xho=plot(xx,yy,'w.','Markersize',15);

Ht_efn=text(2.8,Ef,'E_{F}','Color','b','FontSize',12,...
   'FontAngle','italic','FontWeight','bold');

set(Ha_ex,'Visible','off');
hold off;


%E-f diagram
Ha_ef=axes('Xlim',[1e-10 1],'Ylim',[-(aa+ww) (aa+ww)],...
   'Position',[0.4 (0.5-0.1*(aa+ww)) 0.30 0.2*(aa+ww)],...
   'Color',[0.8 0.8 0.8],'FontAngle','italic'); %Axes for E-f diagram
xlabel('PROBABILITY'); ylabel('ENERGY-E_{i} (eV)');
hold on;
xf=[1e-10 1]; 
yfv=[-aa -(aa+0.08) -(aa+0.15) -(aa+0.35)];
yfc=[aa (aa+0.08) (aa+0.15) (aa+0.35)];
for ii=1:4;
   Hp_Evf(ii)=plot(xf,[yfv(ii) yfv(ii)],'g'); 
   Hp_Ecf(ii)=plot(xf,[yfc(ii) yfc(ii)],'r'); 
end
      if ishandle(Hf_ef)==0;
         return; 
      end;
set(Hp_Evf(1),'LineWidth',3'); set(Hp_Ecf(1),'LineWidth',3);

Hl_Eff=plot([1e-10 1],[Ef Ef],'b-.','LineWidth',2);
yef=[-(aa+ww):0.01:(aa+ww)];
xef=1./(1+exp((yef-Ef)./Vt));
Hp_fho=plot((1-xef),yef,'w-','LineWidth',3);
Hp_fel=plot(xef,yef,'k-','LineWidth',3);

Ht_ec=text(0.8,(aa+0.25),'E_{C}','Color','b','FontSize',12,...
   'FontAngle','italic','FontWeight','bold','HorizontalAlignment','right');
Ht_ev=text(0.8,-(aa+0.15),'E_{V}','Color','b','FontSize',12,...
   'FontAngle','italic','FontWeight','bold','HorizontalAlignment','right');
hold off;

%Lin/log button
Hr_lin=uicontrol(Hf_ef,'Units','normalized','BackgroundColor',[.8 .8 .8],...
   'Position',[0.45 (0.35-0.1*(aa+ww)) 0.15 0.05],...
   'Style','radio','String','LINEAR','CallBack','con_b_usd7=1;','Value',1);
Hr_log=uicontrol(Hf_ef,'Units','normalized','BackgroundColor',[.8 .8 .8],...
   'Position',[0.6 (0.35-0.1*(aa+ww)) 0.07 0.05],...
   'Style','radio','String','LOG','CallBack','con_b_usd7=2;','Value',0);

%The function text
Ha_pos=axes('Xlim',[0 1],'Ylim',[0 1],'Visible','off',...
                 'Position',[0 0 1 1]); %Axes for equation position
Ht_f=text(0.4,(0.63+0.1*(aa+ww)),...
       'f=1/ [1+exp((E-E_{F} )/kT)]','Color','k','FontSize',12,...
       'FontAngle','italic');
Ht_fh=text(0.4,(0.57+0.1*(aa+ww)),...
       'f_{h}=1-f','Color','w','FontSize',12,...
       'FontAngle','italic');

%MAIN LOOP
while con_b_ex_usd7<1 % Main loop
   
   if con_b_usd7>0; %Select lin/log
      if con_b_usd7==1;
         if ishandle(Hf_ef)==0;
            return; 
         end;
         set(Hr_lin,'Value',1); set(Hr_log,'Value',0);
         set(Ha_ef,'XScale','linear');
         set(Ht_ec,'Position',[0.8 (aa+0.25)]);
         set(Ht_ev,'Position',[0.8 -(aa+0.15)]);
      else
         if ishandle(Hf_ef)==0;
            return; 
         end;
         set(Hr_lin,'Value',0); set(Hr_log,'Value',1);
         set(Ha_ef,'XScale','log');
         set(Ht_ec,'Position',[0.01 (aa+0.25)]);
         set(Ht_ev,'Position',[0.01 -(aa+0.15)]);
      end;   
   end; %if con_b_usd7>0
   con_b_usd7=0;
   
   if con_sl_usd7>0; %Set the slider values
      switch con_sl_usd7;
      case 1; %Eg slider
         if ishandle(Hf_ef)==0;
            return; 
         end;
         set(Ht_egd,'String',...
            ['Eg= ' num2str(round(100*get(Hsl_eg,'Value'))/100) ' eV']);
         Eg=round(100*get(Hsl_eg,'Value'))/100;
      case 2; %N-type doping slider
         if ishandle(Hf_ef)==0;
            return; 
         end;
         set(Ht_dnd, 'String',...
            [num2str(round(10^(get(Hsl_dn,'Value')))) ' cm-3']);
         Nd=round(10^(get(Hsl_dn,'Value')));
      case 3; %P-type doping slider
         if ishandle(Hf_ef)==0;
            return; 
         end;
         set(Ht_dpd, 'String',...
            [num2str(round(10^(get(Hsl_dp,'Value')))) ' cm-3']);
         Na=round(10^(get(Hsl_dp,'Value')));
      case 4; %Temperature slider
         if ishandle(Hf_ef)==0;
            return; 
         end;
         set(Ht_td,'String',...
            [num2str(round(get(Hsl_t,'Value'))) ' K']);
         T=(round(get(Hsl_t,'Value')));
      end; %switch   
      
      %Perform the calculations
      Vt=k*T; ni=Ncv*exp(-Eg/(2*Vt));
      if (Nd-Na)>0;
         n=0.5*(Nd-Na+sqrt((Nd-Na)*(Nd-Na)+4*ni*ni));
         p=ni^2/n;
         Ef=0.5*Eg-Vt*log(Ncv/n);
      else if (Nd-Na)<0;
         p=0.5*(Na-Nd+sqrt((Na-Nd)*(Na-Nd)+4*ni*ni));
         n=ni^2/p;
         Ef=-0.5*Eg+Vt*log(Ncv/p);
      else
         n=ni; p=ni; Ef=0;
      end; end   
   
      %Re-plotting the graphs
      %% E-x diagram
      aa=Eg/2;
      if ishandle(Hf_ef)==0;
         return; 
      end;
      set(Hp_Evx,'YData',[-aa -aa -aa-ww -aa-ww]);
      set(Hp_Ecx,'YData',[aa aa aa+ww aa+ww]);
      set(Hl_Efx,'YData',[Ef Ef]);
      yxv=[-aa -(aa+0.08) -(aa+0.15) -(aa+0.35)];
      yxc=[aa (aa+0.08) (aa+0.15) (aa+0.35)];
      for ii=1:4;
         set(Hl_Evx(ii),'YData',[yxv(ii) yxv(ii)]); 
         set(Hl_Ecx(ii),'YData',[yxc(ii) yxc(ii)]); 
      end
      set(Ht_efn,'Position',[2.8 Ef]);
      
      %% E-f diagram
      set(Ha_ef,'Ylim',[-(aa+ww) (aa+ww)],...
        'Position',[0.4 (0.5-0.1*(aa+ww)) 0.30 0.2*(aa+ww)]); 
     yfv=[-aa -(aa+0.08) -(aa+0.15) -(aa+0.35)];
     yfc=[aa (aa+0.08) (aa+0.15) (aa+0.35)];
     for ii=1:4;
        set(Hp_Evf(ii),'YData',[yfv(ii) yfv(ii)]); 
        set(Hp_Ecf(ii),'YData',[yfc(ii) yfc(ii)]); 
     end
     
     if get(Hr_lin,'Value')<1;
     set(Ht_ec,'Position',[0.01 (aa+0.25)]);
     set(Ht_ev,'Position',[0.01 -(aa+0.15)]);
     else
     set(Ht_ec,'Position',[0.8 (aa+0.25)]);
     set(Ht_ev,'Position',[0.8 -(aa+0.15)]);
     end
     
      set(Hl_Eff,'YData',[Ef Ef]);
      yef=[-(aa+ww):0.01:(aa+ww)];
      xef=1./(1+exp((yef-Ef)./Vt));
      set(Hp_fho,'XData',(1-xef),'YData',yef);
      set(Hp_fel,'XData',xef,'YData',yef);
                 
      %% Log/Lin buttons
      set(Hr_lin,'Position',[0.45 (0.35-0.1*(aa+ww)) 0.15 0.05]);
      set(Hr_log,'Position',[0.60 (0.35-0.1*(aa+ww)) 0.07 0.05]);
      
      %%Function text
      set(Ht_f,'Position',[0.4 (0.63+0.1*(aa+ww))]);
      set(Ht_fh,'Position',[0.4 (0.57+0.1*(aa+ww))]);
      
   end; %if con_sl_usd7>0
   con_sl_usd7=0;
   if ishandle(Hf_ef)==0;
      return; 
   end;

   drawnow;
   
   %Electron/hole random generation
   if n>2e7;
      Nel=round((log10(n)-10)^3/7)+3;
   else
      Nel=0;
   end
   if Nel>1;
      Ekin=-5*Vt*log(rand(1,Nel));
      xx=0; yy=0;
      xx=rand(1,Nel)*2.8; yy=aa+Ekin;
      if ishandle(Hf_ef)==0;
         return; 
      end;
      set(Hp_xel,'XData',xx,'YData',yy','Visible','on');
   else
      if ishandle(Hf_ef)==0;
         return; 
      end;
      set(Hp_xel,'Visible','off');
      
   end
   if p>2e7;
      Nho=round((log10(p)-10)^3/7)+3;
   else
      Nho=0;
   end
   if Nho>1;
      Ekin=-5*Vt*log(rand(1,Nho));
      xx=0; yy=0;
      xx=rand(1,Nho)*2.8; yy=-aa-Ekin;
      if ishandle(Hf_ef)==0;
         return; 
      end;
      set(Hp_xho,'XData',xx,'YData',yy','Visible','on');
   else
      if ishandle(Hf_ef)==0;
         return; 
      end;
      set(Hp_xho,'Visible','off');
   end
   
   if ishandle(Hf_ef)==0;
      return; 
   end;
   tic;
   while toc<(300/T);
      pause(0.1);
   end   
   if ishandle(Hf_ef)==0;
      return; 
   end;
  
end % End main loop
close
