docs menu

Module Ymir.Direction

Directions in a frame.

A direction is a frame and float64 vectors [...; 3], each read as the ray it spans. A unit norm cannot be a property of the type: jvp tangents, grad cotangents and Nx.Ptree.map build directions of any norm. So the contract has two halves.

Constructors return unit vectors. Direction.lonlat and Direction.of_xyz give |xyz| = 1 within 3 ulp, and Direction.rotate preserves the norm within 3 ulp.
Readers read the ray. Direction.lon, Direction.lat, Direction.separation and Direction.position_angle read every finite row as the direction it points in. A positive multiple of a row reads the same up to the rounding of the multiple, and bit for bit when the multiple is a power of two that rounds no component, subnormal ones included. A NaN component gives NaN, and an infinite one raises Invalid_argument. Direction.rotate and Direction.xyz are linear.

So a sum of two directions read as a direction is their chordal midpoint.

Derivatives. Every function has its analytic derivative wherever it is differentiable. Where it has none it gives a stated value with derivative 0, never NaN:

Direction.lon at a pole (x = y = 0): 0.
Direction.lat at a pole: ±π/2; at a zero row: 0.
Direction.separation where a × b = 0, that is coincident or antipodal rows or a zero row: π if a · b < 0, else 0.
Direction.position_angle at the same points: 0. With respect to a at a pole, its derivative is 0.

At lat's poles and separation's 0 and π the function behaves like |x|, and 0 is a subgradient; the other points are discontinuities. Near a pole the derivative of lon grows as 1/ρ, ρ the distance from the axis.

Every angle in or out is a float64 quantity. Indexing and reshaping go through Direction.ptree:

ocaml

let near =
  Nx.Ptree.map (Direction.ptree ())
    (fun _ x -> Nx.take ~axis:0 indices x)
    stars

#type ('f, 'd) t

The type for directions in frame 'f: float64 vectors [...; 3], each read as the ray it spans.

val lonlat : 
  'f Frame.t ->
  lon:(float, Nx.float64_elt, 'd) Nx.t Quantity.t ->
  lat:(float, Nx.float64_elt, 'd) Nx.t Quantity.t ->
  ('f, 'd) t

lonlat f ~lon ~lat is the unit vectors (cos b cos l, cos b sin l, sin b) at longitude l and latitude b, their batch axes broadcast. Every finite angle names a direction: latitude 100° at longitude 0° is latitude 80° at longitude 180°. A NaN angle gives a NaN row.

Raises Invalid_argument if lon or lat is not an angle, or where an angle is infinite, as in "Direction.lonlat: lat at [3] is infinite".

#val of_xyz : 'f Frame.t -> (float, Nx.float64_elt, 'd) Nx.t -> ('f, 'd) t

of_xyz f v is v / |v| along v's last axis. A NaN row stays NaN. Its derivative is the projection onto the tangent plane.

Raises Invalid_argument if the last axis does not have 3 elements, or where a vector is infinite or zero, as in "Direction.of_xyz: the vector at [17] is zero and names no direction".

#val frame : ('f, 'd) t -> 'f Frame.t

frame d is d's frame.

#val xyz : ('f, 'd) t -> (float, Nx.float64_elt, 'd) Nx.t

xyz d is d's vectors, [...; 3]: unit vectors for values the constructors build.

#val lon : ('f, 'd) t -> (float, Nx.float64_elt, 'd) Nx.t Quantity.t

lon d is the longitude in [+0, 2π) radians, atan2 y x. It is 0 at a pole, so a function that needs a meridian at a pole takes the meridian of longitude 0.

#val lat : ('f, 'd) t -> (float, Nx.float64_elt, 'd) Nx.t Quantity.t

lat d is the latitude in [−π/2, π/2] radians, atan2 z ρ with ρ = √(x² + y²).

val rotate : 
  ('b Frame.fixed as 'g) Frame.t ->
  ('a Frame.fixed, 'd) t ->
  ('g, 'd) t

rotate g d is d in frame g: Frame.matrix (frame d) g times each vector. It is linear in d, within 1.5 · 2⁻⁵² · |v| of the exact product per component, and rounds the same eagerly and compiled.

#val convert : 'd Earth.t -> 'g Frame.t -> ('f, 'd) t -> ('g, 'd) t

convert e g d is d in frame g at e's instants, their batch axes broadcast. Within a family it is Frame.matrix_at times each vector. From a fixed frame to a geocentric one it reads d as the barycentric direction to a source at infinite distance and gives the geocentric one as iauAtciq does without proper motion or parallax: deflected by the Sun (iauLdsun), aberrated (iauAb), then turned onto g's axes. Within about 5′ of the Sun's centre the deflection falls linearly to 0 at the centre; the planets' is left out. The other way it inverts that model: convert e f (convert e g d) is d within 1e-15.

Between the families it reads rows as rays and returns unit vectors, within 2e-15 of SOFA's; a zero row gives NaN.

Raises Invalid_argument where a row has an infinite component, as in "Direction.convert: the vector at [3] is infinite".

val separation : 
  ('f, 'd) t ->
  ('f, 'd) t ->
  (float, Nx.float64_elt, 'd) Nx.t Quantity.t

separation a b is the angle between the rays, in [0, π] radians. Leading axes broadcast. It is atan2 |a × b| (a · b), each component of a × b within 2 ulp through a fused multiply-add and lifted by a power of two where its square would underflow, so it is within 1.3e-15 relative at any separation above float64's least normal, whatever the rows' norms.

val position_angle : 
  ('f, 'd) t ->
  ('f, 'd) t ->
  (float, Nx.float64_elt, 'd) Nx.t Quantity.t

position_angle a b is b's bearing from a, east of north, in [+0, 2π) radians. Leading axes broadcast. At a pole, north is the limit along longitude 0: at the north pole it points toward longitude 180° and east toward 90°. It reads a × b as separation does, east |a| (a × b)_z and north ((a × b) × a)_z, so it is within 2.1e-15 radians at any separation above float64's least normal, whatever the rows' norms.

#val walk : ('a, 'b, 'd, 'e) Nx.Ptree.Walk.cursor -> ('f, 'd) t -> ('f, 'e) t

walk c d walks d: the frame as Frame.walk reports it, then the vectors as a tensor of a fixed type, at c's path. Nx.Ptree.cast keeps them float64, Nx.Ptree.place moves them, and Nx.Ptree.map and rune's transformations treat them as any tensor.

#val ptree : unit -> (('f, 'd) t, 'd) Nx.Ptree.t

ptree () is Nx.Ptree.v walk, the structure of directions.