path: root/main.cpp

#include <cmath>
#include <fstream>
#include <iostream>
#include <sstream>
#include <sys/stat.h>
#include <sys/types.h>
#include <vector>

#include "gnuplot_i.hpp"

#include "vec3.h"
#include "curve.h"

using namespace std;

/* A current or charge distribution */
struct Entity {
    enum { CHARGE, CURRENT } 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;
    //cout << "Magn: " << ds.magnitude() << endl;

    return 0;
}

void dump_path(ostream &out, Curve *c)
{
    dump_ofs = &out;
    c->integrate(dump, D);
}

void dump_paths(ostream &out, vector<Entity> &e)
{
    for(int i = 0; i < e.size(); i++)
    {
        dump_path(out, e[i].path);

        /* two blank lines mark an index in gnuplot */
        out << endl << endl;
    }
}

/* 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]);
}

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, min_radius, maj_angle, pitch);
    }
    else throw "unknown curve type (must be line, arc, spiral, or toroid)";
}

int main(int argc, char *argv[])
{
    Gnuplot gp;
    gp << "set view equal xyz";

    while(cin)
    {
        cout << "fieldviz> " << flush;
        string line;
        getline(cin, line);

        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: " << typeid(*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 == "plot")
            {
                string type;

                vec3 lower, upper;
                scalar delta;

                if(!(ss >> type >> lower >> upper >> delta))
                    throw "plot requires <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;
            }
        } 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
}