From 291bd26fd8920831181e8207e1fcdf544cd6cd6f Mon Sep 17 00:00:00 2001 From: Franklin Wei Date: Mon, 11 Feb 2019 12:52:45 -0500 Subject: Reorganize, use readline --- src/main.cpp | 423 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 423 insertions(+) create mode 100644 src/main.cpp (limited to 'src/main.cpp') 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 +#include +#include +#include +#include +#include +#include +#include + +#include +#include + +#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 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 &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 ( is a 3-tuple specifying a vector):" << endl; + cout << " line " << endl; + cout << " arc
angle" << endl; + cout << " solenoid
angle pitch" << endl; + cout << " toroid
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] 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 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 +} -- cgit v1.1