gpx.c 21 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588
  1. #ifdef PC_BUILD
  2. #include "../gps-test-tool/main.h"
  3. #else
  4. #include "main.h"
  5. #endif
  6. #define KALMAN_Q 8.5e-6
  7. #define KALMAN_R 4e-5
  8. #define KALMAN_ERR_MAX 6e-4
  9. __flash const char xml_header[] = "<?xml version=\"1.0\" encoding=\"UTF-8\"?>\n"
  10. "<gpx xmlns:xsi=\"http://www.w3.org/2001/XMLSchema-instance\" xmlns=\"http://www.topografix.com/GPX/1/1\" xsi:schemaLocation=\"http://www.topografix.com/GPX/1/1 http://www.topografix.com/GPX/1/1/gpx.xsd\" version=\"1.1\" creator=\"k4be\">\n";
  11. __flash const char xml_trk_start[] = "\t<trk>\n";
  12. __flash const char xml_trkseg_end[] = "\t\t</trkseg>\n";
  13. __flash const char xml_trkseg_start[] = "\t\t<trkseg>\n";
  14. FIL gpx_file;
  15. static char buf[sizeof(xml_header)+sizeof(xml_trk_start)+2];
  16. struct kalman_s {
  17. unsigned char initialized;
  18. float x_est_last;
  19. float P_last;
  20. float Q;
  21. float R;
  22. float K;
  23. };
  24. #define PREV_POINTS_LENGTH 4
  25. #define AVG_COUNT 3
  26. #define MIN_DIST_DELTA 2.0
  27. /* Elevation dead-band + smoothing, applied to the filtered ascent/descent totals */
  28. #define ELEV_SMOOTH_TAU 30.0 /* seconds */
  29. #define ELEV_DEADBAND 5.0 /* meters */
  30. /* Large-jump detection thresholds. POS_JUMP_THRESHOLD is informational
  31. * logging only - the distance-diff spike check above already rejects
  32. * inconsistent position jumps regardless of size. ALT_JUMP_THRESHOLD both
  33. * logs and (combined with the same spike shape) rejects an inconsistent
  34. * altitude jump - unlike position there was previously no altitude
  35. * consistency check at all, letting a poor-geometry fix's bogus altitude
  36. * (right sats/HDOP, wrong solution - e.g. just after reacquiring signal)
  37. * straight into the accepted track. */
  38. #define POS_JUMP_THRESHOLD 100.0 /* meters between consecutive raw fixes */
  39. #define ALT_JUMP_THRESHOLD 50.0 /* meters between consecutive raw fixes */
  40. /* Post-gap settling: a freshly reacquired fix's position/altitude solution
  41. * can keep drifting for several seconds after reacquisition even while
  42. * reported quality (satellite count, HDOP) already looks good throughout -
  43. * seen on a real recorded track exiting an underground space: HDOP was
  44. * already 0.81/17 sats just 7s after reacquiring, yet elevation kept
  45. * changing by several meters per sample for another ~10s after that before
  46. * settling to normal (sub-meter) noise. Sat count/HDOP were tried as a
  47. * "has it converged" signal and rejected: both fluctuate (sat count can
  48. * drop as easily as rise while the receiver's tracking set changes) in a
  49. * way that gave false "stable" readings well before the position had
  50. * actually settled. A fixed wall-clock window matches what the real data
  51. * showed converging by, without depending on a quality signal that proved
  52. * unreliable for this specific purpose in practice.
  53. * Withhold distance/elevation/GPX credit for points within this window
  54. * after a gap, while still updating the Kalman/last-position state normally
  55. * so a genuine tunnel-with-real-displacement case isn't penalized beyond
  56. * losing GAP_SETTLE_SECONDS worth of distance right at the exit. */
  57. #define GAP_SETTLE_SECONDS 18 /* real data showed convergence taking ~15-17s after reacquisition */
  58. struct prev_points_s {
  59. struct location_s data[PREV_POINTS_LENGTH];
  60. unsigned char start;
  61. unsigned char count;
  62. };
  63. struct avg_store_s {
  64. float lat;
  65. float lon;
  66. time_t time;
  67. };
  68. struct elevation_s {
  69. unsigned char initialized;
  70. float smoothed_alt;
  71. float baseline;
  72. time_t last_time;
  73. };
  74. static struct gpx_s {
  75. struct prev_points_s prev_points;
  76. unsigned char avg_count;
  77. unsigned char paused;
  78. unsigned char point_count;
  79. struct avg_store_s avg_store;
  80. struct location_s last_saved;
  81. struct location_s last_distance_point; /* Last accepted point for distance calculation */
  82. struct kalman_s kalman[2];
  83. struct elevation_s elevation;
  84. time_t gap_settle_until; /* points before this time are held back post-gap; 0 = not settling */
  85. } gpx;
  86. float kalman_predict(struct kalman_s *k, float data);
  87. void kalman_init(struct kalman_s *k);
  88. float elevation_process(float alt, time_t time);
  89. void add_elevation_filtered(float amount, unsigned char is_gain);
  90. void prev_points_append(struct location_s *new){
  91. gpx.prev_points.data[(gpx.prev_points.start + gpx.prev_points.count)%PREV_POINTS_LENGTH] = *new;
  92. if(++gpx.prev_points.count > PREV_POINTS_LENGTH){
  93. gpx.prev_points.count--;
  94. gpx.prev_points.start++;
  95. gpx.prev_points.start %= PREV_POINTS_LENGTH;
  96. }
  97. }
  98. struct location_s *prev_points_get(unsigned char index){
  99. unsigned char i, addr = gpx.prev_points.start;
  100. for(i=0; i<index; i++){
  101. addr++;
  102. }
  103. addr %= PREV_POINTS_LENGTH;
  104. return &gpx.prev_points.data[addr];
  105. }
  106. unsigned char gpx_init(FIL *file) {
  107. unsigned int bw;
  108. kalman_init(&gpx.kalman[0]);
  109. kalman_init(&gpx.kalman[1]);
  110. gpx.prev_points.count = 0;
  111. gpx.avg_count = 0;
  112. gpx.last_saved.lon = 0;
  113. gpx.last_saved.lat = 0;
  114. gpx.last_saved.time = 0;
  115. gpx.last_distance_point.lon = 0;
  116. gpx.last_distance_point.lat = NAN; /* sentinel: no prior point yet, distinct from a real lat=0 (equator) fix */
  117. gpx.last_distance_point.time = 0;
  118. gpx.elevation.initialized = 0;
  119. gpx.point_count = 0;
  120. gpx.gap_settle_until = 0;
  121. gpx.paused = 1; /* make it add a <trkseg> tag */
  122. strcpy_P(buf, xml_header);
  123. strcat_P(buf, xml_trk_start);
  124. return f_write(file, buf, strlen(buf), &bw);
  125. }
  126. void gpx_save_single_point(struct location_s *loc) {
  127. FIL gpx;
  128. UINT bw;
  129. unsigned char err = 0;
  130. char *time = get_iso_time(loc->time, 1);
  131. iso_time_to_filename(time);
  132. xsprintf(buf, PSTR("%s-POINT.GPX"), time);
  133. xprintf(PSTR("Writing single point in %s\r\n"), buf);
  134. if ((err = f_open(&gpx, buf, FA_WRITE | FA_OPEN_ALWAYS))) {
  135. f_close(&gpx);
  136. // System.status = STATUS_FILE_OPEN_ERROR;
  137. xputs_P(PSTR("File open error\r\n"));
  138. return; /* Failed to open file */
  139. }
  140. strcpy_P(buf, xml_header);
  141. err |= f_write(&gpx, buf, strlen(buf), &bw);
  142. xsprintf(buf, PSTR("\t<wpt lat=\"%.8f\" lon=\"%.8f\"></wpt>\n</gpx>\n"), loc->lat, loc->lon);
  143. err |= f_write(&gpx, buf, strlen(buf), &bw);
  144. err |= f_close(&gpx);
  145. if (err) {
  146. /* TODO */
  147. }
  148. }
  149. unsigned char is_paused(void) {
  150. return tracking_is_paused();
  151. }
  152. unsigned char gpx_write(struct location_s *loc, FIL *file) {
  153. unsigned int bw;
  154. const char *time;
  155. unsigned char paused = tracking_is_paused();
  156. if (paused) {
  157. if (!gpx.paused) {
  158. strcpy_P(buf, xml_trkseg_end);
  159. gpx.paused = 1;
  160. gpx.point_count = 0;
  161. } else {
  162. return 0; /* nothing to store */
  163. }
  164. } else {
  165. if (gpx.paused) {
  166. strcpy_P(buf, xml_trkseg_start);
  167. f_write(file, buf, strlen(buf), &bw);
  168. gpx.paused = 0;
  169. }
  170. time = get_iso_time(loc->time, 0);
  171. xsprintf(buf, PSTR("\t\t\t<trkpt lat=\"%.8f\" lon=\"%.8f\">\n\t\t\t\t<ele>%.2f</ele>\n\t\t\t\t<time>%s</time>\n"), loc->lat, loc->lon, loc->alt, time);
  172. strcat_P(buf, PSTR("\t\t\t</trkpt>\n"));
  173. }
  174. {
  175. unsigned char ret;
  176. unsigned long int t0 = get_uptime_ms();
  177. ret = f_write(file, buf, strlen(buf), &bw);
  178. io_mark(t0);
  179. return ret;
  180. }
  181. }
  182. unsigned char gpx_close(FIL *file) {
  183. unsigned int bw;
  184. buf[0] = '\0';
  185. if (!gpx.paused)
  186. strcpy_P(buf, xml_trkseg_end);
  187. strcat_P(buf, PSTR("\t</trk>\n</gpx>\n"));
  188. f_write(file, buf, strlen(buf), &bw);
  189. return f_close(file);
  190. }
  191. /* Call once, whenever fix_gap_detected() reports a gap on the point about to
  192. * be passed to gpx_process_point() - see fix_gap_detected()'s comment for why
  193. * a gap (regardless of cause: fix loss, a 2D fix, poor HDOP, or the receiver
  194. * simply not producing lines) needs handling at all. Reset the Kalman/window
  195. * state that would otherwise splice the old and new fixes together, and
  196. * re-run the existing skip_points warm-up: a freshly reacquired fix can carry
  197. * a wildly wrong position/altitude for a few seconds while it converges (seen
  198. * after a real signal loss - a nominally-3D, in-range-HDOP fix still a couple
  199. * hundred meters off in altitude until more satellites came in), the same way
  200. * a fix has never been trusted right after power-on. */
  201. void gpx_reset_gap(void) {
  202. if (get_flag(CONFFLAG_VERBOSE_LOG))
  203. xputs_P(PSTR("FIX GAP: resetting warm-up/filter state\r\n"));
  204. kalman_init(&gpx.kalman[0]);
  205. kalman_init(&gpx.kalman[1]);
  206. gpx.prev_points.count = 0;
  207. gpx.avg_count = 0;
  208. gpx.avg_store.lat = 0;
  209. gpx.avg_store.lon = 0;
  210. gpx.avg_store.time = 0;
  211. gpx.point_count = 0;
  212. /* Without this, the first point(s) accepted after the gap (still inside
  213. * the prev_points warm-up above, so exempt from the distdiff/altdiff
  214. * spike checks until that window refills) would add one uncontrolled
  215. * distance/elevation segment straight from the pre-gap position - which
  216. * may be stale by an arbitrary, unknown amount of real movement during
  217. * the gap - with none of the usual spike protection. Clearing it forces
  218. * distance/elevation accumulation to (re)start fresh from the first
  219. * post-gap point actually reached once warm-up completes, the same way
  220. * it starts fresh after power-on. */
  221. gpx.last_distance_point.lat = NAN;
  222. /* -1 (all-ones): "settling, but the deadline isn't computed yet" - the
  223. * real deadline needs the first post-gap point's own timestamp, which
  224. * isn't available here (this runs before gpx_process_point() sees that
  225. * point), so gpx_process_point() computes gap_settle_until = that
  226. * point's time + GAP_SETTLE_SECONDS on its first point after this. */
  227. gpx.gap_settle_until = (time_t)-1;
  228. }
  229. void gpx_process_point(struct location_s *loc, FIL *file){
  230. float lon_est, lon_err, lat_est, lat_err, dist;
  231. struct location_s *ptr;
  232. static struct location_s nloc;
  233. static struct location_s filtered_loc;
  234. if (tracking_is_paused())
  235. System.points_paused++;
  236. if (gpx.point_count < System.conf.skip_points) { /* Skipping initial points */
  237. gpx.point_count++;
  238. System.points_skipped++;
  239. return;
  240. }
  241. if (gpx.gap_settle_until) {
  242. if (gpx.gap_settle_until == (time_t)-1)
  243. gpx.gap_settle_until = loc->time + GAP_SETTLE_SECONDS;
  244. if (loc->time < gpx.gap_settle_until) {
  245. if (get_flag(CONFFLAG_VERBOSE_LOG))
  246. xputs_P(PSTR("GAP SETTLE: waiting for fix to stabilize\r\n"));
  247. System.points_skipped++;
  248. return;
  249. }
  250. gpx.gap_settle_until = 0; /* window elapsed - resume normal processing below */
  251. }
  252. /* Always apply Kalman filtering for distance calculation */
  253. lat_est = kalman_predict(&gpx.kalman[0], loc->lat);
  254. lon_est = kalman_predict(&gpx.kalman[1], loc->lon);
  255. filtered_loc.lat = lat_est;
  256. filtered_loc.lon = lon_est;
  257. filtered_loc.time = loc->time;
  258. filtered_loc.alt = loc->alt;
  259. if (get_flag(CONFFLAG_DISABLE_FILTERS)) {
  260. /* Write unfiltered data to GPX, but always calculate distance/elevation from filtered data */
  261. if (get_flag(CONFFLAG_VERBOSE_LOG))
  262. xputs_P(PSTR("Write with filters disabled\r\n"));
  263. gpx_write(loc, file);
  264. System.points_written++;
  265. System.points_accepted++;
  266. /* Calculate distance and elevation from filtered points */
  267. if (!isnan(gpx.last_distance_point.lat)) {
  268. float ele_change;
  269. dist = distance(&gpx.last_distance_point, &filtered_loc);
  270. if (dist > POS_JUMP_THRESHOLD)
  271. log_jump(0, dist);
  272. add_distance(dist);
  273. ele_change = filtered_loc.alt - gpx.last_distance_point.alt;
  274. if (fabs(ele_change) > ALT_JUMP_THRESHOLD)
  275. log_jump(1, ele_change);
  276. add_elevation(ele_change);
  277. }
  278. elevation_process(filtered_loc.alt, filtered_loc.time);
  279. gpx.last_distance_point = filtered_loc;
  280. } else {
  281. /* Apply Kalman error check. KALMAN_ERR_MAX is a threshold in degrees
  282. * of latitude - a degree of longitude covers less ground the further
  283. * from the equator, by a factor of cos(latitude), so lon_err is
  284. * scaled up here to the equivalent latitude-degree error before the
  285. * same threshold is applied to both. Without this, the effective
  286. * east-west tolerance in meters would shrink at higher latitudes
  287. * (e.g. roughly half as permissive at 60 degrees as at the equator),
  288. * causing spurious rejects that have nothing to do with fix quality. */
  289. lat_err = fabs(loc->lat - lat_est);
  290. lon_err = fabs(loc->lon - lon_est) * cosf(loc->lat * M_PI / 180.0);
  291. // xprintf(PSTR("lat_err: %e, lon_err: %e, limit: %e\r\n"), lat_err, lon_err, (float)KALMAN_ERR_MAX);
  292. if(lat_err > KALMAN_ERR_MAX || lon_err > KALMAN_ERR_MAX){
  293. if (get_flag(CONFFLAG_VERBOSE_LOG))
  294. xputs_P(PSTR("KALMAN REJECT\r\n"));
  295. log_reject(REJECT_REASON_KALMAN);
  296. return;
  297. }
  298. prev_points_append(&filtered_loc);
  299. if(gpx.prev_points.count == PREV_POINTS_LENGTH){
  300. float dist12 = distance(prev_points_get(0), prev_points_get(1));
  301. float dist34 = distance(prev_points_get(2), prev_points_get(3));
  302. float dist32 = distance(prev_points_get(2), prev_points_get(1));
  303. /* Same spike test as the distance check below, applied to altitude:
  304. * a poor-geometry fix (e.g. right after reacquiring signal) can pass
  305. * the satellite-count/HDOP gate in nmea.c yet still carry a wildly
  306. * wrong altitude, so catch it here the same way a position spike is
  307. * caught - a jump surrounded by two much smaller ones on both sides.
  308. * TODO: GSA's VDOP would be a more direct altitude-quality gate than
  309. * this shape heuristic, but VDOP is only in GSA, which is parsed
  310. * after GGA within the epoch (see nmea.c's gp_gga_parse HDOP
  311. * comment) - using it here would need buffering a point across the
  312. * epoch boundary until VDOP arrives. Left as a future refinement. */
  313. float alt12 = fabs(prev_points_get(0)->alt - prev_points_get(1)->alt);
  314. float alt34 = fabs(prev_points_get(2)->alt - prev_points_get(3)->alt);
  315. float alt32 = fabs(prev_points_get(2)->alt - prev_points_get(1)->alt);
  316. if (get_flag(CONFFLAG_VERBOSE_LOG))
  317. xprintf(PSTR("New distance: %fm\r\n"), dist32);
  318. /* A real spike must stand out on BOTH sides - bigger than the
  319. * older segment AND bigger than the newer one. The previous
  320. * "dist34 > dist12 && dist32 > dist12" shape only checked dist32
  321. * against dist12 twice from different angles, so it stayed true
  322. * for every step of a sustained, monotonic acceleration (dist12 <
  323. * dist32 < dist34, completely normal e.g. when picking up pace
  324. * while hiking) - rejecting real point after real point until the
  325. * trend broke, then crediting all of that accumulated real
  326. * distance to whichever single point finally got through, as one
  327. * artificially fast-looking GPX segment. Comparing dist32 against
  328. * both neighbors instead lets a genuine acceleration/deceleration
  329. * ramp through cleanly while still catching an isolated outlier
  330. * surrounded by two much smaller segments. */
  331. if(dist32 > dist12 && dist32 > dist34){
  332. if (get_flag(CONFFLAG_VERBOSE_LOG))
  333. xputs_P(PSTR("DISTANCE DIFF REJECT\r\n"));
  334. log_reject(REJECT_REASON_DISTDIFF);
  335. return;
  336. }
  337. if(alt32 > ALT_JUMP_THRESHOLD && alt32 > alt12 && alt32 > alt34){
  338. if (get_flag(CONFFLAG_VERBOSE_LOG))
  339. xputs_P(PSTR("ALTITUDE DIFF REJECT\r\n"));
  340. log_reject(REJECT_REASON_ALTDIFF);
  341. return;
  342. }
  343. if (dist32 > POS_JUMP_THRESHOLD)
  344. log_jump(0, dist32);
  345. ptr = prev_points_get(PREV_POINTS_LENGTH - 2);
  346. } else {
  347. if(gpx.prev_points.count >= PREV_POINTS_LENGTH-2){
  348. ptr = prev_points_get(gpx.prev_points.count - 2);
  349. if (get_flag(CONFFLAG_VERBOSE_LOG))
  350. xputs_P(PSTR("NEW\r\n"));
  351. } else {
  352. return;
  353. }
  354. }
  355. if(distance(&gpx.last_saved, ptr) < MIN_DIST_DELTA){
  356. if (get_flag(CONFFLAG_VERBOSE_LOG))
  357. xputs_P(PSTR("Too small position change REJECT\r\n"));
  358. log_reject(REJECT_REASON_MINDIST);
  359. return;
  360. }
  361. if (get_flag(CONFFLAG_VERBOSE_LOG))
  362. xputs_P(PSTR("ACCEPT\r\n"));
  363. log_reject_flush(); /* close out any pending reject burst now that good data has resumed */
  364. System.points_accepted++;
  365. /* Calculate distance and elevation for accepted point */
  366. if (!isnan(gpx.last_distance_point.lat)) {
  367. float ele_change;
  368. dist = distance(&gpx.last_distance_point, ptr);
  369. add_distance(dist);
  370. ele_change = ptr->alt - gpx.last_distance_point.alt;
  371. if (fabs(ele_change) > ALT_JUMP_THRESHOLD)
  372. log_jump(1, ele_change);
  373. add_elevation(ele_change);
  374. }
  375. elevation_process(ptr->alt, ptr->time);
  376. gpx.last_distance_point = *ptr;
  377. gpx.avg_store.lat += ptr->lat;
  378. gpx.avg_store.lon += ptr->lon;
  379. /* Use the newest accepted point's time, not e.g. the midpoint's: the
  380. * averaged lat/lon reflects real movement up through this instant, so
  381. * timestamping it any earlier claims that movement happened faster
  382. * than it did. This matters most exactly when the accepted points
  383. * feeding this average were NOT evenly spaced in real time (e.g. a
  384. * reject streak between two of them) - an earlier timestamp there
  385. * would understate the elapsed time for the position delta and
  386. * inflate the speed an external tool derives from consecutive GPX
  387. * points. */
  388. gpx.avg_store.time = ptr->time;
  389. if(++gpx.avg_count == AVG_COUNT){
  390. nloc.lat = gpx.avg_store.lat / AVG_COUNT;
  391. nloc.lon = gpx.avg_store.lon / AVG_COUNT;
  392. nloc.time = gpx.avg_store.time;
  393. nloc.alt = gpx.elevation.smoothed_alt; /* filtered (smoothed) altitude, for the filtered GPX write */
  394. gpx.avg_count = 0;
  395. gpx.avg_store.lat = 0;
  396. gpx.avg_store.lon = 0;
  397. gpx.avg_store.time = 0;
  398. gpx.last_saved = nloc;
  399. gpx_write(&nloc, file);
  400. System.points_written++;
  401. }
  402. }
  403. if (System.time_start == 0)
  404. System.time_start = utc;
  405. }
  406. void kalman_init(struct kalman_s *k){
  407. k->initialized = 0;
  408. k->P_last = 0;
  409. //the noise in the system
  410. k->Q = KALMAN_Q; // process variance
  411. k->R = KALMAN_R; // measurement variance
  412. k->K = 0;
  413. }
  414. float kalman_predict(struct kalman_s *k, float data){
  415. if(!k->initialized){
  416. //initial values for the kalman filter
  417. k->x_est_last = data;
  418. k->initialized = 1;
  419. return data;
  420. }
  421. //do a prediction
  422. float x_temp_est = k->x_est_last;
  423. float P_temp = k->P_last + k->Q;
  424. //calculate the Kalman gain
  425. k->K = P_temp * (1.0/(P_temp + k->R));
  426. //correct
  427. float x_est = x_temp_est + k->K * (data - x_temp_est);
  428. k->P_last = (1 - k->K) * P_temp;
  429. k->x_est_last = x_est;
  430. return x_est;
  431. }
  432. #define R_EARTH 6371e3 // m
  433. float distance(struct location_s *pos1, struct location_s *pos2){
  434. float lat1 = pos1->lat * M_PI / 180.0;
  435. float lat2 = pos2->lat * M_PI / 180.0;
  436. float dlat = (pos2->lat - pos1->lat) * M_PI / 180.0;
  437. float dlon = (pos2->lon - pos1->lon) * M_PI / 180.0;
  438. float a = sinf(dlat/2.0) * sinf(dlat/2.0) + cosf(lat1) * cosf(lat2) * sinf(dlon/2.0) * sinf(dlon/2.0);
  439. float c = 2 * atan2f(sqrtf(a), sqrtf(1-a));
  440. float ret = R_EARTH * c;
  441. return ret;
  442. }
  443. /* Initial bearing (degrees, 0-360, 0=north) from pos1 to pos2 - used by the
  444. * auto-pause resume-by-distance check (autopause.c) to tell a real, directed
  445. * departure from a smooth-but-directionless position drift (e.g. a satellite
  446. * geometry change reshaping the fix while the receiver doesn't move at all). */
  447. float bearing(struct location_s *pos1, struct location_s *pos2){
  448. float lat1 = pos1->lat * M_PI / 180.0;
  449. float lat2 = pos2->lat * M_PI / 180.0;
  450. float dlon = (pos2->lon - pos1->lon) * M_PI / 180.0;
  451. float y = sinf(dlon) * cosf(lat2);
  452. float x = cosf(lat1) * sinf(lat2) - sinf(lat1) * cosf(lat2) * cosf(dlon);
  453. float deg = atan2f(y, x) * 180.0 / M_PI;
  454. if (deg < 0)
  455. deg += 360.0;
  456. return deg;
  457. }
  458. void add_distance(float dist) {
  459. unsigned char paused = tracking_is_paused();
  460. if (!paused)
  461. System.distance += (dist+0.005)*100.0;
  462. if (get_flag(CONFFLAG_VERBOSE_LOG))
  463. xprintf(PSTR("Distance: %.2f m; sum: %.2f m\r\n"), (double)dist, (double)System.distance/100.0);
  464. }
  465. /* Unfiltered (raw), per-point gain/loss - kept only for comparison against the
  466. * filtered (dead-band+smoothed) totals in the periodic status line/session summary. */
  467. void add_elevation(float ele_change) {
  468. unsigned char paused = tracking_is_paused();
  469. if (!paused) {
  470. if (ele_change > 0) {
  471. System.elevation_gain_raw += (ele_change+0.05)*10.0;
  472. } else if (ele_change < 0) {
  473. System.elevation_loss_raw += (-ele_change+0.05)*10.0;
  474. }
  475. }
  476. if (get_flag(CONFFLAG_VERBOSE_LOG))
  477. xprintf(PSTR("Elevation change: %.1f m; raw gain: %.1f m, raw loss: %.1f m\r\n"),
  478. (double)ele_change, (double)System.elevation_gain_raw/10.0, (double)System.elevation_loss_raw/10.0);
  479. }
  480. /* Filtered (dead-band + smoothed) gain/loss, used for display, GPX ascent stats
  481. * and the periodic status line/session summary. */
  482. void add_elevation_filtered(float amount, unsigned char is_gain) {
  483. unsigned long int dm = (unsigned long int)(amount*10.0 + 0.5);
  484. if (is_gain)
  485. System.elevation_gain += dm;
  486. else
  487. System.elevation_loss += dm;
  488. }
  489. /* Exponential smoothing (~30s time constant) followed by a 5m dead-band on the
  490. * result, so that a step in the smoothed altitude only counts once it clears
  491. * the dead-band, and only the amount past the dead-band edge is credited. */
  492. float elevation_process(float alt, time_t time) {
  493. unsigned char paused = tracking_is_paused();
  494. unsigned long int alt_dm;
  495. float delta;
  496. if (!gpx.elevation.initialized) {
  497. gpx.elevation.smoothed_alt = alt;
  498. gpx.elevation.baseline = alt;
  499. gpx.elevation.last_time = time;
  500. gpx.elevation.initialized = 1;
  501. } else {
  502. float dt = (float)(time - gpx.elevation.last_time);
  503. float alpha;
  504. if (dt <= 0)
  505. dt = 1.0;
  506. gpx.elevation.last_time = time;
  507. alpha = dt / (ELEV_SMOOTH_TAU + dt);
  508. gpx.elevation.smoothed_alt += alpha * (alt - gpx.elevation.smoothed_alt);
  509. }
  510. if (gpx.elevation.smoothed_alt > 0) {
  511. alt_dm = (unsigned long int)(gpx.elevation.smoothed_alt*10.0 + 0.5);
  512. if (alt_dm > System.alt_max)
  513. System.alt_max = alt_dm;
  514. }
  515. if (paused) {
  516. /* Discard the dead-band reference drift accumulated while paused,
  517. * the same way distance/raw elevation data is discarded when paused. */
  518. gpx.elevation.baseline = gpx.elevation.smoothed_alt;
  519. } else {
  520. delta = gpx.elevation.smoothed_alt - gpx.elevation.baseline;
  521. if (delta > ELEV_DEADBAND) {
  522. add_elevation_filtered(delta - ELEV_DEADBAND, 1);
  523. gpx.elevation.baseline = gpx.elevation.smoothed_alt - ELEV_DEADBAND;
  524. } else if (delta < -ELEV_DEADBAND) {
  525. add_elevation_filtered(-delta - ELEV_DEADBAND, 0);
  526. gpx.elevation.baseline = gpx.elevation.smoothed_alt + ELEV_DEADBAND;
  527. }
  528. }
  529. return gpx.elevation.smoothed_alt;
  530. }