function X = Mxrnd(m, n, M, T)

% Implements a random number generator for the Maxwell distribution
% using an acceptance-rejection method.

% m and n are the dimensions of the output array.
% M is the molar mass of the gas in kg.
% T is the absolute temperature in degrees Kelvin.

R = 8.3*10^3; % Perfect gas constant in MKS.

% Maxwell pdf f(x):
A = 4*pi*((M/(2*pi*R*T))^(3/2));
B = M/(2*R*T);
f = inline('A*(x.^2).*exp(-B*(x.^2))', 'x', 'A', 'B');
  
p = 1/sqrt(B); % x value at distribution peak.
P = f(p, A, B); % y value at distribution peak.

% Comparison pdf g(x) (uniform):
fcross = inline('A*(x.^2).*exp(-B*(x.^2))-eps', 'x', 'A', 'B'); % f - eps.
q = fzero(fcross, 10^4, [], A, B); % x value at numerical right end of f(x) distribution.
Q = 1/q;  % y value (height) of uniform distribution g(x) on the interval [0, q].

c = P/Q;  % Scale factor to guarantee c*g(x) >= f(x) for all x.

X = zeros(m, n);

for i = 1:m*n
    accept = 0;
    while accept == 0
       r = q*rand(1);
       u = rand;
       if c*u <= f(r, A, B)/Q
           X(i) = r;
           accept = 1;
       end
    end
end