diff options
Diffstat (limited to 'ao-tools/lib')
-rw-r--r-- | ao-tools/lib/Makefile.am | 24 | ||||
-rw-r--r-- | ao-tools/lib/cc-analyse.c | 267 | ||||
-rw-r--r-- | ao-tools/lib/cc-convert.c | 284 | ||||
-rw-r--r-- | ao-tools/lib/cc-dsp.c | 145 | ||||
-rw-r--r-- | ao-tools/lib/cc-integrate.c | 81 | ||||
-rw-r--r-- | ao-tools/lib/cc-log.c | 121 | ||||
-rw-r--r-- | ao-tools/lib/cc-logfile.c | 336 | ||||
-rw-r--r-- | ao-tools/lib/cc-period.c | 56 | ||||
-rw-r--r-- | ao-tools/lib/cc-process.c | 166 | ||||
-rw-r--r-- | ao-tools/lib/cc-telem.c | 212 | ||||
-rw-r--r-- | ao-tools/lib/cc-telemetry.c | 62 | ||||
-rw-r--r-- | ao-tools/lib/cc-telemetry.h | 243 | ||||
-rw-r--r-- | ao-tools/lib/cc-usb.c | 195 | ||||
-rw-r--r-- | ao-tools/lib/cc-usb.h | 12 | ||||
-rw-r--r-- | ao-tools/lib/cc-usbdev.c | 317 | ||||
-rw-r--r-- | ao-tools/lib/cc-util.c | 80 | ||||
-rw-r--r-- | ao-tools/lib/cc.h | 373 | ||||
-rw-r--r-- | ao-tools/lib/cephes.h | 122 | ||||
-rw-r--r-- | ao-tools/lib/chbevl.c | 81 | ||||
-rw-r--r-- | ao-tools/lib/cmath.h | 179 | ||||
-rw-r--r-- | ao-tools/lib/i0.c | 414 | ||||
-rw-r--r-- | ao-tools/lib/mconf.h | 211 |
22 files changed, 3921 insertions, 60 deletions
diff --git a/ao-tools/lib/Makefile.am b/ao-tools/lib/Makefile.am index 9584e216..1f8f2e42 100644 --- a/ao-tools/lib/Makefile.am +++ b/ao-tools/lib/Makefile.am @@ -1,6 +1,9 @@ noinst_LIBRARIES = libao-tools.a -AM_CFLAGS=$(WARN_CFLAGS) $(LIBUSB_CFLAGS) +AM_CFLAGS=$(WARN_CFLAGS) $(LIBUSB_CFLAGS) $(GNOME_CFLAGS) + +libao_tools_a_uneeded = \ + cc-log.c libao_tools_a_SOURCES = \ ccdbg-command.c \ @@ -14,9 +17,26 @@ libao_tools_a_SOURCES = \ ccdbg-memory.c \ ccdbg-rom.c \ ccdbg-state.c \ + cc-analyse.c \ + cc-convert.c \ + cc-dsp.c \ + cc-integrate.c \ + cc-period.c \ + cc-process.c \ cc-usb.c \ cc-usb.h \ + cc.h \ + cc-usbdev.c \ + cc-util.c \ cc-bitbang.c \ cc-bitbang.h \ + cc-logfile.c \ + cc-telem.c \ + cc-telemetry.c \ + cc-telemetry.h \ cp-usb-async.c \ - cp-usb-async.h + cp-usb-async.h \ + i0.c \ + chbevl.c \ + mconf.h \ + cephes.h diff --git a/ao-tools/lib/cc-analyse.c b/ao-tools/lib/cc-analyse.c new file mode 100644 index 00000000..27c416a6 --- /dev/null +++ b/ao-tools/lib/cc-analyse.c @@ -0,0 +1,267 @@ +/* + * Copyright © 2009 Keith Packard <keithp@keithp.com> + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; version 2 of the License. + * + * This program is distributed in the hope that it will be useful, but + * WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * General Public License for more details. + * + * You should have received a copy of the GNU General Public License along + * with this program; if not, write to the Free Software Foundation, Inc., + * 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA. + */ + +#include "cc.h" +#include <math.h> + +void +cc_timedata_limits(struct cc_timedata *d, double min_time, double max_time, int *start, int *stop) +{ + int i; + + *start = -1; + for (i = 0; i < d->num; i++) { + if (*start < 0 && min_time <= d->data[i].time) + *start = i; + if (d->data[i].time <= max_time) + *stop = i; + } +} + +int +cc_timedata_min(struct cc_timedata *d, double min_time, double max_time) +{ + int i; + int set = 0; + int min_i = -1; + double min; + + if (d->num == 0) + return -1; + for (i = 0; i < d->num; i++) + if (min_time <= d->data[i].time && d->data[i].time <= max_time) + if (!set || d->data[i].value < min) { + min_i = i; + min = d->data[i].value; + set = 1; + } + return min_i; +} + +int +cc_timedata_min_mag(struct cc_timedata *d, double min_time, double max_time) +{ + int i; + int set = 0; + int min_i = -1; + double min; + + if (d->num == 0) + return -1; + for (i = 0; i < d->num; i++) + if (min_time <= d->data[i].time && d->data[i].time <= max_time) + if (!set || fabs(d->data[i].value) < min) { + min_i = i; + min = fabs(d->data[i].value); + set = 1; + } + return min_i; +} + +int +cc_timedata_max(struct cc_timedata *d, double min_time, double max_time) +{ + int i; + double max; + int max_i = -1; + int set = 0; + + if (d->num == 0) + return -1; + for (i = 0; i < d->num; i++) + if (min_time <= d->data[i].time && d->data[i].time <= max_time) + if (!set || d->data[i].value > max) { + max_i = i; + max = d->data[i].value; + set = 1; + } + return max_i; +} + +int +cc_timedata_max_mag(struct cc_timedata *d, double min_time, double max_time) +{ + int i; + double max; + int max_i = -1; + int set = 0; + + if (d->num == 0) + return -1; + for (i = 0; i < d->num; i++) + if (min_time <= d->data[i].time && d->data[i].time <= max_time) + if (!set || fabs(d->data[i].value) > max) { + max_i = i; + max = fabs(d->data[i].value); + set = 1; + } + return max_i; +} + +double +cc_timedata_average(struct cc_timedata *td, double start_time, double stop_time) +{ + int i; + double prev_time; + double next_time; + double interval; + double sum = 0.0; + double period = 0.0; + + prev_time = start_time; + for (i = 0; i < td->num; i++) { + if (start_time <= td->data[i].time && td->data[i].time <= stop_time) { + if (i < td->num - 1 && td->data[i+1].time < stop_time) + next_time = (td->data[i].time + td->data[i+1].time) / 2.0; + else + next_time = stop_time; + interval = next_time - prev_time; + sum += td->data[i].value * interval; + period += interval; + prev_time = next_time; + } + } + return sum / period; +} + +int +cc_perioddata_limits(struct cc_perioddata *d, double min_time, double max_time, int *start, int *stop) +{ + double start_d, stop_d; + + if (d->num == 0) + return 0; + start_d = ceil((min_time - d->start) / d->step); + if (start_d < 0) + start_d = 0; + stop_d = floor((max_time - d->start) / d->step); + if (stop_d >= d->num) + stop_d = d->num - 1; + if (stop_d < start_d) + return 0; + *start = (int) start_d; + *stop = (int) stop_d; + return 1; +} + +int +cc_perioddata_min(struct cc_perioddata *d, double min_time, double max_time) +{ + int i; + double min; + int min_i; + int start, stop; + + if (!cc_perioddata_limits(d, min_time, max_time, &start, &stop)) + return -1; + min = d->data[start]; + min_i = start; + for (i = start + 1; i <= stop; i++) + if (d->data[i] < min) { + min = d->data[i]; + min_i = i; + } + return min_i; +} + +int +cc_perioddata_min_mag(struct cc_perioddata *d, double min_time, double max_time) +{ + int start, stop; + int i; + double min; + int min_i; + + if (!cc_perioddata_limits(d, min_time, max_time, &start, &stop)) + return -1; + min = d->data[start]; + min_i = start; + for (i = start + 1; i <= stop; i++) + if (fabs(d->data[i]) < min) { + min = fabs(d->data[i]); + min_i = i; + } + return min_i; +} + +int +cc_perioddata_max(struct cc_perioddata *d, double min_time, double max_time) +{ + int start, stop; + int i; + double max; + int max_i; + + if (!cc_perioddata_limits(d, min_time, max_time, &start, &stop)) + return -1; + max = d->data[start]; + max_i = start; + for (i = start + 1; i <= stop; i++) + if (d->data[i] > max) { + max = d->data[i]; + max_i = i; + } + return max_i; +} + +int +cc_perioddata_max_mag(struct cc_perioddata *d, double min_time, double max_time) +{ + int start, stop; + int i; + double max; + int max_i; + + if (!cc_perioddata_limits(d, min_time, max_time, &start, &stop)) + return -1; + max = d->data[start]; + max_i = start; + for (i = start + 1; i <= stop; i++) + if (fabs(d->data[i]) > max) { + max = fabs(d->data[i]); + max_i = i; + } + return max_i; +} + +double +cc_perioddata_average(struct cc_perioddata *d, double min_time, double max_time) +{ + int start, stop; + int i; + double sum = 0.0; + + if (!cc_perioddata_limits(d, min_time, max_time, &start, &stop)) + return 0.0; + for (i = start; i <= stop; i++) + sum += d->data[i]; + return sum / (stop - start + 1); +} + +double +cc_perioddata_average_mag(struct cc_perioddata *d, double min_time, double max_time) +{ + int start, stop; + int i; + double sum = 0.0; + + if (!cc_perioddata_limits(d, min_time, max_time, &start, &stop)) + return 0.0; + for (i = start; i <= stop; i++) + sum += fabs(d->data[i]); + return sum / (stop - start + 1); +} diff --git a/ao-tools/lib/cc-convert.c b/ao-tools/lib/cc-convert.c new file mode 100644 index 00000000..8d6876a0 --- /dev/null +++ b/ao-tools/lib/cc-convert.c @@ -0,0 +1,284 @@ +/* + * Copyright © 2009 Keith Packard <keithp@keithp.com> + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; version 2 of the License. + * + * This program is distributed in the hope that it will be useful, but + * WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * General Public License for more details. + * + * You should have received a copy of the GNU General Public License along + * with this program; if not, write to the Free Software Foundation, Inc., + * 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA. + */ + +#include "cc.h" +#include <math.h> + +/* + * Pressure Sensor Model, version 1.1 + * + * written by Holly Grimes + * + * Uses the International Standard Atmosphere as described in + * "A Quick Derivation relating altitude to air pressure" (version 1.03) + * from the Portland State Aerospace Society, except that the atmosphere + * is divided into layers with each layer having a different lapse rate. + * + * Lapse rate data for each layer was obtained from Wikipedia on Sept. 1, 2007 + * at site <http://en.wikipedia.org/wiki/International_Standard_Atmosphere + * + * Height measurements use the local tangent plane. The postive z-direction is up. + * + * All measurements are given in SI units (Kelvin, Pascal, meter, meters/second^2). + * The lapse rate is given in Kelvin/meter, the gas constant for air is given + * in Joules/(kilogram-Kelvin). + */ + +#define GRAVITATIONAL_ACCELERATION -9.80665 +#define AIR_GAS_CONSTANT 287.053 +#define NUMBER_OF_LAYERS 7 +#define MAXIMUM_ALTITUDE 84852.0 +#define MINIMUM_PRESSURE 0.3734 +#define LAYER0_BASE_TEMPERATURE 288.15 +#define LAYER0_BASE_PRESSURE 101325 + +/* lapse rate and base altitude for each layer in the atmosphere */ +static const double lapse_rate[NUMBER_OF_LAYERS] = { + -0.0065, 0.0, 0.001, 0.0028, 0.0, -0.0028, -0.002 +}; + +static const int base_altitude[NUMBER_OF_LAYERS] = { + 0, 11000, 20000, 32000, 47000, 51000, 71000 +}; + +/* outputs atmospheric pressure associated with the given altitude. altitudes + are measured with respect to the mean sea level */ +double +cc_altitude_to_pressure(double altitude) +{ + + double base_temperature = LAYER0_BASE_TEMPERATURE; + double base_pressure = LAYER0_BASE_PRESSURE; + + double pressure; + double base; /* base for function to determine pressure */ + double exponent; /* exponent for function to determine pressure */ + int layer_number; /* identifies layer in the atmosphere */ + int delta_z; /* difference between two altitudes */ + + if (altitude > MAXIMUM_ALTITUDE) /* FIX ME: use sensor data to improve model */ + return 0; + + /* calculate the base temperature and pressure for the atmospheric layer + associated with the inputted altitude */ + for(layer_number = 0; layer_number < NUMBER_OF_LAYERS - 1 && altitude > base_altitude[layer_number + 1]; layer_number++) { + delta_z = base_altitude[layer_number + 1] - base_altitude[layer_number]; + if (lapse_rate[layer_number] == 0.0) { + exponent = GRAVITATIONAL_ACCELERATION * delta_z + / AIR_GAS_CONSTANT / base_temperature; + base_pressure *= exp(exponent); + } + else { + base = (lapse_rate[layer_number] * delta_z / base_temperature) + 1.0; + exponent = GRAVITATIONAL_ACCELERATION / + (AIR_GAS_CONSTANT * lapse_rate[layer_number]); + base_pressure *= pow(base, exponent); + } + base_temperature += delta_z * lapse_rate[layer_number]; + } + + /* calculate the pressure at the inputted altitude */ + delta_z = altitude - base_altitude[layer_number]; + if (lapse_rate[layer_number] == 0.0) { + exponent = GRAVITATIONAL_ACCELERATION * delta_z + / AIR_GAS_CONSTANT / base_temperature; + pressure = base_pressure * exp(exponent); + } + else { + base = (lapse_rate[layer_number] * delta_z / base_temperature) + 1.0; + exponent = GRAVITATIONAL_ACCELERATION / + (AIR_GAS_CONSTANT * lapse_rate[layer_number]); + pressure = base_pressure * pow(base, exponent); + } + + return pressure; +} + + +/* outputs the altitude associated with the given pressure. the altitude + returned is measured with respect to the mean sea level */ +double +cc_pressure_to_altitude(double pressure) +{ + + double next_base_temperature = LAYER0_BASE_TEMPERATURE; + double next_base_pressure = LAYER0_BASE_PRESSURE; + + double altitude; + double base_pressure; + double base_temperature; + double base; /* base for function to determine base pressure of next layer */ + double exponent; /* exponent for function to determine base pressure + of next layer */ + double coefficient; + int layer_number; /* identifies layer in the atmosphere */ + int delta_z; /* difference between two altitudes */ + + if (pressure < 0) /* illegal pressure */ + return -1; + if (pressure < MINIMUM_PRESSURE) /* FIX ME: use sensor data to improve model */ + return MAXIMUM_ALTITUDE; + + /* calculate the base temperature and pressure for the atmospheric layer + associated with the inputted pressure. */ + layer_number = -1; + do { + layer_number++; + base_pressure = next_base_pressure; + base_temperature = next_base_temperature; + delta_z = base_altitude[layer_number + 1] - base_altitude[layer_number]; + if (lapse_rate[layer_number] == 0.0) { + exponent = GRAVITATIONAL_ACCELERATION * delta_z + / AIR_GAS_CONSTANT / base_temperature; + next_base_pressure *= exp(exponent); + } + else { + base = (lapse_rate[layer_number] * delta_z / base_temperature) + 1.0; + exponent = GRAVITATIONAL_ACCELERATION / + (AIR_GAS_CONSTANT * lapse_rate[layer_number]); + next_base_pressure *= pow(base, exponent); + } + next_base_temperature += delta_z * lapse_rate[layer_number]; + } + while(layer_number < NUMBER_OF_LAYERS - 1 && pressure < next_base_pressure); + + /* calculate the altitude associated with the inputted pressure */ + if (lapse_rate[layer_number] == 0.0) { + coefficient = (AIR_GAS_CONSTANT / GRAVITATIONAL_ACCELERATION) + * base_temperature; + altitude = base_altitude[layer_number] + + coefficient * log(pressure / base_pressure); + } + else { + base = pressure / base_pressure; + exponent = AIR_GAS_CONSTANT * lapse_rate[layer_number] + / GRAVITATIONAL_ACCELERATION; + coefficient = base_temperature / lapse_rate[layer_number]; + altitude = base_altitude[layer_number] + + coefficient * (pow(base, exponent) - 1); + } + + return altitude; +} + +/* + * Values for our MP3H6115A pressure sensor + * + * From the data sheet: + * + * Pressure range: 15-115 kPa + * Voltage at 115kPa: 2.82 + * Output scale: 27mV/kPa + * + * + * 27 mV/kPa * 2047 / 3300 counts/mV = 16.75 counts/kPa + * 2.82V * 2047 / 3.3 counts/V = 1749 counts/115 kPa + */ + +static const double counts_per_kPa = 27 * 2047 / 3300; +static const double counts_at_101_3kPa = 1674.0; + +double +cc_barometer_to_pressure(double count) +{ + return ((count / 16.0) / 2047.0 + 0.095) / 0.009 * 1000.0; +} + +double +cc_barometer_to_altitude(double baro) +{ + double Pa = cc_barometer_to_pressure(baro); + return cc_pressure_to_altitude(Pa); +} + +static const double count_per_mss = 27.0; + +double +cc_accelerometer_to_acceleration(double accel, double ground_accel) +{ + return (ground_accel - accel) / count_per_mss; +} + +/* Value for the CC1111 built-in temperature sensor + * Output voltage at 0°C = 0.755V + * Coefficient = 0.00247V/°C + * Reference voltage = 1.25V + * + * temp = ((value / 32767) * 1.25 - 0.755) / 0.00247 + * = (value - 19791.268) / 32768 * 1.25 / 0.00247 + */ + +double +cc_thermometer_to_temperature(double thermo) +{ + return (thermo - 19791.268) / 32728.0 * 1.25 / 0.00247; +} + +double +cc_battery_to_voltage(double battery) +{ + return battery / 32767.0 * 5.0; +} + +double +cc_ignitor_to_voltage(double ignite) +{ + return ignite / 32767 * 15.0; +} + +static inline double sqr(double a) { return a * a; } + +void +cc_great_circle (double start_lat, double start_lon, + double end_lat, double end_lon, + double *dist, double *bearing) +{ + const double rad = M_PI / 180; + const double earth_radius = 6371.2 * 1000; /* in meters */ + double lat1 = rad * start_lat; + double lon1 = rad * -start_lon; + double lat2 = rad * end_lat; + double lon2 = rad * -end_lon; + +// double d_lat = lat2 - lat1; + double d_lon = lon2 - lon1; + + /* From http://en.wikipedia.org/wiki/Great-circle_distance */ + double vdn = sqrt(sqr(cos(lat2) * sin(d_lon)) + + sqr(cos(lat1) * sin(lat2) - + sin(lat1) * cos(lat2) * cos(d_lon))); + double vdd = sin(lat1) * sin(lat2) + cos(lat1) * cos(lat2) * cos(d_lon); + double d = atan2(vdn,vdd); + double course; + + if (cos(lat1) < 1e-20) { + if (lat1 > 0) + course = M_PI; + else + course = -M_PI; + } else { + if (d < 1e-10) + course = 0; + else + course = acos((sin(lat2)-sin(lat1)*cos(d)) / + (sin(d)*cos(lat1))); + if (sin(lon2-lon1) > 0) + course = 2 * M_PI-course; + } + *dist = d * earth_radius; + *bearing = course * 180/M_PI; +} diff --git a/ao-tools/lib/cc-dsp.c b/ao-tools/lib/cc-dsp.c new file mode 100644 index 00000000..518c1a68 --- /dev/null +++ b/ao-tools/lib/cc-dsp.c @@ -0,0 +1,145 @@ +/* + * Copyright © 2009 Keith Packard <keithp@keithp.com> + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; version 2 of the License. + * + * This program is distributed in the hope that it will be useful, but + * WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * General Public License for more details. + * + * You should have received a copy of the GNU General Public License along + * with this program; if not, write to the Free Software Foundation, Inc., + * 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA. + */ + +#include "cc.h" +#include "cephes.h" +#include <math.h> +#include <stdlib.h> + +static inline double sqr (double x) { return x * x; } + +/* + * Kaiser Window digital filter + */ + +#if 0 +/* not used in this program */ +static double highpass (double n, double m, double wc) +{ + double alpha = m/2; + double dist; + + dist = n - alpha; + if (dist == 0) + return (M_PI/2 - wc) / M_PI; + return -sin(dist * (M_PI/2-wc)) / (M_PI * dist); +} +#endif + +static double lowpass (double n, double m, double wc) +{ + double alpha = m/2; + double dist; + dist = n - alpha; + if (dist == 0) + return wc / M_PI; + return sin (wc * dist) / (M_PI * dist); +} + +static double kaiser (double n, double m, double beta) +{ + double alpha = m / 2; + return i0 (beta * sqrt (1 - sqr((n - alpha) / alpha))) / i0(beta); +} + +static double beta (double A) +{ + if (A > 50) + return 0.1102 * (A - 8.7); + else if (A >= 21) + return 0.5842 * pow((A - 21), 0.4) + 0.07886 * (A - 21); + else + return 0.0; +} + +static int M (double A, double delta_omega) +{ + if (A > 21) + return ceil ((A - 7.95) / (2.285 * delta_omega)); + else + return ceil(5.79 / delta_omega); +} + +struct filter_param { + double omega_pass; + double delta_omega; + double A; + double beta; + int M; +} filter_param_t; + +static struct filter_param +filter (double omega_pass, double omega_stop, double error) +{ + struct filter_param p; + p.omega_pass = omega_pass; + p.delta_omega = omega_stop - omega_pass; + p.A = -20 * log10 (error); + p.beta = beta (p.A); + p.M = M (p.A, p.delta_omega); + if ((p.M & 1) == 1) + p.M++; + return p; +} + +static double * +make_low_pass_filter(double omega_pass, double omega_stop, double error, int *length_p) +{ + struct filter_param p = filter(omega_pass, omega_stop, error); + int length; + int n; + double *lpf; + + length = p.M + 1; + lpf = calloc (length, sizeof(double)); + for (n = 0; n < length; n++) + lpf[n] = lowpass(n, p.M, omega_pass) * kaiser(n, p.M, p.beta); + *length_p = length; + return lpf; +} + +static double * +convolve(double *d, int d_len, double *e, int e_len) +{ + int w = (e_len - 1) / 2; + int n; + double *con = calloc (d_len, sizeof (double)); + + for (n = 0; n < d_len; n++) { + double v = 0; + int o; + for (o = -w; o <= w; o++) { + int p = n + o; + double sample = p < 0 ? d[0] : p >= d_len ? d[d_len-1] : d[p]; + v += sample * e[o + w]; + } + con[n] = v; + } + return con; +} + +double * +cc_low_pass(double *data, int data_len, double omega_pass, double omega_stop, double error) +{ + int fir_len; + double *fir = make_low_pass_filter(omega_pass, omega_stop, error, &fir_len); + double *result; + + result = convolve(data, data_len, fir, fir_len); + free(fir); + return result; +} diff --git a/ao-tools/lib/cc-integrate.c b/ao-tools/lib/cc-integrate.c new file mode 100644 index 00000000..ba50761b --- /dev/null +++ b/ao-tools/lib/cc-integrate.c @@ -0,0 +1,81 @@ +/* + * Copyright © 2009 Keith Packard <keithp@keithp.com> + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; version 2 of the License. + * + * This program is distributed in the hope that it will be useful, but + * WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * General Public License for more details. + * + * You should have received a copy of the GNU General Public License along + * with this program; if not, write to the Free Software Foundation, Inc., + * 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA. + */ + +#include "cc.h" +#include <stdlib.h> + +struct cc_timedata * +cc_timedata_convert(struct cc_timedata *d, double (*f)(double v, double a), double a) +{ + struct cc_timedata *r; + int n; + + r = calloc (1, sizeof (struct cc_timedata)); + r->num = d->num; + r->size = d->num; + r->data = calloc (r->size, sizeof (struct cc_timedataelt)); + r->time_offset = d->time_offset; + for (n = 0; n < d->num; n++) { + r->data[n].time = d->data[n].time; + r->data[n].value = f(d->data[n].value, a); + } + return r; +} + +struct cc_timedata * +cc_timedata_integrate(struct cc_timedata *d, double min_time, double max_time) +{ + struct cc_timedata *i; + int n, m; + int start, stop; + + cc_timedata_limits(d, min_time, max_time, &start, &stop); + i = calloc (1, sizeof (struct cc_timedata)); + i->num = stop - start + 1; + i->size = i->num; + i->data = calloc (i->size, sizeof (struct cc_timedataelt)); + i->time_offset = d->data[start].time; + for (n = 0; n < i->num; n++) { + m = n + start; + i->data[n].time = d->data[m].time; + if (n == 0) { + i->data[n].value = 0; + } else { + i->data[n].value = i->data[n-1].value + + (d->data[m].value + d->data[m-1].value) / 2 * + ((d->data[m].time - d->data[m-1].time) / 100.0); + } + } + return i; +} + +struct cc_perioddata * +cc_perioddata_differentiate(struct cc_perioddata *i) +{ + struct cc_perioddata *d; + int n; + + d = calloc (1, sizeof (struct cc_perioddata)); + d->num = i->num; + d->start = i->start; + d->step = i->step; + d->data = calloc (d->num, sizeof(double)); + for (n = 1; n < d->num; n++) + d->data[n] = (i->data[n] - i->data[n-1]) / (i->step / 100.0); + d->data[0] = d->data[1]; + return d; +} diff --git a/ao-tools/lib/cc-log.c b/ao-tools/lib/cc-log.c new file mode 100644 index 00000000..ed51f87e --- /dev/null +++ b/ao-tools/lib/cc-log.c @@ -0,0 +1,121 @@ +/* + * Copyright © 2009 Keith Packard <keithp@keithp.com> + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; version 2 of the License. + * + * This program is distributed in the hope that it will be useful, but + * WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * General Public License for more details. + * + * You should have received a copy of the GNU General Public License along + * with this program; if not, write to the Free Software Foundation, Inc., + * 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA. + */ + +#include <string.h> +#include <stdlib.h> +#include <stdio.h> +#include <unistd.h> +#include <gconf/gconf-client.h> +#include "cc.h" + +static char *cc_file_dir; + +#define ALTOS_DIR_PATH "/apps/aoview/log_dir" +#define DEFAULT_DIR "AltOS" + +static void +cc_file_save_conf(void) +{ + GConfClient *gconf_client; + + g_type_init(); + gconf_client = gconf_client_get_default(); + if (gconf_client) + { + gconf_client_set_string(gconf_client, + ALTOS_DIR_PATH, + cc_file_dir, + NULL); + g_object_unref(G_OBJECT(gconf_client)); + } +} + +static void +cc_file_load_conf(void) +{ + char *file_dir; + GConfClient *gconf_client; + + g_type_init(); + gconf_client = gconf_client_get_default(); + if (gconf_client) + { + file_dir = gconf_client_get_string(gconf_client, + ALTOS_DIR_PATH, + NULL); + g_object_unref(G_OBJECT(gconf_client)); + if (file_dir) + cc_file_dir = strdup(file_dir); + } +} + +void +cc_set_log_dir(char *dir) +{ + cc_file_dir = strdup(dir); + cc_file_save_conf(); +} + +char * +cc_get_log_dir(void) +{ + cc_file_load_conf(); + if (!cc_file_dir) { + cc_file_dir = cc_fullname(getenv("HOME"), DEFAULT_DIR); + cc_file_save_conf(); + } + return cc_file_dir; +} + +char * +cc_make_filename(int serial, int flight, char *ext) +{ + char base[50]; + char seq[20]; + struct tm tm; + time_t now; + char *full; + int r; + int sequence; + + now = time(NULL); + (void) localtime_r(&now, &tm); + cc_mkdir(cc_get_log_dir()); + sequence = 0; + for (;;) { + if (sequence) + snprintf(seq, sizeof(seq), "-seq-%03d", sequence); + else + seq[0] = '\0'; + + snprintf(base, sizeof (base), "%04d-%02d-%02d-serial-%03d-flight-%03d%s.%s", + tm.tm_year + 1900, + tm.tm_mon + 1, + tm.tm_mday, + serial, + flight, + seq, + ext); + full = cc_fullname(cc_get_log_dir(), base); + r = access(full, F_OK); + if (r < 0) + return full; + free(full); + sequence++; + } + +} diff --git a/ao-tools/lib/cc-logfile.c b/ao-tools/lib/cc-logfile.c new file mode 100644 index 00000000..842e5c7c --- /dev/null +++ b/ao-tools/lib/cc-logfile.c @@ -0,0 +1,336 @@ +/* + * Copyright © 2009 Keith Packard <keithp@keithp.com> + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; version 2 of the License. + * + * This program is distributed in the hope that it will be useful, but + * WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * General Public License for more details. + * + * You should have received a copy of the GNU General Public License along + * with this program; if not, write to the Free Software Foundation, Inc., + * 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA. + */ + +#include "cc.h" +#include <stdio.h> +#include <stdlib.h> +#include <string.h> + +static int +timedata_add(struct cc_timedata *data, double time, double value) +{ + struct cc_timedataelt *newdata; + int newsize; + if (data->size == data->num) { + if (data->size == 0) + newdata = malloc((newsize = 256) * sizeof (struct cc_timedataelt)); + else + newdata = realloc (data->data, (newsize = data->size * 2) + * sizeof (struct cc_timedataelt)); + if (!newdata) + return 0; + data->size = newsize; + data->data = newdata; + } + time += data->time_offset; + if (data->num && data->data[data->num-1].time > time) { + data->time_offset += 65536; + time += 65536; + } + data->data[data->num].time = time; + data->data[data->num].value = value; + data->num++; + return 1; +} + +static void +timedata_free(struct cc_timedata *data) +{ + if (data->data) + free(data->data); +} + +static int +gpsdata_add(struct cc_gpsdata *data, struct cc_gpselt *elt) +{ + struct cc_gpselt *newdata; + int newsize; + if (data->size == data->num) { + if (data->size == 0) + newdata = malloc((newsize = 256) * sizeof (struct cc_gpselt)); + else + newdata = realloc (data->data, (newsize = data->size * 2) + * sizeof (struct cc_gpselt)); + if (!newdata) + return 0; + data->size = newsize; + data->data = newdata; + } + elt->time += data->time_offset; + if (data->num && data->data[data->num-1].time > elt->time) { + data->time_offset += 65536; + elt->time += 65536; + } + data->data[data->num] = *elt; + data->num++; + return 1; +} + +static int +gpssat_add(struct cc_gpsdata *data, struct cc_gpssat *sat) +{ + int i, j; + int reuse = 0; + int newsizesats; + struct cc_gpssats *newsats; + + for (i = data->numsats; --i >= 0;) { + if (data->sats[i].sat[0].time == sat->time) { + reuse = 1; + break; + } + if (data->sats[i].sat[0].time < sat->time) + break; + } + if (!reuse) { + if (data->numsats == data->sizesats) { + if (data->sizesats == 0) + newsats = malloc((newsizesats = 256) * sizeof (struct cc_gpssats)); + else + newsats = realloc (data->sats, (newsizesats = data->sizesats * 2) + * sizeof (struct cc_gpssats)); + if (!newsats) + return 0; + data->sats = newsats; + data->sizesats = newsizesats; + } + i = data->numsats++; + data->sats[i].nsat = 0; + } + j = data->sats[i].nsat; + if (j < 12) { + data->sats[i].sat[j] = *sat; + data->sats[i].nsat = j + 1; + } + return 1; +} + +static void +gpsdata_free(struct cc_gpsdata *data) +{ + if (data->data) + free(data->data); +} + +#define AO_LOG_FLIGHT 'F' +#define AO_LOG_SENSOR 'A' +#define AO_LOG_TEMP_VOLT 'T' +#define AO_LOG_DEPLOY 'D' +#define AO_LOG_STATE 'S' +#define AO_LOG_GPS_TIME 'G' +#define AO_LOG_GPS_LAT 'N' +#define AO_LOG_GPS_LON 'W' +#define AO_LOG_GPS_ALT 'H' +#define AO_LOG_GPS_SAT 'V' +#define AO_LOG_GPS_DATE 'Y' + +#define AO_LOG_POS_NONE (~0UL) + +#define GPS_TIME 1 +#define GPS_LAT 2 +#define GPS_LON 4 +#define GPS_ALT 8 + +static int +read_eeprom(const char *line, struct cc_flightraw *f, double *ground_pres, int *ground_pres_count) +{ + char type; + int tick; + int a, b; + static struct cc_gpselt gps; + static int gps_valid; + struct cc_gpssat sat; + int serial; + + if (sscanf(line, "serial-number %u", &serial) == 1) { + f->serial = serial; + return 1; + } + if (sscanf(line, "%c %x %x %x", &type, &tick, &a, &b) != 4) + return 0; + switch (type) { + case AO_LOG_FLIGHT: + f->ground_accel = a; + f->ground_pres = 0; + f->flight = b; + *ground_pres = 0; + *ground_pres_count = 0; + break; + case AO_LOG_SENSOR: + timedata_add(&f->accel, tick, a); + timedata_add(&f->pres, tick, b); + if (*ground_pres_count < 20) { + *ground_pres += b; + (*ground_pres_count)++; + if (*ground_pres_count >= 20) + f->ground_pres = *ground_pres / *ground_pres_count; + } + break; + case AO_LOG_TEMP_VOLT: + timedata_add(&f->temp, tick, a); + timedata_add(&f->volt, tick, b); + break; + case AO_LOG_DEPLOY: + timedata_add(&f->drogue, tick, a); + timedata_add(&f->main, tick, b); + break; + case AO_LOG_STATE: + timedata_add(&f->state, tick, a); + break; + case AO_LOG_GPS_TIME: + /* the flight computer writes TIME first, so reset + * any stale data before adding this record + */ + gps.time = tick; + gps.hour = (a & 0xff); + gps.minute = (a >> 8) & 0xff; + gps.second = (b & 0xff); + gps.flags = (b >> 8) & 0xff; + gps_valid = GPS_TIME; + break; + case AO_LOG_GPS_LAT: + gps.lat = ((int32_t) (a + (b << 16))) / 10000000.0; + gps_valid |= GPS_LAT; + break; + case AO_LOG_GPS_LON: + gps.lon = ((int32_t) (a + (b << 16))) / 10000000.0; + gps_valid |= GPS_LON; + break; + case AO_LOG_GPS_ALT: + gps.alt = (int16_t) a; + gps_valid |= GPS_ALT; + break; + case AO_LOG_GPS_SAT: + sat.time = tick; + sat.svid = a; + sat.c_n = (b >> 8) & 0xff; + gpssat_add(&f->gps, &sat); + break; + case AO_LOG_GPS_DATE: + f->year = 2000 + (a & 0xff); + f->month = (a >> 8) & 0xff; + f->day = (b & 0xff); + break; + default: + return 0; + } + if (gps_valid == 0xf) { + gps_valid = 0; + gpsdata_add(&f->gps, &gps); + } + return 1; +} + +static const char *state_names[] = { + "startup", + "idle", + "pad", + "boost", + "fast", + "coast", + "drogue", + "main", + "landed", + "invalid" +}; + +static enum ao_flight_state +state_name_to_state(char *state_name) +{ + enum ao_flight_state state; + for (state = ao_flight_startup; state < ao_flight_invalid; state++) + if (!strcmp(state_names[state], state_name)) + return state; + return ao_flight_invalid; +} + +static int +read_telem(const char *line, struct cc_flightraw *f) +{ + struct cc_telem telem; + struct cc_gpselt gps; + struct cc_gpssat sat; + int s; + + if (!cc_telem_parse(line, &telem)) + return 0; + f->ground_accel = telem.ground_accel; + f->ground_pres = telem.ground_pres; + f->flight = 0; + timedata_add(&f->accel, telem.tick, telem.flight_accel); + timedata_add(&f->pres, telem.tick, telem.flight_pres); + timedata_add(&f->temp, telem.tick, telem.temp); + timedata_add(&f->volt, telem.tick, telem.batt); + timedata_add(&f->drogue, telem.tick, telem.drogue); + timedata_add(&f->main, telem.tick, telem.main); + timedata_add(&f->state, telem.tick, state_name_to_state(telem.state)); + if (telem.gps.gps_locked) { + f->year = telem.gps.gps_time.year; + f->month = telem.gps.gps_time.month; + f->day = telem.gps.gps_time.day; + gps.time = telem.tick; + gps.lat = telem.gps.lat; + gps.lon = telem.gps.lon; + gps.alt = telem.gps.alt; + gps.hour = telem.gps.gps_time.hour; + gps.minute = telem.gps.gps_time.minute; + gps.second = telem.gps.gps_time.second; + gpsdata_add(&f->gps, &gps); + } + for (s = 0; s < telem.gps_tracking.channels; s++) { + sat.time = telem.tick; + sat.svid = telem.gps_tracking.sats[s].svid; + sat.c_n = telem.gps_tracking.sats[s].c_n0; + gpssat_add(&f->gps, &sat); + } + return 1; +} + +struct cc_flightraw * +cc_log_read(FILE *file) +{ + struct cc_flightraw *f; + char line[8192]; + double ground_pres; + int ground_pres_count; + + f = calloc(1, sizeof (struct cc_flightraw)); + if (!f) + return NULL; + while (fgets(line, sizeof (line), file)) { + if (read_eeprom(line, f, &ground_pres, &ground_pres_count)) + continue; + if (read_telem(line, f)) + continue; + fprintf (stderr, "invalid line: %s", line); + } + return f; +} + +void +cc_flightraw_free(struct cc_flightraw *raw) +{ + timedata_free(&raw->accel); + timedata_free(&raw->pres); + timedata_free(&raw->temp); + timedata_free(&raw->volt); + timedata_free(&raw->main); + timedata_free(&raw->drogue); + timedata_free(&raw->state); + gpsdata_free(&raw->gps); + free(raw); +} diff --git a/ao-tools/lib/cc-period.c b/ao-tools/lib/cc-period.c new file mode 100644 index 00000000..844ac79e --- /dev/null +++ b/ao-tools/lib/cc-period.c @@ -0,0 +1,56 @@ +/* + * Copyright © 2009 Keith Packard <keithp@keithp.com> + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; version 2 of the License. + * + * This program is distributed in the hope that it will be useful, but + * WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * General Public License for more details. + * + * You should have received a copy of the GNU General Public License along + * with this program; if not, write to the Free Software Foundation, Inc., + * 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA. + */ + +#include "cc.h" +#include <stdlib.h> +#include <math.h> + +struct cc_perioddata * +cc_period_make(struct cc_timedata *td, double start_time, double stop_time) +{ + int len = stop_time - start_time + 1; + struct cc_perioddata *pd; + int i, j; + double t; + + pd = calloc(1, sizeof (struct cc_perioddata)); + pd->start = start_time; + pd->step = 1; + pd->num = len; + pd->data = calloc(len, sizeof(double)); + j = 0; + for (i = 0; i < pd->num; i++) { + t = start_time + i * pd->step; + while (j < td->num - 1 && fabs(t - td->data[j].time) >= fabs(t - td->data[j+1].time)) + j++; + pd->data[i] = td->data[j].value; + } + return pd; +} + +struct cc_perioddata * +cc_period_low_pass(struct cc_perioddata *raw, double omega_pass, double omega_stop, double error) +{ + struct cc_perioddata *filtered; + + filtered = calloc (1, sizeof (struct cc_perioddata)); + filtered->start = raw->start; + filtered->step = raw->step; + filtered->num = raw->num; + filtered->data = cc_low_pass(raw->data, raw->num, omega_pass, omega_stop, error); + return filtered; +} diff --git a/ao-tools/lib/cc-process.c b/ao-tools/lib/cc-process.c new file mode 100644 index 00000000..c756b570 --- /dev/null +++ b/ao-tools/lib/cc-process.c @@ -0,0 +1,166 @@ +/* + * Copyright © 2009 Keith Packard <keithp@keithp.com> + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; version 2 of the License. + * + * This program is distributed in the hope that it will be useful, but + * WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * General Public License for more details. + * + * You should have received a copy of the GNU General Public License along + * with this program; if not, write to the Free Software Foundation, Inc., + * 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA. + */ + +#include "cc.h" +#include <stdlib.h> +#include <math.h> +#include <string.h> + +static void +cook_timed(struct cc_timedata *td, struct cc_perioddata *pd, + double start_time, double stop_time, + double omega_pass, double omega_stop, double error) +{ + struct cc_perioddata *unfiltered, *filtered; + + unfiltered = cc_period_make(td, start_time, stop_time); + filtered = cc_period_low_pass (unfiltered, omega_pass, omega_stop, error); + *pd = *filtered; + free (filtered); + free (unfiltered->data); + free (unfiltered); +} + +static double +barometer_to_altitude(double b, double pad_alt) +{ + return cc_barometer_to_altitude(b) - pad_alt; +} + +struct cc_flightcooked * +cc_flight_cook(struct cc_flightraw *raw) +{ + struct cc_flightcooked *cooked; + double flight_start; + double flight_stop; + int start_set = 0; + int stop_set = 0; + int i; + struct cc_timedata *accel; + struct cc_timedata *accel_speed; + struct cc_timedata *accel_pos; + struct cc_timedata *pres; + struct cc_perioddata *pres_speed; + struct cc_perioddata *pres_accel; + + if (raw->accel.num == 0) + return NULL; + + cooked = calloc (1, sizeof (struct cc_flightcooked)); + + /* + * Find flight start and stop times by looking at + * state transitions. The stop time is set to the time + * of landing, which may be long after it landed (due to radio + * issues). Refine this value by looking through the sensor data + */ + for (i = 0; i < raw->state.num; i++) { + if (!start_set && raw->state.data[i].value > ao_flight_pad) { + flight_start = raw->state.data[i].time - 10; + start_set = 1; + } + if (!stop_set && raw->state.data[i].value > ao_flight_main) { + flight_stop = raw->state.data[i].time; + stop_set = 1; + } + } + + if (!start_set || flight_start < raw->accel.data[0].time) + flight_start = raw->accel.data[0].time; + if (stop_set) { + for (i = 0; i < raw->accel.num - 1; i++) { + if (raw->accel.data[i+1].time >= flight_stop) { + flight_stop = raw->accel.data[i].time; + break; + } + } + } else { + flight_stop = raw->accel.data[raw->accel.num-1].time; + } + cooked->flight_start = flight_start; + cooked->flight_stop = flight_stop; + + /* Integrate the accelerometer data to get speed and position */ + accel = cc_timedata_convert(&raw->accel, cc_accelerometer_to_acceleration, raw->ground_accel); + cooked->accel = *accel; + free(accel); + accel_speed = cc_timedata_integrate(&cooked->accel, flight_start - 10, flight_stop); + accel_pos = cc_timedata_integrate(accel_speed, flight_start - 10, flight_stop); + +#define ACCEL_OMEGA_PASS (2 * M_PI * 20 / 100) +#define ACCEL_OMEGA_STOP (2 * M_PI * 30 / 100) +#define BARO_OMEGA_PASS (2 * M_PI * .5 / 100) +#define BARO_OMEGA_STOP (2 * M_PI * 1 / 100) +#define FILTER_ERROR (1e-8) + + cook_timed(&cooked->accel, &cooked->accel_accel, + flight_start, flight_stop, + ACCEL_OMEGA_PASS, ACCEL_OMEGA_STOP, FILTER_ERROR); + cook_timed(accel_speed, &cooked->accel_speed, + flight_start, flight_stop, + ACCEL_OMEGA_PASS, ACCEL_OMEGA_STOP, FILTER_ERROR); + free(accel_speed->data); free(accel_speed); + cook_timed(accel_pos, &cooked->accel_pos, + flight_start, flight_stop, + ACCEL_OMEGA_PASS, ACCEL_OMEGA_STOP, FILTER_ERROR); + free(accel_pos->data); free(accel_pos); + + /* Filter the pressure data */ + pres = cc_timedata_convert(&raw->pres, barometer_to_altitude, + cc_barometer_to_altitude(raw->ground_pres)); + cooked->pres = *pres; + free(pres); + cook_timed(&cooked->pres, &cooked->pres_pos, + flight_start, flight_stop, + BARO_OMEGA_PASS, BARO_OMEGA_STOP, FILTER_ERROR); + /* differentiate twice to get to acceleration */ + pres_speed = cc_perioddata_differentiate(&cooked->pres_pos); + pres_accel = cc_perioddata_differentiate(pres_speed); + + cooked->pres_speed = *pres_speed; + free(pres_speed); + cooked->pres_accel = *pres_accel; + free(pres_accel); + + /* copy state */ + cooked->state.num = raw->state.num; + cooked->state.size = raw->state.num; + cooked->state.data = calloc(cooked->state.num, sizeof (struct cc_timedataelt)); + memcpy(cooked->state.data, raw->state.data, cooked->state.num * sizeof (struct cc_timedataelt)); + cooked->state.time_offset = raw->state.time_offset; + return cooked; +} + +#define if_free(x) ((x) ? free(x) : (void) 0) + +void +cc_flightcooked_free(struct cc_flightcooked *cooked) +{ + if_free(cooked->accel_accel.data); + if_free(cooked->accel_speed.data); + if_free(cooked->accel_pos.data); + if_free(cooked->pres_pos.data); + if_free(cooked->pres_speed.data); + if_free(cooked->pres_accel.data); + if_free(cooked->gps_lat.data); + if_free(cooked->gps_lon.data); + if_free(cooked->gps_alt.data); + if_free(cooked->state.data); + if_free(cooked->accel.data); + if_free(cooked->pres.data); + free(cooked); +} diff --git a/ao-tools/lib/cc-telem.c b/ao-tools/lib/cc-telem.c new file mode 100644 index 00000000..aa52b7c5 --- /dev/null +++ b/ao-tools/lib/cc-telem.c @@ -0,0 +1,212 @@ +/* + * Copyright © 2009 Keith Packard <keithp@keithp.com> + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; version 2 of the License. + * + * This program is distributed in the hope that it will be useful, but + * WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * General Public License for more details. + * + * You should have received a copy of the GNU General Public License along + * with this program; if not, write to the Free Software Foundation, Inc., + * 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA. + */ + +#include "cc.h" +#include <string.h> +#include <stdlib.h> + +static void +cc_parse_string(char *target, int len, char *source) +{ + strncpy(target, source, len-1); + target[len-1] = '\0'; +} + +static void +cc_parse_int(int *target, char *source) +{ + *target = strtol(source, NULL, 0); +} + +static void +cc_parse_hex(int *target, char *source) +{ + *target = strtol(source, NULL, 16); +} + +static void +cc_parse_pos(double *target, char *source) +{ + int deg; + double min; + char dir; + double r; + + if (sscanf(source, "%d°%lf'%c", °, &min, &dir) != 3) { + *target = 0; + return; + } + r = deg + min / 60.0; + if (dir == 'S' || dir == 'W') + r = -r; + *target = r; +} + +#define PARSE_MAX_WORDS 512 + +int +cc_telem_parse(const char *input_line, struct cc_telem *telem) +{ + char *saveptr; + char *raw_words[PARSE_MAX_WORDS]; + char **words; + int version = 0; + int nword; + char line_buf[8192], *line; + int tracking_pos; + + /* avoid smashing our input parameter */ + strncpy (line_buf, input_line, sizeof (line_buf)-1); + line_buf[sizeof(line_buf) - 1] = '\0'; + line = line_buf; + for (nword = 0; nword < PARSE_MAX_WORDS; nword++) { + raw_words[nword] = strtok_r(line, " \t\n", &saveptr); + line = NULL; + if (raw_words[nword] == NULL) + break; + } + if (nword < 36) + return FALSE; + words = raw_words; + if (strcmp(words[0], "VERSION") == 0) { + cc_parse_int(&version, words[1]); + words += 2; + nword -= 2; + } + + if (strcmp(words[0], "CALL") != 0) + return FALSE; + cc_parse_string(telem->callsign, sizeof (telem->callsign), words[1]); + cc_parse_int(&telem->serial, words[3]); + + if (version >= 2) { + cc_parse_int(&telem->flight, words[5]); + words += 2; + nword -= 2; + } else + telem->flight = 0; + + cc_parse_int(&telem->rssi, words[5]); + if (version <= 2) { + /* Older telemetry versions mis-computed the rssi value */ + telem->rssi = (telem->rssi + 74) / 2 - 74; + } + cc_parse_string(telem->state, sizeof (telem->state), words[9]); + cc_parse_int(&telem->tick, words[10]); + cc_parse_int(&telem->accel, words[12]); + cc_parse_int(&telem->pres, words[14]); + cc_parse_int(&telem->temp, words[16]); + cc_parse_int(&telem->batt, words[18]); + cc_parse_int(&telem->drogue, words[20]); + cc_parse_int(&telem->main, words[22]); + cc_parse_int(&telem->flight_accel, words[24]); + cc_parse_int(&telem->ground_accel, words[26]); + cc_parse_int(&telem->flight_vel, words[28]); + cc_parse_int(&telem->flight_pres, words[30]); + cc_parse_int(&telem->ground_pres, words[32]); + if (version >= 1) { + cc_parse_int(&telem->accel_plus_g, words[34]); + cc_parse_int(&telem->accel_minus_g, words[36]); + words += 4; + nword -= 4; + } else { + telem->accel_plus_g = telem->ground_accel; + telem->accel_minus_g = telem->ground_accel + 530; + } + cc_parse_int(&telem->gps.nsat, words[34]); + if (strcmp (words[36], "unlocked") == 0) { + telem->gps.gps_connected = 1; + telem->gps.gps_locked = 0; + telem->gps.gps_time.year = telem->gps.gps_time.month = telem->gps.gps_time.day = 0; + telem->gps.gps_time.hour = telem->gps.gps_time.minute = telem->gps.gps_time.second = 0; + telem->gps.lat = telem->gps.lon = 0; + telem->gps.alt = 0; + tracking_pos = 37; + } else if (nword >= 40) { + telem->gps.gps_locked = 1; + telem->gps.gps_connected = 1; + if (version >= 2) { + sscanf(words[36], "%d-%d-%d", + &telem->gps.gps_time.year, + &telem->gps.gps_time.month, + &telem->gps.gps_time.day); + words += 1; + nword -= 1; + } else { + telem->gps.gps_time.year = telem->gps.gps_time.month = telem->gps.gps_time.day = 0; + } + sscanf(words[36], "%d:%d:%d", &telem->gps.gps_time.hour, &telem->gps.gps_time.minute, &telem->gps.gps_time.second); + cc_parse_pos(&telem->gps.lat, words[37]); + cc_parse_pos(&telem->gps.lon, words[38]); + sscanf(words[39], "%dm", &telem->gps.alt); + tracking_pos = 46; + } else { + telem->gps.gps_connected = 0; + telem->gps.gps_locked = 0; + telem->gps.gps_time.year = telem->gps.gps_time.month = telem->gps.gps_time.day = 0; + telem->gps.gps_time.hour = telem->gps.gps_time.minute = telem->gps.gps_time.second = 0; + telem->gps.lat = telem->gps.lon = 0; + telem->gps.alt = 0; + tracking_pos = -1; + } + if (nword >= 46) { + telem->gps.gps_extended = 1; + sscanf(words[40], "%lfm/s", &telem->gps.ground_speed); + sscanf(words[41], "%d", &telem->gps.course); + sscanf(words[42], "%lfm/s", &telem->gps.climb_rate); + sscanf(words[43], "%lf", &telem->gps.hdop); + sscanf(words[44], "%d", &telem->gps.h_error); + sscanf(words[45], "%d", &telem->gps.v_error); + } else { + telem->gps.gps_extended = 0; + telem->gps.ground_speed = 0; + telem->gps.course = 0; + telem->gps.climb_rate = 0; + telem->gps.hdop = 0; + telem->gps.h_error = 0; + telem->gps.v_error = 0; + } + if (tracking_pos >= 0 && nword >= tracking_pos + 2 && strcmp(words[tracking_pos], "SAT") == 0) { + int c, n, pos; + int per_sat; + int state; + + if (version >= 2) + per_sat = 2; + else + per_sat = 3; + cc_parse_int(&n, words[tracking_pos + 1]); + pos = tracking_pos + 2; + if (nword >= pos + n * per_sat) { + telem->gps_tracking.channels = n; + for (c = 0; c < n; c++) { + cc_parse_int(&telem->gps_tracking.sats[c].svid, + words[pos + 0]); + if (version < 2) + cc_parse_hex(&state, words[pos + 1]); + cc_parse_int(&telem->gps_tracking.sats[c].c_n0, + words[pos + per_sat - 1]); + pos += per_sat; + } + } else { + telem->gps_tracking.channels = 0; + } + } else { + telem->gps_tracking.channels = 0; + } + return TRUE; +} diff --git a/ao-tools/lib/cc-telemetry.c b/ao-tools/lib/cc-telemetry.c new file mode 100644 index 00000000..99da2680 --- /dev/null +++ b/ao-tools/lib/cc-telemetry.c @@ -0,0 +1,62 @@ +/* + * Copyright © 2011 Keith Packard <keithp@keithp.com> + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; version 2 of the License. + * + * This program is distributed in the hope that it will be useful, but + * WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * General Public License for more details. + * + * You should have received a copy of the GNU General Public License along + * with this program; if not, write to the Free Software Foundation, Inc., + * 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA. + */ + +#include "cc.h" +#include <string.h> + +static int +parse_byte(char *data, uint8_t *byte) +{ + char d[3]; + int x; + d[0] = data[0]; + d[1] = data[1]; + d[2] = '\0'; + + if (sscanf(d, "%x", &x) != 1) + return FALSE; + *byte = x; + return TRUE; +} + +int +cc_telemetry_parse(const char *input_line, union ao_telemetry_all *telemetry) +{ + uint8_t *byte; + char *data; + uint8_t hex[35]; + int i; + + if (strncmp(input_line, "TELEM", 5) != 0) + return FALSE; + + data = strchr (input_line, ' '); + if (!data) + return FALSE; + data++; + byte = hex; + for (i = 0; i < 35; i++) { + if (!parse_byte(data, byte)) + return FALSE; + data += 2; + byte++; + } + if (hex[0] != 34) + return FALSE; + memcpy(telemetry, hex+1, 34); + return TRUE; +} diff --git a/ao-tools/lib/cc-telemetry.h b/ao-tools/lib/cc-telemetry.h new file mode 100644 index 00000000..e849cd3b --- /dev/null +++ b/ao-tools/lib/cc-telemetry.h @@ -0,0 +1,243 @@ +/* + * Copyright © 2011 Keith Packard <keithp@keithp.com> + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; version 2 of the License. + * + * This program is distributed in the hope that it will be useful, but + * WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * General Public License for more details. + * + * You should have received a copy of the GNU General Public License along + * with this program; if not, write to the Free Software Foundation, Inc., + * 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA. + */ + +#ifndef _CC_TELEMETRY_H_ +#define _CC_TELEMETRY_H_ +/* + * ao_telemetry.c + */ +#define AO_MAX_CALLSIGN 8 +#define AO_MAX_VERSION 8 +#define AO_MAX_TELEMETRY 128 + +struct ao_telemetry_generic { + uint16_t serial; /* 0 */ + uint16_t tick; /* 2 */ + uint8_t type; /* 4 */ + uint8_t payload[27]; /* 5 */ + uint8_t rssi; /* 32 */ + uint8_t status; /* 33 */ + /* 34 */ +}; + +#define AO_TELEMETRY_SENSOR_TELEMETRUM 0x01 +#define AO_TELEMETRY_SENSOR_TELEMINI 0x02 +#define AO_TELEMETRY_SENSOR_TELENANO 0x03 + +struct ao_telemetry_sensor { + uint16_t serial; /* 0 */ + uint16_t tick; /* 2 */ + uint8_t type; /* 4 */ + + uint8_t state; /* 5 flight state */ + int16_t accel; /* 6 accelerometer (TM only) */ + int16_t pres; /* 8 pressure sensor */ + int16_t temp; /* 10 temperature sensor */ + int16_t v_batt; /* 12 battery voltage */ + int16_t sense_d; /* 14 drogue continuity sense (TM/Tm) */ + int16_t sense_m; /* 16 main continuity sense (TM/Tm) */ + + int16_t acceleration; /* 18 m/s² * 16 */ + int16_t speed; /* 20 m/s * 16 */ + int16_t height; /* 22 m */ + + int16_t ground_pres; /* 24 average pres on pad */ + int16_t ground_accel; /* 26 average accel on pad */ + int16_t accel_plus_g; /* 28 accel calibration at +1g */ + int16_t accel_minus_g; /* 30 accel calibration at -1g */ + /* 32 */ +}; + +#define AO_TELEMETRY_CONFIGURATION 0x04 + +struct ao_telemetry_configuration { + uint16_t serial; /* 0 */ + uint16_t tick; /* 2 */ + uint8_t type; /* 4 */ + + uint8_t device; /* 5 device type */ + uint16_t flight; /* 6 flight number */ + uint8_t config_major; /* 8 Config major version */ + uint8_t config_minor; /* 9 Config minor version */ + uint16_t apogee_delay; /* 10 Apogee deploy delay in seconds */ + uint16_t main_deploy; /* 12 Main deploy alt in meters */ + uint16_t flight_log_max; /* 14 Maximum flight log size in kB */ + char callsign[AO_MAX_CALLSIGN]; /* 16 Radio operator identity */ + char version[AO_MAX_VERSION]; /* 24 Software version */ + /* 32 */ +}; + +#define AO_TELEMETRY_LOCATION 0x05 + +#define AO_GPS_MODE_NOT_VALID 'N' +#define AO_GPS_MODE_AUTONOMOUS 'A' +#define AO_GPS_MODE_DIFFERENTIAL 'D' +#define AO_GPS_MODE_ESTIMATED 'E' +#define AO_GPS_MODE_MANUAL 'M' +#define AO_GPS_MODE_SIMULATED 'S' + +struct ao_telemetry_location { + uint16_t serial; /* 0 */ + uint16_t tick; /* 2 */ + uint8_t type; /* 4 */ + + uint8_t flags; /* 5 Number of sats and other flags */ + int16_t altitude; /* 6 GPS reported altitude (m) */ + int32_t latitude; /* 8 latitude (degrees * 10⁷) */ + int32_t longitude; /* 12 longitude (degrees * 10⁷) */ + uint8_t year; /* 16 (- 2000) */ + uint8_t month; /* 17 (1-12) */ + uint8_t day; /* 18 (1-31) */ + uint8_t hour; /* 19 (0-23) */ + uint8_t minute; /* 20 (0-59) */ + uint8_t second; /* 21 (0-59) */ + uint8_t pdop; /* 22 (m * 5) */ + uint8_t hdop; /* 23 (m * 5) */ + uint8_t vdop; /* 24 (m * 5) */ + uint8_t mode; /* 25 */ + uint16_t ground_speed; /* 26 cm/s */ + int16_t climb_rate; /* 28 cm/s */ + uint8_t course; /* 30 degrees / 2 */ + uint8_t unused[1]; /* 31 */ + /* 32 */ +}; + +#define AO_TELEMETRY_SATELLITE 0x06 + +struct ao_telemetry_satellite_info { + uint8_t svid; + uint8_t c_n_1; +}; + +struct ao_telemetry_satellite { + uint16_t serial; /* 0 */ + uint16_t tick; /* 2 */ + uint8_t type; /* 4 */ + uint8_t channels; /* 5 number of reported sats */ + + struct ao_telemetry_satellite_info sats[12]; /* 6 */ + uint8_t unused[2]; /* 30 */ + /* 32 */ +}; + +#define AO_TELEMETRY_COMPANION 0x07 + +#define AO_COMPANION_MAX_CHANNELS 12 + +struct ao_telemetry_companion { + uint16_t serial; /* 0 */ + uint16_t tick; /* 2 */ + uint8_t type; /* 4 */ + uint8_t board_id; /* 5 */ + + uint8_t update_period; /* 6 */ + uint8_t channels; /* 7 */ + uint16_t companion_data[AO_COMPANION_MAX_CHANNELS]; /* 8 */ + /* 32 */ +}; + +#define AO_TELEMETRY_MEGA_SENSOR 0x08 + +struct ao_telemetry_mega_sensor { + uint16_t serial; /* 0 */ + uint16_t tick; /* 2 */ + uint8_t type; /* 4 */ + + uint8_t pad5; /* 5 */ + int16_t accel; /* 6 Z axis */ + + int32_t pres; /* 8 Pa * 10 */ + int16_t temp; /* 12 °C * 100 */ + + int16_t accel_x; /* 14 */ + int16_t accel_y; /* 16 */ + int16_t accel_z; /* 18 */ + + int16_t gyro_x; /* 20 */ + int16_t gyro_y; /* 22 */ + int16_t gyro_z; /* 24 */ + + int16_t mag_x; /* 26 */ + int16_t mag_y; /* 28 */ + int16_t mag_z; /* 30 */ + /* 32 */ +}; + +#define AO_TELEMETRY_MEGA_DATA 0x09 + +struct ao_telemetry_mega_data { + uint16_t serial; /* 0 */ + uint16_t tick; /* 2 */ + uint8_t type; /* 4 */ + + uint8_t state; /* 5 flight state */ + + int16_t v_batt; /* 6 battery voltage */ + int16_t v_pyro; /* 8 pyro battery voltage */ + int8_t sense[6]; /* 10 continuity sense */ + + int32_t ground_pres; /* 16 average pres on pad */ + int16_t ground_accel; /* 20 average accel on pad */ + int16_t accel_plus_g; /* 22 accel calibration at +1g */ + int16_t accel_minus_g; /* 24 accel calibration at -1g */ + + int16_t acceleration; /* 26 m/s² * 16 */ + int16_t speed; /* 28 m/s * 16 */ + int16_t height; /* 30 m */ + /* 32 */ +}; + +/* #define AO_SEND_ALL_BARO */ + +#define AO_TELEMETRY_BARO 0x80 + +/* + * This packet allows the full sampling rate baro + * data to be captured over the RF link so that the + * flight software can be tested using 'real' data. + * + * Along with this telemetry packet, the flight + * code is modified to send full-rate telemetry all the time + * and never send an RDF tone; this ensure that the full radio + * link is available. + */ +struct ao_telemetry_baro { + uint16_t serial; /* 0 */ + uint16_t tick; /* 2 */ + uint8_t type; /* 4 */ + uint8_t samples; /* 5 number samples */ + + int16_t baro[12]; /* 6 samples */ + /* 32 */ +}; + +union ao_telemetry_all { + struct ao_telemetry_generic generic; + struct ao_telemetry_sensor sensor; + struct ao_telemetry_configuration configuration; + struct ao_telemetry_location location; + struct ao_telemetry_satellite satellite; + struct ao_telemetry_companion companion; + struct ao_telemetry_mega_sensor mega_sensor; + struct ao_telemetry_mega_data mega_data; + struct ao_telemetry_baro baro; +}; + +int +cc_telemetry_parse(const char *input_line, union ao_telemetry_all *telemetry); + +#endif diff --git a/ao-tools/lib/cc-usb.c b/ao-tools/lib/cc-usb.c index 81309983..1580c6d9 100644 --- a/ao-tools/lib/cc-usb.c +++ b/ao-tools/lib/cc-usb.c @@ -30,27 +30,31 @@ #include "cc-usb.h" -#define CC_NUM_READ 16 +#define CC_NUM_HEX_READ 64 /* * AltOS has different buffer sizes for in/out packets */ -#define CC_IN_BUF 256 +#define CC_IN_BUF 65536 #define CC_OUT_BUF 64 #define DEFAULT_TTY "/dev/ttyACM0" -struct cc_read { +struct cc_hex_read { uint8_t *buf; int len; }; struct cc_usb { - int fd; - uint8_t in_buf[CC_IN_BUF]; - int in_count; - uint8_t out_buf[CC_OUT_BUF]; - int out_count; - struct cc_read read_buf[CC_NUM_READ]; - int read_count; + int fd; + uint8_t in_buf[CC_IN_BUF]; + int in_pos; + int in_count; + uint8_t out_buf[CC_OUT_BUF]; + int out_count; + + struct cc_hex_read hex_buf[CC_NUM_HEX_READ]; + int hex_count; + + int remote; }; #define NOT_HEX 0xff @@ -72,61 +76,48 @@ cc_hex_nibble(uint8_t c) * and write them to the waiting buffer */ static void -cc_handle_in(struct cc_usb *cc) +cc_handle_hex_read(struct cc_usb *cc) { uint8_t h, l; - int in_pos; - int read_pos; + int hex_pos; - in_pos = 0; - read_pos = 0; - while (read_pos < cc->read_count && in_pos < cc->in_count) { + hex_pos = 0; + while (hex_pos < cc->hex_count && cc->in_pos < cc->in_count) { /* * Skip to next hex character */ - while (in_pos < cc->in_count && - cc_hex_nibble(cc->in_buf[in_pos]) == NOT_HEX) - in_pos++; + while (cc->in_pos < cc->in_count && + cc_hex_nibble(cc->in_buf[cc->in_pos]) == NOT_HEX) + cc->in_pos++; /* * Make sure we have two characters left */ - if (cc->in_count - in_pos < 2) + if (cc->in_count - cc->in_pos < 2) break; /* * Parse hex number */ - h = cc_hex_nibble(cc->in_buf[in_pos]); - l = cc_hex_nibble(cc->in_buf[in_pos+1]); + h = cc_hex_nibble(cc->in_buf[cc->in_pos]); + l = cc_hex_nibble(cc->in_buf[cc->in_pos+1]); if (h == NOT_HEX || l == NOT_HEX) { fprintf(stderr, "hex read error\n"); break; } - in_pos += 2; + cc->in_pos += 2; /* * Store hex number */ - *cc->read_buf[read_pos].buf++ = (h << 4) | l; - if (--cc->read_buf[read_pos].len <= 0) - read_pos++; - } - - /* Move remaining bytes to the start of the input buffer */ - if (in_pos) { - memmove(cc->in_buf, cc->in_buf + in_pos, - cc->in_count - in_pos); - cc->in_count -= in_pos; + *cc->hex_buf[hex_pos].buf++ = (h << 4) | l; + if (--cc->hex_buf[hex_pos].len <= 0) + hex_pos++; } - /* Move pending reads to the start of the array */ - if (read_pos) { - memmove(cc->read_buf, cc->read_buf + read_pos, - (cc->read_count - read_pos) * sizeof (cc->read_buf[0])); - cc->read_count -= read_pos; + /* Move pending hex reads to the start of the array */ + if (hex_pos) { + memmove(cc->hex_buf, cc->hex_buf + hex_pos, + (cc->hex_count - hex_pos) * sizeof (cc->hex_buf[0])); + cc->hex_count -= hex_pos; } - - /* Once we're done reading, flush any pending input */ - if (cc->read_count == 0) - cc->in_count = 0; } static void @@ -157,8 +148,9 @@ cc_usb_dbg(int indent, uint8_t *bytes, int len) /* * Flush pending writes, fill pending reads */ -void -cc_usb_sync(struct cc_usb *cc) + +static int +_cc_usb_sync(struct cc_usb *cc, int wait_for_input) { int ret; struct pollfd fds; @@ -166,21 +158,33 @@ cc_usb_sync(struct cc_usb *cc) fds.fd = cc->fd; for (;;) { - if (cc->read_count || cc->out_count) - timeout = -1; + if (cc->hex_count || cc->out_count) + timeout = 5000; + else if (wait_for_input && cc->in_pos == cc->in_count) + timeout = wait_for_input; else timeout = 0; fds.events = 0; + /* Move remaining bytes to the start of the input buffer */ + if (cc->in_pos) { + memmove(cc->in_buf, cc->in_buf + cc->in_pos, + cc->in_count - cc->in_pos); + cc->in_count -= cc->in_pos; + cc->in_pos = 0; + } if (cc->in_count < CC_IN_BUF) fds.events |= POLLIN; if (cc->out_count) fds.events |= POLLOUT; ret = poll(&fds, 1, timeout); - if (ret == 0) + if (ret == 0) { + if (timeout) + return -1; break; + } if (ret < 0) { perror("poll"); - break; + return -1; } if (fds.revents & POLLIN) { ret = read(cc->fd, cc->in_buf + cc->in_count, @@ -188,7 +192,8 @@ cc_usb_sync(struct cc_usb *cc) if (ret > 0) { cc_usb_dbg(24, cc->in_buf + cc->in_count, ret); cc->in_count += ret; - cc_handle_in(cc); + if (cc->hex_count) + cc_handle_hex_read(cc); } else if (ret < 0) perror("read"); } @@ -205,6 +210,16 @@ cc_usb_sync(struct cc_usb *cc) perror("write"); } } + return 0; +} + +void +cc_usb_sync(struct cc_usb *cc) +{ + if (_cc_usb_sync(cc, 0) < 0) { + fprintf(stderr, "USB link timeout\n"); + exit(1); + } } void @@ -239,6 +254,38 @@ cc_usb_printf(struct cc_usb *cc, char *format, ...) } int +cc_usb_getchar(struct cc_usb *cc) +{ + while (cc->in_pos == cc->in_count) { + if (_cc_usb_sync(cc, 5000) < 0) { + fprintf(stderr, "USB link timeout\n"); + exit(1); + } + } + return cc->in_buf[cc->in_pos++]; +} + +void +cc_usb_getline(struct cc_usb *cc, char *line, int max) +{ + int c; + + while ((c = cc_usb_getchar(cc)) != '\n') { + switch (c) { + case '\r': + break; + default: + if (max > 1) { + *line++ = c; + max--; + } + break; + } + } + *line++ = '\0'; +} + +int cc_usb_send_bytes(struct cc_usb *cc, uint8_t *bytes, int len) { int this_len; @@ -260,12 +307,18 @@ cc_usb_send_bytes(struct cc_usb *cc, uint8_t *bytes, int len) void cc_queue_read(struct cc_usb *cc, uint8_t *buf, int len) { - struct cc_read *read_buf; - while (cc->read_count >= CC_NUM_READ) + struct cc_hex_read *hex_buf; + + /* At the start of a command sequence, flush any pending input */ + if (cc->hex_count == 0) { + cc_usb_sync(cc); + cc->in_count = 0; + } + while (cc->hex_count >= CC_NUM_HEX_READ) cc_usb_sync(cc); - read_buf = &cc->read_buf[cc->read_count++]; - read_buf->buf = buf; - read_buf->len = len; + hex_buf = &cc->hex_buf[cc->hex_count++]; + hex_buf->buf = buf; + hex_buf->len = len; } int @@ -321,6 +374,29 @@ cc_usb_reset(struct cc_usb *cc) return 1; } +void +cc_usb_open_remote(struct cc_usb *cc, int channel) +{ + if (!cc->remote) { + printf ("channel %d\n", channel); + cc_usb_printf(cc, "\nc r %d\np\nE 0\n", channel); + do { + cc->in_count = cc->in_pos = 0; + _cc_usb_sync(cc, 100); + } while (cc->in_count > 0); + cc->remote = 1; + } +} + +void +cc_usb_close_remote(struct cc_usb *cc) +{ + if (cc->remote) { + cc_usb_printf(cc, "~"); + cc->remote = 0; + } +} + static struct termios save_termios; struct cc_usb * @@ -344,16 +420,19 @@ cc_usb_open(char *tty) save_termios = termios; cfmakeraw(&termios); tcsetattr(cc->fd, TCSAFLUSH, &termios); - cc_usb_printf(cc, "E 0\nm 0\n"); - cc_usb_sync(cc); - sleep(1); - cc_usb_sync(cc); + cc_usb_printf(cc, "\nE 0\nm 0\n"); + do { + cc->in_count = cc->in_pos = 0; + _cc_usb_sync(cc, 100); + } while (cc->in_count > 0); return cc; } void cc_usb_close(struct cc_usb *cc) { + cc_usb_close_remote(cc); + cc_usb_sync(cc); tcsetattr(cc->fd, TCSAFLUSH, &save_termios); close (cc->fd); free (cc); diff --git a/ao-tools/lib/cc-usb.h b/ao-tools/lib/cc-usb.h index d7acfbd2..d3539281 100644 --- a/ao-tools/lib/cc-usb.h +++ b/ao-tools/lib/cc-usb.h @@ -53,7 +53,19 @@ cc_usb_sync(struct cc_usb *cc); void cc_queue_read(struct cc_usb *cc, uint8_t *buf, int len); +int +cc_usb_getchar(struct cc_usb *cc); + +void +cc_usb_getline(struct cc_usb *cc, char *line, int max); + void cc_usb_printf(struct cc_usb *cc, char *format, ...); +void +cc_usb_open_remote(struct cc_usb *cc, int channel); + +void +cc_usb_close_remote(struct cc_usb *cc); + #endif /* _CC_USB_H_ */ diff --git a/ao-tools/lib/cc-usbdev.c b/ao-tools/lib/cc-usbdev.c new file mode 100644 index 00000000..a19e231c --- /dev/null +++ b/ao-tools/lib/cc-usbdev.c @@ -0,0 +1,317 @@ +/* + * Copyright © 2009 Keith Packard <keithp@keithp.com> + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 2 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, but + * WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * General Public License for more details. + * + * You should have received a copy of the GNU General Public License along + * with this program; if not, write to the Free Software Foundation, Inc., + * 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA. + */ + +#define _GNU_SOURCE +#include "cc.h" +#include <ctype.h> +#include <dirent.h> +#include <stdio.h> +#include <stdlib.h> +#include <string.h> + +static char * +load_string(char *dir, char *file) +{ + char *full = cc_fullname(dir, file); + char line[4096]; + char *r; + FILE *f; + int rlen; + + f = fopen(full, "r"); + free(full); + if (!f) + return NULL; + r = fgets(line, sizeof (line), f); + fclose(f); + if (!r) + return NULL; + rlen = strlen(r); + if (r[rlen-1] == '\n') + r[rlen-1] = '\0'; + return strdup(r); +} + +static int +load_hex(char *dir, char *file) +{ + char *line; + char *end; + long i; + + line = load_string(dir, file); + if (!line) + return -1; + i = strtol(line, &end, 16); + free(line); + if (end == line) + return -1; + return i; +} + +static int +load_dec(char *dir, char *file) +{ + char *line; + char *end; + long i; + + line = load_string(dir, file); + if (!line) + return -1; + i = strtol(line, &end, 10); + free(line); + if (end == line) + return -1; + return i; +} + +static int +dir_filter_tty_colon(const struct dirent *d) +{ + return strncmp(d->d_name, "tty:", 4) == 0; +} + +static int +dir_filter_tty(const struct dirent *d) +{ + return strncmp(d->d_name, "tty", 3) == 0; +} + +static char * +usb_tty(char *sys) +{ + char *base; + int num_configs; + int config; + struct dirent **namelist; + int interface; + int num_interfaces; + char endpoint_base[20]; + char *endpoint_full; + char *tty_dir; + int ntty; + char *tty; + + base = cc_basename(sys); + num_configs = load_hex(sys, "bNumConfigurations"); + num_interfaces = load_hex(sys, "bNumInterfaces"); + for (config = 1; config <= num_configs; config++) { + for (interface = 0; interface < num_interfaces; interface++) { + sprintf(endpoint_base, "%s:%d.%d", + base, config, interface); + endpoint_full = cc_fullname(sys, endpoint_base); + + /* Check for tty:ttyACMx style names + */ + ntty = scandir(endpoint_full, &namelist, + dir_filter_tty_colon, + alphasort); + if (ntty > 0) { + free(endpoint_full); + tty = cc_fullname("/dev", namelist[0]->d_name + 4); + free(namelist); + return tty; + } + + /* Check for tty/ttyACMx style names + */ + tty_dir = cc_fullname(endpoint_full, "tty"); + free(endpoint_full); + ntty = scandir(tty_dir, &namelist, + dir_filter_tty, + alphasort); + free (tty_dir); + if (ntty > 0) { + tty = cc_fullname("/dev", namelist[0]->d_name); + free(namelist); + return tty; + } + } + } + return NULL; +} + +static struct cc_usbdev * +usb_scan_device(char *sys) +{ + struct cc_usbdev *usbdev; + + usbdev = calloc(1, sizeof (struct cc_usbdev)); + if (!usbdev) + return NULL; + usbdev->sys = strdup(sys); + usbdev->manufacturer = load_string(sys, "manufacturer"); + usbdev->product = load_string(sys, "product"); + usbdev->serial = load_dec(sys, "serial"); + usbdev->idProduct = load_hex(sys, "idProduct"); + usbdev->idVendor = load_hex(sys, "idVendor"); + usbdev->tty = usb_tty(sys); + return usbdev; +} + +static void +usbdev_free(struct cc_usbdev *usbdev) +{ + free(usbdev->sys); + free(usbdev->manufacturer); + free(usbdev->product); + /* this can get used as a return value */ + if (usbdev->tty) + free(usbdev->tty); + free(usbdev); +} + +#define USB_DEVICES "/sys/bus/usb/devices" + +static int +dir_filter_dev(const struct dirent *d) +{ + const char *n = d->d_name; + char c; + + while ((c = *n++)) { + if (isdigit(c)) + continue; + if (c == '-') + continue; + if (c == '.' && n != d->d_name + 1) + continue; + return 0; + } + return 1; +} + +struct cc_usbdevs * +cc_usbdevs_scan(void) +{ + int e; + struct dirent **ents; + char *dir; + struct cc_usbdev *dev; + struct cc_usbdevs *devs; + int n; + + devs = calloc(1, sizeof (struct cc_usbdevs)); + if (!devs) + return NULL; + + n = scandir (USB_DEVICES, &ents, + dir_filter_dev, + alphasort); + if (!n) + return 0; + for (e = 0; e < n; e++) { + dir = cc_fullname(USB_DEVICES, ents[e]->d_name); + dev = usb_scan_device(dir); + free(dir); + if (dev->idVendor == 0xfffe && dev->tty) { + if (devs->dev) + devs->dev = realloc(devs->dev, + (devs->ndev + 1) * sizeof (struct usbdev *)); + else + devs->dev = malloc (sizeof (struct usbdev *)); + devs->dev[devs->ndev++] = dev; + } + } + free(ents); + return devs; +} + +void +cc_usbdevs_free(struct cc_usbdevs *usbdevs) +{ + int i; + + if (!usbdevs) + return; + for (i = 0; i < usbdevs->ndev; i++) + usbdev_free(usbdevs->dev[i]); + free(usbdevs); +} + +static char * +match_dev(char *product, int serial) +{ + struct cc_usbdevs *devs; + struct cc_usbdev *dev; + int i; + char *tty = NULL; + + devs = cc_usbdevs_scan(); + if (!devs) + return NULL; + for (i = 0; i < devs->ndev; i++) { + dev = devs->dev[i]; + if (product && strncmp (product, dev->product, strlen(product)) != 0) + continue; + if (serial && serial != dev->serial) + continue; + break; + } + if (i < devs->ndev) { + tty = devs->dev[i]->tty; + devs->dev[i]->tty = NULL; + } + cc_usbdevs_free(devs); + return tty; +} + +char * +cc_usbdevs_find_by_arg(char *arg, char *default_product) +{ + char *product; + int serial; + char *end; + char *colon; + char *tty; + + if (arg) + { + /* check for <serial> */ + serial = strtol(arg, &end, 0); + if (end != arg) { + if (*end != '\0') + return NULL; + product = NULL; + } else { + /* check for <product>:<serial> */ + colon = strchr(arg, ':'); + if (colon) { + product = strndup(arg, colon - arg); + serial = strtol(colon + 1, &end, 0); + if (*end != '\0') + return NULL; + } else { + product = arg; + serial = 0; + } + } + } else { + product = NULL; + serial = 0; + } + tty = NULL; + if (!product && default_product) + tty = match_dev(default_product, serial); + if (!tty) + tty = match_dev(product, serial); + if (product && product != arg) + free(product); + return tty; +} diff --git a/ao-tools/lib/cc-util.c b/ao-tools/lib/cc-util.c new file mode 100644 index 00000000..65488ee9 --- /dev/null +++ b/ao-tools/lib/cc-util.c @@ -0,0 +1,80 @@ +/* + * Copyright © 2009 Keith Packard <keithp@keithp.com> + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; version 2 of the License. + * + * This program is distributed in the hope that it will be useful, but + * WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * General Public License for more details. + * + * You should have received a copy of the GNU General Public License along + * with this program; if not, write to the Free Software Foundation, Inc., + * 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA. + */ + +#define _GNU_SOURCE +#include "cc.h" +#include <string.h> +#include <stdlib.h> +#include <unistd.h> +#include <sys/stat.h> +#include <sys/types.h> +#include <errno.h> + +char * +cc_fullname (char *dir, char *file) +{ + char *new; + int dlen = strlen (dir); + int flen = strlen (file); + int slen = 0; + + if (dir[dlen-1] != '/') + slen = 1; + new = malloc (dlen + slen + flen + 1); + if (!new) + return 0; + strcpy(new, dir); + if (slen) + strcat (new, "/"); + strcat(new, file); + return new; +} + +char * +cc_basename(char *file) +{ + char *b; + + b = strrchr(file, '/'); + if (!b) + return file; + return b + 1; +} + +int +cc_mkdir(char *dir) +{ + char *slash; + char *d; + char *part; + + d = dir; + for (;;) { + slash = strchr (d, '/'); + if (!slash) + slash = d + strlen(d); + if (!*slash) + break; + part = strndup(dir, slash - dir); + if (!access(part, F_OK)) + if (mkdir(part, 0777) < 0) + return -errno; + free(part); + d = slash + 1; + } + return 0; +} diff --git a/ao-tools/lib/cc.h b/ao-tools/lib/cc.h new file mode 100644 index 00000000..6257ee44 --- /dev/null +++ b/ao-tools/lib/cc.h @@ -0,0 +1,373 @@ +/* + * Copyright © 2009 Keith Packard <keithp@keithp.com> + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; version 2 of the License. + * + * This program is distributed in the hope that it will be useful, but + * WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * General Public License for more details. + * + * You should have received a copy of the GNU General Public License along + * with this program; if not, write to the Free Software Foundation, Inc., + * 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA. + */ + +#ifndef _CC_H_ +#define _CC_H_ + +#include <stdio.h> +#include <stdint.h> +#include "cc-telemetry.h" + +char * +cc_fullname (char *dir, char *file); + +char * +cc_basename(char *file); + +int +cc_mkdir(char *dir); + +struct cc_usbdev { + char *sys; + char *tty; + char *manufacturer; + char *product; + int serial; /* AltOS always uses simple integer serial numbers */ + int idProduct; + int idVendor; +}; + +struct cc_usbdevs { + struct cc_usbdev **dev; + int ndev; +}; + +void +cc_usbdevs_free(struct cc_usbdevs *usbdevs); + +struct cc_usbdevs * +cc_usbdevs_scan(void); + +char * +cc_usbdevs_find_by_arg(char *arg, char *default_product); + +void +cc_set_log_dir(char *dir); + +char * +cc_get_log_dir(void); + +char * +cc_make_filename(int serial, int flight, char *ext); + +/* + * For sequential data which are not evenly spaced + */ + +struct cc_timedataelt { + double time; + double value; +}; + +struct cc_timedata { + int num; + int size; + struct cc_timedataelt *data; + double time_offset; +}; + + +/* + * For GPS data + */ + +struct cc_gpselt { + double time; + int hour; + int minute; + int second; + int flags; + double lat; + double lon; + double alt; +}; + +#define SIRF_SAT_STATE_ACQUIRED (1 << 0) +#define SIRF_SAT_STATE_CARRIER_PHASE_VALID (1 << 1) +#define SIRF_SAT_BIT_SYNC_COMPLETE (1 << 2) +#define SIRF_SAT_SUBFRAME_SYNC_COMPLETE (1 << 3) +#define SIRF_SAT_CARRIER_PULLIN_COMPLETE (1 << 4) +#define SIRF_SAT_CODE_LOCKED (1 << 5) +#define SIRF_SAT_ACQUISITION_FAILED (1 << 6) +#define SIRF_SAT_EPHEMERIS_AVAILABLE (1 << 7) + +struct cc_gpssat { + double time; + uint16_t svid; + uint8_t c_n; +}; + +struct cc_gpssats { + int nsat; + struct cc_gpssat sat[12]; +}; + +struct cc_gpsdata { + int num; + int size; + struct cc_gpselt *data; + double time_offset; + int numsats; + int sizesats; + struct cc_gpssats *sats; +}; + +/* + * For sequential data which are evenly spaced + */ +struct cc_perioddata { + int num; + double start; + double step; + double *data; +}; + +enum ao_flight_state { + ao_flight_startup = 0, + ao_flight_idle = 1, + ao_flight_pad = 2, + ao_flight_boost = 3, + ao_flight_fast = 4, + ao_flight_coast = 5, + ao_flight_drogue = 6, + ao_flight_main = 7, + ao_flight_landed = 8, + ao_flight_invalid = 9 +}; + +struct cc_flightraw { + int flight; + int serial; + double ground_accel; + double ground_pres; + int year, month, day; + struct cc_timedata accel; + struct cc_timedata pres; + struct cc_timedata temp; + struct cc_timedata volt; + struct cc_timedata main; + struct cc_timedata drogue; + struct cc_timedata state; + struct cc_gpsdata gps; +}; + +struct cc_flightraw * +cc_log_read(FILE *file); + +void +cc_flightraw_free(struct cc_flightraw *raw); + +struct cc_flightcooked { + double flight_start; + double flight_stop; + + struct cc_perioddata accel_accel; + struct cc_perioddata accel_speed; + struct cc_perioddata accel_pos; + struct cc_perioddata pres_pos; + struct cc_perioddata pres_speed; + struct cc_perioddata pres_accel; + struct cc_perioddata gps_lat; + struct cc_perioddata gps_lon; + struct cc_perioddata gps_alt; + + /* unfiltered, but converted */ + struct cc_timedata pres; + struct cc_timedata accel; + struct cc_timedata state; +}; + +/* + * Telemetry data contents + */ + + +struct cc_gps_time { + int year; + int month; + int day; + int hour; + int minute; + int second; +}; + +struct cc_gps { + int nsat; + int gps_locked; + int gps_connected; + struct cc_gps_time gps_time; + double lat; /* degrees (+N -S) */ + double lon; /* degrees (+E -W) */ + int alt; /* m */ + + int gps_extended; /* has extra data */ + double ground_speed; /* m/s */ + int course; /* degrees */ + double climb_rate; /* m/s */ + double hdop; /* unitless? */ + int h_error; /* m */ + int v_error; /* m */ +}; + +#define SIRF_SAT_STATE_ACQUIRED (1 << 0) +#define SIRF_SAT_STATE_CARRIER_PHASE_VALID (1 << 1) +#define SIRF_SAT_BIT_SYNC_COMPLETE (1 << 2) +#define SIRF_SAT_SUBFRAME_SYNC_COMPLETE (1 << 3) +#define SIRF_SAT_CARRIER_PULLIN_COMPLETE (1 << 4) +#define SIRF_SAT_CODE_LOCKED (1 << 5) +#define SIRF_SAT_ACQUISITION_FAILED (1 << 6) +#define SIRF_SAT_EPHEMERIS_AVAILABLE (1 << 7) + +struct cc_gps_sat { + int svid; + int c_n0; +}; + +struct cc_gps_tracking { + int channels; + struct cc_gps_sat sats[12]; +}; + +struct cc_telem { + char callsign[16]; + int serial; + int flight; + int rssi; + char state[16]; + int tick; + int accel; + int pres; + int temp; + int batt; + int drogue; + int main; + int flight_accel; + int ground_accel; + int flight_vel; + int flight_pres; + int ground_pres; + int accel_plus_g; + int accel_minus_g; + struct cc_gps gps; + struct cc_gps_tracking gps_tracking; +}; + +int +cc_telem_parse(const char *input_line, struct cc_telem *telem); + +#ifndef TRUE +#define TRUE 1 +#define FALSE 0 +#endif + +/* Conversion functions */ +double +cc_pressure_to_altitude(double pressure); + +double +cc_altitude_to_pressure(double altitude); + +double +cc_barometer_to_pressure(double baro); + +double +cc_barometer_to_altitude(double baro); + +double +cc_accelerometer_to_acceleration(double accel, double ground_accel); + +double +cc_thermometer_to_temperature(double thermo); + +double +cc_battery_to_voltage(double battery); + +double +cc_ignitor_to_voltage(double ignite); + +void +cc_great_circle (double start_lat, double start_lon, + double end_lat, double end_lon, + double *dist, double *bearing); + +void +cc_timedata_limits(struct cc_timedata *d, double min_time, double max_time, int *start, int *stop); + +int +cc_timedata_min(struct cc_timedata *d, double min_time, double max_time); + +int +cc_timedata_min_mag(struct cc_timedata *d, double min_time, double max_time); + +int +cc_timedata_max(struct cc_timedata *d, double min_time, double max_time); + +int +cc_timedata_max_mag(struct cc_timedata *d, double min_time, double max_time); + +double +cc_timedata_average(struct cc_timedata *d, double min_time, double max_time); + +double +cc_timedata_average_mag(struct cc_timedata *d, double min_time, double max_time); + +int +cc_perioddata_limits(struct cc_perioddata *d, double min_time, double max_time, int *start, int *stop); + +int +cc_perioddata_min(struct cc_perioddata *d, double min_time, double max_time); + +int +cc_perioddata_min_mag(struct cc_perioddata *d, double min_time, double max_time); + +int +cc_perioddata_max(struct cc_perioddata *d, double min_time, double max_time); + +int +cc_perioddata_max_mag(struct cc_perioddata *d, double min_time, double max_time); + +double +cc_perioddata_average(struct cc_perioddata *d, double min_time, double max_time); + +double +cc_perioddata_average_mag(struct cc_perioddata *d, double min_time, double max_time); + +double * +cc_low_pass(double *data, int data_len, double omega_pass, double omega_stop, double error); + +struct cc_perioddata * +cc_period_make(struct cc_timedata *td, double start_time, double stop_time); + +struct cc_perioddata * +cc_period_low_pass(struct cc_perioddata *raw, double omega_pass, double omega_stop, double error); + +struct cc_timedata * +cc_timedata_convert(struct cc_timedata *d, double (*f)(double v, double a), double a); + +struct cc_timedata * +cc_timedata_integrate(struct cc_timedata *d, double min_time, double max_time); + +struct cc_perioddata * +cc_perioddata_differentiate(struct cc_perioddata *i); + +struct cc_flightcooked * +cc_flight_cook(struct cc_flightraw *raw); + +void +cc_flightcooked_free(struct cc_flightcooked *cooked); + +#endif /* _CC_H_ */ diff --git a/ao-tools/lib/cephes.h b/ao-tools/lib/cephes.h new file mode 100644 index 00000000..f8ec264f --- /dev/null +++ b/ao-tools/lib/cephes.h @@ -0,0 +1,122 @@ +/* + * This file comes from the cephes math library, which was + * released under the GPLV2+ license as a part of the Debian labplot + * package (I've included the GPLV2 license reference here to make + * this clear) - Keith Packard <keithp@keithp.com> + * + * Cephes Math Library Release 2.0: April, 1987 + * Copyright 1984, 1987 by Stephen L. Moshier + * Direct inquiries to 30 Frost Street, Cambridge, MA 02140 + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; version 2 of the License. + * + * This program is distributed in the hope that it will be useful, but + * WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * General Public License for more details. + * + * You should have received a copy of the GNU General Public License along + * with this program; if not, write to the Free Software Foundation, Inc., + * 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA. + */ +/* + * Prototypes of Cephes functions + */ + +#ifndef _CEPHES_H_ +#define _CEPHES_H_ + +/* Variable for error reporting. See mtherr.c. */ +extern int merror; + +#if 0 +extern int airy ( double x, double *ai, double *aip, double *bi, double *bip ); +extern double beta ( double a, double b ); +extern double lbeta ( double a, double b ); +extern double chdtrc ( double df, double x ); +extern double chdtr ( double df, double x ); +extern double chdtri ( double df, double y ); +extern double dawsn ( double xx ); +extern double ellie ( double phi, double m ); +extern double ellik ( double phi, double m ); +extern double ellpe ( double x ); +extern double ellpk ( double x ); +extern double expn ( int n, double x ); +extern double fac ( int i ); +extern double fdtrc ( int ia, int ib, double x ); +extern double fdtr ( int ia, int ib, double x ); +extern double fdtri ( int ia, int ib, double y ); +extern double frexp ( double x, int *pw2 ); +extern double ldexp ( double x, int pw2 ); +extern int fresnl ( double xxa, double *ssa, double *cca ); +extern double gdtr ( double a, double b, double x ); +extern double gdtrc ( double a, double b, double x ); +extern double hyp2f0 ( double a, double b, double x, int type, double *err ); +extern double hyp2f1 ( double a, double b, double c, double x ); +extern double hyperg ( double a, double b, double x ); +#endif +extern double i0 ( double x ); +extern double i0e ( double x ); +#if 0 +extern double i1 ( double x ); +extern double i1e ( double x ); +extern double iv ( double v, double x ); +extern double igamc ( double a, double x ); +extern double igam ( double a, double x ); +extern double igami ( double a, double y0_ ); +extern double incbet ( double aa, double bb, double xx ); +extern double incbi ( double aa, double bb, double yy0 ); +extern double jv ( double n, double x ); +extern double k0 ( double x ); +extern double k0e ( double x ); +extern double k1 ( double x ); +extern double k1e ( double x ); +extern double kn ( int nn, double x ); +extern int mtherr ( char *name, int code ); +extern double ndtr ( double a ); +extern double ndtri ( double y0_ ); +extern double pdtrc ( int k, double m ); +extern double pdtr ( int k, double m ); +extern double pdtri ( int k, double y ); +extern double psi ( double x ); +extern void revers ( double y[], double x[], int n ); +extern double true_gamma ( double x ); +extern double rgamma ( double x ); +extern int shichi ( double x, double *si, double *ci ); +extern int sici ( double x, double *si, double *ci ); +extern double spence ( double x ); +extern double stdtr ( int k, double t ); +extern double stdtri ( int k, double p ); +extern double onef2 ( double a, double b, double c, double x, double *err ); +extern double threef0 ( double a, double b, double c, double x, double *err ); +extern double struve ( double v, double x ); +extern double log1p ( double x ); +extern double expm1 ( double x ); +extern double cosm1 ( double x ); +extern double yv ( double v, double x ); +extern double zeta ( double x, double q ); +extern double zetac ( double x ); + +#endif +extern double chbevl ( double x, void *P, int n ); +#if 0 +extern double polevl ( double x, void *P, int n ); +extern double p1evl ( double x, void *P, int n ); + +/* polyn.c */ +extern void polini ( int maxdeg ); +extern void polprt ( double a[], int na, int d ); +extern void polclr ( double *a, int n ); +extern void polmov ( double *a, int na, double *b ); +extern void polmul ( double a[], int na, double b[], int nb, double c[] ); +extern void poladd ( double a[], int na, double b[], int nb, double c[] ); +extern void polsub ( double a[], int na, double b[], int nb, double c[] ); +extern int poldiv ( double a[], int na, double b[], int nb, double c[] ); +extern void polsbt ( double a[], int na, double b[], int nb, double c[] ); +extern double poleva ( double a[], int na, double x ); + +#endif + +#endif /* _CEPHES_H_ */ diff --git a/ao-tools/lib/chbevl.c b/ao-tools/lib/chbevl.c new file mode 100644 index 00000000..22892413 --- /dev/null +++ b/ao-tools/lib/chbevl.c @@ -0,0 +1,81 @@ +/* chbevl.c + * + * Evaluate Chebyshev series + * + * + * + * SYNOPSIS: + * + * int N; + * double x, y, coef[N], chebevl(); + * + * y = chbevl( x, coef, N ); + * + * + * + * DESCRIPTION: + * + * Evaluates the series + * + * N-1 + * - ' + * y = > coef[i] T (x/2) + * - i + * i=0 + * + * of Chebyshev polynomials Ti at argument x/2. + * + * Coefficients are stored in reverse order, i.e. the zero + * order term is last in the array. Note N is the number of + * coefficients, not the order. + * + * If coefficients are for the interval a to b, x must + * have been transformed to x -> 2(2x - b - a)/(b-a) before + * entering the routine. This maps x from (a, b) to (-1, 1), + * over which the Chebyshev polynomials are defined. + * + * If the coefficients are for the inverted interval, in + * which (a, b) is mapped to (1/b, 1/a), the transformation + * required is x -> 2(2ab/x - b - a)/(b-a). If b is infinity, + * this becomes x -> 4a/x - 1. + * + * + * + * SPEED: + * + * Taking advantage of the recurrence properties of the + * Chebyshev polynomials, the routine requires one more + * addition per loop than evaluating a nested polynomial of + * the same degree. + * + */ +/* chbevl.c */ + +/* +Cephes Math Library Release 2.0: April, 1987 +Copyright 1985, 1987 by Stephen L. Moshier +Direct inquiries to 30 Frost Street, Cambridge, MA 02140 +*/ + +#include "cephes.h" + +double chbevl(double x,void* array,int n ) +{ +double b0, b1, b2, *p; +int i; + +p = (double *) array; +b0 = *p++; +b1 = 0.0; +i = n - 1; + +do + { + b2 = b1; + b1 = b0; + b0 = x * b1 - b2 + *p++; + } +while( --i ); + +return( 0.5*(b0-b2) ); +} diff --git a/ao-tools/lib/cmath.h b/ao-tools/lib/cmath.h new file mode 100644 index 00000000..2751aecf --- /dev/null +++ b/ao-tools/lib/cmath.h @@ -0,0 +1,179 @@ +/* + * Grace - GRaphing, Advanced Computation and Exploration of data + * + * Home page: http://plasma-gate.weizmann.ac.il/Grace/ + * + * Copyright (c) 1991-1995 Paul J Turner, Portland, OR + * Copyright (c) 1996-2000 Grace Development Team + * + * Maintained by Evgeny Stambulchik <fnevgeny@plasma-gate.weizmann.ac.il> + * + * + * All Rights Reserved + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 2 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software + * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. + */ + +/* cmath.h - replacement for math.h or missing in libm functions */ + +#if defined(HAVE_MATH_H) +# include <math.h> +#endif +#if defined(HAVE_FLOAT_H) +# include <float.h> +#endif +#if defined(HAVE_IEEEFP_H) +# include <ieeefp.h> +#endif + +#ifndef __GRACE_SOURCE_ + +#ifndef MACHEP +extern double MACHEP; +#endif + +#ifndef UFLOWTHRESH +extern double UFLOWTHRESH; +#endif + +#ifndef MAXNUM +extern double MAXNUM; +#endif + +#endif /* __GRACE_SOURCE_ */ + +#ifndef M_PI +# define M_PI 3.14159265358979323846 +#endif + +#ifndef M_SQRT2 +# define M_SQRT2 1.41421356237309504880 /* sqrt(2) */ +#endif + +#ifndef M_SQRT1_2 +# define M_SQRT1_2 0.70710678118654752440 /* 1/sqrt(2) */ +#endif + +#ifndef M_SQRT1_3 +# define M_SQRT1_3 0.57735026918962576451 /* 1/sqrt(3) */ +#endif + +#ifndef HAVE_HYPOT +# define hypot(x, y) sqrt((x)*(x) + (y)*(y)) +#endif + +extern double round ( double x ); +#ifndef HAVE_RINT +# define rint round +#else +# ifndef HAVE_RINT_DECL +extern double rint ( double x ); +# endif +#endif + +#ifndef HAVE_CBRT_DECL +extern double cbrt ( double x ); +#endif + +/* Cygnus gnuwin32 has the log2 macro */ +#ifdef log2 +# undef log2 +#endif + +#ifndef HAVE_LOG2_DECL +extern double log2 ( double x ); +#endif + +#ifndef HAVE_LGAMMA +extern int sgngam; +# define lgamma lgam +# define signgam sgngam +extern double lgam ( double x ); +#else +# ifndef HAVE_LGAMMA_DECL +extern double lgamma ( double x ); +# endif +# ifndef HAVE_SIGNGAM_DECL +extern int signgam; +# endif +# define lgam lgamma +# define sgngam signgam +#endif + +#ifndef HAVE_ACOSH_DECL +extern double acosh ( double x ); +#endif + +#ifndef HAVE_ASINH_DECL +extern double asinh ( double x ); +#endif + +#ifndef HAVE_ATANH_DECL +extern double atanh ( double x ); +#endif + +#ifndef HAVE_ERF_DECL +extern double erf ( double x ); +#endif + +#ifndef HAVE_ERFC_DECL +extern double erfc ( double x ); +#endif + +#ifndef HAVE_Y0_DECL +extern double y0 ( double x ); +#endif +#ifndef HAVE_Y1_DECL +extern double y1 ( double x ); +#endif +#ifndef HAVE_YN_DECL +extern double yn ( int n, double x ); +#endif +#ifndef HAVE_J0_DECL +extern double j0 ( double x ); +#endif +#ifndef HAVE_J1_DECL +extern double j1 ( double x ); +#endif +#ifndef HAVE_JN_DECL +extern double jn ( int n, double x ); +#endif + +/* isfinite is a macro */ +#ifdef isfinite +# define HAVE_ISFINITE_MACRO +#endif + +#ifndef HAVE_FINITE +# define finite isfinite +# if !defined(HAVE_ISFINITE_DECL) && !defined(HAVE_ISFINITE_MACRO) +extern int isfinite ( double x ); +# endif +#else +# ifndef HAVE_FINITE_DECL +extern int finite ( double x ); +# endif +#endif + +#ifndef HAVE_ISNAN_DECL +#ifdef __FreeBSD__ +# include <sys/param.h> +# if __FreeBSD_version < 500100 +extern int isnan ( double x ); +# endif +#endif +#else +extern int isnan ( double x ); +#endif diff --git a/ao-tools/lib/i0.c b/ao-tools/lib/i0.c new file mode 100644 index 00000000..6f7b5a57 --- /dev/null +++ b/ao-tools/lib/i0.c @@ -0,0 +1,414 @@ +/* + * This file comes from the cephes math library, which was + * released under the GPLV2+ license as a part of the Debian labplot + * package (I've included the GPLV2 license reference here to make + * this clear) - Keith Packard <keithp@keithp.com> + * + * Cephes Math Library Release 2.0: April, 1987 + * Copyright 1984, 1987 by Stephen L. Moshier + * Direct inquiries to 30 Frost Street, Cambridge, MA 02140 + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; version 2 of the License. + * + * This program is distributed in the hope that it will be useful, but + * WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * General Public License for more details. + * + * You should have received a copy of the GNU General Public License along + * with this program; if not, write to the Free Software Foundation, Inc., + * 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA. + */ + +/* i0.c + * + * Modified Bessel function of order zero + * + * + * + * SYNOPSIS: + * + * double x, y, i0(); + * + * y = i0( x ); + * + * + * + * DESCRIPTION: + * + * Returns modified Bessel function of order zero of the + * argument. + * + * The function is defined as i0(x) = j0( ix ). + * + * The range is partitioned into the two intervals [0,8] and + * (8, infinity). Chebyshev polynomial expansions are employed + * in each interval. + * + * + * + * ACCURACY: + * + * Relative error: + * arithmetic domain # trials peak rms + * DEC 0,30 6000 8.2e-17 1.9e-17 + * IEEE 0,30 30000 5.8e-16 1.4e-16 + * + */ +/* i0e.c + * + * Modified Bessel function of order zero, + * exponentially scaled + * + * + * + * SYNOPSIS: + * + * double x, y, i0e(); + * + * y = i0e( x ); + * + * + * + * DESCRIPTION: + * + * Returns exponentially scaled modified Bessel function + * of order zero of the argument. + * + * The function is defined as i0e(x) = exp(-|x|) j0( ix ). + * + * + * + * ACCURACY: + * + * Relative error: + * arithmetic domain # trials peak rms + * IEEE 0,30 30000 5.4e-16 1.2e-16 + * See i0(). + * + */ + +/* i0.c */ + + +/* +Cephes Math Library Release 2.0: April, 1987 +Copyright 1984, 1987 by Stephen L. Moshier +Direct inquiries to 30 Frost Street, Cambridge, MA 02140 +*/ + +#include <math.h> +#include "mconf.h" +#include "cephes.h" + +/* Chebyshev coefficients for exp(-x) I0(x) + * in the interval [0,8]. + * + * lim(x->0){ exp(-x) I0(x) } = 1. + */ + +#ifdef UNK +static double A[] = +{ +-4.41534164647933937950E-18, + 3.33079451882223809783E-17, +-2.43127984654795469359E-16, + 1.71539128555513303061E-15, +-1.16853328779934516808E-14, + 7.67618549860493561688E-14, +-4.85644678311192946090E-13, + 2.95505266312963983461E-12, +-1.72682629144155570723E-11, + 9.67580903537323691224E-11, +-5.18979560163526290666E-10, + 2.65982372468238665035E-9, +-1.30002500998624804212E-8, + 6.04699502254191894932E-8, +-2.67079385394061173391E-7, + 1.11738753912010371815E-6, +-4.41673835845875056359E-6, + 1.64484480707288970893E-5, +-5.75419501008210370398E-5, + 1.88502885095841655729E-4, +-5.76375574538582365885E-4, + 1.63947561694133579842E-3, +-4.32430999505057594430E-3, + 1.05464603945949983183E-2, +-2.37374148058994688156E-2, + 4.93052842396707084878E-2, +-9.49010970480476444210E-2, + 1.71620901522208775349E-1, +-3.04682672343198398683E-1, + 6.76795274409476084995E-1 +}; +#endif + +#ifdef DEC +static unsigned short A[] = { +0121642,0162671,0004646,0103567, +0022431,0115424,0135755,0026104, +0123214,0023533,0110365,0156635, +0023767,0033304,0117662,0172716, +0124522,0100426,0012277,0157531, +0025254,0155062,0054461,0030465, +0126010,0131143,0013560,0153604, +0026517,0170577,0006336,0114437, +0127227,0162253,0152243,0052734, +0027724,0142766,0061641,0160200, +0130416,0123760,0116564,0125262, +0031066,0144035,0021246,0054641, +0131537,0053664,0060131,0102530, +0032201,0155664,0165153,0020652, +0132617,0061434,0074423,0176145, +0033225,0174444,0136147,0122542, +0133624,0031576,0056453,0020470, +0034211,0175305,0172321,0041314, +0134561,0054462,0147040,0165315, +0035105,0124333,0120203,0162532, +0135427,0013750,0174257,0055221, +0035726,0161654,0050220,0100162, +0136215,0131361,0000325,0041110, +0036454,0145417,0117357,0017352, +0136702,0072367,0104415,0133574, +0037111,0172126,0072505,0014544, +0137302,0055601,0120550,0033523, +0037457,0136543,0136544,0043002, +0137633,0177536,0001276,0066150, +0040055,0041164,0100655,0010521 +}; +#endif + +#ifdef IBMPC +static unsigned short A[] = { +0xd0ef,0x2134,0x5cb7,0xbc54, +0xa589,0x977d,0x3362,0x3c83, +0xbbb4,0x721e,0x84eb,0xbcb1, +0x5eba,0x93f6,0xe6d8,0x3cde, +0xfbeb,0xc297,0x5022,0xbd0a, +0x2627,0x4b26,0x9b46,0x3d35, +0x1af0,0x62ee,0x164c,0xbd61, +0xd324,0xe19b,0xfe2f,0x3d89, +0x6abc,0x7a94,0xfc95,0xbdb2, +0x3c10,0xcc74,0x98be,0x3dda, +0x9556,0x13ae,0xd4fe,0xbe01, +0xcb34,0xa454,0xd903,0x3e26, +0x30ab,0x8c0b,0xeaf6,0xbe4b, +0x6435,0x9d4d,0x3b76,0x3e70, +0x7f8d,0x8f22,0xec63,0xbe91, +0xf4ac,0x978c,0xbf24,0x3eb2, +0x6427,0xcba5,0x866f,0xbed2, +0x2859,0xbe9a,0x3f58,0x3ef1, +0x1d5a,0x59c4,0x2b26,0xbf0e, +0x7cab,0x7410,0xb51b,0x3f28, +0xeb52,0x1f15,0xe2fd,0xbf42, +0x100e,0x8a12,0xdc75,0x3f5a, +0xa849,0x201a,0xb65e,0xbf71, +0xe3dd,0xf3dd,0x9961,0x3f85, +0xb6f0,0xf121,0x4e9e,0xbf98, +0xa32d,0xcea8,0x3e8a,0x3fa9, +0x06ea,0x342d,0x4b70,0xbfb8, +0x88c0,0x77ac,0xf7ac,0x3fc5, +0xcd8d,0xc057,0x7feb,0xbfd3, +0xa22a,0x9035,0xa84e,0x3fe5, +}; +#endif + +#ifdef MIEEE +static unsigned short A[] = { +0xbc54,0x5cb7,0x2134,0xd0ef, +0x3c83,0x3362,0x977d,0xa589, +0xbcb1,0x84eb,0x721e,0xbbb4, +0x3cde,0xe6d8,0x93f6,0x5eba, +0xbd0a,0x5022,0xc297,0xfbeb, +0x3d35,0x9b46,0x4b26,0x2627, +0xbd61,0x164c,0x62ee,0x1af0, +0x3d89,0xfe2f,0xe19b,0xd324, +0xbdb2,0xfc95,0x7a94,0x6abc, +0x3dda,0x98be,0xcc74,0x3c10, +0xbe01,0xd4fe,0x13ae,0x9556, +0x3e26,0xd903,0xa454,0xcb34, +0xbe4b,0xeaf6,0x8c0b,0x30ab, +0x3e70,0x3b76,0x9d4d,0x6435, +0xbe91,0xec63,0x8f22,0x7f8d, +0x3eb2,0xbf24,0x978c,0xf4ac, +0xbed2,0x866f,0xcba5,0x6427, +0x3ef1,0x3f58,0xbe9a,0x2859, +0xbf0e,0x2b26,0x59c4,0x1d5a, +0x3f28,0xb51b,0x7410,0x7cab, +0xbf42,0xe2fd,0x1f15,0xeb52, +0x3f5a,0xdc75,0x8a12,0x100e, +0xbf71,0xb65e,0x201a,0xa849, +0x3f85,0x9961,0xf3dd,0xe3dd, +0xbf98,0x4e9e,0xf121,0xb6f0, +0x3fa9,0x3e8a,0xcea8,0xa32d, +0xbfb8,0x4b70,0x342d,0x06ea, +0x3fc5,0xf7ac,0x77ac,0x88c0, +0xbfd3,0x7feb,0xc057,0xcd8d, +0x3fe5,0xa84e,0x9035,0xa22a +}; +#endif + + +/* Chebyshev coefficients for exp(-x) sqrt(x) I0(x) + * in the inverted interval [8,infinity]. + * + * lim(x->inf){ exp(-x) sqrt(x) I0(x) } = 1/sqrt(2pi). + */ + +#ifdef UNK +static double B[] = +{ +-7.23318048787475395456E-18, +-4.83050448594418207126E-18, + 4.46562142029675999901E-17, + 3.46122286769746109310E-17, +-2.82762398051658348494E-16, +-3.42548561967721913462E-16, + 1.77256013305652638360E-15, + 3.81168066935262242075E-15, +-9.55484669882830764870E-15, +-4.15056934728722208663E-14, + 1.54008621752140982691E-14, + 3.85277838274214270114E-13, + 7.18012445138366623367E-13, +-1.79417853150680611778E-12, +-1.32158118404477131188E-11, +-3.14991652796324136454E-11, + 1.18891471078464383424E-11, + 4.94060238822496958910E-10, + 3.39623202570838634515E-9, + 2.26666899049817806459E-8, + 2.04891858946906374183E-7, + 2.89137052083475648297E-6, + 6.88975834691682398426E-5, + 3.36911647825569408990E-3, + 8.04490411014108831608E-1 +}; +#endif + +#ifdef DEC +static unsigned short B[] = { +0122005,0066672,0123124,0054311, +0121662,0033323,0030214,0104602, +0022515,0170300,0113314,0020413, +0022437,0117350,0035402,0007146, +0123243,0000135,0057220,0177435, +0123305,0073476,0144106,0170702, +0023777,0071755,0017527,0154373, +0024211,0052214,0102247,0033270, +0124454,0017763,0171453,0012322, +0125072,0166316,0075505,0154616, +0024612,0133770,0065376,0025045, +0025730,0162143,0056036,0001632, +0026112,0015077,0150464,0063542, +0126374,0101030,0014274,0065457, +0127150,0077271,0125763,0157617, +0127412,0104350,0040713,0120445, +0027121,0023765,0057500,0001165, +0030407,0147146,0003643,0075644, +0031151,0061445,0044422,0156065, +0031702,0132224,0003266,0125551, +0032534,0000076,0147153,0005555, +0033502,0004536,0004016,0026055, +0034620,0076433,0142314,0171215, +0036134,0146145,0013454,0101104, +0040115,0171425,0062500,0047133 +}; +#endif + +#ifdef IBMPC +static unsigned short B[] = { +0x8b19,0x54ca,0xadb7,0xbc60, +0x9130,0x6611,0x46da,0xbc56, +0x8421,0x12d9,0xbe18,0x3c89, +0x41cd,0x0760,0xf3dd,0x3c83, +0x1fe4,0xabd2,0x600b,0xbcb4, +0xde38,0xd908,0xaee7,0xbcb8, +0xfb1f,0xa3ea,0xee7d,0x3cdf, +0xe6d7,0x9094,0x2a91,0x3cf1, +0x629a,0x7e65,0x83fe,0xbd05, +0xbb32,0xcf68,0x5d99,0xbd27, +0xc545,0x0d5f,0x56ff,0x3d11, +0xc073,0x6b83,0x1c8c,0x3d5b, +0x8cec,0xfa26,0x4347,0x3d69, +0x8d66,0x0317,0x9043,0xbd7f, +0x7bf2,0x357e,0x0fd7,0xbdad, +0x7425,0x0839,0x511d,0xbdc1, +0x004f,0xabe8,0x24fe,0x3daa, +0x6f75,0xc0f4,0xf9cc,0x3e00, +0x5b87,0xa922,0x2c64,0x3e2d, +0xd56d,0x80d6,0x5692,0x3e58, +0x616e,0xd9cd,0x8007,0x3e8b, +0xc586,0xc101,0x412b,0x3ec8, +0x9e52,0x7899,0x0fa3,0x3f12, +0x9049,0xa2e5,0x998c,0x3f6b, +0x09cb,0xaca8,0xbe62,0x3fe9 +}; +#endif + +#ifdef MIEEE +static unsigned short B[] = { +0xbc60,0xadb7,0x54ca,0x8b19, +0xbc56,0x46da,0x6611,0x9130, +0x3c89,0xbe18,0x12d9,0x8421, +0x3c83,0xf3dd,0x0760,0x41cd, +0xbcb4,0x600b,0xabd2,0x1fe4, +0xbcb8,0xaee7,0xd908,0xde38, +0x3cdf,0xee7d,0xa3ea,0xfb1f, +0x3cf1,0x2a91,0x9094,0xe6d7, +0xbd05,0x83fe,0x7e65,0x629a, +0xbd27,0x5d99,0xcf68,0xbb32, +0x3d11,0x56ff,0x0d5f,0xc545, +0x3d5b,0x1c8c,0x6b83,0xc073, +0x3d69,0x4347,0xfa26,0x8cec, +0xbd7f,0x9043,0x0317,0x8d66, +0xbdad,0x0fd7,0x357e,0x7bf2, +0xbdc1,0x511d,0x0839,0x7425, +0x3daa,0x24fe,0xabe8,0x004f, +0x3e00,0xf9cc,0xc0f4,0x6f75, +0x3e2d,0x2c64,0xa922,0x5b87, +0x3e58,0x5692,0x80d6,0xd56d, +0x3e8b,0x8007,0xd9cd,0x616e, +0x3ec8,0x412b,0xc101,0xc586, +0x3f12,0x0fa3,0x7899,0x9e52, +0x3f6b,0x998c,0xa2e5,0x9049, +0x3fe9,0xbe62,0xaca8,0x09cb +}; +#endif + +double i0(double x) +{ +double y; + +if( x < 0 ) + x = -x; +if( x <= 8.0 ) + { + y = (x/2.0) - 2.0; + return( exp(x) * chbevl( y, A, 30 ) ); + } + +return( exp(x) * chbevl( 32.0/x - 2.0, B, 25 ) / sqrt(x) ); + +} + + + + +double i0e(double x ) +{ +double y; + +if( x < 0 ) + x = -x; +if( x <= 8.0 ) + { + y = (x/2.0) - 2.0; + return( chbevl( y, A, 30 ) ); + } + +return( chbevl( 32.0/x - 2.0, B, 25 ) / sqrt(x) ); + +} diff --git a/ao-tools/lib/mconf.h b/ao-tools/lib/mconf.h new file mode 100644 index 00000000..af1ebb51 --- /dev/null +++ b/ao-tools/lib/mconf.h @@ -0,0 +1,211 @@ +/* + * This file comes from the cephes math library, which was + * released under the GPLV2+ license as a part of the Debian labplot + * package (I've included the GPLV2 license reference here to make + * this clear) - Keith Packard <keithp@keithp.com> + * + * Cephes Math Library Release 2.0: April, 1987 + * Copyright 1984, 1987 by Stephen L. Moshier + * Direct inquiries to 30 Frost Street, Cambridge, MA 02140 + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; version 2 of the License. + * + * This program is distributed in the hope that it will be useful, but + * WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * General Public License for more details. + * + * You should have received a copy of the GNU General Public License along + * with this program; if not, write to the Free Software Foundation, Inc., + * 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA. + */ +/* mconf.h + * + * Common include file for math routines + * + * + * + * SYNOPSIS: + * + * #include "mconf.h" + * + * + * + * DESCRIPTION: + * + * This file contains definitions for error codes that are + * passed to the common error handling routine mtherr() + * (which see). + * + * The file also includes a conditional assembly definition + * for the type of computer arithmetic (IEEE, DEC, Motorola + * IEEE, or UNKnown). + * + * For Digital Equipment PDP-11 and VAX computers, certain + * IBM systems, and others that use numbers with a 56-bit + * significand, the symbol DEC should be defined. In this + * mode, most floating point constants are given as arrays + * of octal integers to eliminate decimal to binary conversion + * errors that might be introduced by the compiler. + * + * For little-endian computers, such as IBM PC, that follow the + * IEEE Standard for Binary Floating Point Arithmetic (ANSI/IEEE + * Std 754-1985), the symbol IBMPC should be defined. These + * numbers have 53-bit significands. In this mode, constants + * are provided as arrays of hexadecimal 16 bit integers. + * + * Big-endian IEEE format is denoted MIEEE. On some RISC + * systems such as Sun SPARC, double precision constants + * must be stored on 8-byte address boundaries. Since integer + * arrays may be aligned differently, the MIEEE configuration + * may fail on such machines. + * + * To accommodate other types of computer arithmetic, all + * constants are also provided in a normal decimal radix + * which one can hope are correctly converted to a suitable + * format by the available C language compiler. To invoke + * this mode, define the symbol UNK. + * + * An important difference among these modes is a predefined + * set of machine arithmetic constants for each. The numbers + * MACHEP (the machine roundoff error), MAXNUM (largest number + * represented), and several other parameters are preset by + * the configuration symbol. Check the file const.c to + * ensure that these values are correct for your computer. + * + * Configurations NANS, INFINITIES, MINUSZERO, and DENORMAL + * may fail on many systems. Verify that they are supposed + * to work on your computer. + */ + +/* +Cephes Math Library Release 2.3: June, 1995 +Copyright 1984, 1987, 1989, 1995 by Stephen L. Moshier + +Adjusted for use with ACE/gr by Evgeny Stambulchik, October 1997 +*/ + +#define __GRACE_SOURCE_ + +#include "cmath.h" + +/* Type of computer arithmetic */ +/* In ACE/gr, defined as a compiler directive - no need to define here */ + +/* PDP-11, Pro350, VAX: + */ +#if defined(HAVE_DEC_FPU) +# define DEC 1 +#endif + +/* Intel IEEE, low order words come first: + */ +#if defined(HAVE_LIEEE_FPU) +# define IBMPC 1 +#endif + +/* Motorola IEEE, high order words come first + * (Sun 680x0 workstation): + */ +#if defined(HAVE_BIEEE_FPU) +# define MIEEE 1 +#endif + +/* UNKnown arithmetic, invokes coefficients given in + * normal decimal format. Beware of range boundary + * problems (MACHEP, MAXLOG, etc. in const.c) and + * roundoff problems in pow.c: + * (Sun SPARCstation) + */ + +#if (!defined (DEC) && !defined (IBMPC) && !defined (MIEEE)) +# define UNK 1 +#endif + +/* Define this `volatile' if your compiler thinks + * that floating point arithmetic obeys the associative + * and distributive laws. It will defeat some optimizations + * (but probably not enough of them). + * + * #define VOLATILE volatile + */ + +#ifndef VOLATILE +# define VOLATILE +#endif + +#ifdef PI +# undef PI +#endif + +#ifdef NAN +# undef NAN +#endif + +#ifdef INFINITY +# undef INFINITY +#endif + +/* Constant definitions for math error conditions + */ + +#if defined(DOMAIN) +# undef DOMAIN +#endif +#define DOMAIN 1 /* argument domain error */ + +#if defined(SING) +# undef SING +#endif +#define SING 2 /* argument singularity */ + +#if defined(OVERFLOW) +# undef OVERFLOW +#endif +#define OVERFLOW 3 /* overflow range error */ + +#if defined(UNDERFLOW) +# undef UNDERFLOW +#endif +#define UNDERFLOW 4 /* underflow range error */ + +#if defined(TLOSS) +# undef TLOSS +#endif +#define TLOSS 5 /* total loss of precision */ + +#if defined(PLOSS) +# undef PLOSS +#endif +#define PLOSS 6 /* partial loss of precision */ + +#if defined(EDOM) +# undef EDOM +#endif +#define EDOM 33 + +#if defined(ERANGE) +# undef ERANGE +#endif +#define ERANGE 34 + +#if !defined (UNK) + /* Define to support tiny denormal numbers, else undefine. */ +# define DENORMAL 1 + + /* Define to ask for infinity support, else undefine. */ +# define INFINITIES 1 + + /* Define to ask for support of numbers that are Not-a-Number, + else undefine. This may automatically define INFINITIES in some files. */ +# define NANS 1 + + /* Define to distinguish between -0.0 and +0.0. */ +# define MINUSZERO 1 +#endif + +/* Define 1 for ANSI C atan2() function + See atan.c and clog.c. */ +#define ANSIC 1 |