%
%   MAXWELL_2B.M   Plots Maxwellian Distribution for Several Temperatures and 
%                illustrates some of Matlab's 3-D surface plotting capabilities
%
%   This demo illustrates several aspects of programming within the Matlab 
%   environment.  The goal here is to evaluate and plot the Maxwellian 
%   Distribution for several temperatures and to illustrate some 3-D visualization
%   capability in Matlab.  The function of interest is
%        f(E,T) = (2*pi/(pi*kT)^1.5)*sqrt(E)*exp(-E/kT)
%   and this is evaluated on linear grids for both energy and temperature.
%
%   Related files, MAXWELL_1.M and MAXWELL_2A.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_2A.M does a similar evaluation as the current 
%   file, with a focus on plotting f(E,T) as a family of curves in a simple 2-D plot.
%
%   File prepared by J. R. White, UMass-Lowell (last update: September 2017)
%

      clear all,  close all,  nfig = 0;
%
%   define the Boltzmann constant, desired temperatures and energies
      k = 8.6170e-5;                           % Boltzmann constant (eV/K)
      Tc = 0:10:300;                           % desired temperature vector (C)
      T = Tc+273;                              % desired temperature vector (K)
      Emax = 0.2;  E = linspace(0,Emax,101);   % desired energy vector (eV)
%
%   define 2-D arrays containing Tj,Ei grid points
      [EE,TT] = meshgrid(E,T);    TTc = TT-273;                            
%
%   now evaluate Maxwellian at all grid points (be careful with dot arithmetic)
      kT = k*TT;   c = 2*pi./(pi*kT).^1.5;
      F = c.*sqrt(EE).*exp(-EE./kT);
%
%   Plot family of curves for f(E) for every 5th temperature on single 2-D plot 
      Ncm = 6;   scm = ['r- ';'g--';'b-.';'m: ';'c- ';'k--']; % set color/marker code 
      nfig = nfig+1;   figure(nfig),   hold on
      n = 0;  nn = 0;  rr = 1:5:length(Tc);  st = cell(length(rr),1);
      for j = rr                                % j is temp profile of interest
        n = n+1;                                % counter for number of curves
        nn = nn+1;  if nn > Ncm,  nn = 1;  end  % counter for line style/color
        plot(E,F(j,:),scm(nn,:),'LineWidth',2)
        tt = sprintf('T = %3g C',Tc(j));   st(n) = cellstr(tt);
      end
      title('Maxwell\_2b:  Maxwellian Distribution for Different Temps')
      xlabel('Energy (eV)'),ylabel('Probability per eV, f(E,T)')
      grid on,  legend(st),  hold off
%
%   Let's try a mesh plot of f(E,T) --- using mesh
      nfig = nfig+1;   figure(nfig)
      mesh(EE,TTc,F)
      title('Maxwell\_2b:  Maxwellian Distribution (mesh command)')
      xlabel('Energy (eV)'),ylabel('Temperature (C)'),
      zlabel('Probability per eV, f(E,T)')
%
%   Let's try a surface plot of f(E,T) --- using surf
      nfig = nfig+1;   figure(nfig)
      surf(EE,TTc,F)
      title('Maxwell\_2b:  Maxwellian Distribution (surf command)')
      xlabel('Energy (eV)'),ylabel('Temperature (C)'),
      zlabel('Probability per eV, f(E,T)')
%
%   Let's add a color bar to the surface plot --- using surf and colorbar
      nfig = nfig+1;   figure(nfig)
      surf(EE,TTc,F),  colorbar
      title('Maxwell\_2b:  Maxwellian Distribution (surf and colorbar)')
      xlabel('Energy (eV)'),ylabel('Temperature (C)'),
      zlabel('Probability per eV, f(E,T)')
%
%   Let's redo the same 3D plot with a different view --- using view
      nfig = nfig+1;   figure(nfig)
      surf(EE,TTc,F),  colorbar,  view(20,30)
      title('Maxwell\_2b:  Maxwellian Distribution (surf, colorbar, & view)')
      xlabel('Energy (eV)'),ylabel('Temperature (C)'),
      zlabel('Probability per eV, f(E,T)')
%
%   One more time - this time with "interpolated shading" --- using shading
      nfig = nfig+1;   figure(nfig)
      surf(EE,TTc,F),  colorbar,  shading interp
      title('Maxwell\_2b:  Maxwellian Distribution (surf, colorbar, & shading)')
      xlabel('Energy (eV)'),ylabel('Temperature (C)'),
      zlabel('Probability per eV, f(E,T)')
%
%   And, as a final plot, how about a contour plot with a colorbar?
      nfig = nfig+1;   figure(nfig)
      [cs,h] = contour(EE,TTc,F);  clabel(cs),  colorbar,  colormap(jet),  grid
      set(h,'LineWidth',2)
      title('Maxwellian Dist. --  Lines of Constant Prob. per eV')
      xlabel('Energy (eV)'),ylabel('Temperature (C)'),
%
%   Well, one final final plot, how about a contour plot where we set the levels?
      nfig = nfig+1;   figure(nfig)
      [cs,h] = contour(EE,TTc,F,[0.5 4 8 12 16 20]);  
      clabel(cs),  colorbar,  colormap(jet),  grid
      set(h,'LineWidth',2)
      title('Maxwellian Dist. --  Lines of Constant Prob. per eV')
      xlabel('Energy (eV)'),ylabel('Temperature (C)')
%
%   enough already...  end of program
