Vector and Matrix Operations

a2af(ndp::Integer, angle::Real)

Decompose radians into degrees, arcminutes, arcseconds, and fraction.

Input

• ndp – number of useful digits
• angle – angle in radians

Output

• dms – angle in sign, degrees, minutes, seconds, and fraction

Examples

julia> a2af(4, 2.345)
(sign = '+', degree = 134, minute = 21, second = 30, fraction = 9706)

Note

1. The argument ndp is interpreted as follows:

   ndp resolution : ...0000 00 00 -7 1000 00 00 -6 100 00 00 -5 10 00 00 -4 1 00 00 -3 0 10 00 -2 0 01 00 -1 0 00 10 0 0 00 01 1 0 00 00.1 2 0 00 00.01 3 0 00 00.001 : 0 00 00.000...

2. The largest positive useful value for ndp is determined by the size of angle, the format of Float64 on the target platform, and the risk of overflowing dms[3]. On a typical platform, for angle up to 2pi, the available floating-point precision might correspond to ndp=12. However, the practical limit is typically ndp=9, set by the capacity of a 32-bit int, or ndp=4 if int is only 16 bits.

3. The absolute value of angle may exceed 2pi. In cases where it does not, it is up to the caller to test for and handle the case where angle is very nearly 2pi and rounds up to 360 degrees, by testing for dms[0]=360 and setting dms[0-3] to zero.

a2tf(ndp::Integer, angle::Real)

Decompose radians into hours, minutes, seconds, and fraction.

Input

• ndp – number of useful digits
• angle – angle in radians

Output

• hms – angle in sign, hour, minutes, seconds, and fraction

Note

1. The argument ndp is interpreted as follows:

   ndp resolution : ...0000 00 00 -7 1000 00 00 -6 100 00 00 -5 10 00 00 -4 1 00 00 -3 0 10 00 -2 0 01 00 -1 0 00 10 0 0 00 01 1 0 00 00.1 2 0 00 00.01 3 0 00 00.001 : 0 00 00.000...

2. The largest positive useful value for ndp is determined by the size of angle, the format of Float64 on the target platform, and the risk of overflowing hms[3]. On a typical platform, for angle up to 2pi, the available floating-point precision might correspond to ndp=12. However, the practical limit is typically ndp=9, set by the capacity of a 32-bit int, or ndp=4 if int is only 16 bits.

3. The absolute value of angle may exceed 2pi. In cases where it does not, it is up to the caller to test for and handle the case where angle is very nearly 2pi and rounds up to 24 hours, by testing for ihmsf[0]=24 and setting ihmsf[0-3] to zero.

af2a(sign::Char, degree::Integer, minute::Integer, second::Real)

Convert degrees, arcminutes, arcseconds to radians.

Input

• sign – sign of arc
• degree – degrees of arc
• minute – minutes of arc
• second – seconds of arc

Output

• angle – angle in radians

Examples

julia> af2a('-', 45, 13, 27.2)
-0.7893115794313645

Note

1. Unlike the SOFA C function, which computes a result even when the range checks fail (flagging a warning status), this function throws an AssertionError if any argument is out of range.

anp(angle::Real)

Normalize angle into the range 0 <= a < 2pi.

Input

• angle – angle in radians

Output

• angle – angle in radians in range 0-2pi

Examples

julia> anp(-0.1)
6.183185307179587

anpm(angle::Real)

Normalize angle into the range -pi <= a < +pi

Input

• angle – angle in radians

Output

• angle – angle in radians in range +/-pi

c2s(pos::AbstractVector{<:Real})

P-vector to spherical coordinates.

Input

• p – p-vector

Output

• θ – longitude angle (radians)
• ϕ – latitude angle (radians)

Examples

julia> c2s([100.0, -50.0, 25.0])
(θ = -0.4636476090008061, ϕ = 0.21998797739545944)

Note

1. The vector p can have any magnitude; only its direction is used.

2. If p is null, zero θ and ϕ are returned.

3. At either pole, zero θ is returned.

cp(p::AbstractVector{<:Real})

Copy a p-vector.

Input

• p – p-vector to be copied

Output

• c – copy

cpv(pv::AbstractVector{<:AbstractVector{<:Real}})

Copy a position/velocity vector.

Input

• pv – pv-vector to be copied

Output

• c – copy

