Finish up toroid

Still isn't perfect (pitch parameter doesn't work well)

4 files changed, 48 insertions, 17 deletions

diff --git a/anim.sh b/anim.sh
new file mode 100755
index 0000000..7561db9
--- /dev/null
+++ b/anim.sh
@@ -0,0 +1,12 @@
+#!/bin/bash
+rm field -rf
+mkdir field
+pitches="0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40 0.45 0.50 0.55 0.60 0.65 0.70 0.75 0.80 0.85 0.90 0.95 1.00"
+for f in $pitches
+do
+    ./a.out $f > field/$f
+done
+for f in field/*
+do
+    cat $f | gnuplot -e "set terminal png;set xrange[-2:2];set yrange[-2:2];set zrange[-.2:.2];sp '-' w lines;" > $f.png
+done
diff --git a/curve.cpp b/curve.cpp
index c71a8c3..19ce1a8 100644
--- a/curve.cpp
+++ b/curve.cpp
@@ -120,8 +120,10 @@ vec3 Toroid::integrate(vec3 (*integrand)(vec3 s, vec3 ds), scalar dl) const
     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;
-    
+    scalar pitch_length = this->pitch * major_r;
+
+    scalar dpar = this->pitch * dl / (2 * M_PI * sqrt(minor_r * minor_r + pitch_length * pitch_length)) * major_r;
+
     /* how far each dl rotates the major vector */
     scalar major_dtheta = dpar / major_r;
 
@@ -144,6 +146,12 @@ vec3 Toroid::integrate(vec3 (*integrand)(vec3 s, vec3 ds), scalar dl) const
     quat major_radius = this->major_radius;
     quat minor_radius = this->major_radius.normalize() * this->minor_r;
 
+    //cout << "minor_dtheta = " << minor_dtheta << endl;
+    //cout << "minor_dper = " << dper << endl;
+    //cout << "dl = " << dl << endl;
+    //cout << "dpar = " << dpar << endl;
+    //cout << "pitch = " << this->pitch << endl;
+
     for(scalar theta = 0; theta < this->major_angle; theta += major_dtheta)
     {
         vec3 dpar_vec = minor_normal * dpar;
@@ -157,11 +165,16 @@ vec3 Toroid::integrate(vec3 (*integrand)(vec3 s, vec3 ds), scalar dl) const
          * 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_rot = quat::from_angleaxis(minor_dtheta, minor_normal);
         minor_crot = minor_rot.conjugate();
+        //minor_rot = major_rot * minor_rot;
+        //minor_crot = minor_rot.conjugate();
+
+        //cout << "Minor radius is " << minor_radius << endl;
+        //cout << "Minor rot quat is " << minor_rot << endl;
     }
 
 #if 0
diff --git a/curve.h b/curve.h
index 52c7623..19c41c7 100644
--- a/curve.h
+++ b/curve.h
@@ -75,6 +75,7 @@ private:
      * radius vector) */
     scalar pitch;
 public:
+    Toroid() {};
     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;
diff --git a/test.cpp b/test.cpp
index 0a42f40..f9d2cdf 100644
--- a/test.cpp
+++ b/test.cpp
@@ -34,15 +34,15 @@ 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(-.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);
+Toroid loop;
 
 ostream *dump_ofs = NULL;
 
 vec3 dump(vec3 s, vec3 ds)
 {
-    *dump_ofs << s << endl;
+    *dump_ofs << s << " " << ds << endl;
     //cout << "Magn: " << ds.magnitude() << endl;
-        
+
     return 0;
 }
 
@@ -57,9 +57,8 @@ void dump_path(ostream &out, const Curve *c)
 
 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)
+void dump_field(scalar x0, scalar y0, scalar z0, scalar x1, scalar y1, scalar z1, scalar delta)
 {
     for(scalar z = z0; z <= z1; z += delta)
         for(scalar y = y0; y <= y1; y += delta)
@@ -93,8 +92,13 @@ void dump_fieldline(ostream &out, vec3 x, scalar len)
     }
 }
 
-int main()
+int main(int argc, char *argv[])
 {
+    if(argc != 2)
+        return 1;
+    loop = Toroid(vec3(0, 0, 0), vec3(1, 0, 0), vec3(0, 0, 1), 1 * M_PI, .1, atof(argv[1]));
+    
+    
     //LineSegment wire(vec3(0, -100, 0), vec3(0, 100, 0));
 
     //std::cout << "length = " << loop.integrate(integrand, 1e-2) << endl;
@@ -108,17 +112,18 @@ int main()
     //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(-2, -2, 0,
-                12, 2, 0);
+    //dump_field(-2, -2, 0,
+//             2, 2, 0,
+    //.1);
     //dump_field(0,0,0,0,0,0);
 
-    stringstream ss;
-    ss << "curve.fld";
-    ofstream ofs(ss.str());
+    //stringstream ss;
+    //ss << "curve.fld";
+    //ofstream ofs(ss.str());
 
-    dump_path(ofs, (Curve*)&loop);
+    dump_path(cout, (Curve*)&loop);
     
-    ofs.close();
+    //ofs.close();
     
     
 #if 0