Work on adding toroid

7 files changed, 238 insertions, 70 deletions

diff --git a/curve.cpp b/curve.cpp
index 73f40c1..c71a8c3 100644
--- a/curve.cpp
+++ b/curve.cpp
@@ -4,15 +4,15 @@
 
 using namespace std;
 
-vec3 LineSegment::integrate(vec3 (*integrand)(vec3 s, vec3 ds), double dl)
+vec3 LineSegment::integrate(vec3 (*integrand)(vec3 s, vec3 ds), scalar dl) const
 {
     vec3 diff = this->b - this->a, sum = 0;
 
-    double len = diff.magnitude();
+    scalar len = diff.magnitude();
 
     vec3 diffnorm = diff / len, s = this->a, ds = diffnorm * dl;
 
-    double l;
+    scalar l;
 
     for(l = 0; l < len; l += dl, s += ds)
         sum += integrand(s, ds);
@@ -23,11 +23,11 @@ vec3 LineSegment::integrate(vec3 (*integrand)(vec3 s, vec3 ds), double dl)
     return sum;
 }
 
-vec3 Arc::integrate(vec3 (*integrand)(vec3 s, vec3 ds), double dl)
+vec3 Arc::integrate(vec3 (*integrand)(vec3 s, vec3 ds), scalar dl) const
 {
     //cout << "Integrating loop of length " << this->angle << " radians" << endl;
 
-    double r = this->radius.magnitude(), dtheta = dl / r;
+    scalar r = this->radius.magnitude(), dtheta = dl / r;
 
     //cout << "R = " << r << ", dtheta = " << dtheta << endl;
 
@@ -35,7 +35,7 @@ vec3 Arc::integrate(vec3 (*integrand)(vec3 s, vec3 ds), double dl)
 
     //cout << "rot = " << rot << ", crot = " << crot << endl;
 
-    double len = this->angle * r, l;
+    scalar len = this->angle * r, l;
 
     vec3 sum = 0;
 
@@ -63,11 +63,11 @@ vec3 Arc::integrate(vec3 (*integrand)(vec3 s, vec3 ds), double dl)
     return sum;
 }
 
-vec3 Spiral::integrate(vec3 (*integrand)(vec3 s, vec3 ds), double dl)
+vec3 Spiral::integrate(vec3 (*integrand)(vec3 s, vec3 ds), scalar dl) const
 {
     //cout << "Integrating loop of length " << this->angle << " radians" << endl;
 
-    double r = this->radius.magnitude(), dtheta = dl / r;
+    scalar r = this->radius.magnitude(), dtheta = dl / r;
 
     //cout << "R = " << r << ", dtheta = " << dtheta << endl;
 
@@ -75,7 +75,8 @@ vec3 Spiral::integrate(vec3 (*integrand)(vec3 s, vec3 ds), double dl)
 
     //cout << "rot = " << rot << ", crot = " << crot << endl;
 
-    double len = this->angle * r, l;
+    /* only flat (doesn't include stretched length) */
+    scalar len = this->angle * r, l;
 
     vec3 sum = 0;
 
@@ -110,3 +111,73 @@ vec3 Spiral::integrate(vec3 (*integrand)(vec3 s, vec3 ds), double dl)
 
     return sum;
 }
