ceph/src/boost/libs/math/example/HSO3SO4.cpp

   1 // test file for HSO3.hpp and HSO4.hpp
   2
   3 //  (C) Copyright Hubert Holin 2001.
   4 //  Distributed under the Boost Software License, Version 1.0. (See
   5 //  accompanying file LICENSE_1_0.txt or copy at
   6 //  http://www.boost.org/LICENSE_1_0.txt)
   7
   8
   9 #include <iostream>
  10
  11 #include <boost/math/quaternion.hpp>
  12
  13 #include "HSO3.hpp"
  14 #include "HSO4.hpp"
  15
  16
  17 const int    number_of_intervals = 5;
  18
  19 const float    pi = ::std::atan(1.0f)*4;
  20
  21
  22
  23 void    test_SO3();
  24
  25 void    test_SO4();
  26
  27
  28 int    main()
  29
  30 {
  31     test_SO3();
  32
  33     test_SO4();
  34
  35     ::std::cout << "That's all folks!" << ::std::endl;
  36 }
  37
  38
  39 //
  40 //    Test of quaternion and R^3 rotation relationship
  41 //
  42
  43 void    test_SO3_spherical()
  44 {
  45     ::std::cout << "Testing spherical:" << ::std::endl;
  46     ::std::cout << ::std::endl;
  47
  48     const float    rho = 1.0f;
  49
  50     float        theta;
  51     float        phi1;
  52     float        phi2;
  53
  54     for    (int idxphi2 = 0; idxphi2 <= number_of_intervals; idxphi2++)
  55     {
  56         phi2 = (-pi/2)+(idxphi2*pi)/number_of_intervals;
  57
  58         for    (int idxphi1 = 0; idxphi1 <= number_of_intervals; idxphi1++)
  59         {
  60             phi1 = (-pi/2)+(idxphi1*pi)/number_of_intervals;
  61
  62             for    (int idxtheta = 0; idxtheta <= number_of_intervals; idxtheta++)
  63             {
  64                 theta = -pi+(idxtheta*(2*pi))/number_of_intervals;
  65
  66                 //::std::cout << "theta = " << theta << " ; ";
  67                 //::std::cout << "phi1 = " << phi1 << " ; ";
  68                 //::std::cout << "phi2 = " << phi2;
  69                 //::std::cout << ::std::endl;
  70
  71                 ::boost::math::quaternion<float>    q = ::boost::math::spherical(rho, theta, phi1, phi2);
  72
  73                 //::std::cout << "q = " << q << ::std::endl;
  74
  75                 R3_matrix<float>                    rot = quaternion_to_R3_rotation(q);
  76
  77                 //::std::cout << "rot = ";
  78                 //::std::cout << "\t" << rot.a11 << "\t" << rot.a12 << "\t" << rot.a13 << ::std::endl;
  79                 //::std::cout << "\t";
  80                 //::std::cout << "\t" << rot.a21 << "\t" << rot.a22 << "\t" << rot.a23 << ::std::endl;
  81                 //::std::cout << "\t";
  82                 //::std::cout << "\t" << rot.a31 << "\t" << rot.a32 << "\t" << rot.a33 << ::std::endl;
  83
  84                 ::boost::math::quaternion<float>    p = R3_rotation_to_quaternion(rot, &q);
  85
  86                 //::std::cout << "p = " << p << ::std::endl;
  87
  88                 //::std::cout << "round trip discrepancy: " << ::boost::math::abs(q-p) << ::std::endl;
  89
  90                 //::std::cout << ::std::endl;
  91             }
  92         }
  93     }
  94
  95     ::std::cout << ::std::endl;
  96 }
  97
  98
  99 void    test_SO3_semipolar()
 100 {
 101     ::std::cout << "Testing semipolar:" << ::std::endl;
 102     ::std::cout << ::std::endl;
 103
 104     const float    rho = 1.0f;
 105
 106     float        alpha;
 107     float        theta1;
 108     float        theta2;
 109
 110     for    (int idxalpha = 0; idxalpha <= number_of_intervals; idxalpha++)
 111     {
 112         alpha = (idxalpha*(pi/2))/number_of_intervals;
 113
 114         for    (int idxtheta1 = 0; idxtheta1 <= number_of_intervals; idxtheta1++)
 115         {
 116             theta1 = -pi+(idxtheta1*(2*pi))/number_of_intervals;
 117
 118             for    (int idxtheta2 = 0; idxtheta2 <= number_of_intervals; idxtheta2++)
 119             {
 120                 theta2 = -pi+(idxtheta2*(2*pi))/number_of_intervals;
 121
 122                 //::std::cout << "alpha = " << alpha << " ; ";
 123                 //::std::cout << "theta1 = " << theta1 << " ; ";
 124                 //::std::cout << "theta2 = " << theta2;
 125                 //::std::cout << ::std::endl;
 126
 127                 ::boost::math::quaternion<float>    q = ::boost::math::semipolar(rho, alpha, theta1, theta2);
 128
 129                 //::std::cout << "q = " << q << ::std::endl;
 130
 131                 R3_matrix<float>                    rot = quaternion_to_R3_rotation(q);
 132
 133                 //::std::cout << "rot = ";
 134                 //::std::cout << "\t" << rot.a11 << "\t" << rot.a12 << "\t" << rot.a13 << ::std::endl;
 135                 //::std::cout << "\t";
 136                 //::std::cout << "\t" << rot.a21 << "\t" << rot.a22 << "\t" << rot.a23 << ::std::endl;
 137                 //::std::cout << "\t";
 138                 //::std::cout << "\t" << rot.a31 << "\t" << rot.a32 << "\t" << rot.a33 << ::std::endl;
 139
 140                 ::boost::math::quaternion<float>    p = R3_rotation_to_quaternion(rot, &q);
 141
 142                 //::std::cout << "p = " << p << ::std::endl;
 143
 144                 //::std::cout << "round trip discrepancy: " << ::boost::math::abs(q-p) << ::std::endl;
 145
 146                 //::std::cout << ::std::endl;
 147             }
 148         }
 149     }
 150
 151     ::std::cout << ::std::endl;
 152 }
 153
 154
 155 void    test_SO3_multipolar()
 156 {
 157     ::std::cout << "Testing multipolar:" << ::std::endl;
 158     ::std::cout << ::std::endl;
 159
 160     float    rho1;
 161     float    rho2;
 162
 163     float    theta1;
 164     float    theta2;
 165
 166     for    (int idxrho = 0; idxrho <= number_of_intervals; idxrho++)
 167     {
 168         rho1 = (idxrho*1.0f)/number_of_intervals;
 169         rho2 = ::std::sqrt(1.0f-rho1*rho1);
 170
 171         for    (int idxtheta1 = 0; idxtheta1 <= number_of_intervals; idxtheta1++)
 172         {
 173             theta1 = -pi+(idxtheta1*(2*pi))/number_of_intervals;
 174
 175             for    (int idxtheta2 = 0; idxtheta2 <= number_of_intervals; idxtheta2++)
 176             {
 177                 theta2 = -pi+(idxtheta2*(2*pi))/number_of_intervals;
 178
 179                 //::std::cout << "rho1 = " << rho1 << " ; ";
 180                 //::std::cout << "theta1 = " << theta1 << " ; ";
 181                 //::std::cout << "theta2 = " << theta2;
 182                 //::std::cout << ::std::endl;
 183
 184                 ::boost::math::quaternion<float>    q = ::boost::math::multipolar(rho1, theta1, rho2, theta2);
 185
 186                 //::std::cout << "q = " << q << ::std::endl;
 187
 188                 R3_matrix<float>                    rot = quaternion_to_R3_rotation(q);
 189
 190                 //::std::cout << "rot = ";
 191                 //::std::cout << "\t" << rot.a11 << "\t" << rot.a12 << "\t" << rot.a13 << ::std::endl;
 192                 //::std::cout << "\t";
 193                 //::std::cout << "\t" << rot.a21 << "\t" << rot.a22 << "\t" << rot.a23 << ::std::endl;
 194                 //::std::cout << "\t";
 195                 //::std::cout << "\t" << rot.a31 << "\t" << rot.a32 << "\t" << rot.a33 << ::std::endl;
 196
 197                 ::boost::math::quaternion<float>    p = R3_rotation_to_quaternion(rot, &q);
 198
 199                 //::std::cout << "p = " << p << ::std::endl;
 200
 201                 //::std::cout << "round trip discrepancy: " << ::boost::math::abs(q-p) << ::std::endl;
 202
 203                 //::std::cout << ::std::endl;
 204             }
 205         }
 206     }
 207
 208     ::std::cout << ::std::endl;
 209 }
 210
 211
 212 void    test_SO3_cylindrospherical()
 213 {
 214     ::std::cout << "Testing cylindrospherical:" << ::std::endl;
 215     ::std::cout << ::std::endl;
 216
 217     float    t;
 218
 219     float    radius;
 220     float    longitude;
 221     float    latitude;
 222
 223     for    (int idxt = 0; idxt <= number_of_intervals; idxt++)
 224     {
 225         t = -1.0f+(idxt*2.0f)/number_of_intervals;
 226         radius = ::std::sqrt(1.0f-t*t);
 227
 228         for    (int idxlatitude = 0; idxlatitude <= number_of_intervals; idxlatitude++)
 229         {
 230             latitude = (-pi/2)+(idxlatitude*pi)/number_of_intervals;
 231
 232             for    (int idxlongitude = 0; idxlongitude <= number_of_intervals; idxlongitude++)
 233             {
 234                 longitude = -pi+(idxlongitude*(2*pi))/number_of_intervals;
 235
 236                 //::std::cout << "t = " << t << " ; ";
 237                 //::std::cout << "longitude = " << longitude;
 238                 //::std::cout << "latitude = " << latitude;
 239                 //::std::cout << ::std::endl;
 240
 241                 ::boost::math::quaternion<float>    q = ::boost::math::cylindrospherical(t, radius, longitude, latitude);
 242
 243                 //::std::cout << "q = " << q << ::std::endl;
 244
 245                 R3_matrix<float>                    rot = quaternion_to_R3_rotation(q);
 246
 247                 //::std::cout << "rot = ";
 248                 //::std::cout << "\t" << rot.a11 << "\t" << rot.a12 << "\t" << rot.a13 << ::std::endl;
 249                 //::std::cout << "\t";
 250                 //::std::cout << "\t" << rot.a21 << "\t" << rot.a22 << "\t" << rot.a23 << ::std::endl;
 251                 //::std::cout << "\t";
 252                 //::std::cout << "\t" << rot.a31 << "\t" << rot.a32 << "\t" << rot.a33 << ::std::endl;
 253
 254                 ::boost::math::quaternion<float>    p = R3_rotation_to_quaternion(rot, &q);
 255
 256                 //::std::cout << "p = " << p << ::std::endl;
 257
 258                 //::std::cout << "round trip discrepancy: " << ::boost::math::abs(q-p) << ::std::endl;
 259
 260                 //::std::cout << ::std::endl;
 261             }
 262         }
 263     }
 264
 265     ::std::cout << ::std::endl;
 266 }
 267
 268
 269 void    test_SO3_cylindrical()
 270 {
 271     ::std::cout << "Testing cylindrical:" << ::std::endl;
 272     ::std::cout << ::std::endl;
 273
 274     float    r;
 275     float    angle;
 276
 277     float    h1;
 278     float    h2;
 279
 280     for    (int idxh2 = 0; idxh2 <= number_of_intervals; idxh2++)
 281     {
 282         h2 = -1.0f+(idxh2*2.0f)/number_of_intervals;
 283
 284         for    (int idxh1 = 0; idxh1 <= number_of_intervals; idxh1++)
 285         {
 286             h1 = ::std::sqrt(1.0f-h2*h2)*(-1.0f+(idxh2*2.0f)/number_of_intervals);
 287             r = ::std::sqrt(1.0f-h1*h1-h2*h2);
 288
 289             for    (int idxangle = 0; idxangle <= number_of_intervals; idxangle++)
 290             {
 291                 angle = -pi+(idxangle*(2*pi))/number_of_intervals;
 292
 293                 //::std::cout << "angle = " << angle << " ; ";
 294                 //::std::cout << "h1 = " << h1;
 295                 //::std::cout << "h2 = " << h2;
 296                 //::std::cout << ::std::endl;
 297
 298                 ::boost::math::quaternion<float>    q = ::boost::math::cylindrical(r, angle, h1, h2);
 299
 300                 //::std::cout << "q = " << q << ::std::endl;
 301
 302                 R3_matrix<float>                    rot = quaternion_to_R3_rotation(q);
 303
 304                 //::std::cout << "rot = ";
 305                 //::std::cout << "\t" << rot.a11 << "\t" << rot.a12 << "\t" << rot.a13 << ::std::endl;
 306                 //::std::cout << "\t";
 307                 //::std::cout << "\t" << rot.a21 << "\t" << rot.a22 << "\t" << rot.a23 << ::std::endl;
 308                 //::std::cout << "\t";
 309                 //::std::cout << "\t" << rot.a31 << "\t" << rot.a32 << "\t" << rot.a33 << ::std::endl;
 310
 311                 ::boost::math::quaternion<float>    p = R3_rotation_to_quaternion(rot, &q);
 312
 313                 //::std::cout << "p = " << p << ::std::endl;
 314
 315                 //::std::cout << "round trip discrepancy: " << ::boost::math::abs(q-p) << ::std::endl;
 316
 317                 //::std::cout << ::std::endl;
 318             }
 319         }
 320     }
 321
 322     ::std::cout << ::std::endl;
 323 }
 324
 325
 326 void    test_SO3()
 327 {
 328     ::std::cout << "Testing SO3:" << ::std::endl;
 329     ::std::cout << ::std::endl;
 330
 331     test_SO3_spherical();
 332
 333     test_SO3_semipolar();
 334
 335     test_SO3_multipolar();
 336
 337     test_SO3_cylindrospherical();
 338
 339     test_SO3_cylindrical();
 340 }
 341
 342
 343 //
 344 //    Test of quaternion and R^4 rotation relationship
 345 //
 346
 347 void    test_SO4_spherical()
 348 {
 349     ::std::cout << "Testing spherical:" << ::std::endl;
 350     ::std::cout << ::std::endl;
 351
 352     const float    rho1 = 1.0f;
 353     const float    rho2 = 1.0f;
 354
 355     float        theta1;
 356     float        phi11;
 357     float        phi21;
 358
 359     float        theta2;
 360     float        phi12;
 361     float        phi22;
 362
 363     for    (int idxphi21 = 0; idxphi21 <= number_of_intervals; idxphi21++)
 364     {
 365         phi21 = (-pi/2)+(idxphi21*pi)/number_of_intervals;
 366
 367         for    (int idxphi22 = 0; idxphi22 <= number_of_intervals; idxphi22++)
 368         {
 369             phi22 = (-pi/2)+(idxphi22*pi)/number_of_intervals;
 370
 371             for    (int idxphi11 = 0; idxphi11 <= number_of_intervals; idxphi11++)
 372             {
 373                 phi11 = (-pi/2)+(idxphi11*pi)/number_of_intervals;
 374
 375                 for    (int idxphi12 = 0; idxphi12 <= number_of_intervals; idxphi12++)
 376                 {
 377                     phi12 = (-pi/2)+(idxphi12*pi)/number_of_intervals;
 378
 379                     for    (int idxtheta1 = 0; idxtheta1 <= number_of_intervals; idxtheta1++)
 380                     {
 381                         theta1 = -pi+(idxtheta1*(2*pi))/number_of_intervals;
 382
 383                         for    (int idxtheta2 = 0; idxtheta2 <= number_of_intervals; idxtheta2++)
 384                         {
 385                             theta2 = -pi+(idxtheta2*(2*pi))/number_of_intervals;
 386
 387                             //::std::cout << "theta1 = " << theta1 << " ; ";
 388                             //::std::cout << "phi11 = " << phi11 << " ; ";
 389                             //::std::cout << "phi21 = " << phi21;
 390                             //::std::cout << "theta2 = " << theta2 << " ; ";
 391                             //::std::cout << "phi12 = " << phi12 << " ; ";
 392                             //::std::cout << "phi22 = " << phi22;
 393                             //::std::cout << ::std::endl;
 394
 395                             ::boost::math::quaternion<float>    p1 = ::boost::math::spherical(rho1, theta1, phi11, phi21);
 396
 397                             //::std::cout << "p1 = " << p1 << ::std::endl;
 398
 399                             ::boost::math::quaternion<float>    q1 = ::boost::math::spherical(rho2, theta2, phi12, phi22);
 400
 401                             //::std::cout << "q1 = " << q1 << ::std::endl;
 402
 403                             ::std::pair< ::boost::math::quaternion<float> , ::boost::math::quaternion<float> >    pq1 =
 404                                 ::std::make_pair(p1,q1);
 405
 406                             R4_matrix<float>                    rot = quaternions_to_R4_rotation(pq1);
 407
 408                             //::std::cout << "rot = ";
 409                             //::std::cout << "\t" << rot.a11 << "\t" << rot.a12 << "\t" << rot.a13 << "\t" << rot.a14 << ::std::endl;
 410                             //::std::cout << "\t";
 411                             //::std::cout << "\t" << rot.a21 << "\t" << rot.a22 << "\t" << rot.a23 << "\t" << rot.a24 << ::std::endl;
 412                             //::std::cout << "\t";
 413                             //::std::cout << "\t" << rot.a31 << "\t" << rot.a32 << "\t" << rot.a33 << "\t" << rot.a34 << ::std::endl;
 414                             //::std::cout << "\t";
 415                             //::std::cout << "\t" << rot.a41 << "\t" << rot.a42 << "\t" << rot.a43 << "\t" << rot.a44 << ::std::endl;
 416
 417                             ::std::pair< ::boost::math::quaternion<float> , ::boost::math::quaternion<float> >    pq2 =
 418                                 R4_rotation_to_quaternions(rot, &pq1);
 419
 420                             //::std::cout << "p1 = " << pq.first << ::std::endl;
 421                             //::std::cout << "p2 = " << pq.second << ::std::endl;
 422
 423                             //::std::cout << "round trip discrepancy: " << ::std::sqrt(::boost::math::norm(pq1.first-pq2.first)+::boost::math::norm(pq1.second-pq2.second)) << ::std::endl;
 424
 425                             //::std::cout << ::std::endl;
 426                         }
 427                     }
 428                 }
 429             }
 430         }
 431     }
 432
 433     ::std::cout << ::std::endl;
 434 }
 435
 436
 437 void    test_SO4()
 438 {
 439     ::std::cout << "Testing SO4:" << ::std::endl;
 440     ::std::cout << ::std::endl;
 441
 442     test_SO4_spherical();
 443 }
 444
 445