Add circular arc and spiral curves

5 files changed, 188 insertions, 3 deletions

diff --git a/curve.cpp b/curve.cpp
index 9b63e8d..73f40c1 100644
--- a/curve.cpp
+++ b/curve.cpp
@@ -1,3 +1,4 @@
+#include <iostream>
 #include <cmath>
 #include "curve.h"
 
@@ -21,3 +22,91 @@ vec3 LineSegment::integrate(vec3 (*integrand)(vec3 s, vec3 ds), double dl)
 
     return sum;
 }
+
+vec3 Arc::integrate(vec3 (*integrand)(vec3 s, vec3 ds), double dl)
+{
+    //cout << "Integrating loop of length " << this->angle << " radians" << endl;
+
+    double r = this->radius.magnitude(), dtheta = dl / r;
+
+    //cout << "R = " << r << ", dtheta = " << dtheta << endl;
+
+    quat rot = quat::from_angleaxis(dtheta, normal), crot = rot.conjugate();
+
+    //cout << "rot = " << rot << ", crot = " << crot << endl;
+
+    double len = this->angle * r, l;
+
+    vec3 sum = 0;
+
+    quat radius = this->radius;
+
+    for(l = 0; l < len; l += dl)
+    {
+        vec3 ds = this->normal.cross(radius).normalize() * dl;
+        sum += integrand(this->center + radius, ds);
+
+        radius = rot * radius * crot;
+    }
+
+    if(l < len)
+    {
+        dl = len - l;
+        dtheta = dl / r;
+        rot = quat::from_angleaxis(dtheta, normal);
+        crot = rot.conjugate();
+
+        vec3 ds = this->normal.cross(radius).normalize() * dl;
+        sum += integrand(this->center + radius, ds);
+    }
+
+    return sum;
+}
+
+vec3 Spiral::integrate(vec3 (*integrand)(vec3 s, vec3 ds), double dl)
+{
+    //cout << "Integrating loop of length " << this->angle << " radians" << endl;
+
+    double r = this->radius.magnitude(), dtheta = dl / r;
+
+    //cout << "R = " << r << ", dtheta = " << dtheta << endl;
+
+    quat rot = quat::from_angleaxis(dtheta, normal), crot = rot.conjugate();
+
+    //cout << "rot = " << rot << ", crot = " << crot << endl;
+
+    double len = this->angle * r, l;
+
+    vec3 sum = 0;
+
+    quat radius = this->radius;
+
+    /* how far along the axis we've moved */
+    vec3 axis = 0;
+
+    /* we add this axial component each iteration */
+    vec3 dp = this->normal * dtheta * pitch / (2 * M_PI);
+
+    for(l = 0; l < len; l += dl, axis += dp)
+    {
+        vec3 ds = this->normal.cross(radius).normalize() * dl + dp;
+        sum += integrand(this->origin + axis + radius, ds);
+
+        /* rotate by dtheta about normal */
+        radius = rot * radius * crot;
+    }
+
+    if(l < len)
+    {
+        dl = len - l;
+        dtheta = dl / r;
+        rot = quat::from_angleaxis(dtheta, normal);
+        crot = rot.conjugate();
+        dp = this->normal * dtheta * pitch / (2 * M_PI);
+
+        vec3 ds = this->normal.cross(radius).normalize() * dl + dp;
+        sum += integrand(this->origin + axis + radius, ds);
+    }
+
+    return sum;
+}
diff --git a/curve.h b/curve.h
index c896a23..e9e9642 100644
--- a/curve.h
+++ b/curve.h
@@ -1,5 +1,10 @@
+#ifndef CURVE_H
+#define CURVE_H
+
 #include <cmath>
 #include "vec3.h"