+
+vec3 Toroid::integrate(vec3 (*integrand)(vec3 s, vec3 ds), scalar dl) const
+{
+    /* this applies to the minor rotation */
+    scalar minor_r = this->minor_r;
+
+    scalar major_r = this->major_radius.magnitude();
+
+    /* how far parallel to the major median each dl takes us */
+    scalar dpar = this->pitch * dl / (2 * M_PI * minor_r) * major_r;
+    
+    /* how far each dl rotates the major vector */
+    scalar major_dtheta = dpar / major_r;
+
+    /* how long perpendicular to the major median each dl takes us (from
+     * Pythagoras) */
+    scalar dper = sqrt(dl * dl - dpar * dpar);
+
+    scalar minor_dtheta = dper / minor_r;
+
+    /* minor normal is binormal to major normal and radius */
+    vec3 minor_normal = this->major_normal.cross(this->major_radius).normalize();
+
+    /* quaternion describing rotation of minor vector */
+    quat minor_rot = quat::from_angleaxis(minor_dtheta, minor_normal), minor_crot = minor_rot.conjugate();
+
+    quat major_rot = quat::from_angleaxis(major_dtheta, major_normal), major_crot = major_rot.conjugate();
+
+    vec3 sum = 0;
+
+    quat major_radius = this->major_radius;
+    quat minor_radius = this->major_radius.normalize() * this->minor_r;
+
+    for(scalar theta = 0; theta < this->major_angle; theta += major_dtheta)
+    {
+        vec3 dpar_vec = minor_normal * dpar;
+        vec3 ds = minor_normal.cross(minor_radius).normalize() * dper + dpar_vec;
+        sum += integrand(this->origin + major_radius + minor_radius, ds);
+
+        /* we need to rotate some vectors: first the major radius and
+         * minor normal about the major normal, and then the minor
+         * radius */
+        /* we also need to rotate the minor radius rotation quaternion
+         * about the origin */
+        major_radius = major_rot * major_radius * major_crot;
+        minor_normal = major_rot * minor_normal * major_crot;
+            
+        minor_radius = minor_rot * minor_radius * minor_crot;
+
+        minor_rot = minor_rot * major_rot;
+        minor_crot = minor_rot.conjugate();
+    }
+
+#if 0
+    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);
+    }
+#endif
+
+    return sum;
+}
+
diff --git a/curve.h b/curve.h
index e9e9642..52c7623 100644
--- a/curve.h
+++ b/curve.h
@@ -9,7 +9,7 @@ using namespace std;
 
 class Curve {
 public:
-    virtual vec3 integrate(vec3 (*integrand)(vec3 s, vec3 ds), double delta) = 0;
+    virtual vec3 integrate(vec3 (*integrand)(vec3 s, vec3 ds), scalar delta) const = 0;
 };
 
 class LineSegment : Curve {
@@ -18,7 +18,7 @@ private:
 public:
     LineSegment(vec3 a_, vec3 b_) : a(a_), b(b_) {};
 
-    vec3 integrate(vec3 (*integrand)(vec3 s, vec3 ds), double delta);
+    vec3 integrate(vec3 (*integrand)(vec3 s, vec3 ds), scalar delta) const;
 };
 
 class Arc : Curve {
@@ -30,11 +30,11 @@ private:
     vec3 radius, normal;
 
     /* how many radians the arc extends for (can be greater than 2pi) */
-    double angle;
+    scalar angle;
 public:
-    Arc(vec3 c_, vec3 r_, vec3 n_, double th) : center(c_), radius(r_), normal(n_), angle(th) {};
+    Arc(vec3 c_, vec3 r_, vec3 n_, scalar th) : center(c_), radius(r_), normal(n_), angle(th) {};
 
-    vec3 integrate(vec3 (*integrand)(vec3 s, vec3 ds), double delta);
+    vec3 integrate(vec3 (*integrand)(vec3 s, vec3 ds), scalar delta) const;
 };
 
 class Spiral : Curve {
@@ -46,13 +46,37 @@ private:
     vec3 radius, normal;
 
     /* how many radians the arc extends for (can be greater than 2pi) */
-    double angle;
+    scalar angle;
 
-    /* space between turns (2pi) */
-    double pitch;
+    /* linear distance between turns (2pi) */
+    scalar pitch;
 public:
-    Spiral(vec3 c_, vec3 r_, vec3 n_, double th, double p) : origin(c_), radius(r_), normal(n_), angle(th), pitch(p) {};
+    Spiral(vec3 c_, vec3 r_, vec3 n_, scalar th, scalar p) : origin(c_), radius(r_), normal(n_), angle(th), pitch(p) {};
 
-    vec3 integrate(vec3 (*integrand)(vec3 s, vec3 ds), double delta);
+    vec3 integrate(vec3 (*integrand)(vec3 s, vec3 ds), scalar delta) const;
+};
+
+class Toroid : Curve {
+private:
+    vec3 origin;
+
+    /* these are relative to the center (direction will be determined
+     * by RHR of normal), and should be orthonormal */
+    vec3 major_radius, major_normal;
+
+    /* "thickness" of toroid */
+    scalar minor_r;
+    
+    /* how many radians (about the center) the toroid extends for
+     * (can be greater than 2pi) */
+    scalar major_angle;
+
+    /* central angle between successive turns (2pi rotation of small
+     * radius vector) */
+    scalar pitch;
+public:
+    Toroid(vec3 o, vec3 maj_r, vec3 maj_n, scalar ang, scalar min_r, scalar p) : origin(o), major_radius(maj_r), major_normal(maj_n), major_angle(ang), minor_r(min_r), pitch(p) {};
+
+    vec3 integrate(vec3 (*integrand)(vec3 s, vec3 ds), scalar delta) const;
 };
 #endif
