Adjusted code so asteroids should correctly follow geodesics.

src/lib.rs

@@ -30,7 +30,7 @@ use cga::CGA;
 
 const PI: f64 = f64::consts::PI;
 const TWO_PI: f64 = 2.*PI;
-const ERROR: f64 = 0.0001;
+const ERROR: f64 = 0.000001;
 
 const WIDTH: f64 = 600.;
 const HEIGHT: f64 = 600.;
@@ -47,9 +47,9 @@ const BLUE: &str       = "#0000ff";
 const SCALAR: usize = 0; 
 const E1: usize = 1;
 const E2: usize = 2;
-//const E3: usize = 3;
+const E3: usize = 3;
-//const E4: usize = 4;
+const E4: usize = 4;
-//const E5: usize = 5;
+const E5: usize = 5;
 const E12: usize = 6;
 const E13: usize = 7;
 const E14: usize = 8;
@@ -60,46 +60,42 @@ const E25: usize = 12;
 const E34: usize = 13;
 const E35: usize = 14;
 const E45: usize = 15;
-//const E123: usize = 16;
+const E123: usize = 16;
-//const E124: usize = 17;
+const E124: usize = 17;
-//const E125: usize = 18;
+const E125: usize = 18;
-//const E134: usize = 19;
+const E134: usize = 19;
-//const E135: usize = 20;
+const E135: usize = 20;
-//const E145: usize = 21;
+const E145: usize = 21;
-//const E234: usize = 22;
+const E234: usize = 22;
-//const E235: usize = 23;
+const E235: usize = 23;
-//const E245: usize = 24;
+const E245: usize = 24;
-//const E345: usize = 25;
+const E345: usize = 25;
-//const E1234: usize = 26;
+const E1234: usize = 26;
-//const E1235: usize = 27;
+const E1235: usize = 27;
-//const E1245: usize = 28;
+const E1245: usize = 28;
-//const E1345: usize = 29;
+const E1345: usize = 29;
-//const E2345: usize = 30;
+const E2345: usize = 30;
-//const E12345: usize = 31;
+const E12345: usize = 31;
 
 const DRAG: f64 = 1.;
 const ACCEL_STR: f64 = 1.;
 
 struct Ship {
     com_rotor: CGA,
-    verts: (CGA, CGA, CGA),
+    verts: Vec::<CGA>,
     vel: CGA,
     orientation: CGA,
 }
 
 impl Ship {
-    fn update(&mut self, dt_s: f64, keys: &HashSet<String>) {
+    fn update(&mut self, dt_s: f64, keys: &HashSet<String>, boundary: &CGA) {
         let vel = self.vel.clone(); //Velocity at beginning of frame
         let mut vel_2 = self.vel.clone();
 
         let mut a_vec = - DRAG * &vel_2;
 
         if keys.contains("w") {
-            let orient = &self.orientation; //Transforms body frame to local frame
+            let accel_world = self.fixed_to_world(&(ACCEL_STR * CGA::e15()));
-            let pos = &self.com_rotor; //Transforms local frame to world frame
-            let accel_body = ACCEL_STR * CGA::e15(); //The acceleration in the body frame.
-            let accel_local = orient * accel_body * orient.Reverse(); //The acceleration in the local frame.
-            let accel_world = pos * accel_local * pos.Reverse(); //The acceleration in the world frame.
             a_vec = a_vec - accel_world; //Sum of drag and thrust.
         }
 
@@ -111,50 +107,133 @@ impl Ship {
             self.orientation = bivector_exponential(&(-dt_s * CGA::e12())) * &self.orientation;
         }
 
-
         vel_2 = vel_2 + a_vec*dt_s; //Velocity at end of frame
         let vel_3 = 0.5*(vel+&vel_2);
         let delta = vel_3 * dt_s;
         let pos = bivector_exponential(&delta)*&self.com_rotor; //Update position.
-        self.vel = vel_2;
+        let old = self.com_rotor.clone();
         self.com_rotor = pos;
-
+        self.vel = vel_2;
+        if self.intersects_circle(boundary) {
+            let pos = self.fixed_to_world(&gen_hyperbolic_point(&CGA::zero()));
+            let x = pos[E1];
+            let y = pos[E2];
+            let r_sqr = x*x+y*y;
+            let r = r_sqr.sqrt();
+            let clamp = x/r *CGA::e15() + y/r*CGA::e25();
+            let rotor = bivector_exponential(&(clamp*0.01));
+            self.com_rotor = rotor * old;
+            self.vel = CGA::zero();
+        }
     }
 
     fn draw(&self, context: &CanvasRenderingContext2d) {
         let origin = gen_hyperbolic_point(&CGA::zero()); //Origin in the body frame.
-        let orient_rotor = &self.orientation; //Transforms body frame to local frame.
+        let com = self.fixed_to_world(&origin);
-        let com = &self.com_rotor; //Transforms local frame to world frame.
+        let p1 = self.fixed_to_world(&self.verts[0]);
-        let v0_rotor = com*orient_rotor*bivector_exponential(&self.verts.0); //Combined transformation to locate a vertex of the ship in the world frame.
+        let p2 = self.fixed_to_world(&self.verts[1]);
-        let v1_rotor = com*orient_rotor*bivector_exponential(&self.verts.1); //Combined transformation ...
+        let p3 = self.fixed_to_world(&self.verts[2]);
-        let v2_rotor = com*orient_rotor*bivector_exponential(&self.verts.2); //Combined transformation ...
-        let p0 = &v0_rotor*&origin*&v0_rotor.Reverse(); //Produce a vertex of the ship in the world frame.
-        let p1 = &v1_rotor*&origin*&v1_rotor.Reverse(); //...
-        let p2 = &v2_rotor*&origin*&v2_rotor.Reverse(); //...
 
-        draw_point(context, &(com*&origin*com.Reverse()));
+        draw_point(context, &com);
-        draw_point(context, &p0);
         draw_point(context, &p1);
         draw_point(context, &p2);
-        draw_line_between(context, &p0, &p1);
+        draw_point(context, &p3);
         draw_line_between(context, &p1, &p2);
-        draw_line_between(context, &p2, &p0);
+        draw_line_between(context, &p2, &p3);
+        draw_line_between(context, &p3, &p1);
 
         //Construct and draw line pointing straight ahead, to help orient the player.
-        let com_p = com * gen_hyperbolic_point(&CGA::zero()) * com.Reverse();
+        let delta_gen = self.fixed_to_world(&CGA::e14());
-        let delta_gen = com * orient_rotor * CGA::e14() * &orient_rotor.Reverse() * com.Reverse();
+        let dot = delta_gen | &com;
-        let dot = delta_gen | &com_p;
+        let line = CGA::e4() ^ com ^ dot;
-        let line = CGA::e4() ^ com_p ^ dot ^ CGA::e3();
         draw_line(context, &line);
+    }
 
+    fn fixed_to_world(&self, vec: &CGA) -> CGA {
+        &self.com_rotor * &self.orientation * vec * &self.orientation.Reverse() * &self.com_rotor.Reverse()
+    }
 
+    fn intersects_circle(&self, circle: &CGA) -> bool {
+        //Explicitly select only tri-vector part.
+        let circle = grade_selection(circle, 3);
+        //Find vertices in world space.
+        let p1 = self.fixed_to_world(&self.verts[0]);
+        let p2 = self.fixed_to_world(&self.verts[1]);
+        let p3 = self.fixed_to_world(&self.verts[2]);
+        //Construct sides of ship in world space.
+        let e4 = CGA::e4();
+        let side_1 = &e4 ^ &p1 ^ &p2;
+        let side_2 = &e4 ^ &p1 ^ &p3;
+        let side_3 = &e4 ^ &p2 ^ &p3;
+        //Find intersections between any side and the given circle.
+        let int_1 = circle_intersection(&side_1, &circle);
+        let int_2 = circle_intersection(&side_2, &circle);
+        let int_3 = circle_intersection(&side_3, &circle);
+        //For each side, test if the intersection points are between the two vertices defining the
+        //side.
+        match int_1 {
+            CircleIntersection::TwoPoints(a, b) => { 
+                if between_points(&a, &p1, &p2) {
+                    return true;
+                }
+                if between_points(&b, &p1, &p2) {
+                    return true;
+                }
+            },
+            CircleIntersection::None => {},
+        }
+        match int_2 {
+            CircleIntersection::TwoPoints(a, b) => { 
+                if between_points(&a, &p1, &p3) {
+                    return true;
+                }
+                if between_points(&b, &p1, &p3) {
+                    return true;
+                }
+            },
+            CircleIntersection::None => {},
+        }
+        match int_3 {
+            CircleIntersection::TwoPoints(a, b) => { 
+                if between_points(&a, &p2, &p3) {
+                    return true;
+                }
+                if between_points(&b, &p2, &p3) {
+                    return true;
+                }
+            },
+            CircleIntersection::None => {},
+        }
+        false
+    }
+}
+
+fn hyperbolic_distance(point_1:  &CGA, point_2: &CGA) -> f64 {
+    let a = grade_selection(point_1, 1);
+    let b = grade_selection(point_2, 1);
+    let a = hyperbolic_normalize(&a);
+    let b = hyperbolic_normalize(&b);
+    let dot = (&a | &b)[SCALAR];
+    let rad = -0.5*dot;
+    let root = rad.sqrt();
+    let dist = root.asinh();
+    dist
+}
+
+fn between_points(a: &CGA, p1: &CGA, p2: &CGA) -> bool {
+    let total = hyperbolic_distance(p1, p2);
+    let p1_a = hyperbolic_distance(p1, a);
+    let a_p2 = hyperbolic_distance(a, p2);
+    if p1_a + a_p2 - total < ERROR {
+        true
+    } else {
+        false
     }
 }
 
 enum Geometry {
     Circle(f64, f64, f64), //x,y,r
-    Line(f64, f64, f64), //x,r,theta
+    Line(f64, f64, f64), //x,y,theta
-    Null,
 }
 
 /* Translation of a Euclidean Point. Need to use hyperbolic*/
@@ -199,7 +278,6 @@ fn draw_point(context: &CanvasRenderingContext2d, point: &CGA) {
     let (x,y) = point_to_screen_space(x,y);
 
     context.begin_path();
-    context.set_fill_style_str(BLUE);
     context.arc(x, y, 2., 0., TWO_PI).unwrap();
     context.fill();
 
@@ -234,6 +312,11 @@ fn get_hyperbolic_point(point: &CGA) -> (f64, f64) {
     (r*x,r*y)
 }
 
+fn hyperbolic_normalize(point: &CGA) -> CGA {
+    let norm = -(point | CGA::e4())[SCALAR];
+    point * (1. / norm)
+}
+
 fn get_point(point: &CGA) -> (f64, f64) {
     let x = point[E1];
     let y = point[E2];
@@ -242,40 +325,36 @@ fn get_point(point: &CGA) -> (f64, f64) {
     (x/norm, y/norm)
 }
 
-fn get_geometry(circle: &CGA) -> Geometry {
+fn get_geometry(trivector: &CGA) -> Geometry {
-    //Naively assume we are given a valid line/circle
+    let trivector = grade_selection(trivector, 3);
+    let i4 = CGA::e1() ^ CGA::e2() ^ CGA::e4() ^ CGA::e5(); //Left out e3 intentionally
+    let i2 = CGA::e1() ^ CGA::e2();
     let n_bar = CGA::e4() - CGA::e5();
     let n_vec = CGA::e4() + CGA::e5();
-    let dual = (&circle).Dual();
-    let l = (&dual | &n_vec)[SCALAR];
-    let c = (&dual | &n_bar)[SCALAR];
-    let a = (&dual | CGA::e1())[SCALAR];
-    let b = (&dual | CGA::e2())[SCALAR];
 
-    if l.abs() < ERROR {
+    let wedge = &trivector ^ &n_vec;
-        //LINE
+    let mag = magnitude(&wedge);
-        let theta = (-a).atan2(b);
+    let dual = &trivector * &i4;
-        let norm_sqr = a*a + b*b;
+    if mag < ERROR {
-        if norm_sqr < ERROR*ERROR {
+        //Line
-            return Geometry::Null;
+        let d = 0.5*(&dual | &n_bar)[SCALAR];
-        }
+        let m = -1.*(&dual - 0.5*d*&n_vec)*&i2;
-        let x = a*c/norm_sqr/2.;
+        let theta = (m[E2]).atan2(m[E1]);
-        let y = b*c/norm_sqr/2.;
+        return Geometry::Line(0., -d/m[E1], theta);
-        return Geometry::Line(x,y,theta);
     } else {
         //Circle
-        let c_x = a/l;
+        let normalize = -(&dual | &n_vec)[SCALAR];
-        let c_y = b/l;
+        let dual = dual * (1.0/normalize);
-        let r_sqr = c_x*c_x + c_y*c_y + c/l;
+        let center = &trivector * &n_vec * &trivector;
-        let r = r_sqr.sqrt();
+        let (x,y) = get_point(&center);
-        return Geometry::Circle(-c_x, -c_y, r);
+        let sqr = (&dual * &dual)[SCALAR];
+        let r = sqr.sqrt();
+        return Geometry::Circle(x, y, r);
     }
 }
 
 fn draw_line(context: &CanvasRenderingContext2d, line: &CGA) {
     context.begin_path();
-    context.set_stroke_style_str(GREEN);
-
     let geometry = get_geometry(line);
     match geometry {
         Geometry::Circle(x,y,r) => {
@@ -292,16 +371,14 @@ fn draw_line(context: &CanvasRenderingContext2d, line: &CGA) {
                 context.move_to(canvas_x1, canvas_y1);
                 context.line_to(canvas_x2, canvas_y2);
             },
-            Geometry::Null => {},
     }
     context.stroke();
 }
 
 fn draw_line_between(context: &CanvasRenderingContext2d, a: &CGA, b: &CGA) {
-    let line = CGA::e4() ^ a ^ b ^ CGA::e3();
+    let line = CGA::e4() ^ a ^ b;
     let geometry = get_geometry(&line);
     context.begin_path();
-    context.set_stroke_style_str(GREEN);
     match geometry {
         Geometry::Circle(x,y,r) => {
             let (a_x, a_y) = get_point(a);
@@ -352,6 +429,102 @@ fn _get_hyperbolic_translation(vec: &CGA) -> CGA {
     t
 }
 
+enum CircleIntersection {
+    TwoPoints(CGA, CGA),
+    //OnePoint(CGA),
+    None,
+}
+
+fn extract_point_pair(pair: &CGA) -> (CGA, CGA) {
+    let n_vec = CGA::e4() + CGA::e5();
+    let pair = grade_selection(&pair, 2); //Examine only the bivector part
+    let f = pair.normalized();
+    let p = CGA::new(0.5, SCALAR) + 0.5*&f; //Projection Operator
+    let p_rev = CGA::new(0.5, SCALAR) - 0.5*&f; //Projection Operator
+    let diff = &pair | &n_vec; //point A - point B (or a scalar multiple)
+    let a = -1. * &p_rev * &diff;
+    let b = 1. * &p * &diff;
+    (a,b)
+}
+
+fn circle_intersection(circle_1: &CGA, circle_2: &CGA) -> CircleIntersection {
+    //Assumes circles intersect in two places.
+    let circle_1 = grade_selection(circle_1, 3); //Explicitly check only trivector part
+    let circle_2 = grade_selection(circle_2, 3);
+    let i4 = CGA::e1245(); //xy plane in 3D CGA.
+    let ab = &circle_1 * &circle_2;
+    let ba = &circle_2 * &circle_1;
+    let comm = &ab - &ba; //Twice the commutator of circle_1 and circle_2. Unnormalized, but that isn't important.
+    let inter = &i4 * &comm; //The bivector of the pair of intersection points.
+    let sqr = (&inter * &inter)[0];
+    if sqr > ERROR {
+        let (a,b) = extract_point_pair(&inter);
+        CircleIntersection::TwoPoints(a,b)
+    } else {
+        CircleIntersection::None //Technically could have one intersection, but that's not tested yet
+    }
+}
+
+fn _line_intersection(line_1: &CGA, line_2: &CGA) -> bool {
+    let line_1 = grade_selection(line_1, 4);
+    let line_2 = grade_selection(line_2, 4);
+    let meet = line_1 & line_2;
+    let mag = magnitude(&meet);
+    if mag < ERROR {
+        true
+    } else {
+        false
+    }
+}
+
+fn _get_line_intersection(line_1: &CGA, line_2: &CGA) -> Option<CGA> {
+    if !_line_intersection(line_1, line_2) {
+        return None;
+    }
+    let n = CGA::e4() + CGA::e5();
+    let x = CGA::up(0., 0., 0.); //Pick an arbitrary point
+    let l_1_tick = line_2 * line_1 * line_2; //reflect line_1 by line_2
+    let l_1_2tick = line_1 - l_1_tick; //bisector of line_1 and its reflection
+    let x_tick = &l_1_2tick * &x * & l_1_2tick; //reflection of x in the bisector
+    let x_2tick = 0.5*x + 0.5*&x_tick; //midpoint, on the bisector
+    let x_3tick = line_2*&x_2tick*line_2; //Reflection of midpoint in line_2
+    let p_tick = 0.5*&x_2tick + 0.5 * x_3tick; //Intersection point of lines, but with extra n component.
+    let dot = &p_tick | &n;
+    let dot_sqr = dot[SCALAR]*dot[SCALAR];
+    let p = -1.*&p_tick*&n*&p_tick*(0.5/dot_sqr); //Get the point with the correct n component.
+    Some(p)
+}
+
+fn grade_selection(multi: &CGA, grade: usize) -> CGA {
+    match grade {
+        0 => { CGA::new(multi[SCALAR], SCALAR) },
+        1 => { multi[E1]*CGA::e1() + multi[E2]*CGA::e2() + multi[E3]*CGA::e3() + multi[E4]*CGA::e4() + multi[E5]*CGA::e5() },
+        2 => { multi[E12]*CGA::e12() + multi[E13]*CGA::e13() + multi[E14]*CGA::e14() + multi[E15]*CGA::e15()
+            + multi[E23]*CGA::e23() + multi[E24]*CGA::e24() + multi[E25]*CGA::e25()
+            + multi[E34]*CGA::e34() + multi[E35]*CGA::e35() + multi[E45]*CGA::e45() },
+        3 => { multi[E123]*CGA::e123() + multi[E124]*CGA::e124() + multi[E125]*CGA::e125() + multi[E134]*CGA::e134()
+            + multi[E135]*CGA::e135() + multi[E145]*CGA::e145() + multi[E234]*CGA::e234() + multi[E235]*CGA::e235()
+            + multi[E245]*CGA::e245() + multi[E345]*CGA::e345() },
+        4 => { multi[E1234]*CGA::e1234() + multi[E1235]*CGA::e1235() + multi[E1245]*CGA::e1245() + multi[E1345]*CGA::e1345() 
+            + multi[E2345]*CGA::e2345() },
+        5 => { multi[E12345]*CGA::e12345() }
+        _ => { CGA::zero() }
+    }
+}
+
+fn magnitude(multi: &CGA) -> f64 {
+    multi[SCALAR]*multi[SCALAR]
+    + multi[E1]*multi[E1] + multi[E2]*multi[E2] + multi[E3]*multi[E3] + multi[E4]*multi[E4] + multi[E5]*multi[E5]
+    + multi[E12]*multi[E12] + multi[E13]*multi[E13] + multi[E14]*multi[E14] + multi[E15]*multi[E15]
+    + multi[E23]*multi[E23] + multi[E24]*multi[E24] + multi[E25]*multi[E25]
+    + multi[E34]*multi[E34] + multi[E35]*multi[E35] + multi[E45]*multi[E45]
+    + multi[E123]*multi[E123] + multi[E124]*multi[E124] + multi[E125]*multi[E125]
+    + multi[E134]*multi[E134] + multi[E135]*multi[E135] + multi[E145]*multi[E145]
+    + multi[E234]*multi[E234] + multi[E235]*multi[E235] + multi[E245]*multi[E245] + multi[E345]*multi[E345]
+    + multi[E1234]*multi[E1234] + multi[E1235]*multi[E1235] + multi[E1245]*multi[E1245] + multi[E1345]*multi[E1345] + multi[E2345]*multi[E2345]
+    + multi[E12345]*multi[E12345]
+}
+
 fn gen_hyperbolic_point(vec: &CGA) -> CGA {
     let vec_sqr = (vec * vec)[SCALAR];
     if vec_sqr >= 1. {
@@ -365,20 +538,7 @@ fn gen_hyperbolic_point(vec: &CGA) -> CGA {
 fn bivector_exponential(bivector: &CGA) -> CGA {
     //Explicitly take only the bivector part.
     //The library doesn't provide a grade 2 selection.
-    let a = bivector[E12];
+    let bivector = grade_selection(bivector, 2);
-    let b = bivector[E13];
-    let c = bivector[E14];
-    let d = bivector[E15];
-    let e = bivector[E23];
-    let f = bivector[E24];
-    let g = bivector[E25];
-    let h = bivector[E34];
-    let i = bivector[E35];
-    let j = bivector[E45];
-
-    let bivector = a*CGA::e12() + b*CGA::e13() + c*CGA::e14() 
-        + d*CGA::e15()  + e*CGA::e23() + f*CGA::e24() + g*CGA::e25()
-        + h*CGA::e34() + i *CGA::e35() + j*CGA::e45();
     let sqr = (&bivector*&bivector)[SCALAR];
     if sqr.abs() < ERROR {
         CGA::new(1., SCALAR) + bivector
@@ -392,39 +552,124 @@ fn bivector_exponential(bivector: &CGA) -> CGA {
 }
 
 struct Asteroid {
-    circle: CGA,
+    com_rotor: CGA,
+    verts: Vec::<CGA>,
     vel: CGA,
 }
 
+fn circle_from_coords(x: f64, y: f64, r: f64) -> CGA {
+    let (a_x, a_y) = (x + r * 0.0_f64.cos(), y + r * 0.0_f64.sin());
+    let (b_x, b_y) = (x + r * 1.0_f64.cos(), y + r * 1.0_f64.sin());
+    let (c_x, c_y) = (x + r * 2.0_f64.cos(), y + r * 2.0_f64.sin());
+    let a = point_to_cga(a_x, a_y);
+    let b = point_to_cga(b_x, b_y);
+    let c = point_to_cga(c_x, c_y);
+    let circle = a^b^c;
+    circle
+}
+
 impl Asteroid {
     fn new() -> Asteroid {
-        let v = 0.1;
+        let vec = 0.8*CGA::e1() + 0.0*CGA::e2();
+        let vec_sqr = (&vec*&vec)[SCALAR];
+        let scale = 1. / (1.-vec_sqr).sqrt();
+        let com_rotor = scale * (CGA::new(1., SCALAR) + CGA::e5()*vec);
+        let tx_gen = CGA::e15();
+        let ty_gen = CGA::e25();
+        let v1_gen = -0.2*&tx_gen;
+        let v2_gen = 0.2*&tx_gen;
+        let v3_gen = 0.2*&ty_gen;
+        let v1_rot = bivector_exponential(&v1_gen);
+        let v2_rot = bivector_exponential(&v2_gen);
+        let v3_rot = bivector_exponential(&v3_gen);
+
+        let origin = gen_hyperbolic_point(&CGA::zero());
+        let v1 = &v1_rot * &origin * &v1_rot.Reverse();
+        let v2 = &v2_rot * &origin * &v2_rot.Reverse();
+        let v3 = &v3_rot * &origin * &v3_rot.Reverse();
+        let verts = vec![v1, v2, v3];
+
+        let v = 1.;
         let (v_x, v_y) = (v * 1.0_f64.cos(), v*1.0_f64.sin());
-        let vel = v_x*CGA::e15() + v_y*CGA::e25();
+        let vel = v_x*CGA::e1() + v_y*CGA::e2();
-        Self::new_from_coords(0.0_f64, 0.0_f64, 0.1_f64, vel)
+        Asteroid { com_rotor, vel, verts }
     }
-
+/*
-    fn new_from_coords(x: f64, y: f64, r: f64, vel: CGA) -> Asteroid {
+    fn _new_from_coords(x: f64, y: f64, r: f64, vel: CGA) -> Asteroid {
-        let (a_x, a_y) = (x + r * 0.0_f64.cos(), y + r * 0.0_f64.sin());
+        let (circle, verts) = circle_from_coords(x, y, r);
-        let (b_x, b_y) = (x + r * 1.0_f64.cos(), y + r * 1.0_f64.sin());
+        Asteroid { circle, vel, verts }
-        let (c_x, c_y) = (x + r * 2.0_f64.cos(), y + r * 2.0_f64.sin());
-        let a = point_to_cga(a_x, a_y);
-        let b = point_to_cga(b_x, b_y);
-        let c = point_to_cga(c_x, c_y);
-        let circle = a^b^c^CGA::e3();
-        Asteroid { circle, vel }
     }
+*/
 
-    fn update(&mut self, dt_s: f64) {
+    fn update(&mut self, dt_s: f64, boundary: &CGA) {
-        let vel = &self.vel;
+        let vel = &self.vel*CGA::e5();
-        let delta = vel * dt_s;
+        let delta = &vel * dt_s;
+        let delta = self.fixed_to_world(&delta);//&self.com_rotor * delta * &self.com_rotor.Reverse();
         let vel_rotor = bivector_exponential(&delta);
-        let circle = &vel_rotor*&self.circle*&vel_rotor.Reverse();
+        //let vel_rotor = &self.com_rotor.Reverse() * vel_rotor * &self.com_rotor;
-        self.circle = circle;
+        let new_pos = &vel_rotor * &self.com_rotor;
+        let p1 = &new_pos * &self.verts[0] * &new_pos.Reverse();
+        let p2 = &new_pos * &self.verts[1] * &new_pos.Reverse();
+        let p3 = &new_pos * &self.verts[2] * &new_pos.Reverse();
+        let circle = p1 ^ p2 ^ p3;
+        let intersect = circle_intersection(&circle, boundary);
+        match intersect {
+            CircleIntersection::None => { self.com_rotor = new_pos; }
+            CircleIntersection::TwoPoints(_,_) => {
+                self.reflect_off_boundary();
+            }
+        }
     }
 
     fn draw(&self, context: &CanvasRenderingContext2d) {
-        draw_line(context, &self.circle);
+        let p1 = self.fixed_to_world(&self.verts[0]);
+        let p2 = self.fixed_to_world(&self.verts[1]);
+        let p3 = self.fixed_to_world(&self.verts[2]);
+        draw_point(context, &p1);
+        draw_point(context, &p2);
+        draw_point(context, &p3);
+        let circle = &p1 ^ &p2 ^ &p3;
+        draw_line(context, &circle);
+        //Construct and draw line pointing straight ahead, to help orient the player.
+        let origin = gen_hyperbolic_point(&CGA::zero());
+        let com = self.fixed_to_world(&origin);
+        let delta_gen = self.fixed_to_world(&(&self.vel*CGA::e5()));
+        let dot = delta_gen | &com;
+        let line = CGA::e4() ^ com ^ dot;
+        draw_line(context, &line);
+    }
+
+    fn fixed_to_world(&self, vec: &CGA) -> CGA {
+        &self.com_rotor*vec*&self.com_rotor.Reverse()
+    }
+
+    fn get_circle(&self) -> CGA {
+        let p1 = self.fixed_to_world(&self.verts[0]);
+        let p2 = self.fixed_to_world(&self.verts[1]);
+        let p3 = self.fixed_to_world(&self.verts[2]);
+        p1 ^ p2 ^ p3
+    }
+
+    fn reflect_off_boundary(&mut self) {
+        let origin = gen_hyperbolic_point(&CGA::zero());
+        let mut pos = self.fixed_to_world(&origin);
+        pos[E4] = 0.;
+        pos[E5] = 0.;
+        let norm = pos.norm();
+        pos = pos * (1./ norm);
+        let norm_old = self.vel.norm();
+        let vel = self.fixed_to_world(&self.vel);
+        let dot = &vel | &pos;
+        let vel_para_radial = dot * &pos;
+        let vel_perp_radial = &vel - &vel_para_radial;
+        let vel_tick = vel_perp_radial - vel_para_radial;
+        let mut vel_fixed = &self.com_rotor.Reverse()*vel_tick*&self.com_rotor;
+        vel_fixed[E4] = 0.;
+        vel_fixed[E5] = 0.;
+        let norm = vel_fixed.norm();
+        vel_fixed = (1. / norm)*vel_fixed;
+        vel_fixed = norm_old * vel_fixed;
+        self.vel = vel_fixed;
     }
 }
 
@@ -433,7 +678,9 @@ pub struct Game {
     context: CanvasRenderingContext2d,
     ship: Ship,
     asteroids: Vec::<Asteroid>,
+    boundary: CGA,
     keys: HashSet<String>,
+    dt: f64
 }
 
 #[wasm_bindgen]
@@ -456,18 +703,33 @@ impl Game {
         let tx_gen = CGA::e15();
         let ty_gen = CGA::e25();
 
+        let v1_gen = -0.1*&tx_gen;
+        let v2_gen = 0.05*&tx_gen+0.02*&ty_gen;
+        let v3_gen = 0.05*&tx_gen-0.02*&ty_gen;
+        let v1_rot = bivector_exponential(&v1_gen);
+        let v2_rot = bivector_exponential(&v2_gen);
+        let v3_rot = bivector_exponential(&v3_gen);
+
+        let origin = gen_hyperbolic_point(&CGA::zero());
+        let v1 = &v1_rot * &origin * &v1_rot.Reverse();
+        let v2 = &v2_rot * &origin * &v2_rot.Reverse();
+        let v3 = &v3_rot * &origin * &v3_rot.Reverse();
+
         let ship = Ship {
             com_rotor: CGA::new(1., SCALAR),
-            verts: (-0.1*&tx_gen, 0.05*&tx_gen+0.02*&ty_gen, 0.05*&tx_gen-0.02*ty_gen),
+            verts: vec![ v1, v2, v3],
             vel: CGA::zero(),
             orientation: CGA::new(1., SCALAR),
         };
 
+        let boundary = circle_from_coords(0., 0., 0.9);
         Game { 
             context: context,
             ship,
+            boundary,
             asteroids: vec![Asteroid::new()],
             keys: HashSet::new(),
+            dt: 100.
         }
     }
 
@@ -479,10 +741,32 @@ impl Game {
         self.context.set_stroke_style_str(BLACK);
         self.context.arc(CENTER_X, CENTER_Y, RADIUS, 0., TWO_PI).unwrap();
         self.context.stroke();
-        self.ship.draw(&self.context);
+
+        let fps = 1000. / self.dt;
+        self.context.set_fill_style_str(BLACK);
+        self.context.set_font("20px Arial");
+        let message = format!("FPS: {}", fps);
+        self.context.fill_text(&message, 10., 30.).unwrap();
+
+        draw_line(&self.context, &self.boundary);
+
+        let mut intersects = false;
+
+        self.context.set_stroke_style_str(RED);
         for asteroid in &self.asteroids {
+            let circle = asteroid.get_circle();
+            if self.ship.intersects_circle(&circle) {
+                intersects = true;
+            }
             asteroid.draw(&self.context);
         }
+
+        if intersects {
+            self.context.set_stroke_style_str(RED);
+        } else {
+            self.context.set_stroke_style_str(GREEN);
+        }
+        self.ship.draw(&self.context);
     }
 
     pub fn update(&mut self, dt_m: f64) {
@@ -492,11 +776,13 @@ impl Game {
             dt_m / 1000.
         };
 
+        self.dt = dt_m;
+
         for asteroid in &mut self.asteroids {
-            asteroid.update(dt_s);
+            asteroid.update(dt_s, &self.boundary);
         }
 
-        self.ship.update(dt_s, &self.keys);
+        self.ship.update(dt_s, &self.keys, &self.boundary);
     }
 
     pub fn key_down(&mut self, key: String) {