 # Boost :

Subject: Re: [boost] [units] - learning to use, continued :)
From: Matthias Schabel (boost_at_[hidden])
Date: 2011-11-30 12:18:30

> But still I am quite happy today, because I managed to get a
> correctly working plots of this harmonic oscillator, but without
> units, unfortunately. One small problem is that I couldn't get to
> higher energies because double precision has limited factorials.
> I would need a working factorial of 100! which currently seems
> unlikely. And my currently maximum is only 25!

Using Stirling's approximation will help here; see below code. However, you will likely run into other problems with the Hermite polynomials at large n.

> So, for your pleasure I am presenting you with a harmonic oscillator code :)
>
> Unfortunately without units :) But if you will try to convert this to
> units, I might learn even more about using this library. Especially
> because dimension of wavefunction in 1 dimensional space is
> 1/sqrt(meter), and I have no idea how to declare this
> derived_dimension, this didn't work:

Here's a fully-unitized code for the 1D harmonic oscillator:

#include <complex>
#include <iostream>
#include <limits>

using namespace std;

#include <boost/mpl/list.hpp>

// use math constants
#include <boost/math/constants/constants.hpp>

// use math factorial
#include <boost/math/special_functions/factorials.hpp>
#include <boost/math/special_functions/hermite.hpp>
#include <boost/typeof/std/complex.hpp>

#include <boost/units/cmath.hpp>

#include <boost/units/io.hpp>
#include <boost/units/pow.hpp>
#include <boost/units/systems/si.hpp>
#include <boost/units/systems/si/prefixes.hpp>
#include <boost/units/systems/si/codata/electron_constants.hpp>
#include <boost/units/systems/si/codata/physico-chemical_constants.hpp>
#include <boost/units/systems/si/codata/universal_constants.hpp>

using namespace boost::units;
using namespace boost::units::si;

// MCS - need codata namespace to find hbar
using namespace boost::units::si::constants::codata;

// MCS - derived_dimension only works with integer powers; use make_dimension_list here instead
typedef make_dimension_list<boost::mpl::list<dim<length_base_dimension,static_rational<-1,2> > > >::type quantum_wavefunction_1D_dimension;
typedef unit<quantum_wavefunction_1D_dimension,si::system> quantum_wavefunction_1D;

double factorial_to_minus_one_half(int n)
{
if (n <= 10)
{
return 1./sqrt(boost::math::factorial<double>(n));
}
else
{
// use Stirling's approximation for n>10
const double pi(boost::math::constants::pi<double>()),
logfac(log(2*pi*n)/2. + n*(log(n)-1));

return exp(-0.5*logfac);
}
}

// this function returns 1/sqrt(meter)
quantity<quantum_wavefunction_1D>
quantum_harmonic_oscillator_wavefunction_1d(int n,
const quantity<mass>& m,
const quantity<frequency>& omega,
const quantity<length>& x)
{
const double pi(boost::math::constants::pi<double>()),
prefac(pow(2.,-n/2.)*factorial_to_minus_one_half(n));

const quantity<energy_time> h_(hbar);

const double Hn(boost::math::hermite(n,root<2>(m*omega/h_)*x)),
gaussian(exp(-m*omega*pow<2>(x)/(2.*h_)));

return prefac*root<4>(m*omega/(pi*h_))*gaussian*Hn;
}

int main()
{
// MCS - change -1 -> +1 here otherwise loop doesn't run
quantity<length> s(15*nano*meter),
ds(0.1*nano*meter);

for (quantity<length> x=-s ; x<=s ; x+=ds)
{
const quantity<temperature> T(273.15*kelvin);
const quantity<frequency> omega(k_B*T/hbar.value());

std::cout << x.value() << "\t";

for (int n=0;n<25;++n)
{
const quantity<quantum_wavefunction_1D> psi(quantum_harmonic_oscillator_wavefunction_1d(n,m_e,omega,x));

std::cout << psi.value() << "\t";
}

std::cout << std::endl;
}

return 0;
}

and PDFs for n=0..25 :

Matthias