1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
|
#include <cmath>
#include <iostream>
#include <fstream>
#include <sys/stat.h>
#include <sys/types.h>
#include <sstream>
#include "vec3.h"
#include "curve.h"
using namespace std;
vec3 integrand(vec3 s, vec3 ds)
{
cout << "point " << s << "ds = " << ds << endl;
return ds;
}
vec3 point;
/* dl x r / (|r| ^ 2) */
vec3 dB(vec3 s, vec3 ds)
{
vec3 r = s - point;
scalar r2 = r.magnitudeSquared();
vec3 rnorm = r / std::sqrt(r2);
return ds.cross(rnorm) / r2;
}
//Arc loop(vec3(0, 0, 0), vec3(0, 1, 0), vec3(1, 0, 0), M_PI * 2);
//Spiral loop(vec3(0, 0, 0), vec3(0, 1, 0), vec3(1, 0, 0), M_PI * 2 * 10, 1);
//LineSegment loop(vec3(-.1, .1, 0), vec3(.1, .1, 0));
Toroid loop(vec3(0, 0, 0), vec3(1, 0, 0), vec3(0, 0, 1), M_PI * 2, .1, 2*M_PI / 10);
ostream *dump_ofs = NULL;
vec3 dump(vec3 s, vec3 ds)
{
*dump_ofs << s << endl;
//cout << "Magn: " << ds.magnitude() << endl;
return 0;
}
void dump_path(ostream &out, const Curve *c)
{
dump_ofs = &out;
c->integrate(dump, 1e-2);
}
/* dump the field (gnuplot format) at z = 0 */
/* requires x0 < x1, y0 < y1 */
const scalar U0 = 4e-7 * M_PI;
const scalar I = 1;
const scalar delta = .2;
void dump_field(scalar x0, scalar y0, scalar z0, scalar x1, scalar y1, scalar z1)
{
for(scalar z = z0; z <= z1; z += delta)
for(scalar y = y0; y <= y1; y += delta)
for(scalar x = x0; x <= x1; x += delta)
{
point = vec3(x, y, z);
vec3 B = loop.integrate(dB, 1e-1) * U0 * I;
if(B.magnitude() > 1e-8)
{
B=B.normalize() / 10;
cout << point[0] << " " << point[1] << " " << point[2] << " ";
cout << B[0] << " " << B[1] << " " << B[2] << endl;
}
}
}
void dump_fieldline(ostream &out, vec3 x, scalar len)
{
point = x;
scalar delta = .1;
while(len > 0)
{
out << point << endl;
vec3 B = loop.integrate(dB, 1e-1) * U0 * I;
point += delta * B;
len -= delta;
}
}
int main()
{
//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;
dump_field(-2, -2, 0,
12, 2, 0);
//dump_field(0,0,0,0,0,0);
stringstream ss;
ss << "curve.fld";
ofstream ofs(ss.str());
dump_path(ofs, (Curve*)&loop);
ofs.close();
#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
}
|
