diff options
Diffstat (limited to 'test.cpp')
| -rw-r--r-- | test.cpp | 231 |
1 files changed, 0 insertions, 231 deletions
diff --git a/test.cpp b/test.cpp deleted file mode 100644 index eb5c72a..0000000 --- a/test.cpp +++ /dev/null @@ -1,231 +0,0 @@ -#include <cmath> -#include <fstream> -#include <iostream> -#include <sstream> -#include <sys/stat.h> -#include <sys/types.h> -#include <vector> - -#include "vec3.h" -#include "curve.h" - -using namespace std; - -/* A current or charge distribution */ -struct Entity { - enum { CHARGE, CURRENT } type; - union { - double Q_density; /* linear charge density */ - double 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; - -void add_current(scalar I, Curve *path) -{ - Entity n = { Entity::CURRENT, .I = I, .path = path }; - entities.push_back(n); -} - -void add_charge(scalar Q_density, Curve *path) -{ - Entity n = { Entity::CHARGE, .Q_density = Q_density, .path = path }; - entities.push_back(n); -} - -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 */ -void dump_field(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; - cout << 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; - } -} - -int main(int argc, char *argv[]) -{ - if(argc != 2) - return 1; - - Toroid toroid(vec3(0, 0, 0), vec3(1, 0, 0), vec3(0, 0, 1), 2 * M_PI, .1, atof(argv[1])); - Spiral solenoid(vec3(-1,0,0), vec3(0,1,0), vec3(1,0,0), 2*M_PI * 10, .2); - //LineSegment wire(vec3(0, -10, 0), vec3(0, 10, 0)); - - add_charge(1, (Curve*) &toroid); - add_current(1, (Curve*) &solenoid); - //add_current(1, (Curve*) &wire); - - dump_field(FieldType::E, - vec3(-3, -3, -1), - vec3(3, 3, 1), - .5); - - stringstream ss; - ss << "curve.fld"; - ofstream ofs(ss.str()); - - dump_paths(ofs, entities); - - ofs.close(); - - - //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 -} |