Create a gist now

Instantly share code, notes, and snippets.

@wchristian /marp.pl Secret
Last active Dec 21, 2015

What would you like to do?
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;
},
);
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment