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
|
#include <cmath>
#include <iostream>
#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;
double 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(0, 0, 0), vec3(10, 0, 0));
/* dump the field (gnuplot format) at z = 0 */
/* requires x0 < x1, y0 < y1 */
void dump_field(double x0, double y0, double z0, double x1, double y1, double z1)
{
const double delta = .2;
for(double z = z0; z <= z1; z += delta)
for(double y = y0; y <= y1; y += delta)
for(double x = x0; x <= x1; x += delta)
{
point = vec3(x, y, z);
const double U0 = 4e-7 * M_PI;
const double I = 1;
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;
}
}
}
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;
//double 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(-1, -1.5, -1.5,
11, 1.5, 1.5);
}
|
