Skip to content

Commit

Permalink
Fix cpr_nl for high latitude.
Browse files Browse the repository at this point in the history 
cpr_nl was returning 0 in some cases, which could lead to panics. This
change switches cpr_nl from a trigonometry-based implementation
to (essentially) a table lookup.

I added a couple unit tests, and a criterion benchmark for cpr_nl to
validate the table-based speedup.

On my machine this version is, in the worst case situation, 40%
faster. In the best case it's 89% faster.

```
cpr_nl_high_lat         time:   [10.764 ns 10.795 ns 10.830 ns]
                        change: [-40.716% -39.922% -39.255%] (p = 0.00 < 0.05)
                        Performance has improved.

cpr_nl_low_lat          time:   [2.2473 ns 2.2617 ns 2.2767 ns]
                        change: [-89.576% -89.479% -89.386%] (p = 0.00 < 0.05)
                        Performance has improved.
```

Fixes issue asmarques#3.
  • Loading branch information
@wiseman
wiseman committed Jan 22, 2021
1 parent 0eb41ce commit 0e24334
Show file tree
Hide file tree
Showing 2 changed files with 226 additions and 6 deletions.
12 changes: 11 additions & 1 deletion benches/parser.rs
View file
Original file line number Diff line number Diff line change
@@ -1,4 +1,4 @@
use adsb::parse_binary;
use adsb::{cpr::cpr_nl, parse_binary};
use criterion::{black_box, criterion_group, criterion_main, Criterion};

fn criterion_benchmark(c: &mut Criterion) {
Expand Down Expand Up @@ -30,6 +30,16 @@ fn criterion_benchmark(c: &mut Criterion) {
))
})
});
c.bench_function("cpr_nl_high_lat", |b| {
b.iter(|| {
cpr_nl(black_box(89.0))
})
});
c.bench_function("cpr_nl_low_lat", |b| {
b.iter(|| {
cpr_nl(black_box(0.0))
})
});
}

criterion_group!(benches, criterion_benchmark);
Expand Down
220 changes: 215 additions & 5 deletions src/cpr.rs
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,17 +2,199 @@

use crate::types::{CPRFrame, Parity, Position};
use std::cmp;
use std::f64::consts::PI;

const NZ: f64 = 15.0;
const D_LAT_EVEN: f64 = 360.0 / (4.0 * NZ);
const D_LAT_ODD: f64 = 360.0 / (4.0 * NZ - 1.0);
const CPR_MAX: f64 = 131_072.0;

fn cpr_nl(lat: f64) -> u64 {
let x = 1.0 - (PI / (2.0 * NZ)).cos();
let y = ((PI / 180.0) * lat).cos().powi(2);
((2.0 * PI) / (1.0 - (x / y)).acos()).floor() as u64
// The NL function uses the precomputed table from 1090-WP-9-14
// This code is translated from https://github.com/wiedehopf/readsb/blob/dev/cpr.c

pub fn cpr_nl(lat: f64) -> u64 {
let mut lat = lat;
if lat < 0.0 {
// Table is symmetric about the equator
lat = -lat;
}
if lat < 29.91135686 {
if lat < 10.47047130 {
return 59;
}
if lat < 14.82817437 {
return 58;
}
if lat < 18.18626357 {
return 57;
}
if lat < 21.02939493 {
return 56;
}
if lat < 23.54504487 {
return 55;
}
if lat < 25.82924707 {
return 54;
}
if lat < 27.93898710 {
return 53;
}
// < 29.91135686
return 52;
}
if lat < 44.19454951 {
if lat < 31.77209708 {
return 51;
}
if lat < 33.53993436 {
return 50;
}
if lat < 35.22899598 {
return 49;
}
if lat < 36.85025108 {
return 48;
}
if lat < 38.41241892 {
return 47;
}
if lat < 39.92256684 {
return 46;
}
if lat < 41.38651832 {
return 45;
}
if lat < 42.80914012 {
return 44;
}
// < 44.19454951
return 43;
}
if lat < 59.95459277 {
if lat < 45.54626723 {
return 42;
}
if lat < 46.86733252 {
return 41;
}
if lat < 48.16039128 {
return 40;
}
if lat < 49.42776439 {
return 39;
}
if lat < 50.67150166 {
return 38;
}
if lat < 51.89342469 {
return 37;
}
if lat < 53.09516153 {
return 36;
}
if lat < 54.27817472 {
return 35;
}
if lat < 55.44378444 {
return 34;
}
if lat < 56.59318756 {
return 33;
}
if lat < 57.72747354 {
return 32;
}
if lat < 58.84763776 {
return 31;
}
// < 59.95459277
return 30;
}
if lat < 61.04917774 {
return 29;
}
if lat < 62.13216659 {
return 28;
}
if lat < 63.20427479 {
return 27;
}
if lat < 64.26616523 {
return 26;
}
if lat < 65.31845310 {
return 25;
}
if lat < 66.36171008 {
return 24;
}
if lat < 67.39646774 {
return 23;
}
if lat < 68.42322022 {
return 22;
}
if lat < 69.44242631 {
return 21;
}
if lat < 70.45451075 {
return 20;
}
if lat < 71.45986473 {
return 19;
}
if lat < 72.45884545 {
return 18;
}
if lat < 73.45177442 {
return 17;
}
if lat < 74.43893416 {
return 16;
}
if lat < 75.42056257 {
return 15;
}
if lat < 76.39684391 {
return 14;
}
if lat < 77.36789461 {
return 13;
}
if lat < 78.33374083 {
return 12;
}
if lat < 79.29428225 {
return 11;
}
if lat < 80.24923213 {
return 10;
}
if lat < 81.19801349 {
return 9;
}
if lat < 82.13956981 {
return 8;
}
if lat < 83.07199445 {
return 7;
}
if lat < 83.99173563 {
return 6;
}
if lat < 84.89166191 {
return 5;
}
if lat < 85.75541621 {
return 4;
}
if lat < 86.53536998 {
return 3;
}
if lat < 87.00000000 {
return 2;
}
return 1;
}

/// Calculates a globally unambiguous position based on a pair of frames containing position information
Expand Down Expand Up @@ -100,6 +282,13 @@ mod tests {
use super::*;
use assert_approx_eq::assert_approx_eq;

#[test]
fn test_cpr_nl() {
assert_eq!(cpr_nl(89.9), 1);
assert_eq!(cpr_nl(-89.9), 1);
assert_eq!(cpr_nl(86.9), 2);
}

#[test]
fn cpr_calculate_position() {
let odd = CPRFrame {
Expand All @@ -122,4 +311,25 @@ mod tests {
assert_approx_eq!(position.latitude, 52.2572021484375);
assert_approx_eq!(position.longitude, 3.91937255859375);
}

#[test]
fn cpr_calculate_position_high_lat() {
let even = CPRFrame {
position: Position {
latitude: 108011.0,
longitude: 110088.0,
},
parity: Parity::Even,
};
let odd = CPRFrame {
position: Position {
latitude: 75050.0,
longitude: 36777.0,
},
parity: Parity::Odd,
};
let position = get_position((&even, &odd)).unwrap();
assert_approx_eq!(position.latitude, 88.91747426178496);
assert_approx_eq!(position.longitude, 101.01104736328125);
}
}

0 comments on commit 0e24334

Please sign in to comment.