aboutsummaryrefslogtreecommitdiff
path: root/src
diff options
context:
space:
mode:
Diffstat (limited to 'src')
-rw-r--r--src/curve.cpp196
-rw-r--r--src/curve.h108
-rw-r--r--src/main.cpp423
-rw-r--r--src/quat.cpp38
-rw-r--r--src/quat.h25
-rw-r--r--src/vec3.cpp101
-rw-r--r--src/vec3.h36
7 files changed, 927 insertions, 0 deletions
diff --git a/src/curve.cpp b/src/curve.cpp
new file mode 100644
index 0000000..19ce1a8
--- /dev/null
+++ b/src/curve.cpp
@@ -0,0 +1,196 @@
+#include <iostream>
+#include <cmath>
+#include "curve.h"
+
+using namespace std;
+
+vec3 LineSegment::integrate(vec3 (*integrand)(vec3 s, vec3 ds), scalar dl) const
+{
+ vec3 diff = this->b - this->a, sum = 0;
+
+ scalar len = diff.magnitude();
+
+ vec3 diffnorm = diff / len, s = this->a, ds = diffnorm * dl;
+
+ scalar l;
+
+ for(l = 0; l < len; l += dl, s += ds)
+ sum += integrand(s, ds);
+
+ if(l < len)
+ sum += integrand(s, diffnorm * (len - l));
+
+ return sum;
+}
+
+vec3 Arc::integrate(vec3 (*integrand)(vec3 s, vec3 ds), scalar dl) const
+{
+ //cout << "Integrating loop of length " << this->angle << " radians" << endl;
+
+ scalar r = this->radius.magnitude(), dtheta = dl / r;
+
+ //cout << "R = " << r << ", dtheta = " << dtheta << endl;
+
+ quat rot = quat::from_angleaxis(dtheta, normal), crot = rot.conjugate();
+
+ //cout << "rot = " << rot << ", crot = " << crot << endl;
+
+ scalar len = this->angle * r, l;
+
+ vec3 sum = 0;
+
+ quat radius = this->radius;
+
+ for(l = 0; l < len; l += dl)
+ {
+ vec3 ds = this->normal.cross(radius).normalize() * dl;
+ sum += integrand(this->center + radius, ds);
+
+ radius = rot * radius * crot;
+ }
+
+ if(l < len)
+ {
+ dl = len - l;
+ dtheta = dl / r;
+ rot = quat::from_angleaxis(dtheta, normal);
+ crot = rot.conjugate();
+
+ vec3 ds = this->normal.cross(radius).normalize() * dl;
+ sum += integrand(this->center + radius, ds);
+ }
+
+ return sum;
+}
+
+vec3 Spiral::integrate(vec3 (*integrand)(vec3 s, vec3 ds), scalar dl) const
+{
+ //cout << "Integrating loop of length " << this->angle << " radians" << endl;
+
+ scalar r = this->radius.magnitude(), dtheta = dl / r;
+
+ //cout << "R = " << r << ", dtheta = " << dtheta << endl;
+
+ quat rot = quat::from_angleaxis(dtheta, normal), crot = rot.conjugate();
+
+ //cout << "rot = " << rot << ", crot = " << crot << endl;
+
+ /* only flat (doesn't include stretched length) */
+ scalar len = this->angle * r, l;
+
+ vec3 sum = 0;
+
+ quat radius = this->radius;
+
+ /* how far along the axis we've moved */
+ vec3 axis = 0;
+
+ /* we add this axial component each iteration */
+ vec3 dp = this->normal * dtheta * pitch / (2 * M_PI);
+
+ for(l = 0; l < len; l += dl, axis += dp)
+ {
+ vec3 ds = this->normal.cross(radius).normalize() * dl + dp;
+ sum += integrand(this->origin + axis + radius, ds);
+
+ /* rotate by dtheta about normal */
+ radius = rot * radius * crot;
+ }
+
+ 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);
+ }
+
+ 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 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;
+
+ /* 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;
+
+ //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;
+ 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 = 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
+ 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/src/curve.h b/src/curve.h
new file mode 100644
index 0000000..5f11351
--- /dev/null
+++ b/src/curve.h
@@ -0,0 +1,108 @@
+#ifndef CURVE_H
+#define CURVE_H
+
+#include <cmath>
+#include <iostream>
+#include "vec3.h"
+#include "quat.h"
+
+/* All curves inherit this class */
+class Curve {
+public:
+ virtual vec3 integrate(vec3 (*integrand)(vec3 s, vec3 ds), scalar delta) const = 0;
+ virtual const char *name() const = 0;
+};
+
+class LineSegment : Curve {
+private:
+ vec3 a, b;
+public:
+ LineSegment(vec3 a_, vec3 b_) : a(a_), b(b_) {};
+
+ vec3 integrate(vec3 (*integrand)(vec3 s, vec3 ds), scalar delta) const;
+
+ const char *name() const { return "LineSegment"; }
+
+ friend std::istream operator>>(std::istream &is, LineSegment &ls);
+};
+
+std::istream operator>>(std::istream &is, LineSegment &ls);
+
+class Arc : Curve {
+private:
+ vec3 center;
+
+ /* these are relative to the center (direction will be determined
+ * by RHR of normal), and should be orthonormal */
+ vec3 radius, normal;
+
+ /* how many radians the arc extends for (can be greater than 2pi) */
+ scalar angle;
+public:
+ 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), scalar delta) const;
+
+ const char *name() const { return "Arc"; }
+
+ friend std::istream operator>>(std::istream &is, LineSegment &ls);
+};
+
+std::istream operator>>(std::istream &is, LineSegment &ls);
+
+class Spiral : Curve {
+private:
+ vec3 origin;
+
+ /* these are relative to the center (direction will be determined
+ * by RHR of normal), and should be orthonormal */
+ vec3 radius, normal;
+
+ /* how many radians the arc extends for (can be greater than 2pi) */
+ scalar angle;
+
+ /* linear distance between turns (2pi) */
+ scalar pitch;
+public:
+ 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), scalar delta) const;
+
+ const char *name() const { return "Solenoid"; }
+
+ friend std::istream operator>>(std::istream &is, LineSegment &ls);
+};
+
+std::istream operator>>(std::istream &is, LineSegment &ls);
+
+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() {};
+ 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;
+
+ const char *name() const { return "Toroid"; };
+
+ friend std::istream operator>>(std::istream &is, LineSegment &ls);
+};
+
+std::istream operator>>(std::istream &is, LineSegment &ls);
+#endif
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
+}
diff --git a/src/quat.cpp b/src/quat.cpp
new file mode 100644
index 0000000..df117d0
--- /dev/null
+++ b/src/quat.cpp
@@ -0,0 +1,38 @@
+#include "quat.h"
+#include <cmath>
+
+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) { }
+
+quat::operator vec3()
+{
+ return vec3(this->x, this->y, this->z);
+}
+
+quat operator*(const quat &lhs, const quat &rhs)
+{
+ return quat(lhs.w * rhs.w - lhs.x * rhs.x - lhs.y * rhs.y - lhs.z * rhs.z,
+ lhs.w * rhs.x + lhs.x * rhs.w + lhs.y * rhs.z - lhs.z * rhs.y,
+ lhs.w * rhs.y - lhs.x * rhs.z + lhs.y * rhs.w + lhs.z * rhs.x,
+ lhs.w * rhs.z + lhs.x * rhs.y - lhs.y * rhs.x + lhs.z * rhs.w);
+}
+
+quat quat::conjugate() const
+{
+ return quat(this->w, -this->x, -this->y, -this->z);
+}
+
+quat quat::from_angleaxis(scalar angle, vec3 axis)
+{
+ scalar si = std::sin(angle / 2);
+ scalar co = std::cos(angle / 2);
+ return quat(co, si * axis[0], si * axis[1], si * axis[2]);
+}
+
+std::ostream &operator<<(std::ostream &os, const quat &q)
+{
+ return os << "(" << q.w << ", " << q.x << ", " << q.y << ", " << q.z << ")";
+}
diff --git a/src/quat.h b/src/quat.h
new file mode 100644
index 0000000..5821f81
--- /dev/null
+++ b/src/quat.h
@@ -0,0 +1,25 @@
+#ifndef QUAT_H
+#define QUAT_H
+#include "vec3.h"
+#include <iostream>
+
+class quat {
+public:
+ scalar w, x, y, z;
+public:
+ 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();
+
+ operator vec3();
+
+ quat conjugate() const;
+
+ static quat from_angleaxis(scalar angle, vec3 axis);
+};
+
+quat operator*(const quat &, const quat &);
+std::ostream &operator<<(std::ostream &os, const quat &);
+#endif
diff --git a/src/vec3.cpp b/src/vec3.cpp
new file mode 100644
index 0000000..e894a01
--- /dev/null
+++ b/src/vec3.cpp
@@ -0,0 +1,101 @@
+/* copy-pasted from:
+ * https://www.programming-techniques.com/2013/05/basic-euclidean-vector-operations-in-c.htm
+ */
+
+#include <iostream>
+#include <cmath>
+#include "vec3.h"
+using std::ostream;
+vec3::vec3() {
+ v[0] = 0;
+ v[1] = 0;
+ v[2] = 0;
+}
+vec3::vec3(scalar x) {
+ v[0] = x;
+ v[1] = 0;
+ v[2] = 0;
+}
+vec3::vec3(scalar x, scalar y, scalar z) {
+ v[0] = x;
+ v[1] = y;
+ v[2] = z;
+}
+scalar &vec3::operator[](int index) {
+ return v[index];
+}
+scalar vec3::operator[](int index) const {
+ return v[index];
+}
+vec3 vec3::operator*(scalar scale) const {
+ return vec3(v[0] * scale, v[1] * scale, v[2] * scale);
+}
+vec3 vec3::operator/(scalar scale) const {
+ return vec3(v[0] / scale, v[1] / scale, v[2] / scale);
+}
+vec3 vec3::operator+(const vec3 &other) const{
+ return vec3(v[0] + other.v[0], v[1] + other.v[1], v[2] + other.v[2]);
+}
+vec3 vec3::operator-(const vec3 &other) const {
+ return vec3(v[0] - other.v[0], v[1] - other.v[1], v[2] - other.v[2]);
+}
+vec3 vec3::operator-() const {
+ return vec3(-v[0], -v[1], -v[2]);
+}
+const vec3 &vec3::operator*=(scalar scale) {
+ v[0] *= scale;
+ v[1] *= scale;
+ v[2] *= scale;
+ return *this;
+}
+const vec3 &vec3::operator/=(scalar scale) {
+ v[0] /= scale;
+ v[1] /= scale;
+ v[2] /= scale;
+ return *this;
+}
+const vec3 &vec3::operator+=(const vec3 &other) {
+ v[0] += other.v[0];
+ v[1] += other.v[1];
+ v[2] += other.v[2];
+ return *this;
+}
+const vec3 &vec3::operator-=(const vec3 &other) {
+ v[0] -= other.v[0];
+ v[1] -= other.v[1];
+ v[2] -= other.v[2];
+ return *this;
+}
+scalar vec3::magnitude() const {
+ return sqrt(v[0] * v[0] + v[1] * v[1] + v[2] * v[2]);
+}
+scalar vec3::magnitudeSquared() const {
+ return v[0] * v[0] + v[1] * v[1] + v[2] * v[2];
+}
+vec3 vec3::normalize() const {
+ 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);
+}
+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 {
+ return vec3(v[1] * other.v[2] - v[2] * other.v[1],
+ v[2] * other.v[0] - v[0] * other.v[2],
+ v[0] * other.v[1] - v[1] * other.v[0]);
+}
+std::ostream &operator<<(std::ostream &output, const vec3 &v) {
+ return output << v[0] << " " << v[1] << " " << v[2];
+}
+
+std::istream &operator>>(std::istream &input, vec3 &v)
+{
+ if(!(input >> v[0] >> v[1] >> v[2]))
+ throw "error parsing vector";
+ return input;
+}
+
+vec3 operator*(scalar scale, const vec3 &v)
+{
+ return v * scale;
+}
diff --git a/src/vec3.h b/src/vec3.h
new file mode 100644
index 0000000..df68104
--- /dev/null
+++ b/src/vec3.h
@@ -0,0 +1,36 @@
+#ifndef VEC3_H
+#define VEC3_H
+#include <iostream>
+
+typedef float scalar;
+
+class vec3 {
+ public:
+ scalar v[3];
+ public:
+ vec3();
+ 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*=(scalar scale);
+ const vec3 &operator/=(scalar scale);
+ const vec3 &operator+=(const vec3 &other);
+ const vec3 &operator-=(const vec3 &other);
+ scalar magnitude() const;
+ scalar magnitudeSquared() const;
+ vec3 normalize() 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);
+std::istream &operator>>(std::istream &input, vec3 &v);
+
+#endif
generated by cgit v1.1 at 2026-04-24 23:55:08 +0000