Re: [boost] [units] - learning to use, continued :)

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