did I fix the meshing now?
[antikythera] / src / machine.cc
1 #include <stdlib.h>
2 #include <string.h>
3 #include <math.h>
4 #include <float.h>
5 #include <assert.h>
6 #include "machine.h"
7
8 static float delta_angle(float a, float b);
9
10 Machine::Machine()
11 {
12         meshing = 0;
13         meshing_valid = false;
14         visited = 0;
15 }
16
17 Machine::~Machine()
18 {
19         int ngears = (int)gears.size();
20         for(int i=0; i<ngears; i++) {
21                 delete gears[i];
22         }
23
24         if(meshing) {
25                 delete [] meshing[0];
26                 delete [] meshing;
27         }
28         delete [] visited;
29 }
30
31 void Machine::add_gear(Gear *g)
32 {
33         gears.push_back(g);
34         meshing_valid = false;
35 }
36
37 void Machine::add_motor(int gearidx, float speed_hz)
38 {
39         Motor m;
40         m.drive = gearidx;
41         m.speed = speed_hz;
42         motors.push_back(m);
43 }
44
45 void Machine::invalidate_meshing()
46 {
47         meshing_valid = false;
48 }
49
50 void Machine::calc_meshing()
51 {
52         int ngears = (int)gears.size();
53
54         if(!meshing) {
55                 meshing = new bool*[ngears];
56                 meshing[0] = new bool[ngears * ngears];
57
58                 for(int i=1; i<ngears; i++) {
59                         meshing[i] = meshing[i - 1] + ngears;
60                 }
61         }
62
63         if(!visited) {
64                 visited = new bool[ngears];
65         }
66
67         // we're going to need the planar position of each gear on its plane, so let's cache it
68         Vec3 *ppos = (Vec3*)alloca(ngears * sizeof *ppos);
69         for(int i=0; i<ngears; i++) {
70                 ppos[i] = gears[i]->get_planar_position();
71         }
72
73         for(int i=0; i<ngears; i++) {
74                 for(int j=i; j<ngears; j++) {
75                         meshing[i][j] = meshing[j][i] = false;
76
77                         if(i == j) continue;
78
79                         if(1.0 - fabs(dot(gears[i]->axis, gears[j]->axis)) < 1e-5) {
80                                 // co-planar, just check Z range after inverse-transforming to the XY plane
81                                 if(fabs(ppos[i].z - ppos[j].z) > (gears[i]->thickness + gears[j]->thickness) / 2.0) {
82                                         continue;
83                                 }
84                                 // Z interval match, check distance
85                                 float dsq = length_sq(ppos[i].xy() - ppos[j].xy());
86
87                                 float outer_rad_sum = gears[i]->radius + gears[j]->radius;
88                                 float inner_rad_sum = outer_rad_sum - gears[i]->teeth_length - gears[j]->teeth_length;
89
90                                 if(dsq <= outer_rad_sum * outer_rad_sum && dsq >= inner_rad_sum * inner_rad_sum) {
91                                         //printf("connecting co-planar gears %d - %d\n", i, j);
92                                         meshing[i][j] = meshing[j][i] = true;
93                                 }
94
95                         } else {
96                                 /* TODO: not co-planar
97                                  * - calc line of intersection between the two planes
98                                  * - find distance of each gear to that line
99                                  * - profit...
100                                  */
101                         }
102                 }
103         }
104
105         // fix the initial angles so that teeth mesh as best as possible
106         // should work in one pass as long as the gear train is not impossible
107         for(int i=0; i<ngears; i++) {
108                 float rnd = gears[i]->angle + gears[i]->get_angular_pitch() / 2.0;
109                 float snap = rnd - fmod(rnd, gears[i]->get_angular_pitch());
110                 gears[i]->set_angle(snap);
111                 gears[i]->set_angular_offset(0);
112         }
113
114         for(int i=0; i<ngears; i++) {
115                 for(int j=i; j<ngears; j++) {
116                         if(meshing[i][j]) {
117                                 assert(i != j);
118
119                                 Vec2 dir = normalize(ppos[j].xy() - ppos[i].xy());
120                                 float rel_angle = atan2(dir.y, dir.x);
121
122                                 float frac_i = fmod((gears[i]->init_angle + rel_angle) / gears[i]->get_angular_pitch() + 100.0, 1.0);
123                                 float frac_j = fmod((gears[j]->init_angle - rel_angle) / gears[j]->get_angular_pitch() + 100.0, 1.0);
124                                 assert(frac_i >= 0.0 && frac_j >= 0.0);
125                                 float delta = frac_j - frac_i;
126
127                                 float correction = 0.5 - delta;
128                                 float prev_offs = gears[j]->get_angular_offset();
129                                 gears[j]->set_angular_offset(prev_offs + correction * gears[j]->get_angular_pitch());
130                         }
131                 }
132         }
133 }
134
135 void Machine::update_gear(int idx, float angle)
136 {
137         if(visited[idx]) {
138                 if(delta_angle(angle, gears[idx]->angle) > 0.25 / gears[idx]->nteeth) {
139                         fprintf(stderr, "warning: trying to transmit different values to gear %s (%d)\n",
140                                         gears[idx]->name.c_str(), idx);
141                         gears[idx]->angle = 0;
142                 }
143                 return;
144         }
145
146         gears[idx]->set_angle(angle);
147         visited[idx] = true;
148
149         int ngears = (int)gears.size();
150         for(int i=0; i<ngears; i++) {
151                 if(!meshing[idx][i]) continue;
152                 assert(idx != i);
153
154                 float ratio = -(float)gears[idx]->nteeth / (float)gears[i]->nteeth;
155                 update_gear(i, angle * ratio);
156         }
157 }
158
159 void Machine::update(float dt)
160 {
161         int ngears = (int)gears.size();
162
163         if(!meshing_valid) {
164                 calc_meshing();
165                 meshing_valid = true;
166         }
167
168         memset(visited, 0, ngears * sizeof *visited);
169         for(size_t i=0; i<motors.size(); i++) {
170                 int gidx = motors[i].drive;
171                 if(gidx < 0) continue;
172
173                 update_gear(gidx, gears[gidx]->angle + dt * motors[i].speed);
174         }
175 }
176
177 void Machine::draw() const
178 {
179         for(size_t i=0; i<gears.size(); i++) {
180                 gears[i]->draw();
181         }
182 }
183
184 Gear *Machine::intersect_gear(const Ray &ray, HitPoint *hitp) const
185 {
186         Gear *res = 0;
187         HitPoint nearest;
188         nearest.dist = FLT_MAX;
189
190         for(size_t i=0; i<gears.size(); i++) {
191                 Vec3 pos = gears[i]->get_global_position();
192                 float rad = gears[i]->radius;
193
194                 Plane plane = Plane(pos, gears[i]->axis);
195
196                 HitPoint hit;
197                 if(plane.intersect(ray, &hit) && hit.dist < nearest.dist &&
198                                 length_sq(hit.pos - pos) <= rad * rad) {
199                         nearest = hit;
200                         res = gears[i];
201                 }
202         }
203
204         if(hitp) *hitp = nearest;
205         return res;
206 }
207
208 static float delta_angle(float a, float b)
209 {
210         float api = fmod(a + M_PI, 2.0 * M_PI);
211         float bpi = fmod(b + M_PI, 2.0 * M_PI);
212         return std::min(fabs(a - b), fabs(api - bpi));
213 }