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
use std::f64::consts::FRAC_PI_2;
use super::proj::*;
use crate::{xy_geom::ellipse::*, Customf64};
/// TODO: Look at the Java code for a more robust approach!!
#[derive(Debug)]
pub struct EllipticalCone {
center: ProjSIN, // the center of the projection is the center of the cone
ellipse: Ellipse,
// Ellipse parameters
a: f64, // semi-major axis, in radians
b: f64, // semi-minor axis, in radians
theta_sin_cos: (f64, f64), // (pi/2 - position angle).sin_cos()
}
impl EllipticalCone {
/// #Inputs
/// - `lon` longitude of the ellipse center, in radians
/// - `lat` latitude of the ellipse center, in radians
/// - `a` semi-major axis of the ellipse, in radians
/// - `b` semi-minor axis of the ellipse, in radians
/// - `pa` position angle (east of north), in radians
pub fn new(lon: f64, lat: f64, a: f64, b: f64, pa: f64) -> EllipticalCone {
let center = ProjSIN::new(lon, lat);
let theta_sin_cos = (FRAC_PI_2 - pa).sin_cos();
let ellipse = Ellipse::from_oriented(a.sin(), b.sin(), theta_sin_cos);
EllipticalCone {
center,
ellipse,
a,
b,
theta_sin_cos,
}
}
/// Returns `true` if the given point on the unit sphere is inside the elliptical cone.
/// # Inputs
/// - `lon` longitude of the point we want to test, in radians
/// - `lat` latitude of the point we want to test, in radians
pub fn contains(&self, lon: f64, lat: f64) -> bool {
match self.center.proj(lon, lat) {
Some((x, y)) => self.ellipse.contains(x, y),
None => false,
}
}
/// Returns `true` if the given cone overlap the elliptical cone.
/// # Inputs
/// - `lon` longitude of the center of the cone, in radians
/// - `lat` latitude of the center of the cone, in radians
/// - `radius` cone radius, in radians
///
/// # Info
/// To properly address the question, one should try to see if the following equations system has
/// at least one solution:
/// ```math
/// \frac{x^2}{\tan^2a} + \frac{y^2}{\tan^2b} = z^2
/// (x - x_0)^2 + (y - y_0)^2 + (z - z_0)^2 = 4\sin^2\frac{\theta}{2}
/// x^2 + y^2 + z^2 = 0
/// ```
/// with `a` and `b` which are angles and
/// ```math
/// x = \cos\delta\cos\alpha
/// y = \cos\delta\sin\alpha
/// z = \sin\delta
/// ```
/// This system can be replaced by (the second equation being the angular distance formula):
/// ```math
/// \frac{\cos^2\alpha}{\tan^2a} + \frac{\sin^2\alpha}{\tan^2b} = \tan^2\delta
/// \sin\delta\sin\delta_0 + \cos\delta\cos\delta_0\cos(\alpha-\alpha_0) = \cos\theta
/// ```
/// - The first equation is the elliptical cone equation
/// - The second is the classical cone equation
/// - $(x_ 0, y_ 0, z_ 0)$ are obtained by rotation of the original cone center such that the
/// elliptical cone is center is $(x = 0, y = 0, z = 1)$, and the x-axis is along the semi-major
/// axis.
///
/// # Info 2
/// Contrary to the above Info, another (the best?) approach is to compute the coordinates of the
/// two ellipse foci F0 and F1. The sum of the distances f0 and f1 to both foci is constant with:
/// ```math
/// f0 + f1 = 2a
/// ```
///
pub fn overlap_cone(&self, lon: f64, lat: f64, radius: f64) -> bool {
assert!(radius > 0.0);
let ((x, y), ang_dist, _top_hemisphere) = self.center.forced_proj_and_distance(lon, lat);
if self.a + radius < ang_dist {
// Quick rejection test
return false;
}
// The projection of a cone on a plane is an ellipse
let proj_d_min = (ang_dist - radius).sin();
let proj_d_max = (ang_dist + radius).sin();
let proj_a = radius.sin(); // projected ellipse semi-major axis
let proj_b = 0.5 * (proj_d_max - proj_d_min).abs(); // projected ellipse semi-minor axis
let one_over_norm = 1.0 / (x.pow2() + y.pow2()).sqrt(); // distance from (0, 0) to the cone projected center (!= ellipse center)
if !one_over_norm.is_finite() {
// Special case: the center of the HEALPix cell is the center of the ellipse
return radius <= self.b;
}
// phi = angle of the (x, y) position in the euclidean plane
let cos_phi = x * one_over_norm;
let sin_phi = y * one_over_norm;
let sin_cos_angle = if sin_phi >= 0.0 {
(-cos_phi, sin_phi)
} else {
(cos_phi, -sin_phi)
};
let proj_ell = Ellipse::from_oriented(proj_a, proj_b, sin_cos_angle);
let proj_ell_dist = 0.5 * (proj_d_max + proj_d_min); // distance from (0, 0) to ellipse center
let proj_ell_x = proj_ell_dist * cos_phi;
let proj_ell_y = proj_ell_dist * sin_phi;
/*if top_hemisphere {
self.ellipse.overlap(proj_ell_x, proj_ell_y, &proj_ell)
} else {
// The projection of the center of the cone is in the back hemisphere
// so the area to be considered looks like a crescent moon.
// So far we apply the same algo, knowing that at deeper depth cells will have less changes
// to spuriously overlap
self.ellipse.overlap(proj_ell_x, proj_ell_y, &proj_ell)
}*/
self.ellipse.overlap(proj_ell_x, proj_ell_y, &proj_ell)
}
/// Returns `true` if the given cone is fully inside the elliptical cone.
/// # Inputs
/// - `lon` longitude of the center of the cone, in radians
/// - `lat` latitude of the center of the cone, in radians
/// - `radius` cone radius, in radians
pub fn contains_cone(&self, lon: f64, lat: f64, radius: f64) -> bool {
if radius >= self.b {
false
} else {
match self.center.proj(lon, lat) {
Some((x, y)) => {
let eucl_a = (self.a - radius).sin();
let eucl_b = (self.b - radius).sin();
// WARNING: not sure this is a 100% reliable for large distances!
// A 100% reliable solution would be to compute the distance to both foci:
// (f0 + f1)/2 <= a - r
Ellipse::from_oriented(eucl_a, eucl_b, self.theta_sin_cos).contains(x, y)
}
None => false,
}
}
}
#[allow(dead_code)]
pub fn path_along_edge(
lon: f64,
lat: f64,
a: f64,
b: f64,
theta: f64,
half_num_points: usize,
) -> Box<[(f64, f64)]> {
let center = ProjSIN::new(lon, lat);
let mut ar = Ellipse::path_along_edge(a.sin(), b.sin(), FRAC_PI_2 - theta, half_num_points);
for coo in ar.iter_mut() {
*coo = match center.unproj(coo.0, coo.1) {
Some(lonlat) => lonlat,
None => (f64::NAN, f64::NAN),
};
}
ar
}
}