cr(r::AbstractMatrix{<:Real})

Copy an r-matrix.

Input

• r – r-matrix to be copied

Output

• c – copy

d2tf(ndp::Integer, day::Real)

Decompose days to sign, hours, minutes, seconds, fraction.

Input

• ndp – number of useful digits
• day – interval in days

Output

• hms – hms in sign, hours, minutes, seconds, and fraction

Note

1. The argument ndp is interpreted as follows:

   ndp resolution : ...0000 00 00 -7 1000 00 00 -6 100 00 00 -5 10 00 00 -4 1 00 00 -3 0 10 00 -2 0 01 00 -1 0 00 10 0 0 00 01 1 0 00 00.1 2 0 00 00.01 3 0 00 00.001 : 0 00 00.000...

2. The largest positive useful value for ndp is determined by the size of days, the format of Float64 on the target platform, and the risk of overflowing hms[3]. On a typical platform, for days up to 1.0, the available floating-point precision might correspond to ndp=12. However, the practical limit is typically ndp=9, set by the capacity of a 32-bit int, or ndp=4 if int is only 16 bits.

3. The absolute value of days may exceed 1.0. In cases where it does not, it is up to the caller to test for and handle the case where days is very nearly 1.0 and rounds up to 24 hours, by testing for ihmsf[0]=24 and setting ihmsf[0-3] to zero.

ir()

Initialize an r-matrix to the identity matrix.

Output

• r – r-matrix

p2pv(p::AbstractVector{<:Real})

Extend a p-vector to a pv-vector by appending a zero velocity.

Input

• p – p-vector

Output

• pv – pv-vector

p2s(pos::AbstractVector{<:Real})

P-vector to spherical polar coordinates.

Input

• p – p-vector

Output

• θ – longitude angle (radians)
• ϕ – latitude angle (radians)
• r – radial distance

Note

1. If P is null, zero θ, ϕ and r are returned.

2. At either pole, zero θ is returned.

pap(a::AbstractVector{<:Real}, b::AbstractVector{<:Real})

Position-angle from two p-vectors.

Input

• a – direction of reference point
• b – direction of point whose PA is required

Output

• θ – position angle of b with respect to a (radians)

Note

1. The result is the position angle, in radians, of direction b with respect to direction a. It is in the range -pi to +pi. The sense is such that if b is a small distance "north" of a the position angle is approximately zero, and if b is a small distance "east" of a the position angle is approximately +pi/2.

2. The vectors a and b need not be of unit length.

3. Zero is returned if the two directions are the same or if either vector is null.

4. If vector a is at a pole, the result is ill-defined.

pas(λa::Real, ϕa::Real, λb::Real, ϕb::Real)

Position-angle from spherical coordinates.

Input

• λa – longitude of point A (e.g. RA) in radians
• ϕa – latitude of point A (e.g. Dec) in radians
• λb – longitude of point B
• ϕb – latitude of point B

Output

• θ – position angle of B with respect to A

Note

1. The result is the bearing (position angle), in radians, of point B with respect to point A. It is in the range -pi to +pi. The sense is such that if B is a small distance "east" of point A, the bearing is approximately +pi/2.

2. Zero is returned if the two points are coincident.

pdp(a::AbstractVector{<:Real}, b::AbstractVector{<:Real})

P-vector inner (=scalar=dot) product.

Input

• a – first p-vector
• b – second p-vector

Output

• r – a . b

Examples

julia> pdp([2.0, 2.0, 3.0], [1.0, 3.0, 4.0])
20.0

pm(p::AbstractVector{<:Real}) = norm(p)

Modulus of p-vector.

Input

• p – p-vector

Output

• r – modulus

pmp(a::AbstractVector{<:Real}, b::AbstractVector{<:Real})

P-vector subtraction.

Input

• a – first p-vector
• b – second p-vector

Output

• amb – a - b

pn(p::AbstractVector{<:Real})

Convert a p-vector into modulus and unit vector.

Input

• p – p-vector

Output

• r – modulus
• u – unit vector

Note

1. If p is null, the result is null. Otherwise the result is a unit vector.

ppp(a::AbstractVector{<:Real}, b::AbstractVector{<:Real})

P-vector addition.

Input

• a – first p-vector
• b – second p-vector

Output

• apb – a + b

ppsp(a::AbstractVector{<:Real}, s::Real, b::AbstractVector{<:Real})

