1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
use crate::coords::LonLatT;
use crate::{
FRAC_PI_2, FRAC_PI_4, LonLat, SQRT_6, TRANSITION_LATITUDE, TRANSITION_Z, proj, unchecked,
};
impl super::Layer {
/// Returns the cell number (hash value) associated with the given position on the unit sphere
/// # Inputs
/// - `lon`: longitude in radians, support reasonably large positive and negative values
/// producing accurate results with a naive range reduction like modulo 2*pi
/// (i.e. without having to resort on Cody-Waite or Payne Hanek range reduction).
/// - `lat`: latitude in radians, must be in `[-pi/2, pi/2]`
/// # Output
/// - the cell number (hash value) associated with the given position on the unit sphere,
/// in `[0, 12*nside^2[`
/// # Panics
/// If `lat` **not in** `[-pi/2, pi/2]`, this method panics.
/// # Examples
/// ```rust
/// use healpix::checked::nside;
/// use healpix::coords::Degrees;
/// use healpix::get;
///
/// let depth = 12_u8;
/// let nside = nside(depth) as u64;
/// let nested12 = get(depth);
/// assert_eq!(nside * nside - 1, nested12.hash(Degrees(12.5_f64, 89.99999_f64)));
/// ```
pub fn hash(&self, coords: impl LonLatT) -> u64 {
self.hash_impl(coords.lon(), coords.lat())
}
fn hash_impl(&self, lon: f64, lat: f64) -> u64 {
crate::check_lat(lat);
let (d0h, l_in_d0c, h_in_d0c) = Self::d0h_lh_in_d0c(lon, lat);
// Coords inside the base cell
// - ok to cast on u32 since small negative values due to numerical inaccuracies (like -1e-15), are rounded to 0
let i =
f64::from_bits((self.time_half_nside + (h_in_d0c + l_in_d0c).to_bits() as i64) as u64)
as u32;
let j =
f64::from_bits((self.time_half_nside + (h_in_d0c - l_in_d0c).to_bits() as i64) as u64)
as u32;
// - deals with numerical inaccuracies, rare so branch miss-prediction negligible
let i = if i == self.nside {
self.nside_minus_1
} else {
i
};
let j = if j == self.nside {
self.nside_minus_1
} else {
j
};
self.build_hash_from_parts(d0h, i, j)
}
pub fn try_hash(&self, coords: impl LonLatT) -> Result<u64, f64> {
let lat = coords.lat();
if lat.abs() > FRAC_PI_2 {
return Err(lat);
}
Ok(self.hash_impl(coords.lon(), lat))
}
#[inline]
fn d0h_lh_in_d0c(lon: f64, lat: f64) -> (u8, f64, f64) {
let (x_pm1, q) = Self::xpm1_and_q(lon);
if lat > TRANSITION_LATITUDE {
// North polar cap, Collignon projection.
// - set the origin to (PI/4, 0)
let sqrt_3_one_min_z = SQRT_6 * (0.5 * lat + FRAC_PI_4).cos();
let (x_proj, y_proj) = (x_pm1 * sqrt_3_one_min_z, 2.0 - sqrt_3_one_min_z);
let d0h = q;
(d0h, x_proj, y_proj)
} else if lat < -TRANSITION_LATITUDE {
// South polar cap, Collignon projection
// - set the origin to (PI/4, -PI/2)
let sqrt_3_one_min_z = SQRT_6 * (0.5 * lat - FRAC_PI_4).cos(); // cos(-x) = cos(x)
let (x_proj, y_proj) = (x_pm1 * sqrt_3_one_min_z, sqrt_3_one_min_z);
let d0h = q + 8;
(d0h, x_proj, y_proj)
} else {
// Equatorial region, Cylindrical equal area projection
// - set the origin to (PI/4, 0) if q = 2
// - set the origin to (PI/4, -PI/2) if q = 0
// - set the origin to (0, -TRANSITION_LAT) if q = 3
// - set the origin to (PI/2, -TRANSITION_LAT) if q = 1
// let zero_or_one = (x_cea as u8) & 1;
let y_pm1 = lat.sin() / TRANSITION_Z;
// Inequalities have been carefully chosen so that S->E and S->W axis are part of the cell,
// and not E->N and W->N
// Version with branch
// |\3/|
// .2X1.
// |/0\|
/*let q13 = (x_pm1 >= -y_pm1) as u8; /* 0\1 */ debug_assert!(q12 == 0 || q12 == 1);
let q23 = (x_pm1 <= y_pm1) as u8; /* 1/0 */ debug_assert!(q23 == 0 || q23 == 1);
match q13 | (q23 << 1) {
0 => ( q , x_pm1 , y_pm1 + 2.0),
1 => ((q + 5) & 7, x_pm1 - 1.0, y_pm1 + 1.0), // (q + 5) & 7 <=> (q + 1) | 4
2 => ( q + 4 , x_pm1 + 1.0, y_pm1 + 1.0),
3 => ( q + 8 , x_pm1 , y_pm1),
_ => unreachable!(),
}*/
// Branch free version
// |\2/|
// .3X1.
// |/0\|
let q01 = (x_pm1 > y_pm1) as u8; /* 0/1 */
debug_assert!(q01 == 0 || q01 == 1);
let q12 = (x_pm1 >= -y_pm1) as u8; /* 0\1 */
debug_assert!(q12 == 0 || q12 == 1);
let q1 = q01 & q12; /* = 1 if q1, 0 else */
debug_assert!(q1 == 0 || q1 == 1);
let q013 = q01 + (1 - q12); // = q01 + q03; /* 0/1 + 1\0 + */
// x: x_pm1 + 1 if q3 | x_pm1 - 1 if q1 | x_pm1 if q0 or q2
let x_proj = x_pm1 - ((q01 + q12) as i8 - 1) as f64;
// y: y_pm1 + 0 if q2 | y_pm1 + 1 if q1 or q3 | y_pm1 + 2 if q0
let y_proj = y_pm1 + q013 as f64;
// d0h: +8 if q0 | +4 if q3 | +5 if q1
let d0h = (q013 << 2) + ((q + q1) & 3);
(d0h, x_proj, y_proj)
}
}
/// Returns the cell number (hash value) associated with the given position on the unit sphere,
/// together with the offset `(dx, dy)` on the Euclidean plane of the projected position with
/// respect to the origin of the cell (South vertex).
/// # Inputs
/// - `lon`: longitude in radians, support reasonably large positive and negative values
/// producing accurate results with a naive range reduction like modulo 2*pi
/// (i.e. without having to resort on Cody-Waite or Payne Hanek range reduction).
/// - `lat`: latitude in radians, must be in `[-pi/2, pi/2]`
/// # Output
/// - the cell number (hash value) associated with the given position on the unit sphere,
/// in `[0, 12*nside^2[`
/// - `dx`: the positional offset $\in [0, 1[$ along the south-to-east axis
/// - `dy`: the positional offset $\in [0, 1[$ along the south-to-west axis
/// # Panics
/// If `lat` **not in** `[-pi/2, pi/2]`, this method panics.
/// # Examples
/// ```rust
/// use healpix::checked::nside;
/// use healpix::get;
///
/// let depth = 12_u8;
/// let nside = nside(depth) as u64;
/// let nested12 = get(depth);
/// let h_org = nside * nside - 1;
/// let [h_ra, h_dec] = nested12.center(h_org).as_f64s();
/// let (h, dx, dy) = nested12.hash_with_dxdy([h_ra, h_dec]);
/// assert_eq!(h_org, h);
/// // A precision of 1e-12 in a cell of depth 12 (side of ~51.5 arcsec)
/// // leads to an absolute precision of ~0.05 nanoarcsec
/// assert!((dx - 0.5).abs() < 1e-12_f64);
/// assert!((dy - 0.5).abs() < 1e-12_f64);
/// ```
pub fn hash_with_dxdy(&self, coords: impl LonLatT) -> (u64, f64, f64) {
self.hash_with_dxdy_impl(coords.lon(), coords.lat())
}
fn hash_with_dxdy_impl(&self, lon: f64, lat: f64) -> (u64, f64, f64) {
crate::check_lat(lat);
let (d0h, l_in_d0c, h_in_d0c) = Self::d0h_lh_in_d0c(lon, lat);
let x =
f64::from_bits((self.time_half_nside + (h_in_d0c + l_in_d0c).to_bits() as i64) as u64);
debug_assert!(
-1e-14 < x || x < self.nside as f64 * (1.0_f64 + 1e-14),
"x: {}, x_proj: {}; y_proj: {}",
&x,
&h_in_d0c,
&l_in_d0c
);
let y =
f64::from_bits((self.time_half_nside + (h_in_d0c - l_in_d0c).to_bits() as i64) as u64);
debug_assert!(
-1e-14 < y || y < self.nside as f64 * (1.0_f64 + 1e-14),
"y: {}, x_proj: {}; y_proj: {}",
&y,
&h_in_d0c,
&l_in_d0c
);
// - ok to cast on u32 since small negative values due to numerical inaccuracies (like -1e-15), are rounded to 0
let i = x as u32;
let j = y as u32;
// - deals with numerical inaccuracies, rare so branch miss-prediction negligible
let i = if i == self.nside {
self.nside_minus_1
} else {
i
};
let j = if j == self.nside {
self.nside_minus_1
} else {
j
};
(
self.build_hash_from_parts(d0h, i, j),
(x - (i as f64)),
(y - (j as f64)),
)
}
/// Compute the position on the unit sphere of the center (in the Euclidean projection plane)
/// of the cell associated to the given hash value.
///
/// # Input
/// - `hash`: the hash value of the cell we look for the unprojected center
///
/// # Output
/// - `(lon, lat)` in radians, the unprojected position (on the unit sphere) of the center of
/// the cell in the Euclidean plane
/// - `lon`, longitude in `[0, 2pi]` radians;
/// - `lat`, latitude in `[-pi/2, pi/2]` radians.
///
/// # Panics
/// If the given `hash` value is not in `[0, 12*nside^2[`, this method panics.
///
/// # Example
/// ```rust
/// use std::f64::consts::PI;
/// use healpix::geo::distance;
/// use healpix::{TRANSITION_LATITUDE, get};
///
/// let depth = 0u8;
/// let nested0 = get(depth);
///
/// assert!(distance([PI / 4f64, TRANSITION_LATITUDE], nested0.center(0u64)) < 1e-15);
/// ```
///
#[inline]
pub fn center(&self, hash: u64) -> LonLat {
let (x, y) = self.center_of_projected_cell(hash);
proj::unproj(x, y)
}
/// Compute the position on the unit sphere of the position '(dx, dy)' from the south vertex of
/// the HEALPix cell associated to the given hash value.
/// The x-axis is the South-East axis while the y-axis is the south-west axis.
///
/// # Input
/// - `hash`: the hash value of the cell in which are defined `dx` and `dy`
/// - `dx`: the positional offset $\in [0, 1[$ along the south-to-east axis
/// - `dy`: the positional offset $\in [0, 1[$ along the south-to-west axis
///
/// # Output
/// - `(lon, lat)` in radians, the unprojected position (on the unit sphere) of the given position
/// inside the given cell in the Euclidean plane
/// - `lon`, longitude in `[0, 2pi]` radians;
/// - `lat`, latitude in `[-pi/2, pi/2]` radians.
///
/// # Panics
/// This method panics if either:
/// - the given `hash` value is not in `[0, 12*nside^2[`,
/// - `dx` or `dy` is not $\in [0, 1[$
///
/// # Example
/// ```rust
/// use healpix::geo::distance;
/// use healpix::get;
///
/// let depth = 0u8;
/// let nested0 = get(depth);
///
/// assert!(distance(nested0.sph_coo(0, 0.5, 0.5), nested0.center(0)) < 1e-15);
/// ```
///
pub fn sph_coo(&self, hash: u64, dx: f64, dy: f64) -> [f64; 2] {
assert!((0.0..1.0).contains(&dx));
assert!((0.0..1.0).contains(&dy));
let (mut x, mut y) = self.center_of_projected_cell(hash);
x += (dx - dy) * self.one_over_nside;
y += (dx + dy - 1.0) * self.one_over_nside;
proj::unproj(x.rem_euclid(8.0), y).as_f64s()
}
#[inline]
pub(crate) fn build_hash_from_parts(&self, d0h: u8, i: u32, j: u32) -> u64 {
self.build_hash((d0h as u64) << self.twice_depth, i, j)
}
#[inline]
pub(crate) fn build_hash(&self, d0h_bits: u64, i: u32, j: u32) -> u64 {
use crate::zoc::ZOrderCurve;
debug_assert!(
i < self.nside && j < self.nside,
"nside: {}; i: {}, j: {}",
self.nside,
i,
j
);
d0h_bits | self.z_order_curve.ij2h(i, j)
}
/// Conveniency function simply calling the [hash](fn.to_uniq.html) method with this layer *depth*.
pub fn to_uniq(&self, hash: u64) -> u64 {
// depth already tested, so we call the unsafe method
unchecked::to_uniq(self.depth, hash)
}
/// Conveniency function simply calling the [hash](fn.to_uniq_ivoa.html) method with this layer *depth*.
pub fn to_uniq_ivoa(&self, hash: u64) -> u64 {
// depth already tested, so we call the unsafe method
unchecked::to_uniq_ivoa(self.depth, hash)
}
}