aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
authorFranklin Wei <me@fwei.tk>2019-05-30 23:03:17 -0400
committerFranklin Wei <me@fwei.tk>2019-05-30 23:04:12 -0400
commitbc6dcafc3868d55d2653081d27f1eaf771c2d532 (patch)
tree3fa44332d6981538247a85e0ad83d10c069ae431
parentcdfd5b37012935f7b0fb0a41ea8ca119ef8313b6 (diff)
downloadfieldviz-bc6dcafc3868d55d2653081d27f1eaf771c2d532.zip
fieldviz-bc6dcafc3868d55d2653081d27f1eaf771c2d532.tar.gz
fieldviz-bc6dcafc3868d55d2653081d27f1eaf771c2d532.tar.bz2
fieldviz-bc6dcafc3868d55d2653081d27f1eaf771c2d532.tar.xz
Generalize to 2-manifolds, refactor, improve
Moves some stuff to libfml
-rw-r--r--CMakeLists.txt2
-rw-r--r--src/curve.cpp24
-rw-r--r--src/curve.h19
-rw-r--r--src/main.cpp195
-rw-r--r--src/manifold.h17
-rw-r--r--src/surface.cpp149
-rw-r--r--src/surface.h121
7 files changed, 465 insertions, 62 deletions
diff --git a/CMakeLists.txt b/CMakeLists.txt
index cd5836a..a99d693 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -1,6 +1,6 @@
cmake_minimum_required (VERSION 2.6)
project (fieldviz)
-add_executable(fieldviz src/main.cpp src/curve.cpp)
+add_executable(fieldviz src/main.cpp src/curve.cpp src/surface.cpp)
add_definitions(-std=c++14 -O2 -g)
diff --git a/src/curve.cpp b/src/curve.cpp
index 19ce1a8..83bec06 100644
--- a/src/curve.cpp
+++ b/src/curve.cpp
@@ -31,7 +31,7 @@ vec3 Arc::integrate(vec3 (*integrand)(vec3 s, vec3 ds), scalar dl) const
//cout << "R = " << r << ", dtheta = " << dtheta << endl;
- quat rot = quat::from_angleaxis(dtheta, normal), crot = rot.conjugate();
+ quat rot = quat::from_angleaxis(dtheta, normal);
//cout << "rot = " << rot << ", crot = " << crot << endl;
@@ -46,15 +46,13 @@ vec3 Arc::integrate(vec3 (*integrand)(vec3 s, vec3 ds), scalar dl) const
vec3 ds = this->normal.cross(radius).normalize() * dl;
sum += integrand(this->center + radius, ds);
- radius = rot * radius * crot;
+ radius = radius.rotateby(rot);
}
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);
@@ -71,7 +69,7 @@ vec3 Spiral::integrate(vec3 (*integrand)(vec3 s, vec3 ds), scalar dl) const
//cout << "R = " << r << ", dtheta = " << dtheta << endl;
- quat rot = quat::from_angleaxis(dtheta, normal), crot = rot.conjugate();
+ quat rot = quat::from_angleaxis(dtheta, normal);
//cout << "rot = " << rot << ", crot = " << crot << endl;
@@ -94,15 +92,13 @@ vec3 Spiral::integrate(vec3 (*integrand)(vec3 s, vec3 ds), scalar dl) const
sum += integrand(this->origin + axis + radius, ds);
/* rotate by dtheta about normal */
- radius = rot * radius * crot;
+ radius = radius.rotateby(rot);
}
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;
@@ -137,9 +133,9 @@ vec3 Toroid::integrate(vec3 (*integrand)(vec3 s, vec3 ds), scalar dl) const
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 minor_rot = quat::from_angleaxis(minor_dtheta, minor_normal);
- quat major_rot = quat::from_angleaxis(major_dtheta, major_normal), major_crot = major_rot.conjugate();
+ quat major_rot = quat::from_angleaxis(major_dtheta, major_normal);
vec3 sum = 0;
@@ -163,13 +159,12 @@ vec3 Toroid::integrate(vec3 (*integrand)(vec3 s, vec3 ds), scalar dl) const
* 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;
+ major_radius = major_radius.rotateby(major_rot);
+ minor_normal = minor_normal.rotateby(major_rot);
- minor_radius = minor_rot * minor_radius * minor_crot;
+ minor_radius = minor_radius.rotateby(minor_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();
@@ -193,4 +188,3 @@ vec3 Toroid::integrate(vec3 (*integrand)(vec3 s, vec3 ds), scalar dl) const
return sum;
}
-
diff --git a/src/curve.h b/src/curve.h
index 891a6af..a3d743a 100644
--- a/src/curve.h
+++ b/src/curve.h
@@ -4,16 +4,19 @@
#include <cmath>
#include <iostream>
#include <fml/fml.h>
+
+#include "manifold.h"
+
using namespace fml;
-/* All curves inherit this class */
-class Curve {
+/* All curves inherit this class (which is empty because Manifold has
+ * everything we need). */
+class Curve : public Manifold {
public:
- virtual vec3 integrate(vec3 (*integrand)(vec3 s, vec3 ds), scalar delta) const = 0;
- virtual const char *name() const = 0;
+ const int dimension() const { return 1; }
};
-class LineSegment : Curve {
+class LineSegment : public Curve {
private:
vec3 a, b;
public:
@@ -28,7 +31,7 @@ public:
std::istream operator>>(std::istream &is, LineSegment &ls);
-class Arc : Curve {
+class Arc : public Curve {
private:
vec3 center;
@@ -50,7 +53,7 @@ public:
std::istream operator>>(std::istream &is, LineSegment &ls);
-class Spiral : Curve {
+class Spiral : public Curve {
private:
vec3 origin;
@@ -75,7 +78,7 @@ public:
std::istream operator>>(std::istream &is, LineSegment &ls);
-class Toroid : Curve {
+class Toroid : public Curve {
private:
vec3 origin;
diff --git a/src/main.cpp b/src/main.cpp
index b6b6c10..f56209f 100644
--- a/src/main.cpp
+++ b/src/main.cpp
@@ -1,11 +1,12 @@
#include <cmath>
+#include <csignal>
#include <cstdlib>
#include <fstream>
#include <iostream>
+#include <map>
#include <sstream>
#include <sys/stat.h>
#include <sys/types.h>
-#include <vector>
#include <readline/readline.h>
#include <readline/history.h>
@@ -13,8 +14,12 @@
#include "gnuplot_i.hpp"
#include "curve.h"
+#include "surface.h"
#include <fml/fml.h>
+// under $HOME
+#define HISTORY_FILE ".fieldviz_history"
+
using namespace fml;
using namespace std;
@@ -27,7 +32,7 @@ struct Entity {
scalar I; /* current */
};
- Curve *path;
+ Manifold *path;
};
vec3 point;
@@ -56,30 +61,31 @@ vec3 dE(vec3 s, vec3 ds)
return rnorm * ds.magnitude() / r2;
}
-vector<Entity> entities;
+int ent_counter = 0;
+map<int, Entity> entities;
-int add_current(scalar I, Curve *path)
+int add_entity(Entity e)
{
- Entity n = { Entity::CURRENT, I, path };
- entities.push_back(n);
-
- /* index */
- return entities.size() - 1;
+ entities[ent_counter] = e;
+ return ent_counter++;
}
-int add_charge(scalar Q_density, Curve *path)
+int add_current(scalar I, Manifold *path)
{
- Entity n = { Entity::CHARGE, Q_density, path };
- entities.push_back(n);
+ return add_entity((Entity){ Entity::CURRENT, I, path });
+}
- return entities.size() - 1;
+int add_charge(scalar Q_density, Manifold *path)
+{
+ return add_entity((Entity){ Entity::CHARGE, Q_density, path });
}
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-2;
+const scalar DEFAULT_D = 1e-1;
+scalar D = DEFAULT_D;
vec3 calc_Bfield(vec3 x)
{
@@ -87,10 +93,11 @@ vec3 calc_Bfield(vec3 x)
vec3 B = 0;
- for(int i = 0; i < entities.size(); i++)
+ for(map<int, Entity>::iterator i = entities.begin(); i != entities.end(); i++)
{
- if(entities[i].type == Entity::CURRENT)
- B += entities[i].path->integrate(dB, D) * U0 * entities[i].I;
+ Entity &e = i->second;
+ if(e.type == Entity::CURRENT)
+ B += e.path->integrate(dB, D) * U0 * e.I;
}
return B;
@@ -102,10 +109,11 @@ vec3 calc_Efield(vec3 x)
vec3 E = 0;
- for(int i = 0; i < entities.size(); i++)
+ for(map<int, Entity>::iterator i = entities.begin(); i != entities.end(); i++)
{
- if(entities[i].type == Entity::CHARGE)
- E += entities[i].path->integrate(dE, D) * K_E * entities[i].Q_density;
+ Entity &e = i->second;
+ if(e.type == Entity::CHARGE)
+ E += e.path->integrate(dE, D) * K_E * e.Q_density;
}
return E;
@@ -120,20 +128,21 @@ vec3 dump(vec3 s, vec3 ds)
return 0;
}
-void dump_path(ostream &out, Curve *c)
+void dump_points(ostream &out, Manifold *c)
{
dump_ofs = &out;
c->integrate(dump, D);
}
-int dump_entities(ostream &out, int which, vector<Entity> &e)
+int dump_entities(ostream &out, int which, map<int, Entity> &en)
{
int count = 0;
- for(int i = 0; i < e.size(); i++)
+ for(map<int, Entity>::iterator i = en.begin(); i != en.end(); i++)
{
- if(which & e[i].type)
+ Entity &e = i->second;
+ if(which & e.type)
{
- dump_path(out, e[i].path);
+ dump_points(out, e.path);
/* two blank lines mark an index in gnuplot */
out << endl << endl;
@@ -141,6 +150,7 @@ int dump_entities(ostream &out, int which, vector<Entity> &e)
count++;
}
}
+
return count;
}
@@ -189,8 +199,6 @@ void dump_values(vec3 start, vec3 del, int times)
point = start;
while(times--)
{
-
-
point += del;
}
}
@@ -208,7 +216,7 @@ string itoa(int n)
return ss.str();
}
-Curve *parse_curve(stringstream &ss)
+Manifold *parse_curve(stringstream &ss)
{
string type;
ss >> type;
@@ -242,6 +250,41 @@ Curve *parse_curve(stringstream &ss)
return (Curve*)new Toroid(origin, maj_radius, maj_normal, maj_angle, min_radius, pitch);
}
+ else if(type == "plane")
+ {
+ vec3 origin, v1, v2;
+ ss >> origin >> v1 >> v2;
+ return (Surface*)new Plane(origin, v1, v2);
+ }
+ else if(type == "disk")
+ {
+ vec3 center, radius, normal;
+ scalar angle;
+ ss >> center >> radius >> normal >> angle;
+ return (Surface*)new Disk(center, radius, normal, angle);
+ }
+ else if(type == "sphere")
+ {
+ vec3 center;
+ scalar radius;
+ ss >> center >> radius;
+
+ return (Surface*)new Sphere(center, radius);
+ }
+ else if(type == "opencylinder")
+ {
+ vec3 origin, axis;
+ scalar rad;
+ ss >> origin >> axis >> rad;
+ return (Surface*)new OpenCylinder(origin, axis, rad);
+ }
+ else if(type == "closedcylinder")
+ {
+ vec3 origin, axis;
+ scalar rad;
+ ss >> origin >> axis >> rad;
+ return (Surface*)new ClosedCylinder(origin, axis, rad);
+ }
else throw "unknown curve type (must be line, arc, spiral, or toroid)";
}
@@ -252,13 +295,21 @@ void print_help()
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 << " add {I CURRENT|Q DENSITY} MANIFOLD" << endl;
+ cout << " Add an entity of the specified type and the shape MANIFOLD, which can be:" << endl;
cout << " of (<X> is a 3-tuple specifying a vector):" << endl;
+ cout << " 1-manifolds:" << 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 << " 2-manifolds:" << endl;
+ cout << " plane <origin> <vec1> <vec2>" << endl;
+ cout << " disk <center> <radius> <normal> angle" << endl;
+ cout << " sphere <center> radius" << endl;
+ cout << endl;
+ cout << " delete [ID..]" << endl;
+ cout << " Delete an entity by its previously returned identifier." << endl;
cout << endl;
cout << " draw [I|Q] ..." << endl;
cout << " Draw the specified current/charge distributions" << endl;
@@ -266,16 +317,51 @@ void print_help()
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;
+ cout << endl;
+ cout << " newwindow" << endl;
+ cout << " Make future plots go into a new window" << endl;
+ cout << endl;
+ cout << " delta D" << endl;
+ cout << " Set integration fineness to D (smaller is better but slower)" << endl;
+}
+
+vec3 dA(vec3 s, vec3 dA)
+{
+ /* will cast to vec3 */
+ return dA.magnitude();
+}
+
+string hist_path;
+
+void exit_handler()
+{
+ write_history(hist_path.c_str());
+}
+
+void int_handler(int a)
+{
+ (void) a;
+ exit(0);
}
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);
+ Surface *surf = new Sphere(vec3(0, 0, 1), 1);
+ cout << "Area of 10x10 square = " << surf->integrate(dA, D) << endl;
+
+ hist_path = getenv("HOME");
+ hist_path += "/";
+ hist_path += HISTORY_FILE;
+
+ using_history();
+ read_history(hist_path.c_str());
+ atexit(exit_handler);
+ signal(SIGINT, int_handler);
Gnuplot *gp = NULL;
try {
+ Gnuplot::set_terminal_std("qt");
gp = new Gnuplot();
}
catch(GnuplotException e) {
@@ -288,17 +374,20 @@ int main(int argc, char *argv[])
cout << "Welcome to fieldviz!" << endl << endl;
cout << "Type `help' for a command listing." << endl;
+ /* will change to replot on subsequent commands */
+ string plot_cmd = "splot";
+
while(1)
{
char *cs = readline("fieldviz> ");
if(!cs)
return 0;
add_history(cs);
-
+
string line(cs);
free(cs);
-
+
all_lower(line);
/* parse */
@@ -320,9 +409,9 @@ int main(int argc, char *argv[])
double val;
ss >> val;
- Curve *path = parse_curve(ss);
+ Manifold *path = parse_curve(ss);
- cout << "Curve type: " << path->name() << endl;
+ cout << "Manifold type: " << path->name() << endl;
int idx;
if(type == "i")
@@ -333,6 +422,17 @@ int main(int argc, char *argv[])
cout << "Index: " << idx << endl;
}
+ else if(cmd == "delete")
+ {
+ int id;
+ while(ss >> id)
+ {
+ if(entities.erase(id))
+ cout << "Deleted " << id << "." << endl;
+ else
+ cerr << "No entity " << id << "!" << endl;
+ }
+ }
else if(cmd == "field")
{
string type;
@@ -354,8 +454,10 @@ int main(int argc, char *argv[])
out.close();
- string cmd = "splot '" + fname + "' w vectors";
+ string cmd = plot_cmd + " '" + fname + "' w vectors";
*gp << cmd;
+
+ plot_cmd = "replot";
}
else if(cmd == "draw")
{
@@ -384,11 +486,28 @@ int main(int argc, char *argv[])
entities);
out.close();
- string cmd = "splot for[i = 0:" + itoa(n - 1) + "] '" + fname + "' i i w lines";
+ string cmd = plot_cmd + " for[i = 0:" + itoa(n - 1) + "] '" + fname + "' i i w vectors";
*gp << cmd;
+
+ plot_cmd = "replot";
+ }
+ else if(cmd == "delta")
+ {
+ ss >> D;
+ if(D <= 0)
+ {
+ cerr << "D must be positive and non-zero!" << endl;
+ D = DEFAULT_D;
+ }
+ }
+ else if(cmd == "newwindow")
+ {
+ plot_cmd = "splot";
}
else if(cmd == "help")
print_help();
+ else
+ cerr << "invalid" << endl;
} catch(const char *err) {
cerr << "parse error: " << err << endl;
}
diff --git a/src/manifold.h b/src/manifold.h
new file mode 100644
index 0000000..f7fad5b
--- /dev/null
+++ b/src/manifold.h
@@ -0,0 +1,17 @@
+#ifndef MANIFOLD_H
+#define MANIFOLD_H
+
+#include <cmath>
+#include <iostream>
+#include <fml/fml.h>
+using namespace fml;
+
+/* All manifolds inherit this class */
+class Manifold {
+public:
+ virtual vec3 integrate(vec3 (*integrand)(vec3 s, vec3 ds), scalar delta) const = 0;
+ virtual const char *name() const = 0;
+ virtual const int dimension() const = 0; // 0 = point, 1 = curve, 2 = surface, 3 = solid
+};
+
+#endif
diff --git a/src/surface.cpp b/src/surface.cpp
new file mode 100644
index 0000000..8ab52d4
--- /dev/null
+++ b/src/surface.cpp
@@ -0,0 +1,149 @@
+#include <iostream>
+#include <cmath>
+#include "surface.h"
+
+vec3 Plane::integrate(vec3 (*integrand)(vec3 s, vec3 dA), scalar d) const
+{
+ vec3 sum = 0;
+
+ vec3 dA = (d * this->v1).cross(d * this->v2);
+
+ for(scalar s = 0; s < 1; s += d)
+ for(scalar t = 0; t < 1; t += d)
+ {
+ vec3 p = this->p0 + s * this->v1 + t * this->v2;
+ sum += integrand(p, dA);
+ }
+ return sum;
+}
+
+vec3 Disk::integrate(vec3 (*integrand)(vec3 s, vec3 dA), scalar dr) const
+{
+ vec3 sum = 0;
+
+ scalar radius = this->radius.magnitude();
+ vec3 radnorm = this->radius.normalize();
+
+ /* chosen so that the outermost ring will consist of square area
+ * elements */
+ scalar dtheta = dr / radius;
+
+ quat rot = quat::from_angleaxis(dtheta, this->normal);
+
+ for(scalar r = 0; r < radius; r += dr)
+ {
+ vec3 s = this->center + radnorm * r;
+
+ /* area element is constant for given r */
+ vec3 dA = this->normal * r * dr * dtheta;
+
+ for(scalar theta = 0; theta < this->angle; theta += dtheta)
+ {
+ sum += integrand(s, dA);
+
+ s = s.rotateby(rot);
+ }
+ }
+
+ return sum;
+}
+
+vec3 Sphere::integrate(vec3 (*integrand)(vec3 s, vec3 dA), scalar d) const
+{
+ vec3 sum = 0;
+ /*
+ * Coordinate reference (right-handed):
+ *
+ * ^ z
+ * |
+ * | ^ y
+ * | /
+ * | /
+ * | /
+ * |/
+ * O---------> x
+ *
+ */
+
+ /* we will rotate this unit vector clockwise in the X-Z plane by d
+ * each outer loop (scale and offset later) */
+ vec3 rad = vec3(0, 0, 1.0);
+
+ scalar r_sq = this->radius * this->radius;
+
+ quat roty = quat::from_angleaxis(d, vec3(0, 1, 0));
+ quat rotz = quat::from_angleaxis(d, vec3(0, 0, 1));
+
+ for(scalar phi = 0; phi < M_PI; phi += d)
+ {
+ /* operate on a copy to avoid accumulating roundoff error */
+ vec3 rad2 = rad;
+
+ /* from Jacobian: dA = r^2 * sin(phi) * dphi * dtheta */
+ scalar dA = r_sq * sin(phi) * d * d;
+
+ for(scalar theta = 0; theta < 2*M_PI; theta += d)
+ {
+ sum += integrand(this->center + this->radius * rad2, dA * rad2);
+
+ /* rotate rad2 around the z axis */
+ rad2 = rad2.rotateby(rotz);
+ }
+
+ /* rotate radius clockwise around y */
+ rad = rad.rotateby(roty);
+ }
+
+ return sum;
+}
+
+vec3 OpenCylinder::integrate(vec3 (*integrand)(vec3 s, vec3 dA), scalar d) const
+{
+ /*
+ * We loop along the axis length, rotating a vector around it as
+ * we go.
+ *
+ * Offset is trivial.
+ *
+ * rad/v
+ * ^ x
+ * .|x . x
+ * . | x . x
+ * O-------x-------------> axis
+ * x . x . x .
+ * x . x . x .
+ * x x x
+ *
+ */
+
+ /* TODO: normalize d so all integrals run in 1/d^n time? */
+
+ vec3 v = vec3::any_unit_normal(this->axis);
+ quat rot = quat::from_angleaxis(d, this->axis);
+
+ scalar axis_len = this->axis.magnitude();
+ vec3 norm_axis = this->axis.normalize();
+
+ vec3 sum = 0;
+
+ for(scalar l = 0; l < axis_len; l += d)
+ {
+ vec3 rad = v;
+ for(scalar theta = 0; theta < 2*M_PI; theta += d)
+ {
+ sum += integrand(this->origin + l * norm_axis + this->radius * rad, rad);
+ rad = rad.rotateby(rot);
+ }
+ }
+
+ return sum;
+}
+
+vec3 ClosedCylinder::integrate(vec3 (*integrand)(vec3 s, vec3 dA), scalar d) const
+{
+ vec3 sum = OpenCylinder::integrate(integrand, d);
+
+ sum += cap1.integrate(integrand, d) + cap2.integrate(integrand, d);
+
+ return sum;
+}
diff --git a/src/surface.h b/src/surface.h
new file mode 100644
index 0000000..4726521
--- /dev/null
+++ b/src/surface.h
@@ -0,0 +1,121 @@
+#ifndef SURFACE_H
+#define SURFACE_H
+
+#include <fml/fml.h>
+#include "manifold.h"
+
+using namespace fml;
+
+/* All surfaces inherit this class */
+/* The exact meaning of d is surface-dependent (see class
+ * definitions below); the limit of the calculated integral must
+ * approach its true value as d->0. d must be > 0. */
+class Surface : public Manifold {
+public:
+ const int dimension() const { return 2; }
+};
+
+class Plane : public Surface {
+private:
+ /* The surface is specified by all points p = p0 + s v1 + t v2,
+ * such that 0 <= {s, t} < 1
+ *
+ * v1 and v2 must NOT be parallel.
+ *
+ * v1 and v2 should (but do not have to be) be normal (this is
+ * motivated primarily by usability; if the two vectors are indeed
+ * normal, then m1 and m2 have the nice geometric meaning of side
+ * length).
+ *
+ * d = ds = dt.
+ * dA will be in the direction of v1 x v2.
+ */
+ vec3 p0, v1, v2;
+
+public:
+ Plane(vec3 p, vec3 _v1, vec3 _v2) : p0(p), v1(_v1), v2(_v2) {};
+
+ vec3 integrate(vec3 (*integrand)(vec3 s, vec3 dA), scalar d) const;
+ const char *name() const { return "Plane"; }
+};
+
+/* Flat, circular disk (of varying extent) */
+class Disk : public Surface {
+private:
+ /* This represents a (possibly incomplete) circular disk in space,
+ * with a center, radius, and normal vector as specified.
+ *
+ * `angle' specifies the extent of the disk; angle=2pi for a full
+ * circle.
+ *
+ * We use:
+ * d = dr.
+ * d_theta = d / r.
+ *
+ * This makes it so that differential area elements on the edge of
+ * the disk are square.
+ *
+ * dA will be in the direction of normal.
+ */
+ vec3 center, radius, normal;
+ scalar angle;
+
+public:
+ Disk(vec3 c, vec3 r, vec3 n, scalar a) : center(c), radius(r), normal(n), angle(a) {};
+
+ vec3 integrate(vec3 (*integrand)(vec3 s, vec3 dA), scalar d) const;
+ const char *name() const { return "Disk"; }
+};
+
+/* Hollow, spherical shell */
+class Sphere : public Surface {
+private:
+ /* d = dtheta = dphi */
+
+ vec3 center;
+ scalar radius;
+
+public:
+ Sphere(vec3 c, scalar r) : center(c), radius(r) {};
+
+ vec3 integrate(vec3 (*integrand)(vec3 s, vec3 dA), scalar d) const;
+ const char *name() const { return "Sphere"; }
+};
+
+/* Cylinder without end caps */
+class OpenCylinder : public Surface {
+private:
+ /* Like this:
+ *
+ * ___________________
+ * / \ \
+ * origin ---------------> axis
+ * \_/_________________/
+ *
+ */
+
+ vec3 origin, axis;
+ scalar radius;
+
+public:
+ OpenCylinder(vec3 o, vec3 a, scalar r) : origin(o), axis(a), radius(r) {}
+
+ vec3 integrate(vec3 (*integrand)(vec3 s, vec3 dA), scalar d) const;
+ const char *name() const { return "OpenCylinder"; }
+};
+
+/* Capped cylinder */
+class ClosedCylinder : public OpenCylinder {
+private:
+ Disk cap1, cap2;
+
+public:
+ ClosedCylinder(vec3 o, vec3 a, scalar r) : OpenCylinder(o, a, r),
+ cap1(o, vec3::any_unit_normal(a), -a.normalize(), 2*M_PI),
+ cap2(o + a, vec3::any_unit_normal(a), a.normalize(), 2*M_PI) {}
+
+ vec3 integrate(vec3 (*integrand)(vec3 s, vec3 dA), scalar d) const;
+ const char *name() const { return "ClosedCylinder"; }
+};
+
+#endif