P-vector plus scaled p-vector.

Input

• a – first p-vector
• s – scalar (multiplier for b)
• b – second p-vector

Output

• apsb – a + s*b

pv2p(pv::AbstractVector{<:AbstractVector{<:Real}})

Discard velocity component of a pv-vector.

Input

• pv – pv-vector

Output

• p – p-vector

pv2s(pv::AbstractVector{<:AbstractVector{<:Real}})

Convert position/velocity from Cartesian to spherical coordinates.

Input

• posvel – position-velocity-vector

Output

• θ – longitude angle (radians)
• ϕ – latitude angle (radians)
• r – radial distance
• dθ – rate of change of θ
• dϕ – rate of change of ϕ
• dr – rate of change of r

Note

1. If the position part of pv is null, theta, ϕ, td and pd are indeterminate. This is handled by extrapolating the position through unit time by using the velocity part of pv. This moves the origin without changing the direction of the velocity component. If the position and velocity components of pv are both null, zeroes are returned for all six results.

2. If the position is a pole, theta, td and pd are indeterminate. In such cases zeroes are returned for all three.

pvdpv(a::AbstractVector{<:AbstractVector{<:Real}}, b::AbstractVector{<:AbstractVector{<:Real}})

Inner (=scalar=dot) product of two pv-vectors.

Input

• a – first pv-vector
• b – second pv-vector

Output

• adb – a . b (see note)

Note

1. If the position and velocity components of the two pv-vectors are ( ap, av ) and ( bp, bv ), the result, a . b, is the pair of numbers ( ap . bp , ap . bv + av . bp ). The two numbers are the dot-product of the two p-vectors and its derivative.

pvm(pv::AbstractVector{<:AbstractVector{<:Real}})

Modulus of pv-vector.

Input

• pv – pv-vector

Output

• r – modulus of position component
• s – modulus of velocity component

pvmpv(a::AbstractVector{<:AbstractVector{<:Real}}, b::AbstractVector{<:AbstractVector{<:Real}})

Subtract one pv-vector from another.

Input

• a – first pv-vector
• b – second pv-vector

Output

• amb – a - b

pvppv(a::AbstractVector{<:AbstractVector{<:Real}}, b::AbstractVector{<:AbstractVector{<:Real}})

Add one pv-vector to another.

Input

• a – first pv-vector
• b – second pv-vector

Output

• apb – a + b

pvu(dt::Real, pv::AbstractVector{<:AbstractVector{<:Real}})

Update a pv-vector.

Input

• dt – time interval
• pv – pv-vector

Output

• upv – p updated, v unchanged

Note

1. "Update" means "refer the position component of the vector to a new date dt time units from the existing date".

2. The time units of dt must match those of the velocity.

pvup(dt::Real, pv::AbstractVector{<:AbstractVector{<:Real}})

Update a pv-vector, discarding the velocity component.

Input

• dt – time interval
• pv – pv-vector

Output

• p – p-vector

Note

1. "Update" means "refer the position component of the vector to a new date dt time units from the existing date".

2. The time units of dt must match those of the velocity.

pvxpv(a::AbstractVector{<:AbstractVector{<:Real}}, b::AbstractVector{<:AbstractVector{<:Real}})

Outer (=vector=cross) product of two pv-vectors.

Input

• a – first pv-vector
• b – second pv-vector

Output

• axb – a x b

Note

1. If the position and velocity components of the two pv-vectors are ( ap, av ) and ( bp, bv ), the result, a x b, is the pair of vectors ( ap x bp, ap x bv + av x bp ). The two vectors are the cross-product of the two p-vectors and its derivative.

pxp(a::AbstractVector{<:Real}, b::AbstractVector{<:Real})

P-vector outer (=vector=cross) product.

Input

• a – first p-vector
• b – second p-vector

Output

• axb – a x b

rm2v(r::AbstractMatrix{<:Real})

Express an r-matrix as an r-vector.

Input

• r – rotation matrix

Output

• w – rotation vector (Note 1)

Note

1. A rotation matrix describes a rotation through some angle about some arbitrary axis called the Euler axis. The "rotation vector" returned by this function has the same direction as the Euler axis, and its magnitude is the angle in radians. (The magnitude and direction can be separated by means of the function pn.)