+#include "quat.h"
+
 using namespace std;
 
 class Curve {
@@ -19,10 +24,35 @@ public:
 class Arc : Curve {
 private:
     vec3 center;
-    double radius;
-    double angle[2]; /* start and end angles */
+
+    /* these are relative to the center (direction will be determined
+     * by RHR of normal), and should be orthonormal */
+    vec3 radius, normal;
+
+    /* how many radians the arc extends for (can be greater than 2pi) */
+    double angle;
 public:
-    LineSegment(vec3 a_, vec3 b_) : a(a_), b(b_) {};
+    Arc(vec3 c_, vec3 r_, vec3 n_, double th) : center(c_), radius(r_), normal(n_), angle(th) {};
+
+    vec3 integrate(vec3 (*integrand)(vec3 s, vec3 ds), double delta);
+};
+
+class Spiral : Curve {
+private:
+    vec3 origin;
+
+    /* these are relative to the center (direction will be determined
+     * by RHR of normal), and should be orthonormal */
+    vec3 radius, normal;
+
+    /* how many radians the arc extends for (can be greater than 2pi) */
+    double angle;
+
+    /* space between turns (2pi) */
+    double pitch;
+public:
+    Spiral(vec3 c_, vec3 r_, vec3 n_, double th, double p) : origin(c_), radius(r_), normal(n_), angle(th), pitch(p) {};
 
     vec3 integrate(vec3 (*integrand)(vec3 s, vec3 ds), double delta);
 };
+#endif
diff --git a/quat.cpp b/quat.cpp
new file mode 100644
index 0000000..aa53af6
--- /dev/null
+++ b/quat.cpp
@@ -0,0 +1,38 @@
+#include "quat.h"
+#include <cmath>
+
+quat::quat(double w_, double x_, double y_, double z_) : w(w_), x(x_), y(y_), z(z_) { }
+quat::quat(double x_, double y_, double z_) : w(0), x(x_), y(y_), z(z_) { }
+quat::quat(double w_, vec3 vec) : w(w_), x(vec[0]), y(vec[1]), z(vec[2]) { }
+quat::quat(vec3 vec) : w(0), x(vec[0]), y(vec[1]), z(vec[2]) { }
+quat::quat() : w(0), x(0), y(0), z(0) { }
+
+quat::operator vec3()
+{
+    return vec3(this->x, this->y, this->z);
+}
+
+quat operator*(const quat &lhs, const quat &rhs)
+{
+    return quat(lhs.w * rhs.w - lhs.x * rhs.x - lhs.y * rhs.y - lhs.z * rhs.z,
+                lhs.w * rhs.x + lhs.x * rhs.w + lhs.y * rhs.z - lhs.z * rhs.y,
+                lhs.w * rhs.y - lhs.x * rhs.z + lhs.y * rhs.w + lhs.z * rhs.x,
+                lhs.w * rhs.z + lhs.x * rhs.y - lhs.y * rhs.x + lhs.z * rhs.w);
+}
+
+quat quat::conjugate() const
+{
+    return quat(this->w, -this->x, -this->y, -this->z);
+}
+
+quat quat::from_angleaxis(double angle, vec3 axis)
+{
+    double si = std::sin(angle / 2);
+    double co = std::cos(angle / 2);
+    return quat(co, si * axis[0], si * axis[1], si * axis[2]);
+}
+
+std::ostream &operator<<(std::ostream &os, const quat &q)
+{
+    return os << "(" << q.w << ", " << q.x << ", " << q.y << ", " << q.z << ")";
+}
diff --git a/quat.h b/quat.h
new file mode 100644
index 0000000..65a828a
--- /dev/null
+++ b/quat.h
@@ -0,0 +1,25 @@
+#ifndef QUAT_H
+#define QUAT_H
+#include "vec3.h"
+#include <iostream>
+
+class quat {
+public:
+    double w, x, y, z;
+public:
+    quat(double w, double x, double y, double z);
+    quat(double x, double y, double z);
+    quat(double w, vec3 vec);
+    quat(vec3 vec);
+    quat();
+
+    operator vec3();
+
+    quat conjugate() const;
+
+    static quat from_angleaxis(double angle, vec3 axis);
+};
+
+quat operator*(const quat &, const quat &);
+std::ostream &operator<<(std::ostream &os, const quat &);
+#endif
diff --git a/vec3.h b/vec3.h
index 1f01098..abfc6c3 100644
--- a/vec3.h
+++ b/vec3.h
@@ -1,3 +1,5 @@
+#ifndef VEC3_H
+#define VEC3_H
 #include <iostream>
 class vec3 {
     public:
@@ -25,3 +27,4 @@ class vec3 {
 };
 
 std::ostream &operator<<(std::ostream &output, const vec3 &v);
+#endif