hso3so4.cpp

C++的一个好库。。。现在很流行

// test file for HSO3.hpp and HSO4.hpp

//  (C) Copyright Hubert Holin 2001.
//  Distributed under the Boost Software License, Version 1.0. (See
//  accompanying file LICENSE_1_0.txt or copy at
//  http://www.boost.org/LICENSE_1_0.txt)


#include <iostream>

#include <boost/math/quaternion.hpp>

#include "HSO3.hpp"
#include "HSO4.hpp"


const int    number_of_intervals = 5;

const float    pi = ::std::atan(1.0f)*4;



void    test_SO3();
    
void    test_SO4();


int    main()

{
    test_SO3();
    
    test_SO4();
    
    ::std::cout << "That's all folks!" << ::std::endl;
}


//
//    Test of quaternion and R^3 rotation relationship
//

void    test_SO3_spherical()
{
    ::std::cout << "Testing spherical:" << ::std::endl;
    ::std::cout << ::std::endl;
    
    const float    rho = 1.0f;
    
    float        theta;
    float        phi1;
    float        phi2;
    
    for    (int idxphi2 = 0; idxphi2 <= number_of_intervals; idxphi2++)
    {
        phi2 = (-pi/2)+(idxphi2*pi)/number_of_intervals;
        
        for    (int idxphi1 = 0; idxphi1 <= number_of_intervals; idxphi1++)
        {
            phi1 = (-pi/2)+(idxphi1*pi)/number_of_intervals;
            
            for    (int idxtheta = 0; idxtheta <= number_of_intervals; idxtheta++)
            {
                theta = -pi+(idxtheta*(2*pi))/number_of_intervals;
                
                //::std::cout << "theta = " << theta << " ; ";
                //::std::cout << "phi1 = " << phi1 << " ; ";
                //::std::cout << "phi2 = " << phi2;
                //::std::cout << ::std::endl;
                
                ::boost::math::quaternion<float>    q = ::boost::math::spherical(rho, theta, phi1, phi2);
                
                //::std::cout << "q = " << q << ::std::endl;
                
                R3_matrix<float>                    rot = quaternion_to_R3_rotation(q);
                
                //::std::cout << "rot = ";
                //::std::cout << "\t" << rot.a11 << "\t" << rot.a12 << "\t" << rot.a13 << ::std::endl;
                //::std::cout << "\t";
                //::std::cout << "\t" << rot.a21 << "\t" << rot.a22 << "\t" << rot.a23 << ::std::endl;
                //::std::cout << "\t";
                //::std::cout << "\t" << rot.a31 << "\t" << rot.a32 << "\t" << rot.a33 << ::std::endl;
                
                ::boost::math::quaternion<float>    p = R3_rotation_to_quaternion(rot, &q);
                
                //::std::cout << "p = " << p << ::std::endl;
                
                //::std::cout << "round trip discrepancy: " << ::boost::math::abs(q-p) << ::std::endl;
                
                //::std::cout << ::std::endl;
            }
        }
    }
    
    ::std::cout << ::std::endl;
}

    
void    test_SO3_semipolar()
{
    ::std::cout << "Testing semipolar:" << ::std::endl;
    ::std::cout << ::std::endl;
    
    const float    rho = 1.0f;
    
    float        alpha;
    float        theta1;
    float        theta2;
    
    for    (int idxalpha = 0; idxalpha <= number_of_intervals; idxalpha++)
    {
        alpha = (idxalpha*(pi/2))/number_of_intervals;
        
        for    (int idxtheta1 = 0; idxtheta1 <= number_of_intervals; idxtheta1++)
        {
            theta1 = -pi+(idxtheta1*(2*pi))/number_of_intervals;
            
            for    (int idxtheta2 = 0; idxtheta2 <= number_of_intervals; idxtheta2++)
            {
                theta2 = -pi+(idxtheta2*(2*pi))/number_of_intervals;
                
                //::std::cout << "alpha = " << alpha << " ; ";
                //::std::cout << "theta1 = " << theta1 << " ; ";
                //::std::cout << "theta2 = " << theta2;
                //::std::cout << ::std::endl;
                
                ::boost::math::quaternion<float>    q = ::boost::math::semipolar(rho, alpha, theta1, theta2);
                
                //::std::cout << "q = " << q << ::std::endl;
                
                R3_matrix<float>                    rot = quaternion_to_R3_rotation(q);
                
                //::std::cout << "rot = ";
                //::std::cout << "\t" << rot.a11 << "\t" << rot.a12 << "\t" << rot.a13 << ::std::endl;
                //::std::cout << "\t";
                //::std::cout << "\t" << rot.a21 << "\t" << rot.a22 << "\t" << rot.a23 << ::std::endl;
                //::std::cout << "\t";
                //::std::cout << "\t" << rot.a31 << "\t" << rot.a32 << "\t" << rot.a33 << ::std::endl;
                
                ::boost::math::quaternion<float>    p = R3_rotation_to_quaternion(rot, &q);
                
                //::std::cout << "p = " << p << ::std::endl;
                
                //::std::cout << "round trip discrepancy: " << ::boost::math::abs(q-p) << ::std::endl;
                
                //::std::cout << ::std::endl;
            }
        }
    }
    
    ::std::cout << ::std::endl;
}

    
void    test_SO3_multipolar()
{
    ::std::cout << "Testing multipolar:" << ::std::endl;
    ::std::cout << ::std::endl;
    
    float    rho1;
    float    rho2;
    
    float    theta1;
    float    theta2;
    
    for    (int idxrho = 0; idxrho <= number_of_intervals; idxrho++)
    {
        rho1 = (idxrho*1.0f)/number_of_intervals;
        rho2 = ::std::sqrt(1.0f-rho1*rho1);
        
        for    (int idxtheta1 = 0; idxtheta1 <= number_of_intervals; idxtheta1++)
        {
            theta1 = -pi+(idxtheta1*(2*pi))/number_of_intervals;
            
            for    (int idxtheta2 = 0; idxtheta2 <= number_of_intervals; idxtheta2++)
            {
                theta2 = -pi+(idxtheta2*(2*pi))/number_of_intervals;
                
                //::std::cout << "rho1 = " << rho1 << " ; ";
                //::std::cout << "theta1 = " << theta1 << " ; ";
                //::std::cout << "theta2 = " << theta2;
                //::std::cout << ::std::endl;
                
                ::boost::math::quaternion<float>    q = ::boost::math::multipolar(rho1, theta1, rho2, theta2);
                
                //::std::cout << "q = " << q << ::std::endl;
                
                R3_matrix<float>                    rot = quaternion_to_R3_rotation(q);
                
                //::std::cout << "rot = ";
                //::std::cout << "\t" << rot.a11 << "\t" << rot.a12 << "\t" << rot.a13 << ::std::endl;
                //::std::cout << "\t";
                //::std::cout << "\t" << rot.a21 << "\t" << rot.a22 << "\t" << rot.a23 << ::std::endl;
                //::std::cout << "\t";
                //::std::cout << "\t" << rot.a31 << "\t" << rot.a32 << "\t" << rot.a33 << ::std::endl;
                
                ::boost::math::quaternion<float>    p = R3_rotation_to_quaternion(rot, &q);
                
                //::std::cout << "p = " << p << ::std::endl;
                
                //::std::cout << "round trip discrepancy: " << ::boost::math::abs(q-p) << ::std::endl;
                
                //::std::cout << ::std::endl;
            }
        }
    }
    
    ::std::cout << ::std::endl;
}

    
void    test_SO3_cylindrospherical()
{
    ::std::cout << "Testing cylindrospherical:" << ::std::endl;
    ::std::cout << ::std::endl;
    
    float    t;
    
    float    radius;
    float    longitude;
    float    latitude;
    
    for    (int idxt = 0; idxt <= number_of_intervals; idxt++)