2. If r is null, so is the result. If r is not a rotation matrix the result is undefined; r must be proper (i.e. have a positive determinant) and real orthogonal (inverse = transpose).

3. The reference frame rotates clockwise as seen looking along the rotation vector from the origin.

rv2m(w::AbstractVector{<:Real})

Form the r-matrix corresponding to a given r-vector.

Input

• w – rotation vector (Note 1)

Output

• r – rotation matrix

Note

1. A rotation matrix describes a rotation through some angle about some arbitrary axis called the Euler axis. The "rotation vector" supplied to this function has the same direction as the Euler axis, and its magnitude is the angle in radians.

2. If w is null, the identity matrix is returned.

3. The reference frame rotates clockwise as seen looking along the rotation vector from the origin.

rx(ϕ::Real, r::AbstractMatrix{<:Real})

Rotate an r-matrix about the x-axis.

Input

• ϕ – angle (radians)
• r – r-matrix

Output

• r – r-matrix, rotated

Note

1. Calling this function with positive ϕ incorporates in the supplied r-matrix r an additional rotation, about the x-axis, anticlockwise as seen looking towards the origin from positive x.

2. The additional rotation can be represented by this matrix:

   (  1        0            0      )
   (                               )
   (  0   + cos(ϕ)   + sin(ϕ)  )
   (                               )
   (  0   - sin(ϕ)   + cos(ϕ)  )

rxp(r::AbstractMatrix{<:Real}, p::AbstractVector{<:Real})

Multiply a p-vector by an r-matrix.

Input

• r – r-matrix
• p – p-vector

Output

• rp – r * p

Note

1. It is permissible for p and rp to be the same array.

rxpv(r::AbstractMatrix{<:Real}, pv::AbstractVector{<:AbstractVector{<:Real}})

Multiply a pv-vector by an r-matrix.

Input

• r – r-matrix
• pv – pv-vector

Output

• rpv – r * pv

Note

1. The algorithm is for the simple case where the r-matrix r is not a function of time. The case where r is a function of time leads to an additional velocity component equal to the product of the derivative of r and the position vector.

2. It is permissible for pv and rpv to be the same array.

rxr(a::AbstractMatrix{<:Real}, b::AbstractMatrix{<:Real})

Multiply two r-matrices.

Input

• a – first r-matrix
• b – second r-matrix

Output

• atb – a * b

Note

1. It is permissible to re-use the same array for any of the arguments.

ry(θ::Real, r::AbstractMatrix{<:Real})

Rotate an r-matrix about the y-axis.

Input

• θ – angle (radians)
• r – r-matrix

Output

• r – r-matrix, rotated

Note

1. Calling this function with positive theta incorporates in the supplied r-matrix r an additional rotation, about the y-axis, anticlockwise as seen looking towards the origin from positive y.

2. The additional rotation can be represented by this matrix:

   (  + cos(θ)     0      - sin(θ)  )
   (                                        )
   (       0           1           0        )
   (                                        )
   (  + sin(θ)     0      + cos(θ)  )

rz(ψ::Real, r::AbstractMatrix{<:Real})

Rotate an r-matrix about the z-axis.

Input

• ψ – angle (radians)
• r – r-matrix

Output

• r – r-matrix, rotated

Note

1. Calling this function with positive ψ incorporates in the supplied r-matrix r an additional rotation, about the z-axis, anticlockwise as seen looking towards the origin from positive z.

2. The additional rotation can be represented by this matrix:

   (  + cos(ψ)   + sin(ψ)     0  )
   (                              )
   (  - sin(ψ)   + cos(ψ)     0  )
   (                              )
   (       0            0      1  )

s2c(θ::Real, ϕ::Real)

Convert spherical coordinates to Cartesian.

Input

• θ – longitude angle (radians)
• ϕ – latitude angle (radians)

Output

• c – direction cosines

Examples

julia> s2c(3.0123, -0.999)
3-element StaticArraysCore.MVector{3, Float64} with indices SOneTo(3):
 -0.5366267667260525
  0.06977111097651444
 -0.8409302618566215

s2p(θ::Real, ϕ::Real, r::Real)

Convert spherical polar coordinates to p-vector.

Input

• θ – longitude angle (radians)
• ϕ – latitude angle (radians)
• r – radial distance

Output

• p – Cartesian coordinates

s2pv(θ::Real, ϕ::Real, r::Real, dθ::Real, dϕ::Real,
     dr::Real)

