%-----------------------------------------------------------------------
% Find the temperature at which the molecules have the given mean energy
% based on the Maxwell-Boltzmann distribution function
% Programming Note:
%   Call QUAD8/QUADL to integrate
%   Call FZERO to solve f(x)=0
%-----------------------------------------------------------------------
% Instructor: Nam Sun Wang
%-----------------------------------------------------------------------

% Start fresh ----------------------------------------------------------
      clear all

      global xm
      global v_target

% Print out a program header -------------------------------------------
      disp('Find the temperature at which the molecules have the given "mean" energy')
      disp('based on the Maxwell-Boltzmann distribution function')
      disp('             /oo          ')
      disp('  v_mean(T)= |  v*f(T,v)dv')
      disp('             /0           ')
      disp('where v = speed           ')
      disp('      T = temperature     ')
      disp('      f = Maxwell-Boltzmann distribution function')
      disp('      v_mean = mean speed ')
      disp('      "mean" energy = mass*(v_mean)^2/2 ')
      disp(' ')

% Input data -----------------------------------------------------------
      xm     = input ('Enter molecular mass in Da (e.g., 32. for oxygen): ');
      energy = input('Enter "mean" energy in Joule/mole (e.g., 5000.): ');

% Change units to be consistent with the rest --------------------------
      xm = xm/1000.;
      v_target = sqrt(2.*energy/xm);

% Find the temperature by calling a routine to solve f(T)=0 ------------
      guess = 400.;       % Initial guess
      T = fzero('boltfeq0', guess);

% Print out the results ------------------------------------------------
      disp(' ')
      disp([ 'Temperature in K (numerical solution) =', num2str(T) ])

% average energy = xm*v_mean**2/2 = 4*xk*T/pi
%     pi = 3.141593
      xk = 8.319;   % Joule/mole-K
      disp([ 'Temperature in K (analytical results) =', num2str(energy*pi/4./xk) ])

      disp(' ')
      disp('Condition:')
      disp([ '    Molecular mass (Dalton): ', num2str(xm*1000.)   ])
      disp([ '    "Mean" energy (Joules/mole): ', num2str(energy) ])
      disp([ '    Mean speed (m/sec): ', num2str(v_target)        ])