Work on adding toroid

1 files changed, 79 insertions, 14 deletions

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
 }