Convert position/velocity from spherical to Cartesian coordinates.

Input

• θ – longitude angle (radians)
• ϕ – latitude angle (radians)
• r – radial distance
• dθ – rate of change of θ
• dϕ – rate of change of ϕ
• dr – rate of change of r

Output

• pv – pv-vector

s2xpv(s1::Real, s2::Real, pv::AbstractVector{<:AbstractVector{<:Real}})

Multiply a pv-vector by two scalars.

Input

• s1 – scalar to multiply position component by
• s2 – scalar to multiply velocity component by
• pv – pv-vector

Output

• spv – pv-vector: p scaled by s1, v scaled by s2

sepp(a::AbstractVector{<:Real}, b::AbstractVector{<:Real})

Angular separation between two p-vectors.

Input

• a – first p-vector (not necessarily unit length)
• b – second p-vector (not necessarily unit length)

Output

• θ – angular separation (radians, always positive)

Note

1. If either vector is null, a zero result is returned.

2. The angular separation is most simply formulated in terms of scalar product. However, this gives poor accuracy for angles near zero and pi. The present algorithm uses both cross product and dot product, to deliver full accuracy whatever the size of the angle.

seps(λa::Real, ϕa::Real, λb::Real, ϕb::Real)

Angular separation between two sets of spherical coordinates.

Input

• λa – first longitude (radians)
• ϕa – first latitude (radians)
• λb – second longitude (radians)
• ϕb – second latitude (radians)

Output

• θ – angular separation (radians)

Examples

julia> seps(1.0, 0.1, 0.2, -3.0)
2.346722016996999

sxp(s::Real, p::AbstractVector{<:Real})

Multiply a p-vector by a scalar.

Input

• s – scalar
• p – p-vector

Output

• sp – s * p

sxpv(s::Real, pv::AbstractVector{<:AbstractVector{<:Real}})

Multiply a pv-vector by a scalar.

Input

• s – scalar
• pv – pv-vector

Output

• spv – s * pv

tf2a(sign::Char, hour::Integer, minute::Integer, second::Real)

Convert hours, minutes, seconds to radians.

Input

• sign – sign: '-' = negative, otherwise positive
• hour – hours
• minute – minutes
• second – seconds

Output

• angle – angle in radians

Note

1. Unlike the SOFA C function, which computes a result even when the range checks fail (flagging a warning status), this function throws an AssertionError if any argument is out of range.

tf2d(sign::Char, hour::Integer, minute::Integer, second::Real)

Convert hours, minutes, seconds to days.

Input

• sign – sign: '-' = negative, otherwise positive
• hour – hours
• minute – minutes
• second – seconds

Output

• day – interval in days

Note

1. Unlike the SOFA C function, which computes a result even when the range checks fail (flagging a warning status), this function throws an AssertionError if any argument is out of range.

tr(r::AbstractMatrix{<:Real})

Transpose an r-matrix.

Input

• r – r-matrix

Output

• rt – transpose

Note

1. It is permissible for r and rt to be the same array.

2. The name tr follows the SOFA C library (iauTr) but collides with LinearAlgebra.tr (the matrix trace). When both SOFA and LinearAlgebra are loaded, an unqualified tr is ambiguous and must be written as SOFA.tr or LinearAlgebra.tr.

trxp(r::AbstractMatrix{<:Real}, p::AbstractVector{<:Real})

Multiply a p-vector by the transpose of an r-matrix.

Input

• r – r-matrix
• p – p-vector

Output

• trp – r^T * p

Note

1. It is permissible for p and trp to be the same array.

trxpv(r::AbstractMatrix{<:Real}, pv::AbstractVector{<:AbstractVector{<:Real}})

Multiply a pv-vector by the transpose of an r-matrix.

Input

• r – r-matrix
• pv – pv-vector

Output

• trpv – r^T * pv

Note

1. The algorithm is for the simple case where the r-matrix r is not a function of time. The case where r is a function of time leads to an additional velocity component equal to the product of the derivative of the transpose of r and the position vector.

2. It is permissible for pv and rpv to be the same array.

zp()

Zero a p-vector.

Output

• p – zero p-vector

zpv()

Zero a pv-vector.

Output

• pv – zero pv-vector

zr()

Initialize an r-matrix to the null matrix.

Output

• r – r-matrix