initial commit

basic B field integrator

6 files changed, 205 insertions, 0 deletions

diff --git a/curve.cpp b/curve.cpp
new file mode 100644
index 0000000..9b63e8d
--- /dev/null
+++ b/curve.cpp
@@ -0,0 +1,23 @@
+#include <cmath>
+#include "curve.h"
+
+using namespace std;
+
+vec3 LineSegment::integrate(vec3 (*integrand)(vec3 s, vec3 ds), double dl)
+{
+    vec3 diff = this->b - this->a, sum = 0;
+
+    double len = diff.magnitude();
+
+    vec3 diffnorm = diff / len, s = this->a, ds = diffnorm * dl;
+
+    double 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;
+}
diff --git a/curve.h b/curve.h
new file mode 100644
index 0000000..c896a23
--- /dev/null
+++ b/curve.h
@@ -0,0 +1,28 @@
+#include <cmath>
+#include "vec3.h"
+using namespace std;
+
+class Curve {
+public:
+    virtual vec3 integrate(vec3 (*integrand)(vec3 s, vec3 ds), double delta) = 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), double delta);
+};
+
+class Arc : Curve {
+private:
+    vec3 center;
+    double radius;
+    double angle[2]; /* start and end angles */
+public:
+    LineSegment(vec3 a_, vec3 b_) : a(a_), b(b_) {};
+
+    vec3 integrate(vec3 (*integrand)(vec3 s, vec3 ds), double delta);
+};
diff --git a/main.cpp b/main.cpp
new file mode 100644
index 0000000..604782e
--- /dev/null
+++ b/main.cpp
@@ -0,0 +1 @@
+#include <iostream>
diff --git a/test.cpp b/test.cpp
new file mode 100644
index 0000000..8d20829
--- /dev/null
+++ b/test.cpp
@@ -0,0 +1,36 @@
+#include <cmath>
+#include <iostream>
+#include "curve.h"
+
+vec3 integrand(vec3 s, vec3 ds)
+{
+    return ds;
+}
+
+vec3 point;
+
+/* dl x r / (|r| ^ 2) */
+vec3 dB(vec3 s, vec3 ds)
+{
+    vec3 r = s - point;
+
+    double r2 = r.magnitudeSquared();
+
+    vec3 rnorm = r / std::sqrt(r2);
+
+    return ds.cross(rnorm) / r2;
+}
+
+int main()
+{
+    LineSegment wire(vec3(0, -100, 0), vec3(0, 100, 0));
+
+    vec3 dx = .001;
+
+    point = .01;
+
+    double I = 1;
+
+    for(int i = 0; i < 100; i++, point += dx)
+        std::cout << point[0] << " " << wire.integrate(dB, 1e-2)[2] << endl;
+}
diff --git a/vec3.cpp b/vec3.cpp
new file mode 100644
index 0000000..8513916
--- /dev/null
+++ b/vec3.cpp
@@ -0,0 +1,90 @@
+/* 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(double x) {
+    v[0] = x;
+    v[1] = 0;
+    v[2] = 0;
+}
+vec3::vec3(double x, double y, double z) {
+    v[0] = x;
+    v[1] = y;
+    v[2] = z;
+}
+double &vec3::operator[](int index) {
+    return v[index];
+}
+double vec3::operator[](int index) const {
+    return v[index];
+}
+vec3 vec3::operator*(double scale) const {
+    return vec3(v[0] * scale, v[1] * scale, v[2] * scale);
+}
+vec3 vec3::operator/(double 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*=(double scale) {
+    v[0] *= scale;
+    v[1] *= scale;
+    v[2] *= scale;
+    return *this;
+}
+const vec3 &vec3::operator/=(double 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;
+}
+double vec3::magnitude() const {
+    return sqrt(v[0] * v[0] + v[1] * v[1] + v[2] * v[2]);
+}
+double vec3::magnitudeSquared() const {
+    return v[0] * v[0] + v[1] * v[1] + v[2] * v[2];
+}
+vec3 vec3::normalize() const {
+    double 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);
+}
+double 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) {
+    output << '(' << v[0] << ", " << v[1] << ", " << v[2] << ')';
+    return output;
+}
diff --git a/vec3.h b/vec3.h
new file mode 100644
index 0000000..1f01098
--- /dev/null
+++ b/vec3.h
@@ -0,0 +1,27 @@
+#include <iostream>
+class vec3 {
+    public:
+        double v[3];
+    public:
+        vec3();
+        vec3(double x);
+        vec3(double x, double y, double z);
+        double &operator[](int index);
+        double operator[](int index) const;
+        vec3 operator*(double scale) const;
+        vec3 operator/(double scale) const;
+        vec3 operator+(const vec3 &other) const;
+        vec3 operator-(const vec3 &other) const;
+        vec3 operator-() const;
+        const vec3 &operator*=(double scale);
+        const vec3 &operator/=(double scale);
+        const vec3 &operator+=(const vec3 &other);
+        const vec3 &operator-=(const vec3 &other);
+        double magnitude() const;
+        double magnitudeSquared() const;
+        vec3 normalize() const;
+        double dot(const vec3 &other) const;
+        vec3 cross(const vec3 &other) const;
+};
+
+std::ostream &operator<<(std::ostream &output, const vec3 &v);