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
|
#include <iostream>
#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;
}
vec3 Arc::integrate(vec3 (*integrand)(vec3 s, vec3 ds), double dl)
{
//cout << "Integrating loop of length " << this->angle << " radians" << endl;
double r = this->radius.magnitude(), dtheta = dl / r;
//cout << "R = " << r << ", dtheta = " << dtheta << endl;
quat rot = quat::from_angleaxis(dtheta, normal), crot = rot.conjugate();
//cout << "rot = " << rot << ", crot = " << crot << endl;
double len = this->angle * r, l;
vec3 sum = 0;
quat radius = this->radius;
for(l = 0; l < len; l += dl)
{
vec3 ds = this->normal.cross(radius).normalize() * dl;
sum += integrand(this->center + radius, ds);
radius = rot * radius * crot;
}
if(l < len)
{
dl = len - l;
dtheta = dl / r;
rot = quat::from_angleaxis(dtheta, normal);
crot = rot.conjugate();
vec3 ds = this->normal.cross(radius).normalize() * dl;
sum += integrand(this->center + radius, ds);
}
return sum;
}
vec3 Spiral::integrate(vec3 (*integrand)(vec3 s, vec3 ds), double dl)
{
//cout << "Integrating loop of length " << this->angle << " radians" << endl;
double r = this->radius.magnitude(), dtheta = dl / r;
//cout << "R = " << r << ", dtheta = " << dtheta << endl;
quat rot = quat::from_angleaxis(dtheta, normal), crot = rot.conjugate();
//cout << "rot = " << rot << ", crot = " << crot << endl;
double len = this->angle * r, l;
vec3 sum = 0;
quat radius = this->radius;
/* how far along the axis we've moved */
vec3 axis = 0;
/* we add this axial component each iteration */
vec3 dp = this->normal * dtheta * pitch / (2 * M_PI);
for(l = 0; l < len; l += dl, axis += dp)
{
vec3 ds = this->normal.cross(radius).normalize() * dl + dp;
sum += integrand(this->origin + axis + radius, ds);
/* rotate by dtheta about normal */
radius = rot * radius * crot;
}
if(l < len)
{
dl = len - l;
dtheta = dl / r;
rot = quat::from_angleaxis(dtheta, normal);
crot = rot.conjugate();
dp = this->normal * dtheta * pitch / (2 * M_PI);
vec3 ds = this->normal.cross(radius).normalize() * dl + dp;
sum += integrand(this->origin + axis + radius, ds);
}
return sum;
}
|