diff --git a/quat.cpp b/quat.cpp
index aa53af6..df117d0 100644
--- a/quat.cpp
+++ b/quat.cpp
@@ -1,9 +1,9 @@
 #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(scalar w_, scalar x_, scalar y_, scalar z_) : w(w_), x(x_), y(y_), z(z_) { }
+quat::quat(scalar x_, scalar y_, scalar z_) : w(0), x(x_), y(y_), z(z_) { }
+quat::quat(scalar 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) { }
 
@@ -25,10 +25,10 @@ quat quat::conjugate() const
     return quat(this->w, -this->x, -this->y, -this->z);
 }
 
-quat quat::from_angleaxis(double angle, vec3 axis)
+quat quat::from_angleaxis(scalar angle, vec3 axis)
 {
-    double si = std::sin(angle / 2);
-    double co = std::cos(angle / 2);
+    scalar si = std::sin(angle / 2);
+    scalar co = std::cos(angle / 2);
     return quat(co, si * axis[0], si * axis[1], si * axis[2]);
 }
 
diff --git a/quat.h b/quat.h
index 65a828a..5821f81 100644
--- a/quat.h
+++ b/quat.h
@@ -5,11 +5,11 @@
 
 class quat {
 public:
-    double w, x, y, z;
+    scalar 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(scalar w, scalar x, scalar y, scalar z);
+    quat(scalar x, scalar y, scalar z);
+    quat(scalar w, vec3 vec);
     quat(vec3 vec);
     quat();
 
@@ -17,7 +17,7 @@ public:
 
     quat conjugate() const;
 
-    static quat from_angleaxis(double angle, vec3 axis);
+    static quat from_angleaxis(scalar angle, vec3 axis);
 };
 
 quat operator*(const quat &, const quat &);
diff --git a/test.cpp b/test.cpp
index 92c9272..0a42f40 100644
--- a/test.cpp
+++ b/test.cpp
@@ -1,6 +1,13 @@
 #include <cmath>
 #include <iostream>
+#include <fstream>
+#include <sys/stat.h>
+#include <sys/types.h>
+#include <sstream>
+
+#include "vec3.h"
 #include "curve.h"
+
 using namespace std;
 
 vec3 integrand(vec3 s, vec3 ds)
