altos: Compute initial rotation from vertical
[fw/altos] / src / core / ao_sample.c
1 /*
2  * Copyright © 2011 Keith Packard <keithp@keithp.com>
3  *
4  * This program is free software; you can redistribute it and/or modify
5  * it under the terms of the GNU General Public License as published by
6  * the Free Software Foundation; version 2 of the License.
7  *
8  * This program is distributed in the hope that it will be useful, but
9  * WITHOUT ANY WARRANTY; without even the implied warranty of
10  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
11  * General Public License for more details.
12  *
13  * You should have received a copy of the GNU General Public License along
14  * with this program; if not, write to the Free Software Foundation, Inc.,
15  * 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA.
16  */
17
18 #ifndef AO_FLIGHT_TEST
19 #include "ao.h"
20 #include <ao_data.h>
21 #endif
22
23 #if HAS_GYRO
24 #include <ao_quaternion.h>
25 #endif
26
27 /*
28  * Current sensor values
29  */
30
31 #ifndef PRES_TYPE
32 #define PRES_TYPE int32_t
33 #define ALT_TYPE int32_t
34 #define ACCEL_TYPE int16_t
35 #endif
36
37 __pdata uint16_t        ao_sample_tick;         /* time of last data */
38 __pdata pres_t          ao_sample_pres;
39 __pdata alt_t           ao_sample_alt;
40 __pdata alt_t           ao_sample_height;
41 #if HAS_ACCEL
42 __pdata accel_t         ao_sample_accel;
43 #endif
44 #if HAS_GYRO
45 __pdata accel_t         ao_sample_accel_along;
46 __pdata accel_t         ao_sample_accel_across;
47 __pdata accel_t         ao_sample_accel_through;
48 __pdata gyro_t          ao_sample_roll;
49 __pdata gyro_t          ao_sample_pitch;
50 __pdata gyro_t          ao_sample_yaw;
51 __pdata angle_t         ao_sample_orient;
52 #endif
53
54 __data uint8_t          ao_sample_data;
55
56 /*
57  * Sensor calibration values
58  */
59
60 __pdata pres_t          ao_ground_pres;         /* startup pressure */
61 __pdata alt_t           ao_ground_height;       /* MSL of ao_ground_pres */
62
63 #if HAS_ACCEL
64 __pdata accel_t         ao_ground_accel;        /* startup acceleration */
65 __pdata accel_t         ao_accel_2g;            /* factory accel calibration */
66 __pdata int32_t         ao_accel_scale;         /* sensor to m/s² conversion */
67 #endif
68
69 #if HAS_GYRO
70 __pdata accel_t         ao_ground_accel_along;
71 __pdata accel_t         ao_ground_accel_across;
72 __pdata accel_t         ao_ground_accel_through;
73 __pdata int32_t         ao_ground_pitch;
74 __pdata int32_t         ao_ground_yaw;
75 __pdata int32_t         ao_ground_roll;
76 #endif
77
78 static __pdata uint8_t  ao_preflight;           /* in preflight mode */
79
80 static __pdata uint16_t nsamples;
81 __pdata int32_t ao_sample_pres_sum;
82 #if HAS_ACCEL
83 __pdata int32_t ao_sample_accel_sum;
84 #endif
85 #if HAS_GYRO
86 __pdata int32_t ao_sample_accel_along_sum;
87 __pdata int32_t ao_sample_accel_across_sum;
88 __pdata int32_t ao_sample_accel_through_sum;
89 __pdata int32_t ao_sample_pitch_sum;
90 __pdata int32_t ao_sample_yaw_sum;
91 __pdata int32_t ao_sample_roll_sum;
92 static struct ao_quaternion ao_rotation;
93 #endif
94
95 static void
96 ao_sample_preflight_add(void)
97 {
98 #if HAS_ACCEL
99         ao_sample_accel_sum += ao_sample_accel;
100 #endif
101         ao_sample_pres_sum += ao_sample_pres;
102 #if HAS_GYRO
103         ao_sample_accel_along_sum += ao_sample_accel_along;
104         ao_sample_accel_across_sum += ao_sample_accel_across;
105         ao_sample_accel_through_sum += ao_sample_accel_through;
106         ao_sample_pitch_sum += ao_sample_pitch;
107         ao_sample_yaw_sum += ao_sample_yaw;
108         ao_sample_roll_sum += ao_sample_roll;
109 #endif
110         ++nsamples;
111 }
112
113 static void
114 ao_sample_preflight_set(void)
115 {
116 #if HAS_ACCEL
117         ao_ground_accel = ao_sample_accel_sum >> 9;
118         ao_sample_accel_sum = 0;
119 #endif
120         ao_ground_pres = ao_sample_pres_sum >> 9;
121         ao_ground_height = pres_to_altitude(ao_ground_pres);
122         ao_sample_pres_sum = 0;
123 #if HAS_GYRO
124         ao_ground_accel_along = ao_sample_accel_along_sum >> 9;
125         ao_ground_accel_across = ao_sample_accel_across_sum >> 9;
126         ao_ground_accel_through = ao_sample_accel_through_sum >> 9;
127         ao_ground_pitch = ao_sample_pitch_sum;
128         ao_ground_yaw = ao_sample_yaw_sum;
129         ao_ground_roll = ao_sample_roll_sum;
130         ao_sample_accel_along_sum = 0;
131         ao_sample_accel_across_sum = 0;
132         ao_sample_accel_through_sum = 0;
133         ao_sample_pitch_sum = 0;
134         ao_sample_yaw_sum = 0;
135         ao_sample_roll_sum = 0;
136         ao_sample_orient = 0;
137
138         struct ao_quaternion    orient;
139
140         /* Take the pad IMU acceleration values and compute our current direction
141          */
142
143         ao_quaternion_init_vector(&orient,
144                                   (ao_ground_accel_across - ao_config.accel_zero_across),
145                                   (ao_ground_accel_through - ao_config.accel_zero_through),
146                                   (ao_ground_accel_along - ao_config.accel_zero_along));
147
148         ao_quaternion_normalize(&orient,
149                                 &orient);
150
151         /* Here's up */
152
153         struct ao_quaternion    up = { .r = 0, .x = 0, .y = 0, .z = 1 };
154
155         if (ao_config.pad_orientation != AO_PAD_ORIENTATION_ANTENNA_UP)
156                 up.z = -1;
157
158         /* Compute rotation to get from up to our current orientation, set
159          * that as the current rotation vector
160          */
161         ao_quaternion_vectors_to_rotation(&ao_rotation, &up, &orient);
162 #endif  
163         nsamples = 0;
164 }
165
166 #if HAS_GYRO
167 static void
168 ao_sample_rotate(void)
169 {
170 #ifdef AO_FLIGHT_TEST
171         float   dt = (ao_sample_tick - ao_sample_prev_tick) / 100.0;
172 #else
173         static const float dt = 1/100.0;
174 #endif
175         float   x = ao_mpu6000_gyro((float) ((ao_sample_pitch << 9) - ao_ground_pitch) / 512.0f) * dt;
176         float   y = ao_mpu6000_gyro((float) ((ao_sample_yaw << 9) - ao_ground_yaw) / 512.0f) * dt;
177         float   z = ao_mpu6000_gyro((float) ((ao_sample_roll << 9) - ao_ground_roll) / 512.0f) * dt;
178         struct ao_quaternion    rot;
179         struct ao_quaternion    point;
180
181         /* The amount of rotation is just the length of the vector. Now,
182          * here's the trick -- assume that the rotation amount is small. In this case,
183          * sin(x) ≃ x, so we can just make this the sin.
184          */
185
186         n_2 = x*x + y*y + z*z;
187         n = sqrtf(n_2);
188         s = n / 2;
189         if (s > 1)
190                 s = 1;
191         c = sqrtf(1 - s*s);
192
193         /* Make unit vector */
194         if (n > 0) {
195                 x /= n;
196                 y /= n;
197                 z /= n;
198         }
199
200         /* Now compute the unified rotation quaternion */
201
202         ao_quaternion_init_rotation(&rot,
203                                     x, y, z,
204                                     s, c);
205
206         /* Integrate with the previous rotation amount */
207         ao_quaternion_multiply(&ao_rotation, &ao_rotation, &rot);
208
209         /* And normalize to make sure it remains a unit vector */
210         ao_quaternion_normalize(&ao_rotation, &ao_rotation);
211
212         /* Compute pitch angle from vertical by taking the pad
213          * orientation vector and rotating it by the current total
214          * rotation value. That will be a unit vector pointing along
215          * the airframe axis. The Z value will be the cosine of the
216          * change in the angle from vertical since boost.
217          *
218          * rot = ao_rotation * vertical * ao_rotation°
219          * rot = ao_rotation * (0,0,0,1) * ao_rotation°
220          *     = ((a.z, a.y, -a.x, a.r) * (a.r, -a.x, -a.y, -a.z)) .z
221          *
222          *     = (-a.z * -a.z) + (a.y * -a.y) - (-a.x * -a.x) + (a.r * a.r)
223          *     = a.z² - a.y² - a.x² + a.r²
224          *
225          * rot = ao_rotation * (0, 0, 0, -1) * ao_rotation°
226          *     = ((-a.z, -a.y, a.x, -a.r) * (a.r, -a.x, -a.y, -a.z)) .z
227          *
228          *     = (a.z * -a.z) + (-a.y * -a.y) - (a.x * -a.x) + (-a.r * a.r)
229          *     = -a.z² + a.y² + a.x² - a.r²
230          */
231
232         float rotz;
233         rotz = ao_rotation.z * ao_rotation.z - ao_rotation.y * ao_rotation.y - ao_rotation.x * ao_rotation.x + ao_rotation.r * ao_rotation.r;
234
235         ao_sample_orient = acosf(rotz) * (float) (180.0/M_PI);
236 }
237 #endif
238
239 static void
240 ao_sample_preflight(void)
241 {
242         /* startup state:
243          *
244          * Collect 512 samples of acceleration and pressure
245          * data and average them to find the resting values
246          */
247         if (nsamples < 512) {
248                 ao_sample_preflight_add();
249         } else {
250 #if HAS_ACCEL
251                 ao_accel_2g = ao_config.accel_minus_g - ao_config.accel_plus_g;
252                 ao_accel_scale = to_fix32(GRAVITY * 2 * 16) / ao_accel_2g;
253 #endif
254                 ao_sample_preflight_set();
255                 ao_preflight = FALSE;
256         }
257 }
258
259 /*
260  * While in pad mode, constantly update the ground state by
261  * re-averaging the data.  This tracks changes in orientation, which
262  * might be caused by adjustments to the rocket on the pad and
263  * pressure, which might be caused by changes in the weather.
264  */
265
266 static void
267 ao_sample_preflight_update(void)
268 {
269         if (nsamples < 512)
270                 ao_sample_preflight_add();
271         else if (nsamples < 1024)
272                 ++nsamples;
273         else
274                 ao_sample_preflight_set();
275 }
276
277 #if 0
278 #if HAS_GYRO
279 static int32_t  p_filt;
280 static int32_t  y_filt;
281
282 static gyro_t inline ao_gyro(void) {
283         gyro_t  p = ao_sample_pitch - ao_ground_pitch;
284         gyro_t  y = ao_sample_yaw - ao_ground_yaw;
285
286         p_filt = p_filt - (p_filt >> 6) + p;
287         y_filt = y_filt - (y_filt >> 6) + y;
288
289         p = p_filt >> 6;
290         y = y_filt >> 6;
291         return ao_sqrt(p*p + y*y);
292 }
293 #endif
294 #endif
295
296 uint8_t
297 ao_sample(void)
298 {
299         ao_wakeup(DATA_TO_XDATA(&ao_sample_data));
300         ao_sleep((void *) DATA_TO_XDATA(&ao_data_head));
301         while (ao_sample_data != ao_data_head) {
302                 __xdata struct ao_data *ao_data;
303
304                 /* Capture a sample */
305                 ao_data = (struct ao_data *) &ao_data_ring[ao_sample_data];
306                 ao_sample_tick = ao_data->tick;
307
308 #if HAS_BARO
309                 ao_data_pres_cook(ao_data);
310                 ao_sample_pres = ao_data_pres(ao_data);
311                 ao_sample_alt = pres_to_altitude(ao_sample_pres);
312                 ao_sample_height = ao_sample_alt - ao_ground_height;
313 #endif
314
315 #if HAS_ACCEL
316                 ao_sample_accel = ao_data_accel_cook(ao_data);
317                 if (ao_config.pad_orientation != AO_PAD_ORIENTATION_ANTENNA_UP)
318                         ao_sample_accel = ao_data_accel_invert(ao_sample_accel);
319                 ao_data_set_accel(ao_data, ao_sample_accel);
320 #endif
321 #if HAS_GYRO
322                 ao_sample_accel_along = ao_data_along(ao_data);
323                 ao_sample_accel_across = ao_data_across(ao_data);
324                 ao_sample_accel_through = ao_data_through(ao_data);
325                 ao_sample_pitch = ao_data_pitch(ao_data);
326                 ao_sample_yaw = ao_data_yaw(ao_data);
327                 ao_sample_roll = ao_data_roll(ao_data);
328 #endif
329
330                 if (ao_preflight)
331                         ao_sample_preflight();
332                 else {
333                         if (ao_flight_state < ao_flight_boost)
334                                 ao_sample_preflight_update();
335                         ao_kalman();
336 #if HAS_GYRO
337                         ao_sample_rotate();
338 #endif
339                 }
340 #ifdef AO_FLIGHT_TEST
341                 ao_sample_prev_tick = ao_sample_tick;
342 #endif
343                 ao_sample_data = ao_data_ring_next(ao_sample_data);
344         }
345         return !ao_preflight;
346 }
347
348 void
349 ao_sample_init(void)
350 {
351         ao_config_get();
352         nsamples = 0;
353         ao_sample_pres_sum = 0;
354         ao_sample_pres = 0;
355 #if HAS_ACCEL
356         ao_sample_accel_sum = 0;
357         ao_sample_accel = 0;
358 #endif
359 #if HAS_GYRO
360         ao_sample_accel_along_sum = 0;
361         ao_sample_accel_across_sum = 0;
362         ao_sample_accel_through_sum = 0;
363         ao_sample_accel_along = 0;
364         ao_sample_accel_across = 0;
365         ao_sample_accel_through = 0;
366         ao_sample_pitch_sum = 0;
367         ao_sample_yaw_sum = 0;
368         ao_sample_roll_sum = 0;
369         ao_sample_pitch = 0;
370         ao_sample_yaw = 0;
371         ao_sample_roll = 0;
372         ao_sample_orient = 0;
373 #endif
374         ao_sample_data = ao_data_head;
375         ao_preflight = TRUE;
376 }