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:10:22

> 
Here is a corrected version of your code with a few notes:

#include <iostream>
#include <complex>
#include <boost/foreach.hpp>
using namespace std;
// use rational numbers
#include <boost/rational.hpp>
// use math constants
#include <boost/math/constants/constants.hpp>
// use math factorial
#include <boost/math/special_functions/factorials.hpp>
//using namespace boost::math

// use units
#include <boost/typeof/std/complex.hpp>
 
// MCS - need to include this header for make_dimension_list
#include <boost/mpl/list.hpp>
 
#include <boost/units/cmath.hpp>
 
// MCS - need to include this header
#include <boost/units/pow.hpp>

#include <boost/units/io.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;
//use units end

vector<boost::rational<signed long> > hermite_polynomial_coefficients(int order, boost::rational<signed long> lambda_per_alpha)
{
    vector<boost::rational<signed long> > c_even; // can't use static, because lambda_per_alpha changes between calls
    vector<boost::rational<signed long> > c_odd;
    if(c_even.size() == 0)
    {
        c_even.push_back(1);
        c_even.push_back(0);
    }
    if(c_odd.size() == 0)
    {
        c_odd.push_back(0);
        c_odd.push_back(1);
    }
    vector<boost::rational<signed long> >& a( (order % 2 == 0) ?c_even:c_odd );
    
    if(order+1 < a.size())
        return a;

    for(int i = a.size() ; i <= order ; ++i)
    {
        boost::rational<signed long> next = a[i-2]*(2*(i-2)+1-lambda_per_alpha)/(i*(i-1));
        a.push_back(next);
    }
    
//    BOOST_FOREACH(boost::rational<signed long> v,a)
//        std::cout << v << " ";
//    std::cout << "\n";
    
    return a;
}

vector<boost::rational<signed long> > hermite_polynomial_scaled(int order, boost::rational<signed long> lambda_per_alpha)
{
    vector<boost::rational<signed long> > a(hermite_polynomial_coefficients(order,lambda_per_alpha));
    vector<boost::rational<signed long> > b;
    boost::rational<signed long> factor(std::pow(2.,order)/a[order]);
    for(int i = 0 ; i<=order ; ++i)
        b.push_back(a[i]*factor);
    return b;
}

// 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;

 // this function returns 1/sqrt(meter) 
quantity<quantum_wavefunction_1D,complex<double> > 
quantum_oscillator_wavefunction(
       int n					// n - order of wavefunction
     , const quantity<length>		x       // position
     , const quantity<mass>		mass	// mass - mass of particle
     , const quantity<frequency>	omega	// oscillator frequency
 )
{
        // MCS - need to define unit typedef in addition to dimension
        typedef derived_dimension<length_base_dimension,-2>::type	reciprocal_area_dimension;
        typedef unit<reciprocal_area_dimension,si::system>		reciprocal_area;

        quantity<energy_time>		h_(hbar);

        quantity<energy>		E(h_*omega*(n+0.5));			// energy of oscillator
        quantity<reciprocal_area>	lambda(2.0*mass*E/(h_*h_)),
                                        alpha(mass*omega/h_);
        quantity<wavenumber>		alpha_root(root<2>(alpha));
        double				lambda_per_alpha(2*n+1);

        const double	pi(boost::math::constants::pi<double>()),
                        fac(boost::math::factorial<double>(n));

        // N_n should be 1/sqrt(meter) because it has a dimension of quantum_wavefunction_1D
        const quantity<quantum_wavefunction_1D,complex<double> >	N_n(root<2>(root<2>(alpha/pi)/(std::pow(2.,n)*fac)));

        const complex<double> i(0,1);

        quantity<quantum_wavefunction_1D,complex<double> > result(0/root<2>(meter));

        vector<boost::rational<signed long> > a(hermite_polynomial_scaled(n,lambda_per_alpha));

        for(int v=0 ; v<=n ; v++)
        {
                //         result += N_n*a[v]*pow<v>(alpha_root*x)*exp(-0.5*alpha*pow<2>(x));
                // MCS - a[v] is boost::rational -> need to convert to double via rational_cast
                // MCS - pow<v> does not work because v must be known at compile time
                const double		a_v(boost::rational_cast<double>(a[v])),
                                        tmp1(a_v*pow(alpha_root*x,v));
                const complex<double>	tmp2(exp(-0.5*alpha*pow<2>(x))),
                                        tmp3(tmp1*tmp2);

                result += tmp3*N_n;
        }

        return result;
}

int main()
{
        // MCS - change -1 -> +1 here otherwise loop doesn't run
        quantity<length> s(15*nano*meter),
                                         ds(0.01*nano*meter);
        
        std::cout << s << "\t" << ds << std::endl;
        
        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 << " " << quantum_oscillator_wavefunction(25,x,m_e,omega) << "\n";
        }

        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());
                
                const quantity<quantum_wavefunction_1D,complex<double> >	psi(quantum_oscillator_wavefunction(25,x,m_e,omega));
                
                std::cout << x.value() << " " << psi.value().real() << "\t" << psi.value().imag() << "\n";
        }

    return 0;
}

and the n=25 wavefunction:

Cheers,

Matthias