aboutsummaryrefslogtreecommitdiff
path: root/src/main.cpp
diff options
context:
space:
mode:
Diffstat (limited to 'src/main.cpp')
-rw-r--r--src/main.cpp423
1 files changed, 423 insertions, 0 deletions
diff --git a/src/main.cpp b/src/main.cpp
new file mode 100644
index 0000000..26ed520
--- /dev/null
+++ b/src/main.cpp
@@ -0,0 +1,423 @@
+#include <cmath>
+#include <cstdlib>
+#include <fstream>
+#include <iostream>
+#include <sstream>
+#include <sys/stat.h>
+#include <sys/types.h>
+#include <vector>
+
+#include <readline/readline.h>
+#include <readline/history.h>
+
+#include "gnuplot_i.hpp"
+
+#include "vec3.h"
+#include "curve.h"
+
+using namespace std;
+
+/* A current or charge distribution */
+struct Entity {
+ /* can bitwise-OR together */
+ enum { CHARGE = 1 << 0, CURRENT = 1 << 1 } type;
+ union {
+ scalar Q_density; /* linear charge density */
+ scalar I; /* current */
+ };
+
+ Curve *path;
+};
+
+vec3 point;
+
+/* dl x r / (|r| ^ 2) */
+vec3 dB(vec3 s, vec3 ds)
+{
+ vec3 r = point - s;
+
+ scalar r2 = r.magnitudeSquared();
+
+ vec3 rnorm = r / std::sqrt(r2);
+
+ return ds.cross(rnorm) / r2;
+}
+
+/* dl * r / (|r| ^ 2) */
+vec3 dE(vec3 s, vec3 ds)
+{
+ vec3 r = point - s;
+
+ scalar r2 = r.magnitudeSquared();
+
+ vec3 rnorm = r / std::sqrt(r2);
+
+ return rnorm * ds.magnitude() / r2;
+}
+
+vector<Entity> entities;
+
+int add_current(scalar I, Curve *path)
+{
+ Entity n = { Entity::CURRENT, I, path };
+ entities.push_back(n);
+
+ /* index */
+ return entities.size() - 1;
+}
+
+int add_charge(scalar Q_density, Curve *path)
+{
+ Entity n = { Entity::CHARGE, Q_density, path };
+ entities.push_back(n);
+
+ return entities.size() - 1;
+}
+
+const scalar U0 = 4e-7 * M_PI;
+const scalar C = 299792458;
+const scalar E0 = 1 / ( U0 * C * C );
+const scalar K_E = 1 / (4 * M_PI * E0);
+const scalar D = 1e-1;
+
+vec3 calc_Bfield(vec3 x)
+{
+ point = x;
+
+ vec3 B = 0;
+
+ for(int i = 0; i < entities.size(); i++)
+ {
+ if(entities[i].type == Entity::CURRENT)
+ B += entities[i].path->integrate(dB, D) * U0 * entities[i].I;
+ }
+
+ return B;
+}
+
+vec3 calc_Efield(vec3 x)
+{
+ point = x;
+
+ vec3 E = 0;
+
+ for(int i = 0; i < entities.size(); i++)
+ {
+ if(entities[i].type == Entity::CHARGE)
+ E += entities[i].path->integrate(dE, D) * K_E * entities[i].Q_density;
+ }
+
+ return E;
+}
+
+ostream *dump_ofs = NULL;
+
+vec3 dump(vec3 s, vec3 ds)
+{
+ *dump_ofs << s << " " << ds << endl;
+
+ return 0;
+}
+
+void dump_path(ostream &out, Curve *c)
+{
+ dump_ofs = &out;
+ c->integrate(dump, D);
+}
+
+int dump_entities(ostream &out, int which, vector<Entity> &e)
+{
+ int count = 0;
+ for(int i = 0; i < e.size(); i++)
+ {
+ if(which & e[i].type)
+ {
+ dump_path(out, e[i].path);
+
+ /* two blank lines mark an index in gnuplot */
+ out << endl << endl;
+
+ count++;
+ }
+ }
+ return count;
+}
+
+/* dump the field vectors with a spacing of `delta' */
+/* requires x0 < x1, y0 < y1, z0 < z1 */
+
+enum FieldType { E, B };
+
+/* dump field in a region of space to vectors */
+void dump_field(ostream &out,
+ enum FieldType type,
+ vec3 lower_corner, vec3 upper_corner,
+ scalar delta)
+{
+ for(scalar z = lower_corner[2]; z <= upper_corner[2]; z += delta)
+ for(scalar y = lower_corner[1]; y <= upper_corner[1]; y += delta)
+ for(scalar x = lower_corner[0]; x <= upper_corner[0]; x += delta)
+ {
+ vec3 pt(x, y, z);
+ vec3 field = (type == E) ? calc_Efield(pt) : calc_Bfield(pt);
+
+ field = field.normalize() / 10;
+ out << pt << " " << field << endl;
+ }
+}
+
+/* trace a field line */
+void dump_fieldline(ostream &out, vec3 x, scalar len)
+{
+ point = x;
+ scalar delta = .1;
+ while(len > 0)
+ {
+ out << point << endl;
+
+ vec3 B = calc_Bfield(point);
+
+ point += delta * B;
+ len -= delta;
+ }
+}
+
+/* dump field magnitudes along a line */
+void dump_values(vec3 start, vec3 del, int times)
+{
+ point = start;
+ while(times--)
+ {
+
+
+ point += del;
+ }
+}
+
+void all_lower(string &str)
+{
+ for(int i = 0; i < str.length(); i++)
+ str[i] = tolower(str[i]);
+}
+
+string itoa(int n)
+{
+ stringstream ss;
+ ss << n;
+ return ss.str();
+}
+
+Curve *parse_curve(stringstream &ss)
+{
+ string type;
+ ss >> type;
+ if(type == "line" || type == "linesegment")
+ {
+ vec3 a, b;
+ ss >> a >> b;
+ return (Curve*)new LineSegment(a, b);
+ }
+ else if(type == "arc")
+ {
+ vec3 center, radius, normal;
+ scalar angle;
+ ss >> center >> radius >> normal;
+ ss >> angle;
+ return (Curve*)new Arc(center, radius, normal, angle);
+ }
+ else if(type == "spiral" || type == "solenoid")
+ {
+ vec3 origin, radius, normal;
+ scalar angle, pitch;
+ ss >> origin >> radius >> normal >> angle >> pitch;
+ return (Curve*)new Spiral(origin, radius, normal, angle, pitch);
+ }
+ else if(type == "toroid")
+ {
+ vec3 origin, maj_radius, maj_normal;
+ scalar min_radius, maj_angle, pitch;
+ ss >> origin >> maj_radius >> maj_normal;
+ ss >> min_radius >> maj_angle >> pitch;
+
+ return (Curve*)new Toroid(origin, maj_radius, maj_normal, maj_angle, min_radius, pitch);
+ }
+ else throw "unknown curve type (must be line, arc, spiral, or toroid)";
+}
+
+void print_help()
+{
+ cout << endl;
+ cout << "fieldviz 0.1" << endl;
+ cout << "Copyright (C) 2019 Franklin Wei" << endl << endl;
+
+ cout << "Commands:" << endl;
+ cout << " add {I CURRENT|Q DENSITY} CURVE" << endl;
+ cout << " Add an entity of the specified type and the shape CURVE, where CURVE is one" << endl;
+ cout << " of (<X> is a 3-tuple specifying a vector):" << endl;
+ cout << " line <a> <b>" << endl;
+ cout << " arc <center> <radius> <normal> angle" << endl;
+ cout << " solenoid <center> <radius> <normal> angle pitch" << endl;
+ cout << " toroid <center> <radius> <maj_normal> min_radius maj_angle pitch" << endl;
+ cout << endl;
+ cout << " draw [I|Q] ..." << endl;
+ cout << " Draw the specified current/charge distributions" << endl;
+ cout << endl;
+ cout << " field [E|B] <lower_corner> <upper_corner> DELTA" << endl;
+ cout << " Plot the E or B field in the rectangular prism bounded by lower and upper." << endl;
+ cout << " DELTA specifies density." << endl;
+}
+
+int main(int argc, char *argv[])
+{
+ Toroid loop(vec3(0, 0, 0), vec3(1, 0, 0), vec3(0, 0, 1), M_PI * 2, .1, 2*M_PI / 10);
+ //add_current(1, (Curve*)&loop);
+
+ Gnuplot *gp = NULL;
+
+ try {
+ gp = new Gnuplot();
+ }
+ catch(GnuplotException e) {
+ Gnuplot::set_terminal_std("dumb");
+ gp = new Gnuplot();
+ }
+
+ *gp << "set view equal xyz";
+
+ cout << "Welcome to fieldviz!" << endl << endl;
+ cout << "Type `help' for a command listing." << endl;
+
+ while(1)
+ {
+ char *cs = readline("fieldviz> ");
+ if(!cs)
+ return 0;
+ add_history(cs);
+
+ string line(cs);
+
+ free(cs);
+
+ all_lower(line);
+
+ /* parse */
+ stringstream ss(line);
+
+ string cmd;
+ ss >> cmd;
+
+ try {
+ if(cmd == "add")
+ {
+ /* add a current or charge distribution */
+ Entity e;
+
+ string type;
+ ss >> type;
+
+ /* union */
+ double val;
+ ss >> val;
+
+ Curve *path = parse_curve(ss);
+
+ cout << "Curve type: " << path->name() << endl;
+
+ int idx;
+ if(type == "i")
+ idx = add_current(val, path);
+ else if(type == "q")
+ idx = add_charge(val, path);
+ else throw "unknown distribution type (must be I or Q)";
+
+ cout << "Index: " << idx << endl;
+ }
+ else if(cmd == "field")
+ {
+ string type;
+
+ vec3 lower, upper;
+ scalar delta;
+
+ if(!(ss >> type >> lower >> upper >> delta))
+ throw "plot requires <E/B> <lower> <upper> delta";
+
+ FieldType t = (type == "e") ? FieldType::E : FieldType::B;
+
+ ofstream out;
+ string fname = gp->create_tmpfile(out);
+
+ dump_field(out,
+ t,
+ lower, upper, delta);
+
+ out.close();
+
+ string cmd = "splot '" + fname + "' w vectors";
+ *gp << cmd;
+ }
+ else if(cmd == "draw")
+ {
+ int e_types = 0;
+
+ while(ss)
+ {
+ string type_str;
+ if(ss >> type_str)
+ {
+ if(type_str == "i")
+ e_types |= Entity::CURRENT;
+ else if(type_str == "q")
+ e_types |= Entity::CHARGE;
+ else
+ throw "unknown entity type (must be I or Q)";
+ }
+ }
+
+ if(!e_types)
+ e_types |= Entity::CHARGE | Entity::CURRENT;
+
+ ofstream out;
+ string fname = gp->create_tmpfile(out);
+ int n = dump_entities(out, e_types,
+ entities);
+ out.close();
+
+ string cmd = "splot for[i = 0:" + itoa(n - 1) + "] '" + fname + "' i i w lines";
+ *gp << cmd;
+ }
+ else if(cmd == "help")
+ print_help();
+ } catch(const char *err) {
+ cerr << "parse error: " << err << endl;
+ }
+ }
+
+ //LineSegment wire(vec3(0, -100, 0), vec3(0, 100, 0));
+
+ //std::cout << "length = " << loop.integrate(integrand, 1e-2) << endl;
+
+ //vec3 dx = .01;
+
+ //point = 0;
+
+ //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;
+
+#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
+}
generated by cgit v1.1 at 2026-04-24 19:44:24 +0000