From 722bc71ed8d4e1ae5616c5c8284fbffe21c4ffa4 Mon Sep 17 00:00:00 2001
From: Jaron Kent-Dobias <jaron@kent-dobias.com>
Date: Tue, 17 Jul 2018 23:27:10 -0400
Subject: working generator of rotations near the current state

---
 lib/orthogonal.h | 22 +++++++++++++++++-----
 1 file changed, 17 insertions(+), 5 deletions(-)

(limited to 'lib/orthogonal.h')

diff --git a/lib/orthogonal.h b/lib/orthogonal.h
index 6166223..d8ad33d 100644
--- a/lib/orthogonal.h
+++ b/lib/orthogonal.h
@@ -167,9 +167,9 @@ orthogonal_t <q, double> generate_rotation_uniform (gsl_rng *r, const state_t <o
 
 template <q_t q>
 orthogonal_t <q, double> generate_rotation_perturbation (gsl_rng *r, const state_t <orthogonal_t <q, double>, vector_t <q, double>> *s, double epsilon) {
-  orthogonal_t <q, double> ptr;
+  orthogonal_t <q, double> m;
   vector_t <q, double> tmp_v;
-  ptr.is_reflection = true;
+  m.is_reflection = true;
 
   tmp_v.x = (double *)malloc(q * sizeof(double));
 
@@ -196,10 +196,22 @@ orthogonal_t <q, double> generate_rotation_perturbation (gsl_rng *r, const state
     tmp_v.x[i] /= mag_v;
   }
 
-  vector_t <q, double> v = act(s->R, tmp_v);
+  vector_t <q, double> v = act_inverse(s->R, tmp_v);
   free(tmp_v.x);
-  ptr.x = v.x;
+  m.x = v.x;
 
-  return ptr;
+  v2 = 0;
+
+  for (q_t i = 0; i < q; i++) {
+    v2 += m.x[i] * m.x[i];
+  }
+
+  mag_v = sqrt(v2);
+
+  for (q_t i = 0; i < q; i++) {
+    m.x[i] /= mag_v;
+  }
+
+  return m;
 }
 
-- 
cgit v1.2.3-70-g09d2