Instantly share code, notes, and snippets.

# wchristian/marp.plSecret Last active Dec 21, 2015

 use strict; use warnings; use 5.010; use PDL; # must be called before (!) 'use Inline Pdlpp' calls use IO::All -binary; use Inline 'Pdlpp' => io( "pdlpp.pm" )->all; use Benchmark 'timethis'; my \$v = pdl( 1, 0, 0 ); \$v->rotate( 180, pdl( 0, 1, 0 ) ); say \$v; my \$vec = pdl( 1, 0, 0 ); my \$axis = pdl( 0, 1, 0 ); timethis( 40000, sub { \$vec->rotate( 180, \$axis ) } );
 use strict; pp_def( 'rotate', Pars => 'vec(n); angle(); axis(n)', Code => q{ double PI = 3.141592653589793; double SinHalfAngle = sin( \$angle() / 2.0 / 180.0 * PI ); double CosHalfAngle = cos( \$angle() / 2.0 / 180.0 * PI ); // rotation quaternion double Rx = \$axis(n=>0) * SinHalfAngle; double Ry = \$axis(n=>1) * SinHalfAngle; double Rz = \$axis(n=>2) * SinHalfAngle; double Rw = CosHalfAngle; double Vx = \$vec(n=>0); double Vy = \$vec(n=>1); double Vz = \$vec(n=>2); double Vw = 0; // first step double x1 = ( Rx * Vw ) + ( Rw * Vx ) + ( Ry * Vz ) - ( Rz * Vy ); double y1 = ( Ry * Vw ) + ( Rw * Vy ) + ( Rz * Vx ) - ( Rx * Vz ); double z1 = ( Rz * Vw ) + ( Rw * Vz ) + ( Rx * Vy ) - ( Ry * Vx ); double w1 = ( Rw * Vw ) - ( Rx * Vx ) - ( Ry * Vy ) - ( Rz * Vz ); // conjugated quaternion double Qx = Rx * -1; double Qy = Ry * -1; double Qz = Rz * -1; double Qw = CosHalfAngle; // second step double x2 = ( x1 * Qw ) + ( w1 * Qx ) + ( y1 * Qz ) - ( z1 * Qy ); double y2 = ( y1 * Qw ) + ( w1 * Qy ) + ( z1 * Qx ) - ( x1 * Qz ); double z2 = ( z1 * Qw ) + ( w1 * Qz ) + ( x1 * Qy ) - ( y1 * Qx ); double w2 = ( w1 * Qw ) - ( x1 * Qx ) - ( y1 * Qy ) - ( z1 * Qz ); \$vec(n=>0) = x2; \$vec(n=>1) = y2; \$vec(n=>2) = z2; }, );