use std::sync::Once;
use std::f64::consts::{PI, FRAC_PI_2, FRAC_PI_4};
pub const SQRT6: f64 = 2.449_489_742_783_178_f64;
const ONE_OVER_SQRT6: f64 = 0.408_248_290_463_863_f64;
const HALF: f64 = 0.5_f64;
const EPS_POLE: f64 = 1e-13_f64;
pub const TWICE_PI: f64 = 2.0 * PI;
pub const FOUR_OVER_PI: f64 = 4_f64 / PI;
pub const PI_OVER_FOUR: f64 = 0.25_f64 * PI;
pub const DEPTH_MAX: u8 = 29;
pub const NSIDE_MAX: u32 = 536870912;
pub const TRANSITION_LATITUDE: f64 = 0.729_727_656_226_966_3_f64; pub const TRANSITION_Z: f64 = 2_f64 / 3_f64;
pub const ONE_OVER_TRANSITION_Z: f64 = 1.5_f64;
pub const F64_SIGN_BIT_MASK: u64 = 0x8000000000000000;
pub const F64_BUT_SIGN_BIT_MASK: u64 = 0x7FFFFFFFFFFFFFFF;
static SMALLER_EDGE2OPEDGE_DIST: [f64; 30] = [
0.8410686705685088, 0.37723631722170053, 0.18256386461918295, 0.09000432499034523, 0.04470553761855741, 0.02228115704023076, 0.011122977211214961, 0.005557125022105058, 0.0027774761500209185, 0.0013884670480328143, 6.941658374603201E-4, 3.4706600585087755E-4, 1.7352877579970442E-4, 8.676333125510362E-5, 4.338140148342286E-5, 2.1690634707822447E-5, 1.084530084565172E-5, 5.422646295795749E-6, 2.711322116099695E-6, 1.3556608000873442E-6, 6.778303355805395E-7, 3.389151516386149E-7, 1.69457571754776E-7, 8.472878485272006E-8, 4.236439215502565E-8, 2.1182195982014308E-8, 1.0591097960375205E-8, 5.295548939447981E-9, 2.647774429917369E-9, 1.3238871881399636E-9 ];
pub static LAT_OF_SQUARE_CELL: f64 = 0.399_340_199_478_977_75_f64;
static COS_LAT_OF_SQUARE_CELL: f64 = 0.921_317_731_923_561_3_f64;
static mut CSTS_C2V: [Option<ConstantsC2V>; 30] = [
None, None, None, None, None, None, None, None, None, None, None, None, None, None, None,
None, None, None, None, None, None, None, None, None, None, None, None, None, None, None
];
static CSTS_C2V_INIT: [Once; 30] = [
Once::new(), Once::new(), Once::new(), Once::new(), Once::new(), Once::new(), Once::new(),
Once::new(), Once::new(), Once::new(), Once::new(), Once::new(), Once::new(), Once::new(),
Once::new(), Once::new(), Once::new(), Once::new(), Once::new(), Once::new(), Once::new(),
Once::new(), Once::new(), Once::new(), Once::new(), Once::new(), Once::new(), Once::new(),
Once::new(), Once::new()
];
fn get_or_create(depth: u8) -> &'static ConstantsC2V {
unsafe {
match CSTS_C2V[depth as usize] {
Some(ref v) => return v,
None => {
CSTS_C2V_INIT[depth as usize].call_once(|| {
CSTS_C2V[depth as usize] = Some(ConstantsC2V::new(depth));
});
},
}
match CSTS_C2V[depth as usize] {
Some(ref v) => v,
_ => unreachable!(),
}
}
}
struct ConstantsC2V {
slope_npc: f64,
intercept_npc: f64,
slope_eqr: f64,
intercept_eqr: f64,
coeff_x2_eqr: f64,
coeff_cst_eqr: f64,
}
impl ConstantsC2V {
fn new(depth: u8) -> ConstantsC2V {
let nside = nside_unsafe(depth);
let dist_cw = 1.0_f64 / (nside as f64); let one_min_dist_cw = 1.0_f64 - dist_cw;
let lat_north = f64::asin(1_f64 - (pow2(one_min_dist_cw) / 3_f64));
let d_min_npc = lat_north - TRANSITION_LATITUDE;
let d_max_npc: f64 = sphe_dist(
squared_half_segment(
FRAC_PI_4 * dist_cw, d_min_npc,
lat_north.cos(), TRANSITION_LATITUDE.cos()));
debug_assert!(d_min_npc < d_max_npc);
let slope_npc: f64 = (d_max_npc - d_min_npc) / (FRAC_PI_4 * one_min_dist_cw); let intercept_npc: f64 = d_min_npc; let d_min_eqr_top = PI_OVER_FOUR * dist_cw * COS_LAT_OF_SQUARE_CELL; let d_max_eqr_top = TRANSITION_LATITUDE - f64::asin(one_min_dist_cw * TRANSITION_Z);
debug_assert!(d_min_eqr_top < d_max_eqr_top);
let slope_eqr: f64 = (d_max_eqr_top - d_min_eqr_top) / (TRANSITION_LATITUDE - LAT_OF_SQUARE_CELL); let intercept_eqr: f64 = d_min_eqr_top - slope_eqr * LAT_OF_SQUARE_CELL; let d_max = PI_OVER_FOUR * dist_cw;
let coeff_cst_eqr: f64 = d_max;
let coeff_x2_eqr: f64 = (d_min_eqr_top - d_max) / pow2(LAT_OF_SQUARE_CELL);
ConstantsC2V {
slope_npc,
intercept_npc,
slope_eqr,
intercept_eqr,
coeff_x2_eqr,
coeff_cst_eqr
}
}
}
#[inline]
pub fn haversine_dist(p1_lon: f64, p1_lat: f64, p2_lon: f64, p2_lat: f64) -> f64 {
let shs = squared_half_segment(
p2_lon - p1_lon, p2_lat - p1_lat,
p1_lat.cos(), p2_lat.cos());
sphe_dist(shs)
}
#[inline]
fn sphe_dist(squared_half_segment: f64) -> f64 {
squared_half_segment.sqrt().asin().twice()
}
#[inline]
fn squared_half_segment(dlon: f64, dlat: f64, cos_lat1: f64, cos_lat2: f64) -> f64 {
dlat.half().sin().pow2() + cos_lat1 * cos_lat2 * dlon.half().sin().pow2()
}
#[inline]
fn to_squared_half_segment(spherical_distance: f64) -> f64 {
spherical_distance.half().sin().pow2()
}
#[inline]
fn pow2(x: f64) -> f64 {
x * x
}
pub trait Customf64 {
fn pow2(self) -> f64;
fn twice(self) -> f64;
fn half(self) -> f64;
fn div_eucl(self, rhs: f64) -> f64;
}
impl Customf64 for f64 {
#[inline]
fn pow2(self) -> f64 {
self * self }
#[inline]
fn twice(self) -> f64 {
2.0 * self }
#[inline]
fn half(self) -> f64 {
0.5 * self
}
#[inline]
fn div_eucl(self, rhs: f64) -> f64 {
let q = (self / rhs).trunc();
if self % rhs < 0.0 {
return if rhs > 0.0 { q - 1.0 } else { q + 1.0 }
}
q
}
}
pub fn largest_center_to_vertex_distance(depth: u8, lon: f64, lat: f64) -> f64 {
if depth == 0 {
return FRAC_PI_2 - TRANSITION_LATITUDE;
}
let lat_abs = lat.abs();
if lat_abs >= TRANSITION_LATITUDE {
largest_c2v_dist_in_npc(lon, get_or_create(depth))
} else if lat_abs >= LAT_OF_SQUARE_CELL {
largest_c2v_dist_in_eqr_top(lat_abs, get_or_create(depth))
} else {
largest_c2v_dist_in_eqr_bottom(lat_abs, get_or_create(depth))
}
}
pub fn largest_center_to_vertex_distance_with_radius(depth: u8, lon: f64, lat: f64, radius: f64) -> f64 {
if depth == 0 {
return FRAC_PI_2 - TRANSITION_LATITUDE;
}
let lat_abs = lat.abs();
let lat_max = lat_abs + radius;
let lat_min = lat_abs - radius;
if lat_max >= TRANSITION_LATITUDE {
largest_c2v_dist_in_npc_with_radius(lon, radius, get_or_create(depth))
} else if lat_min >= LAT_OF_SQUARE_CELL {
largest_c2v_dist_in_eqr_top_with_radius(lat_abs, radius, get_or_create(depth))
} else if lat_max <= LAT_OF_SQUARE_CELL {
largest_c2v_dist_in_eqr_bottom_with_radius(lat_abs, radius, get_or_create(depth))
} else {
let csts = get_or_create(depth);
f64::max(
largest_c2v_dist_in_eqr_top_with_radius(lat_abs, radius, csts),
largest_c2v_dist_in_eqr_bottom_with_radius(lat_abs, radius, csts)
)
}
}
pub fn largest_center_to_vertex_distances_with_radius(mut from_depth: u8, to_depth: u8, lon: f64, lat: f64, radius: f64) -> Box<[f64]> {
let mut vec: Vec<f64> = Vec::with_capacity((to_depth - from_depth) as usize);
if from_depth == 0 {
vec.push(FRAC_PI_2 - TRANSITION_LATITUDE);
from_depth = 1_u8;
}
let lat_abs = lat.abs();
let lat_max = lat_abs + radius;
let lat_min = lat_abs - radius;
if lat_max >= TRANSITION_LATITUDE {
let mut lon = (FRAC_PI_4 - (lon % FRAC_PI_2)).abs();
lon = f64::min(lon + radius, FRAC_PI_4);
for depth in from_depth..to_depth {
let csts = get_or_create(depth);
vec.push(linear_approx(lon, csts.slope_npc, csts.intercept_npc));
}
} else if lat_min >= LAT_OF_SQUARE_CELL {
for depth in from_depth..to_depth {
vec.push(largest_c2v_dist_in_eqr_top(lat_max, get_or_create(depth)));
}
} else if lat_max <= LAT_OF_SQUARE_CELL {
let val_min = f64::max(lat_min, 0_f64);
for depth in from_depth..to_depth {
vec.push(largest_c2v_dist_in_eqr_bottom(val_min, get_or_create(depth)));
}
} else {
let val_max = f64::min(lat_max, TRANSITION_LATITUDE);
let val_min = f64::max(lat_min, 0_f64);
for depth in from_depth..to_depth {
let csts = get_or_create(depth);
vec.push(
f64::max(
largest_c2v_dist_in_eqr_top(val_max, csts),
largest_c2v_dist_in_eqr_bottom(val_min, csts)
)
);
}
}
vec.into_boxed_slice()
}
#[inline]
fn largest_c2v_dist_in_npc(lon: f64, csts: &ConstantsC2V) -> f64 {
let lon = (FRAC_PI_4 - (lon % FRAC_PI_2)).abs();
debug_assert!((0_f64..=FRAC_PI_4).contains(&lon));
linear_approx(lon, csts.slope_npc, csts.intercept_npc)
}
#[inline]
fn largest_c2v_dist_in_npc_with_radius(lon: f64, radius: f64, csts: &ConstantsC2V) -> f64 {
debug_assert!(0_f64 < radius);
let mut lon = (FRAC_PI_4 - (lon % FRAC_PI_2)).abs();
debug_assert!((0_f64..=FRAC_PI_4).contains(&lon));
lon = f64::min(lon + radius, FRAC_PI_4);
linear_approx(lon, csts.slope_npc, csts.intercept_npc)
}
#[inline]
fn largest_c2v_dist_in_eqr_top(lat_abs: f64, csts: &ConstantsC2V) -> f64 {
debug_assert!(LAT_OF_SQUARE_CELL <= lat_abs && lat_abs < TRANSITION_LATITUDE);
linear_approx(lat_abs, csts.slope_eqr, csts.intercept_eqr)
}
#[inline]
fn largest_c2v_dist_in_eqr_top_with_radius(lat_abs: f64, radius: f64, csts: &ConstantsC2V) -> f64 {
debug_assert!(0_f64 < radius);
debug_assert!(LAT_OF_SQUARE_CELL <= lat_abs && lat_abs < TRANSITION_LATITUDE);
largest_c2v_dist_in_eqr_top(f64::min(lat_abs + radius, TRANSITION_LATITUDE), csts)
}
#[inline]
fn largest_c2v_dist_in_eqr_bottom(lat_abs: f64, csts: &ConstantsC2V) -> f64 {
debug_assert!(0_f64 <= lat_abs && lat_abs <= LAT_OF_SQUARE_CELL);
csts.coeff_x2_eqr * pow2(lat_abs) + csts.coeff_cst_eqr
}
#[inline]
fn largest_c2v_dist_in_eqr_bottom_with_radius(lat_abs: f64, radius: f64, csts: &ConstantsC2V) -> f64 {
debug_assert!(0_f64 < radius);
debug_assert!(0_f64 <= lat_abs && lat_abs <= LAT_OF_SQUARE_CELL);
largest_c2v_dist_in_eqr_bottom(f64::max(lat_abs - radius, 0_f64), csts)
}
#[inline]
fn linear_approx(x: f64, slope: f64, intercept: f64) -> f64 {
slope * x + intercept
}
#[inline]
pub fn nside(depth: u8) -> u32 {
check_depth(depth);
nside_unsafe(depth)
}
#[inline]
pub const fn nside_unsafe(depth: u8) -> u32 {
1_u32 << depth
}
#[inline]
pub fn nside_square(delta_depth: u8) -> u64 {
check_depth(delta_depth);
nside_square_unsafe(delta_depth)
}
#[inline]
pub const fn nside_square_unsafe(delta_depth: u8) -> u64 {
1_u64 << (delta_depth << 1)
}
#[inline]
fn check_depth(depth: u8) {
assert!(is_depth(depth), "Expected depth in [0, 29]");
}
#[inline]
pub const fn is_depth(depth: u8) -> bool {
depth <= DEPTH_MAX
}
#[inline]
pub fn depth(nside: u32) -> u8 {
check_nside(nside);
depth_unsafe(nside)
}
#[inline]
pub const fn depth_unsafe(nside: u32) -> u8 {
nside.trailing_zeros() as u8
}
#[inline]
fn check_nside(nside: u32) {
assert!(is_nside(nside), "Nside must be a power of 2 in [1-2^29]");
}
#[inline]
pub fn is_nside(nside: u32) -> bool {
is_pow_of_2(nside) && nside > 0 && nside <= NSIDE_MAX
}
#[inline]
fn is_pow_of_2(x: u32) -> bool {
(x & (x - 1)) == 0
}
#[inline]
pub fn n_hash(depth: u8) -> u64 {
check_depth(depth);
n_hash_unsafe(depth)
}
#[inline]
pub const fn n_hash_unsafe(depth: u8) -> u64 { 12u64 << (depth << 1u8) }
#[inline]
pub fn has_best_starting_depth(d_max_rad: f64) -> bool {
d_max_rad < SMALLER_EDGE2OPEDGE_DIST[0]
}
#[inline]
pub fn best_starting_depth(d_max_rad: f64) -> u8 { assert!(d_max_rad < SMALLER_EDGE2OPEDGE_DIST[0],
"Too large value, use first function has_best_starting_depth");
if d_max_rad < SMALLER_EDGE2OPEDGE_DIST[29] {
29
} else if d_max_rad < SMALLER_EDGE2OPEDGE_DIST[15] {
if d_max_rad < SMALLER_EDGE2OPEDGE_DIST[22] {
if d_max_rad < SMALLER_EDGE2OPEDGE_DIST[25] {
if d_max_rad < SMALLER_EDGE2OPEDGE_DIST[27] {
if d_max_rad < SMALLER_EDGE2OPEDGE_DIST[28] {
28
} else {
27
}
} else if d_max_rad < SMALLER_EDGE2OPEDGE_DIST[26] {
26
} else {
25
}
} else if d_max_rad < SMALLER_EDGE2OPEDGE_DIST[24] {
24
} else if d_max_rad < SMALLER_EDGE2OPEDGE_DIST[23] {
23
} else {
22
}
} else if d_max_rad < SMALLER_EDGE2OPEDGE_DIST[18] {
if d_max_rad < SMALLER_EDGE2OPEDGE_DIST[20] {
if d_max_rad < SMALLER_EDGE2OPEDGE_DIST[21] {
21
} else {
20
}
} else if d_max_rad < SMALLER_EDGE2OPEDGE_DIST[19] {
19
} else {
18
}
} else if d_max_rad < SMALLER_EDGE2OPEDGE_DIST[17] {
17
} else if d_max_rad < SMALLER_EDGE2OPEDGE_DIST[16] {
16
} else {
15
}
} else if d_max_rad < SMALLER_EDGE2OPEDGE_DIST[7] {
if d_max_rad < SMALLER_EDGE2OPEDGE_DIST[11] {
if d_max_rad < SMALLER_EDGE2OPEDGE_DIST[13] {
if d_max_rad < SMALLER_EDGE2OPEDGE_DIST[14] {
14
} else {
13
}
} else if d_max_rad < SMALLER_EDGE2OPEDGE_DIST[12] {
12
} else {
11
}
} else if d_max_rad < SMALLER_EDGE2OPEDGE_DIST[9] {
if d_max_rad < SMALLER_EDGE2OPEDGE_DIST[10] {
10
} else {
9
}
} else if d_max_rad < SMALLER_EDGE2OPEDGE_DIST[8] {
8
} else {
7
}
} else if d_max_rad < SMALLER_EDGE2OPEDGE_DIST[3] {
if d_max_rad < SMALLER_EDGE2OPEDGE_DIST[5] {
if d_max_rad < SMALLER_EDGE2OPEDGE_DIST[6] {
6
} else {
5
}
} else if d_max_rad < SMALLER_EDGE2OPEDGE_DIST[4] {
4
} else {
3
}
} else if d_max_rad < SMALLER_EDGE2OPEDGE_DIST[2] {
2
} else if d_max_rad < SMALLER_EDGE2OPEDGE_DIST[1] {
1
} else {
0
}
}
#[inline]
pub fn proj(lon: f64, lat: f64) -> (f64, f64) {
check_lat(lat);
let lon = abs_sign_decompose(lon);
let lat = abs_sign_decompose(lat);
let x = pm1_offset_decompose(lon.abs * FOUR_OVER_PI);
let mut xy = (x.pm1, lat.abs);
if is_in_equatorial_region(lat.abs) {
proj_cea(&mut xy);
} else {
proj_collignon(&mut xy);
}
apply_offset_and_signs(&mut xy, x.offset, lon.sign, lat.sign);
xy
}
pub fn base_cell_from_proj_coo(x: f64, y: f64) -> u8 {
let mut x = 0.5 * ensures_x_is_positive(x);
let mut y = 0.5 * (y + 3.0);
let mut i = x as u8; debug_assert!(i < 4); let mut j = (y as u8) << 1; debug_assert!(j == 0 || j == 2 || j == 4);
x -= i as f64; debug_assert!((0.0..1.0).contains(&x));
y -= (j >> 1) as f64; debug_assert!((0.0..1.0).contains(&y));
let in_northwest = (x <= y) as u8; let in_southeast = (x >= 1.0 - y) as u8; i += in_southeast >> in_northwest; j += in_northwest + in_southeast;
if j == 6 { j = 4; } debug_assert!(j == 2 || j == 3 || j == 4);
((4 - j) << 2) + i
}
#[inline]
pub fn unproj(x: f64, y: f64) -> (f64, f64) {
check_y(y);
let x = abs_sign_decompose(x);
let y = abs_sign_decompose(y);
let lon = pm1_offset_decompose(x.abs);
let mut lonlat= (lon.pm1, y.abs);
if is_in_projected_equatorial_region(y.abs) {
deproj_cea(&mut lonlat);
} else {
deproj_collignon(&mut lonlat);
}
apply_offset_and_signs(&mut lonlat, lon.offset, x.sign, y.sign);
lonlat.0 *= FRAC_PI_4;
lonlat
}
#[inline]
pub(crate) fn ensures_x_is_positive(x: f64) -> f64 {
if x < 0.0 { x + 8.0 } else { x }
}
#[inline]
fn check_lat(lat: f64) {
assert!((-FRAC_PI_2..=FRAC_PI_2).contains(&lat));
}
#[inline]
fn check_lat_res(lat: f64) -> Result<(), String> {
if (-FRAC_PI_2..=FRAC_PI_2).contains(&lat) {
Ok(())
} else {
Err(format!("Wrong latitude. Expected value in [-pi/2, pi/2]. Actual: {}", lat))
}
}
#[inline]
fn check_y(y: f64) { assert!((-2f64..=2f64).contains(&y)); }
#[inline]
pub fn is_in_equatorial_region(abs_lat: f64) -> bool {
abs_lat <= TRANSITION_LATITUDE
}
#[inline]
pub fn is_in_projected_equatorial_region(abs_y: f64) -> bool { abs_y <= 1.0 }
struct AbsAndSign {
abs: f64,
sign: u64,
}
#[inline]
pub(crate) fn abs_sign_decompose(x: f64) -> AbsAndSign {
let bits = f64::to_bits(x);
AbsAndSign {
abs: f64::from_bits(bits & F64_BUT_SIGN_BIT_MASK),
sign: bits & F64_SIGN_BIT_MASK,
}
}
pub(crate) struct OffsetAndPM1 {
offset: u8, pm1: f64, }
#[inline]
pub(crate) fn pm1_offset_decompose(x: f64) -> OffsetAndPM1 {
let floor: u8 = x as u8;
let odd_floor: u8 = floor | 1u8;
OffsetAndPM1 {
offset: odd_floor & 7u8, pm1: x - (odd_floor as f64),
}
}
#[inline]
pub(crate) fn proj_cea(xy: &mut (f64, f64)) {
let (_, ref mut y) = *xy;
*y = f64::sin(*y) * ONE_OVER_TRANSITION_Z;
}
#[inline]
fn deproj_cea(lonlat: &mut (f64, f64)) {
let (_, ref mut lat) = *lonlat;
*lat = f64::asin((*lat) * TRANSITION_Z);
}
#[inline]
pub(crate) fn proj_collignon(xy: &mut (f64, f64)) {
let (ref mut x, ref mut y) = *xy;
*y = SQRT6 * f64::cos(HALF * *y + FRAC_PI_4);
*x *= *y;
*y = 2.0 - *y;
}
#[inline]
fn deproj_collignon(lonlat: &mut (f64, f64)) {
let (ref mut lon, ref mut lat) = *lonlat;
*lat = 2.0 - *lat;
if is_not_near_from_pole(*lat) { *lon /= *lat;
deal_with_numerical_approx_in_edges(lon);
} *lat *= ONE_OVER_SQRT6;
*lat = 2.0 * f64::acos(*lat) - FRAC_PI_2;
}
#[inline]
fn is_not_near_from_pole(sqrt_of_three_time_one_minus_sin_of: f64) -> bool {
sqrt_of_three_time_one_minus_sin_of > EPS_POLE
}
#[inline]
fn deal_with_numerical_approx_in_edges(lon: &mut f64) {
if *lon > 1.0 {
*lon = 1.0;
} else if *lon < -1.0 {
*lon = -1.0;
}
}
#[inline]
pub(crate) fn apply_offset_and_signs(ab: &mut (f64, f64), off: u8, a_sign: u64, b_sign: u64) {
let (ref mut a, ref mut b) = *ab;
*a += off as f64;
*a = f64::from_bits(f64::to_bits(*a) | a_sign);
*b = f64::from_bits(f64::to_bits(*b) | b_sign);
}
pub mod compass_point;
pub mod external_edge;
use crate::compass_point::{MainWind};
use crate::compass_point::MainWind::*;
pub fn neighbour(base_cell: u8, direction: MainWind) -> Option<u8> {
if direction == MainWind::C {
Some(base_cell)
} else {
let d0h_mod_4 = base_cell & 3_u8; match base_cell >> 2 { 0 => npc_neighbour(d0h_mod_4, direction),
1 => eqr_neighbour(d0h_mod_4, direction),
2 => spc_neighbour(d0h_mod_4, direction),
_ => panic!("Base cell must be in [0, 12["),
}
}
}
fn npc_neighbour(d0h_mod_4: u8, direction: MainWind) -> Option<u8> {
match direction {
S => base_cell_opt(iden(d0h_mod_4), 2),
SE => base_cell_opt(next(d0h_mod_4), 1),
SW => base_cell_opt(iden(d0h_mod_4), 1),
NE => base_cell_opt(next(d0h_mod_4), 0),
NW => base_cell_opt(prev(d0h_mod_4), 0),
N => base_cell_opt(oppo(d0h_mod_4), 0),
_ => None,
}
}
fn eqr_neighbour(d0h_mod_4: u8, direction: MainWind) -> Option<u8> {
match direction {
SE => base_cell_opt(iden(d0h_mod_4), 2),
E => base_cell_opt(next(d0h_mod_4), 1),
SW => base_cell_opt(prev(d0h_mod_4), 2),
NE => base_cell_opt(iden(d0h_mod_4), 0),
W => base_cell_opt(prev(d0h_mod_4), 1),
NW => base_cell_opt(prev(d0h_mod_4), 0),
_ => None,
}
}
fn spc_neighbour(d0h_mod_4: u8, direction: MainWind) -> Option<u8> {
match direction {
S => base_cell_opt(oppo(d0h_mod_4), 2),
SE => base_cell_opt(next(d0h_mod_4), 2),
SW => base_cell_opt(prev(d0h_mod_4), 2),
NE => base_cell_opt(next(d0h_mod_4), 1),
NW => base_cell_opt(iden(d0h_mod_4), 1),
N => base_cell_opt(iden(d0h_mod_4), 0),
_ => None,
}
}
pub fn edge_cell_direction_from_neighbour(base_cell: u8, inner_direction: &MainWind, neighbour_direction: &MainWind) -> MainWind {
match base_cell >> 2 { 0 => npc_egde_direction_from_neighbour(inner_direction, neighbour_direction),
1 => eqr_edge_direction_from_neighbour(inner_direction, neighbour_direction),
2 => spc_edge_direction_from_neighbour(inner_direction, neighbour_direction),
_ => panic!("Base cell must be in [0, 12["),
}
}
fn npc_egde_direction_from_neighbour(inner_direction: &MainWind, neighbour_direction: &MainWind) -> MainWind {
match neighbour_direction {
C => panic!("No neighbour in direction {:?}", &neighbour_direction),
E => match inner_direction {
N | NE => N,
E => panic!("No neighbour in direction {:?}", &neighbour_direction),
S | SE => neighbour_direction.opposite(),
_ => unreachable!(),
},
W => match inner_direction {
N | NW => N,
W => panic!("No neighbour in direction {:?}", &neighbour_direction),
S | SW => neighbour_direction.opposite(),
_ => unreachable!(),
},
NE => {
assert!(*inner_direction == N || *inner_direction == E || *inner_direction == NE);
NW
},
NW => {
assert!(*inner_direction == N || *inner_direction == W || *inner_direction == NW);
NE
},
N => match inner_direction {
N => N,
E | NE => W,
W | NW => E,
_ => unreachable!(),
},
_ => neighbour_direction.opposite(),
}
}
fn eqr_edge_direction_from_neighbour(_inner_direction: &MainWind, neighbour_direction: &MainWind) -> MainWind {
neighbour_direction.opposite()
}
fn spc_edge_direction_from_neighbour(inner_direction: &MainWind, neighbour_direction: &MainWind) -> MainWind {
match neighbour_direction {
C => panic!("No neighbour in direction {:?}", &neighbour_direction),
E => match inner_direction {
S | SE => S,
E => panic!("No neighbour in direction {:?}", &neighbour_direction),
N | NE => neighbour_direction.opposite(),
_ => unreachable!(),
},
W => match inner_direction {
S | SW => S,
W => panic!("No neighbour in direction {:?}", &neighbour_direction),
N | NW => neighbour_direction.opposite(),
_ => unreachable!(),
},
SE => {
assert!(*inner_direction == S || *inner_direction == E || *inner_direction == SE);
SW
},
SW => {
assert!(*inner_direction == S || *inner_direction == W || *inner_direction == SW);
SE
},
S => match inner_direction {
S => S,
E | SE => W,
W | SW => E,
_ => unreachable!(),
},
_ => neighbour_direction.opposite(),
}
}
pub fn direction_from_neighbour(base_cell: u8, neighbour_direction: &MainWind) -> MainWind {
match base_cell >> 2 { 0 => npc_direction_from_neighbour(neighbour_direction),
1 => eqr_direction_from_neighbour(neighbour_direction),
2 => spc_direction_from_neighbour(neighbour_direction),
_ => panic!("Base cell must be in [0, 12["),
}
}
fn npc_direction_from_neighbour(neighbour_direction: &MainWind) -> MainWind {
match neighbour_direction {
E | W | C => panic!("No neighbour in direction {:?}", &neighbour_direction),
NE => NW,
NW => NE,
N => N,
_ => neighbour_direction.opposite(),
}
}
fn eqr_direction_from_neighbour(neighbour_direction: &MainWind) -> MainWind {
match neighbour_direction {
S | N | C => panic!("No neighbour in direction {:?}", &neighbour_direction),
_ => neighbour_direction.opposite(),
}
}
fn spc_direction_from_neighbour(neighbour_direction: &MainWind) -> MainWind {
match neighbour_direction {
E | W | C => panic!("No neighbour in direction {:?}", &neighbour_direction),
S => S,
SE => SW,
SW => SE,
_ => neighbour_direction.opposite(),
}
}
#[inline]
fn prev(mod4: u8) -> u8 {
debug_assert!(mod4 < 4);
(((mod4 as i8) - 1) & 3) as u8
}
#[inline]
fn next(mod4: u8) -> u8 {
debug_assert!(mod4 < 4);
(mod4 + 1) & 3
}
#[inline]
fn oppo(mod4: u8) -> u8 {
debug_assert!(mod4 < 4);
(mod4 + 2) & 3
}
#[inline]
fn iden(mod4: u8) -> u8 {
debug_assert!(mod4 < 4);
mod4
}
#[inline]
fn base_cell_opt(i: u8, j: u8) -> Option<u8> {
Some(base_cell(i, j))
}
#[inline]
fn base_cell(i: u8, j: u8) -> u8 {
debug_assert!(i < 4 && j < 3);
(j << 2) + i
}
pub mod nested;
pub mod ring;
mod xy_geom;
pub mod sph_geom;
mod special_points_finder;
#[cfg(test)]
mod tests {
use super::*;
#[test]
fn testok_nside() {
assert_eq!(1, nside(0));
assert_eq!(2, nside(1));
assert_eq!(4, nside(2));
assert_eq!(8, nside(3));
assert_eq!(16, nside(4));
assert_eq!(32, nside(5));
assert_eq!(64, nside(6));
assert_eq!(128, nside(7));
assert_eq!(256, nside(8));
assert_eq!(512, nside(9));
assert_eq!(1024, nside(10));
assert_eq!(2048, nside(11));
assert_eq!(4096, nside(12));
assert_eq!(8192, nside(13));
assert_eq!(16384, nside(14));
assert_eq!(32768, nside(15));
assert_eq!(65536, nside(16));
assert_eq!(131072, nside(17));
assert_eq!(262144, nside(18));
assert_eq!(524288, nside(19));
assert_eq!(1048576, nside(20));
assert_eq!(2097152, nside(21));
assert_eq!(4194304, nside(22));
assert_eq!(8388608, nside(23));
assert_eq!(16777216, nside(24));
assert_eq!(33554432, nside(25));
assert_eq!(67108864, nside(26));
assert_eq!(134217728, nside(27));
assert_eq!(268435456, nside(28));
assert_eq!(536870912, nside(29));
}
#[test]
#[should_panic]
fn testpanic_nside() {
nside(30);
}
#[test]
fn testok_proj() {
fn dist(p1: (f64, f64), p2: (f64, f64)) -> f64 {
f64::sqrt((p2.0 - p1.0) * (p2.0 - p1.0) + (p2.1 - p1.1) * (p2.1 - p1.1))
}
assert!(dist((0.0, 0.0), proj(0.0 * FRAC_PI_2, 0.0)) < 1e-15);
assert!(dist((0.0, 1.0), proj(0.0 * FRAC_PI_2, TRANSITION_LATITUDE)) < 1e-15);
assert!(dist((1.0, 2.0), proj(0.0 * FRAC_PI_2 + FRAC_PI_4, FRAC_PI_2)) < 1e-15);
assert!(dist((2.0, 1.0), proj(1.0 * FRAC_PI_2, TRANSITION_LATITUDE)) < 1e-15);
assert!(dist((3.0, 2.0), proj(1.0 * FRAC_PI_2 + FRAC_PI_4, FRAC_PI_2)) < 1e-15);
assert!(dist((4.0, 1.0), proj(2.0 * FRAC_PI_2, TRANSITION_LATITUDE)) < 1e-15);
assert!(dist((5.0, 2.0), proj(2.0 * FRAC_PI_2 + FRAC_PI_4, FRAC_PI_2)) < 1e-15);
assert!(dist((6.0, 1.0), proj(3.0 * FRAC_PI_2, TRANSITION_LATITUDE)) < 1e-15);
assert!(dist((7.0, 2.0), proj(3.0 * FRAC_PI_2 + FRAC_PI_4, FRAC_PI_2)) < 1e-15);
assert!(dist((0.0, 1.0), proj(4.0 * FRAC_PI_2, TRANSITION_LATITUDE)) < 1e-15);
assert!(dist((0.0, 0.0), proj(4.0 * FRAC_PI_2, 0.0)) < 1e-15);
assert!(dist((0.0, -1.0), proj(0.0 * FRAC_PI_2, -TRANSITION_LATITUDE)) < 1e-15);
assert!(dist((1.0, -2.0), proj(0.0 * FRAC_PI_2 + FRAC_PI_4, -FRAC_PI_2)) < 1e-15);
assert!(dist((2.0, -1.0), proj(1.0 * FRAC_PI_2, -TRANSITION_LATITUDE)) < 1e-15);
assert!(dist((3.0, -2.0), proj(1.0 * FRAC_PI_2 + FRAC_PI_4, -FRAC_PI_2)) < 1e-15);
assert!(dist((4.0, -1.0), proj(2.0 * FRAC_PI_2, -TRANSITION_LATITUDE)) < 1e-15);
assert!(dist((5.0, -2.0), proj(2.0 * FRAC_PI_2 + FRAC_PI_4, -FRAC_PI_2)) < 1e-15);
assert!(dist((6.0, -1.0), proj(3.0 * FRAC_PI_2, -TRANSITION_LATITUDE)) < 1e-15);
assert!(dist((7.0, -2.0), proj(3.0 * FRAC_PI_2 + FRAC_PI_4, -FRAC_PI_2)) < 1e-15);
assert!(dist((0.0, -1.0), proj(4.0 * FRAC_PI_2, -TRANSITION_LATITUDE)) < 1e-15);
assert!(dist((-0.0, 0.0), proj(-0.0 * FRAC_PI_2, 0.0)) < 1e-15);
assert!(dist((-0.0, 1.0), proj(-0.0 * FRAC_PI_2, TRANSITION_LATITUDE)) < 1e-15);
assert!(dist((-1.0, 2.0), proj(-0.0 * FRAC_PI_2 - FRAC_PI_4, FRAC_PI_2)) < 1e-15);
assert!(dist((-2.0, 1.0), proj(-1.0 * FRAC_PI_2, TRANSITION_LATITUDE)) < 1e-15);
assert!(dist((-3.0, 2.0), proj(-1.0 * FRAC_PI_2 - FRAC_PI_4, FRAC_PI_2)) < 1e-15);
assert!(dist((-4.0, 1.0), proj(-2.0 * FRAC_PI_2, TRANSITION_LATITUDE)) < 1e-15);
assert!(dist((-5.0, 2.0), proj(-2.0 * FRAC_PI_2 - FRAC_PI_4, FRAC_PI_2)) < 1e-15);
assert!(dist((-6.0, 1.0), proj(-3.0 * FRAC_PI_2, TRANSITION_LATITUDE)) < 1e-15);
assert!(dist((-7.0, 2.0), proj(-3.0 * FRAC_PI_2 - FRAC_PI_4, FRAC_PI_2)) < 1e-15);
assert!(dist((-0.0, 1.0), proj(-4.0 * FRAC_PI_2, TRANSITION_LATITUDE)) < 1e-15);
assert!(dist((-0.0, 0.0), proj(-4.0 * FRAC_PI_2, 0.0)) < 1e-15);
assert!(dist((-0.0, -1.0), proj(-0.0 * FRAC_PI_2, -TRANSITION_LATITUDE)) < 1e-15);
assert!(dist((-1.0, -2.0), proj(-0.0 * FRAC_PI_2 - FRAC_PI_4, -FRAC_PI_2)) < 1e-15);
assert!(dist((-2.0, -1.0), proj(-1.0 * FRAC_PI_2, -TRANSITION_LATITUDE)) < 1e-15);
assert!(dist((-3.0, -2.0), proj(-1.0 * FRAC_PI_2 - FRAC_PI_4, -FRAC_PI_2)) < 1e-15);
assert!(dist((-4.0, -1.0), proj(-2.0 * FRAC_PI_2, -TRANSITION_LATITUDE)) < 1e-15);
assert!(dist((-5.0, -2.0), proj(-2.0 * FRAC_PI_2 - FRAC_PI_4, -FRAC_PI_2)) < 1e-15);
assert!(dist((-6.0, -1.0), proj(-3.0 * FRAC_PI_2, -TRANSITION_LATITUDE)) < 1e-15);
assert!(dist((-7.0, -2.0), proj(-3.0 * FRAC_PI_2 - FRAC_PI_4, -FRAC_PI_2)) < 1e-15);
assert!(dist((-0.0, -1.0), proj(-4.0 * FRAC_PI_2, -TRANSITION_LATITUDE)) < 1e-15);
}
#[test]
#[should_panic]
fn testpanic_proj_1() {
proj(0.0, -1.58);
}
#[test]
#[should_panic]
fn testpanic_proj_2() {
proj(3.14159, -1.58);
}
#[test]
fn testok_shs() {
let ang_dist = 0.24;
let shs = to_squared_half_segment(ang_dist);
let ang_dist_2 = sphe_dist(shs);
assert!((ang_dist - ang_dist_2).abs() < 1e-4);
}
#[test]
fn testok_largest_center_to_vertex_distance() {
let depth = 8;
let d1_mas = largest_center_to_vertex_distance(depth, 0.0, 0.0)
.to_degrees() * 3600.0;
let d2_mas = largest_center_to_vertex_distance(depth, 0.0, LAT_OF_SQUARE_CELL)
.to_degrees() * 3600.0;
let d3_mas = largest_center_to_vertex_distance(depth, 0.0, TRANSITION_LATITUDE)
.to_degrees() * 3600.0;
let d4_mas = largest_center_to_vertex_distance(depth, PI / 4.0, TRANSITION_LATITUDE)
.to_degrees() * 3600.0;
let d5_mas = largest_center_to_vertex_distance(depth, 0.0, 0.5 * LAT_OF_SQUARE_CELL)
.to_degrees() * 3600.0;
let d6_mas = largest_center_to_vertex_distance(depth, 0.0, 0.5 * TRANSITION_LATITUDE)
.to_degrees() * 3600.0;
println!("d1: {}", d1_mas); println!("d2: {}", d2_mas); println!("d3: {}", d3_mas); println!("d4: {}", d4_mas); println!("d5: {} ~620?", d5_mas); println!("d6: {} ~660?", d6_mas);
assert!((d1_mas - 632.8125000000000).abs() < 1e-12);
assert!((d2_mas - 583.0213772328785).abs() < 1e-12);
assert!((d3_mas - 861.2025252838432).abs() < 1e-12);
assert!((d4_mas - 720.3786531802374).abs() < 1e-12);
assert!((d5_mas - 620.3647193082197).abs() < 1e-12);
assert!((d6_mas - 591.2475292479544).abs() < 1e-12);
}
}