Code covered by the BSD License

### Highlights fromHydrogenic Wavefunction Visulization Tool

from Hydrogenic Wavefunction Visulization Tool by Jos Martin
Visualize hydrogenic wave-functions in 3D and watch animations of atomic state changes

calculateHydrogenicWavefunctionBounds(Z,N,L,threshold)
```function rScale = calculateHydrogenicWavefunctionBounds(Z,N,L,threshold)
% calculateHydrogenicWavefunctionBounds Find wavefunction's bounding box
%   rScale = calculateHydrogenicWavefunctionBounds(Z,N,L,threshold)
%   calculates the smallest radius (in Bohr radii) where the absolute value
%   of the wavefunction is not more than the 'threshold' multiplied by the
%   maximum value of the radial wavefunction. Z,N,L are the central
%   nuclear charge, principal quantum number and angular momentum quantum
%   number respectively.

% Copyright 2009 The MathWorks, Inc.

% To calculate the bounding radius we :
%   * Find an rMax such that the absolute value of the radial wavefunction
%     is definitely below the threshold. I.e. ensure that rScale lies in the
%     range 0 -> rMax.
%   * Search from r = rMax towards r = 0 to find the radius, rScale, where the
%     threshold is first crossed.

%   Find rMax > rScale

% Initial value, this is not critical as long as it is >> 0.
rMax = 100;

% Search for an upper limit on rScale, by increasing rMax until the radial
% equation evaluated at rMax is below the one half of the threshold.
done = false;
while ~done
r = linspace(0,rMax,500);

% Scale such that the maximum value is 1.
R = R/max(R);

% If the wavefunction is well below the threshold at r = rMax.
if R(end) < 0.5*threshold
done = true;
else
rMax = 2*rMax;
end
end

% Now the threshold is definitely crossed in the interval r = 0 -> rMax

% Find where the probability crosses the threshold.
threshCrossing = find( R >= threshold, 1, 'last');

% The value of r at which the threshold is first reached, when approaching
% r=0 from infinity.
rScale = r(threshCrossing);

end

% Function to calculate the (un-normalised) radial analytic solution
% at point r.

% The Laguerre polynomial coefficents
L_coeffs = LaguerreGen(N-L-1, 2*L+1);

% rho := 2Zr/n.a0
rho = 2*Z*r/N;

% R := exp(-rho/2) * rho^L * Laguerre^(2*L+1)_(N-L-1)(rho)
R = polyval( L_coeffs, rho) .* exp(-0.5*rho) .* (rho).^L;
end```