%
%   MAXWELL_1.M   Plots Maxwellian Distribution for Single Temperature
%
%   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 a single temperature.  Both linear and 
%   logarithmic scales are used.  The function of interest is:
%        f(E) = (2*pi/(pi*kT)^1.5)*sqrt(E)*exp(-E/kT)
%  
%   In evaluating this function, we will also try to illustrate some of the
%   vector processing capabilities and simple 2-d plotting functions available
%   in Matlab.  This program also shows our first use of Matlab's input function
%   for communicating interactively with the user.  It also illustrates the use 
%   of the num2str command for use in building plot titles with quantitative 
%   information.
%
%   File prepared by J. R. White, UMass-Lowell (last update: August 2017)
%

      clear all,  close all,  nfig = 0;
%%
%   define the Boltzmann constant 
      k = 8.6170e-5;               % Boltzmann constant (eV/K)
%
%   input desired temperature and evaluate some other equation constants
      Tc = input('Input desired temperature for Maxwellian distribution (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), on a linear energy grid
      E1o = 0;    E1f = 0.5;   NE1 = 251;   E1 = linspace(E1o,E1f,NE1)';
      F1 = c*sqrt(E1).*exp(-E1/kT);   
%
%   Case 2  Evaluate function, f(E), on a logarithmic energy grid
      E2o = -5;   E2f = 0;     NE2 = 251;   E2 = logspace(E2o,E2f,NE2)';
      F2 = c*sqrt(E2).*exp(-E2/kT);   
%%
%   now plot the results (both linear and semilog scales)
      nfig = nfig+1;   figure(nfig)
      subplot(2,1,1),plot(E1,F1,'r-','LineWidth',2), grid
      title(['Maxwell\_1:  Maxwellian Distribution for T = ',num2str(Tc),' C']); 
      xlabel('Energy (eV)'),ylabel('Probability per eV, f(E)')
      v = axis;   v(2) = 0.25;   axis(v);      % resets max E for plots
%
      subplot(2,1,2),semilogx(E2,F2,'r-','LineWidth',2), grid
      xlabel('Energy (eV)'),ylabel('Probability per eV, f(E)')
%
%   end of program
