Add circular arc and spiral curves

1 files changed, 89 insertions, 0 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;
+}