sysrap/tests/stranRotateTest.cc

0001 #include <iostream>
0002
0003 #include "scuda.h"
0004 #include "squad.h"
0005 #include "sqat4.h"
0006 #include "stran.h"
0007
0008 #include "SPath.hh"
0009 #include "NP.hh"
0010
0011
0012 #include <glm/gtc/random.hpp>
0013 #include <glm/gtx/component_wise.hpp>
0014
0015
0016 template<typename T>
0017 struct stranRotateTest
0018 {
0019     static const T MAX_DEVIATION ;
0020     static const char* NAME ;
0021
0022     static T make_rotate_a2b_(const glm::tvec3<T>& a_, const glm::tvec3<T>& b_, bool dump );
0023
0024     static void make_rotate_a2b_0(unsigned num);
0025     static void make_rotate_a2b_1();
0026     static void make_rotate_a2b_2();
0027
0028     static T make_rotate_a2b_3(unsigned num);
0029     static T make_rotate_a2b_4(unsigned num);
0030
0031     static void test_a2b();
0032
0033 };
0034
0035 template<>
0036 const float stranRotateTest<float>::MAX_DEVIATION = 1e-6f ;
0037
0038 template<>
0039 const double stranRotateTest<double>::MAX_DEVIATION = 1e-12f ;
0040
0041
0042 template<>
0043 const char* stranRotateTest<float>::NAME = "stranRotateTest<float>" ;
0044
0045 template<>
0046 const char* stranRotateTest<double>::NAME = "stranRotateTest<double>" ;
0047
0048
0049
0050
0051 template<typename T>
0052 T stranRotateTest<T>::make_rotate_a2b_(const glm::tvec3<T>& a_, const glm::tvec3<T>& b_, bool dump )
0053 {
0054     glm::tvec3<T> a = glm::normalize(a_);
0055     glm::tvec3<T> b = glm::normalize(b_);
0056     bool flip = true ;
0057
0058     const Tran<T>* tr = Tran<T>::make_rotate_a2b( a, b, flip );
0059
0060     T w = 0. ;
0061     glm::tvec4<T> a4(a.x,a.y,a.z,w);
0062
0063     glm::tvec4<T> b4 = tr->t * a4 ;
0064     glm::tvec3<T> b3(b4);
0065     T dmax = glm::compMax( glm::abs( b3 - b ) );
0066
0067     glm::tvec4<T> B4 = a4 * tr->t ;
0068     glm::tvec3<T> B3(B4);
0069     T Dmax = glm::compMax( glm::abs( B3 - b ) );
0070
0071     if(dump)
0072     {
0073         std::cout
0074             << "make_rotate_a2b_ "
0075             << std::endl
0076             << " a " << glm::to_string(a)
0077             << std::endl
0078             << " b " << glm::to_string(b)
0079             << std::endl
0080             << " b4 " << glm::to_string(b4)
0081             << std::endl
0082             << " b3 " << glm::to_string(b3)
0083             << std::endl
0084             << " B3 " << glm::to_string(B3)
0085             << std::endl
0086             << " dmax*1e6 " << std::setw(10) << std::fixed << std::setprecision(10) << dmax*1e6
0087             << std::endl
0088             << " Dmax*1e6 " << std::setw(10) << std::fixed << std::setprecision(10) << Dmax*1e6
0089             << std::endl
0090
0091
0092             ;
0093
0094         std::cout << *tr << std::endl ;
0095
0096         for(int i=0 ; i < 3 ; i++) std::cout
0097                 << std::setw(1) << i << " : "
0098                 << std::setw(10) << std::fixed << std::setprecision(10) << ( b[i] - b4[i] )
0099                 << std::endl
0100                 ;
0101     }
0102
0103
0104
0105     return dmax ;
0106 }
0107
0108
0109
0110 template<typename T>
0111 void stranRotateTest<T>::make_rotate_a2b_0(unsigned num)
0112 {
0113     glm::tvec3<T> px(1,0,0) ;
0114     glm::tvec3<T> py(0,1,0) ;
0115     glm::tvec3<T> pz(0,0,1) ;
0116
0117     glm::tvec3<T> nx(-1, 0, 0) ;
0118     glm::tvec3<T> ny( 0,-1, 0) ;
0119     glm::tvec3<T> nz( 0, 0,-1) ;
0120
0121     if(num > 0)  make_rotate_a2b_(px, py, true);
0122     if(num > 1 ) make_rotate_a2b_(pz, pz, true);
0123     if(num > 2 ) make_rotate_a2b_(pz, nz, true);
0124 }
0125
0126
0127 template<typename T>
0128 void stranRotateTest<T>::make_rotate_a2b_1()
0129 {
0130      glm::tvec3<T> a(1,0,0) ;
0131
0132      for(T p=0. ; p < 2. ; p+=0.1 )
0133      {
0134          T phi = glm::pi<T>()*p  ;
0135
0136          glm::tvec3<T> b(cos(phi),sin(phi),0) ;
0137
0138          std::cout << " a " << glm::to_string(a) << std::endl ;
0139          std::cout << " b " << glm::to_string(b) << std::endl ;
0140
0141          const Tran<T>* t = Tran<T>::make_rotate_a2b( a, b );
0142
0143          std::cout << *t << std::endl ;
0144      }
0145 }
0146
0147 template<typename T>
0148 void stranRotateTest<T>::make_rotate_a2b_2()
0149 {
0150     glm::tvec3<T> a( 1, 1,0) ;
0151     glm::tvec3<T> b(-1,-1,0) ;
0152     make_rotate_a2b_(a, b, true);
0153 }
0154
0155
0156
0157 /**
0158 stranRotateTest::make_rotate_a2b_3
0159 ------------------------
0160
0161 0. Generate two random unit vectors a,b
0162 1. compute a rotation matrix R that transforms a->b
0163 2. use R to transform a
0164 3. compare R a with b
0165
0166 **/
0167
0168 template<typename T>
0169 T stranRotateTest<T>::make_rotate_a2b_3(unsigned num)
0170 {
0171     bool dump = num < 3 ;
0172     //bool dump = false ;
0173     T dmaxx = 0. ;
0174
0175     for(unsigned i=0 ; i < num ; i++)
0176     {
0177         glm::tvec3<T> a = glm::sphericalRand(1.);
0178         glm::tvec3<T> b = glm::sphericalRand(1.);
0179
0180         T dmax = make_rotate_a2b_(a,b, dump);
0181         if( dmax > dmaxx ) dmaxx = dmax ;
0182
0183         if(dump) std::cout
0184             << " i " << std::setw(4) << i
0185             << " a " << std::setw(40) << glm::to_string(a)
0186             << " b " << std::setw(40) << glm::to_string(b)
0187             << " dmax*1e9 " << std::fixed << std::setw(10) << std::setprecision(6) << dmax*1e9
0188             << std::endl
0189             ;
0190     }
0191     std::cout
0192         << " test_make_rotate_a2b_3 "
0193         << " dmaxx " << std::scientific << dmaxx
0194         << " dmaxx*1e9 " << std::fixed << std::setw(10) << std::setprecision(6) << dmaxx*1e9
0195         << std::endl
0196         ;
0197     return dmaxx ;
0198 }
0199
0200
0201 /**
0202 make_rotate_a2b_4 : checking handling of parallel or anti-parallel
0203 --------------------------------------------------------------------------
0204
0205 0. Generate random unit vectors a and obtain b from it by negation
0206 1. compute a rotation matrix R that transforms a->b
0207 2. use R to transform a
0208 3. compare R a with b
0209
0210 **/
0211
0212 template<typename T>
0213 T stranRotateTest<T>::make_rotate_a2b_4(unsigned num )
0214 {
0215     bool dump = num < 3 ;
0216     //bool dump = false ;
0217
0218     T dmaxx = 0. ;
0219
0220     for(unsigned i=0 ; i < num ; i++)
0221     {
0222         glm::tvec3<T> a = glm::sphericalRand(1.);
0223         glm::tvec3<T> b = i % 2 == 0 ? a : -a ;
0224
0225         T dmax = make_rotate_a2b_(a,b, dump);
0226         if( dmax > dmaxx ) dmaxx = dmax ;
0227
0228         if(dump) std::cout
0229             << " i " << std::setw(4) << i
0230             << " a " << std::setw(40) << glm::to_string(a)
0231             << " b " << std::setw(40) << glm::to_string(b)
0232             << " dmax*1e9 " << std::fixed << std::setw(10) << std::setprecision(6) << dmax*1e9
0233             << std::endl
0234             ;
0235     }
0236
0237     std::cout
0238         << " test_make_rotate_a2b_4 "
0239         << " dmaxx " << std::scientific << dmaxx
0240         << " dmaxx*1e9 " << std::fixed << std::setw(10) << std::setprecision(6) << dmaxx*1e9
0241         << std::endl
0242         ;
0243
0244     return dmaxx ;
0245 }
0246
0247
0248 template<typename T>
0249 void stranRotateTest<T>::test_a2b()
0250 {
0251     std::cout << NAME << std::endl ;
0252     //test_make_rotate_a2b_0<T>(1u);
0253
0254     T d3 = make_rotate_a2b_3(1000);
0255     bool d3_expect = d3 < MAX_DEVIATION  ;
0256     assert( d3_expect );
0257     if(!d3_expect) std::raise(SIGINT);
0258
0259     T d4 = make_rotate_a2b_4(1000);
0260     bool d4_expect = d4 < MAX_DEVIATION  ;
0261     assert( d4_expect );
0262     if(!d4_expect) std::raise(SIGINT);
0263
0264 }
0265
0266
0267 int main()
0268 {
0269     stranRotateTest<double>::test_a2b();
0270     stranRotateTest<float>::test_a2b();
0271
0272     return 0 ;
0273 }
0274 // om- ; TEST=stranRotateTest om-t
0275
0276
0277