cdshealpix 0.9.1

Rust implementation of the HEALPix tesselation.
Documentation 
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

use std::f64::consts::FRAC_PI_2;

use super::super::Customf64;
use super::super::TWICE_PI;

/// Represents a spherical projection, i.e. the projection of spherical coordinates
/// (on the unit sphere) on a two dimensional plane.
pub trait Proj: Default {
  /// Set the center of the projection
  /// # Inputs
  /// - `lon` longitude of the centre of the projection, in `[0, 2pi]` radians
  /// - `lat` latitude of the centre of the projection, in `[-pi/2, pi/2]` radians
  /// # info
  ///  If the longitude or the latitude is out of bounds, a time consuming normalization may be performed
  fn set_center(&mut self, lon: f64, lat: f64);

  /// Computes the projected coordinates of the given position on the unit sphere
  /// # Inputs
  /// - `lon` longitude of the coordinate we want to project
  /// - `lat` latitude of the coordinate we want to project
  /// # Output
  /// - `(x, y)` the coordinates in the projection plane (if they exists)
  fn proj(&self, lon: f64, lat: f64) -> Option<(f64, f64)>;

  /// Computes the position in the unit sphere of the given coordinates on the projection plane
  /// # Inputs
  /// - `x` coordinates along the x-axis on the projection plane
  /// - `y` coordinates along the y-axis on the projection plane
  /// # Output
  /// - `(lon, lat)` the position on the unit sphere (in radians), if the input coordinates `(x, y)`
  ///   are valid
  fn unproj(&self, x: f64, y: f64) -> Option<(f64, f64)>;
}

/// Normalize the given longitude and latitude so that their values are in `[0, 2pi]` and  `[-pi/2, pi/2]`.
/// WARNING: the normalization (without a priori) is time consuming and should be avoided as much as possible
fn normalize_lonlat(lon: &mut f64, lat: &mut f64) {
  if *lon < 0.0 || TWICE_PI <= *lon || *lat < -FRAC_PI_2 || FRAC_PI_2 < *lat {
    let (sin_l, cos_l) = (*lon).sin_cos();
    let (sin_b, cos_b) = (*lat).sin_cos();
    let x = cos_b * cos_l;
    let y = cos_b * sin_l;
    let z = sin_b;
    *lon = y.atan2(x);
    if *lon < 0.0_f64 {
      *lon += TWICE_PI;
    }
    *lat = z.atan2((x.pow2() + y.pow2()).sqrt());
  }
}

/// Orthographic projection
#[derive(Debug)]
pub struct ProjSIN {
  center_lon: f64,
  center_lat: f64,
  // derived quantities
  pub cos_center_lat: f64,
  pub sin_center_lat: f64,
}

impl ProjSIN {
  /// # Inputs
  /// - `lon` longitude of the centre of the projection, in `[0, 2pi]` radians
  /// - `lat` latitude of the centre of the projection, in `[-pi/2, pi/2]` radians
  /// # info
  ///  If the longitude or the latitude is out of bounds, a time consuming normalization is performed
  pub fn new(lon: f64, lat: f64) -> ProjSIN {
    let mut proj: ProjSIN = Default::default();
    proj.set_center(lon, lat);
    proj
  }

  /// Returns the (x, y) projected position, even if the source is in the opposite hemisphere.
  /// Also returns the angular distance (computed for large values, not using the Haversine formula).
  /// The boolean tells is the projected point is in the visible hemisphere.
  pub fn forced_proj_and_distance(&self, lon: f64, lat: f64) -> ((f64, f64), f64, bool) {
    let (sin_lat, cos_lat) = lat.sin_cos();
    let dlon = lon - self.center_lon;
    let (sin_dlon, cos_dlon) = dlon.sin_cos();
    (
      (
        cos_lat * sin_dlon,
        self.cos_center_lat * sin_lat - self.sin_center_lat * cos_lat * cos_dlon,
      ),
      (self.sin_center_lat * sin_lat + self.cos_center_lat * cos_lat * cos_dlon).acos(),
      self.sin_center_lat * sin_lat + self.cos_center_lat * cos_lat * cos_dlon > 0.0,
    )
  }
}

impl Default for ProjSIN {
  fn default() -> Self {
    ProjSIN {
      center_lon: 0.0,
      center_lat: 0.0,
      // derived quantities
      cos_center_lat: 1.0,
      sin_center_lat: 0.0,
    }
  }
}

impl Proj for ProjSIN {
  fn set_center(&mut self, lon: f64, lat: f64) {
    // I put assert tests here because the normalization is time consuming and I would like to
    // avoid it as much as possible
    //    debug_assert!(0.0 <= lon && lon < TWICE_PI);
    //    debug_assert!(-FRAC_PI_2 <= lat && lat < FRAC_PI_2);
    self.center_lon = lon;
    self.center_lat = lat;
    normalize_lonlat(&mut self.center_lon, &mut self.center_lat);
    let (sin_lat, cos_lat) = self.center_lat.sin_cos();
    self.cos_center_lat = cos_lat;
    self.sin_center_lat = sin_lat;
  }

  fn proj(&self, lon: f64, lat: f64) -> Option<(f64, f64)> {
    let (sin_lat, cos_lat) = lat.sin_cos();
    let dlon = lon - self.center_lon;
    /*
     let mut dlon = lon - self.center_lon;
    if dlon < -PI {
      dlon += TWICE_PI;
    } else if dlon > PI {
      dlon -= TWICE_PI;
    }*/
    let (sin_dlon, cos_dlon) = dlon.sin_cos();
    if self.sin_center_lat * sin_lat + self.cos_center_lat * cos_lat * cos_dlon > 0.0 {
      Some((
        cos_lat * sin_dlon,
        self.cos_center_lat * sin_lat - self.sin_center_lat * cos_lat * cos_dlon,
      ))
    } else {
      None
    }
  }

  fn unproj(&self, x: f64, y: f64) -> Option<(f64, f64)> {
    let rho2 = x.pow2() + y.pow2();
    if rho2 < 1.0 {
      let z = (1.0 - rho2).sqrt(); // = cos(arcsin(rho)) = rho / tan(arcsin(rho))
      let lat = (self.sin_center_lat * z + y * self.cos_center_lat).asin();
      let lon = self.cos_center_lat * z - y * self.sin_center_lat;
      let lon = self.center_lon + x.atan2(lon);
      Some((
        if lon < 0.0 {
          lon + TWICE_PI
        } else if lon > TWICE_PI {
          lon - TWICE_PI
        } else {
          lon
        },
        lat,
      ))
    } else {
      None
    }
  }
}