%
%   MAXWELL_2A.M   Plots Maxwellian Distribution for Several Temperatures
%
%   This is a demo that illustrates several aspects of programming within the 
%   Matlab environment.  The goal here is simply to evaluate and plot the 
%   the Maxwellian Distribution for several temperatures.  Both linear and  
%   logarithmic energy scales are used.  The function of interest is:
%        f(E,T) = (2*pi/(pi*kT)^1.5)*sqrt(E)*exp(-E/kT)
%  
%   In evaluating this function, we will illustrate some of the vector processing
%   capabilities and simple 2-d plotting functions available in Matlab.
%
%   Related files, MAXWELL_1.M and MAXWELL_2B.M, perform similar operations. 
%   MAXWELL_1.M is a simplified version of this file and only plots f(E) for one 
%   temperature at a time.  MAXWELL_2B.M does a similar evaluation as the current 
%   file, with a focus on some of Matlab's 3-D plotting options. 
%
%   File prepared by J. R. White, UMass-Lowell (last update: September 2017)
%

      clear all,  close all,  nfig = 0;
%
%   define the Boltzmann constant 
      k = 8.6170e-5;               % Boltzmann constant (eV/K)
%
%   define desired temperatures and evaluate some other equation constants
      Tc = [20 150 300];           % desired temperature vector (C)
      T = Tc+273;                  % convert to absolute temperature (K)
      kT = k*T;                    % energy associated with given temperature (eV)
      c = 2*pi./(pi*kT).^1.5;      % normalization constant in equation for f(E)
%
%   Case 1  Evaluate function, f(E,T), on a linear energy grid
      E1o = 0;    E1f = 0.25;   NE1 = 251;   E1 = linspace(E1o,E1f,NE1)';
      NT = length(T);   F1 = zeros(NE1,NT);   st = cell(NT,1);
      for j = 1:NT
        F1(:,j) = c(j)*sqrt(E1).*exp(-E1/kT(j));  % evaluate function for all E & Tj
      end
%
%   Case 2  Evaluate function, f(E,T), on a logarithmic energy grid
      E2o = -5;   E2f = 0;     NE2 = 251;   E2 = logspace(E2o,E2f,NE2)';
      F2 = zeros(NE2,NT);   
      for j = 1:NT
        F2(:,j) = c(j)*sqrt(E2).*exp(-E2/kT(j));  % evaluate function for all E & Tj
      end
%
%   now plot the results (both linear and semilog scales)
      nfig = nfig+1;   figure(nfig)
      subplot(2,1,1),plot(E1,F1,'LineWidth',2), grid
      title('Maxwell\_2a:  Maxwellian Distribution for Different Temperatures'); 
      xlabel('Energy (eV)'),ylabel('Probability per eV, f(E)')
      for j = 1:NT,   tt = sprintf('T = %3g C',Tc(j));   st(j) = cellstr(tt);   end
      legend(st)
%
      subplot(2,1,2),semilogx(E2,F2,'LineWidth',2), grid
      xlabel('Energy (eV)'),ylabel('Probability per eV, f(E)')
      for j = 1:NT;   gtext(['T = ',num2str(Tc(j)),' C']);   end
%
%   plot linear case again -- this time with different line styles (max of 6 lines)
      Ncm = 6;   scm = ['b- ';'g--';'r-.';'m: ';'c- ';'k--']; % set color/marker code 
      nfig = nfig+1;   figure(nfig),   hold on
      for j = 1:NT
        plot(E1,F1(:,j),scm(j,:),'LineWidth',2)
        tt = sprintf('T = %3g C',Tc(j));   st(j) = cellstr(tt); 
      end
      title('Maxwell\_2a:  Maxwellian Distribution for Different Temperatures'); 
      xlabel('Energy (eV)'),ylabel('Probability per eV, f(E)')
      grid on,  legend(st),  hold off
%
%   end of program