@@ -17,7 +24,7 @@ vec3 dB(vec3 s, vec3 ds)
 {
     vec3 r = s - point;
 
-    double r2 = r.magnitudeSquared();
+    scalar r2 = r.magnitudeSquared();
 
     vec3 rnorm = r / std::sqrt(r2);
 
@@ -26,21 +33,39 @@ vec3 dB(vec3 s, vec3 ds)
 
 //Arc loop(vec3(0, 0, 0), vec3(0, 1, 0), vec3(1, 0, 0), M_PI * 2);
 //Spiral loop(vec3(0, 0, 0), vec3(0, 1, 0), vec3(1, 0, 0), M_PI * 2 * 10, 1);
-LineSegment loop(vec3(0, 0, 0), vec3(10, 0, 0));
+//LineSegment loop(vec3(-.1, .1, 0), vec3(.1, .1, 0));
+Toroid loop(vec3(0, 0, 0), vec3(1, 0, 0), vec3(0, 0, 1), M_PI * 2, .1, 2*M_PI / 10);
+
+ostream *dump_ofs = NULL;
+
+vec3 dump(vec3 s, vec3 ds)
+{
+    *dump_ofs << s << endl;
+    //cout << "Magn: " << ds.magnitude() << endl;
+        
+    return 0;
+}
+
+void dump_path(ostream &out, const Curve *c)
+{
+    dump_ofs = &out;
+    c->integrate(dump, 1e-2);
+}
 
 /* dump the field (gnuplot format) at z = 0 */
 /* requires x0 < x1, y0 < y1 */
-void dump_field(double x0, double y0, double z0, double x1, double y1, double z1)
-{
-    const double delta = .2;
 
-    for(double z = z0; z <= z1; z += delta)
-        for(double y = y0; y <= y1; y += delta)
-            for(double x = x0; x <= x1; x += delta)
+const scalar U0 = 4e-7 * M_PI;
+const scalar I = 1;
+const scalar delta = .2;
+
+void dump_field(scalar x0, scalar y0, scalar z0, scalar x1, scalar y1, scalar z1)
+{
+    for(scalar z = z0; z <= z1; z += delta)
+        for(scalar y = y0; y <= y1; y += delta)
+            for(scalar x = x0; x <= x1; x += delta)
             {
                 point = vec3(x, y, z);
-                const double U0 = 4e-7 * M_PI;
-                const double I = 1;
 
                 vec3 B = loop.integrate(dB, 1e-1) * U0 * I;
 
@@ -53,6 +78,21 @@ void dump_field(double x0, double y0, double z0, double x1, double y1, double z1
             }
 }
 
+void dump_fieldline(ostream &out, vec3 x, scalar len)
+{
+    point = x;
+    scalar delta = .1;
+    while(len > 0)
+    {
+        out << point << endl;
+        
+        vec3 B = loop.integrate(dB, 1e-1) * U0 * I;
+
+        point += delta * B;
+        len -= delta;
+    }
+}
+
 int main()
 {
     //LineSegment wire(vec3(0, -100, 0), vec3(0, 100, 0));
@@ -63,11 +103,36 @@ int main()
 
     //point = 0;
 
-    //double I = 1;
+    //scalar I = 1;
 
     //for(int i = 0; i < 1000; i++, point += dx)
     //std::cout << point[0] << " " << U0 / ( 4 * M_PI ) * loop.integrate(dB, 1e-2)[0] << endl;
-
-    dump_field(-1, -1.5, -1.5,
-               11, 1.5, 1.5);
+    
+    dump_field(-2, -2, 0,
+                12, 2, 0);
+    //dump_field(0,0,0,0,0,0);
+
+    stringstream ss;
+    ss << "curve.fld";
+    ofstream ofs(ss.str());
+
+    dump_path(ofs, (Curve*)&loop);
+    
+    ofs.close();
+    
+    
+#if 0
+    mkdir("field", 0755);
+
+    for(scalar y = -1.5; y <= 1.5; y += .1)
+    {
+        stringstream ss;
+        ss << "field/" << y << ".fld";
+        ofstream ofs(ss.str());
+
+        dump_fieldline(ofs, vec3(0, y, 0), 10);
+
+        ofs.close();
+    }
+#endif
 }
diff --git a/vec3.cpp b/vec3.cpp
index 8513916..51d4870 100644
--- a/vec3.cpp
+++ b/vec3.cpp
@@ -11,26 +11,26 @@ vec3::vec3() {
     v[1] = 0;
     v[2] = 0;
 }
-vec3::vec3(double x) {
+vec3::vec3(scalar x) {
     v[0] = x;
     v[1] = 0;
     v[2] = 0;
 }
-vec3::vec3(double x, double y, double z) {
+vec3::vec3(scalar x, scalar y, scalar z) {
     v[0] = x;
     v[1] = y;
     v[2] = z;
 }
-double &vec3::operator[](int index) {
+scalar &vec3::operator[](int index) {
     return v[index];
 }
-double vec3::operator[](int index) const {
+scalar vec3::operator[](int index) const {
     return v[index];
 }
-vec3 vec3::operator*(double scale) const {
+vec3 vec3::operator*(scalar scale) const {
     return vec3(v[0] * scale, v[1] * scale, v[2] * scale);
 }
-vec3 vec3::operator/(double scale) const {
+vec3 vec3::operator/(scalar scale) const {
     return vec3(v[0] / scale, v[1] / scale, v[2] / scale);
 }
 vec3 vec3::operator+(const vec3 &other)  const{
@@ -42,13 +42,13 @@ vec3 vec3::operator-(const vec3 &other) const {
 vec3 vec3::operator-() const {
     return vec3(-v[0], -v[1], -v[2]);
 }
-const vec3 &vec3::operator*=(double scale) {
+const vec3 &vec3::operator*=(scalar scale) {
     v[0] *= scale;
     v[1] *= scale;
     v[2] *= scale;
     return *this;
 }
-const vec3 &vec3::operator/=(double scale) {
+const vec3 &vec3::operator/=(scalar scale) {
     v[0] /= scale;
     v[1] /= scale;
     v[2] /= scale;
@@ -66,17 +66,17 @@ const vec3 &vec3::operator-=(const vec3 &other) {
     v[2] -= other.v[2];
     return *this;
 }
-double vec3::magnitude() const {
+scalar vec3::magnitude() const {
     return sqrt(v[0] * v[0] + v[1] * v[1] + v[2] * v[2]);
 }
-double vec3::magnitudeSquared() const {
+scalar vec3::magnitudeSquared() const {
     return v[0] * v[0] + v[1] * v[1] + v[2] * v[2];
 }
 vec3 vec3::normalize() const {
-    double m = sqrt(v[0] * v[0] + v[1] * v[1] + v[2] * v[2]);
+    scalar m = sqrt(v[0] * v[0] + v[1] * v[1] + v[2] * v[2]);
     return vec3(v[0] / m, v[1] / m, v[2] / m);
 }
-double vec3::dot(const vec3 &other) const {
+scalar vec3::dot(const vec3 &other) const {
     return v[0] * other.v[0] + v[1] * other.v[1] + v[2] * other.v[2];
 }
 vec3 vec3::cross(const vec3 &other) const {
@@ -85,6 +85,10 @@ vec3 vec3::cross(const vec3 &other) const {
                  v[0] * other.v[1] - v[1] * other.v[0]);
 }
 std::ostream &operator<<(std::ostream &output, const vec3 &v) {
-    output << '(' << v[0] << ", " << v[1] << ", " << v[2] << ')';
-    return output;
+    return output << v[0] << " " << v[1] << " " << v[2];
+}
+
+vec3 operator*(scalar scale, const vec3 &v)
+{
+    return v * scale;
 }
diff --git a/vec3.h b/vec3.h
index abfc6c3..cb51f23 100644
--- a/vec3.h
+++ b/vec3.h
@@ -1,30 +1,34 @@
 #ifndef VEC3_H
 #define VEC3_H
 #include <iostream>
+
+typedef float scalar;
+
 class vec3 {
     public:
-        double v[3];
+        scalar v[3];
     public:
         vec3();
-        vec3(double x);
-        vec3(double x, double y, double z);
-        double &operator[](int index);
-        double operator[](int index) const;
-        vec3 operator*(double scale) const;
-        vec3 operator/(double scale) const;
+        vec3(scalar x);
+        vec3(scalar x, scalar y, scalar z);
+        scalar &operator[](int index);
+        scalar operator[](int index) const;
+        vec3 operator*(scalar scale) const;
+        vec3 operator/(scalar scale) const;
         vec3 operator+(const vec3 &other) const;
         vec3 operator-(const vec3 &other) const;
         vec3 operator-() const;
-        const vec3 &operator*=(double scale);
-        const vec3 &operator/=(double scale);
+        const vec3 &operator*=(scalar scale);
+        const vec3 &operator/=(scalar scale);
         const vec3 &operator+=(const vec3 &other);
         const vec3 &operator-=(const vec3 &other);
-        double magnitude() const;
-        double magnitudeSquared() const;
+        scalar magnitude() const;
+        scalar magnitudeSquared() const;
         vec3 normalize() const;
-        double dot(const vec3 &other) const;
+        scalar dot(const vec3 &other) const;
         vec3 cross(const vec3 &other) const;
 };
+vec3 operator*(scalar scale, const vec3 &v);
 
 std::ostream &operator<<(std::ostream &output, const vec3 &v);
 #endif