# Class Gaia

java.lang.Object
uk.ac.starlink.ttools.func.Gaia

public class Gaia extends Object
Functions related to astrometry suitable for use with data from the Gaia astrometry mission.

The methods here are not specific to the Gaia mission, but the parameters of the functions and their units are specified in a form that is convenient for use with Gaia data, in particular the `gaia_source` catalogue available from http://gea.esac.esa.int/archive/ and copies or mirrors.

There are currently three main sets of functions here:

- position and velocity vector calculation and manipulation
- distance estimation from parallaxes
- astrometry propagation to different epochs

Position and velocity vectors

Functions are provided for converting the astrometric parameters contained in the Gaia catalogue to ICRS Cartesian position (XYZ) and velocity (UVW) vectors. Functions are also provided to convert these vectors between ICRS and Galactic or Ecliptic coordinates. The calculations are fairly straightforward, and follow the equations laid out in section 1.5.6 of The Hipparcos and Tycho Catalogues, ESA SP-1200 (1997) and also section 3.1.7 of the Gaia DR2 documentation (2018).

These functions will often be combined; for instance to calculate the position and velocity in galactic coordinates from Gaia catalogue values, the following expressions may be useful:

       xyz_gal = icrsToGal(astromXYZ(ra,dec,parallax))
       uvw_gal = icrsToGal(astromUVW(array(ra,dec,parallax,pmra,pmdec,radial_velocity)))

though note that these particular examples simply invert parallax to provide distance estimates, which is not generally valid. Note also that these functions do not attempt to correct for solar motion. Such adjustments should be carried out by hand on the results of these functions if they are required.

Functions for calculating errors on the Cartesian components based on the error and correlation quantities from the Gaia catalogue are not currently provided. They would require fairly complicated invocations. If there is demand they may be implemented in the future.

Distance estimation

Gaia measures parallaxes, but some scientific use cases require the radial distance instead. While distance in parsec is in principle the reciprocal of parallax in arcsec, in the presence of non-negligable errors on measured parallax, this inversion does not give a good estimate of distance. A thorough discussion of this topic and approaches to estimating distances for Gaia-like data can be found in the papers

- C.A.L.Bailer-Jones, "Estimating distances from parallaxes", PASP 127, p994 (2015) 2015PASP..127..994B
- T.L.Astraatmadja and C.A.L.Bailer-Jones, "Estimating Distances from Parallaxes. II. Performance of Bayesian Distance Estimators on a Gaia-like Catalogue", ApJ 832, a137 (2016) 2016ApJ...832..137A
- X.Luri et al., "Gaia Data Release 2: Using Gaia Parallaxes", A&A in press (2018) arXiv:1804.09376

The functions provided here correspond to calculations from Astraatmadja & Bailer-Jones, "Estimating Distances from Parallaxes. III. Distances of Two Million Stars in the Gaia DR1 Catalogue", ApJ 833, a119 (2016) 2016ApJ...833..119A based on the Exponentially Decreasing Space Density prior defined therein. This implementation was written with reference to the Java implementation by Enrique Utrilla (DPAC).

These functions are parameterised by a length scale L that defines the exponential decay (the mode of the prior PDF is at r=2L). Some value for this length scale, specified in parsec, must be supplied to the functions as the `lpc` parameter.

Note that the values provided by these functions do not match those from the paper Bailer-Jones et al. "Estimating Distances from Parallaxes IV: Distances to 1.33 Billion stars in Gaia Data Release 2", accepted for AJ (2018) arXiv:1804.10121. The calculations of that paper differ from the ones presented here in several ways: it uses a galactic model for the direction-dependent length scale not currently available here, it pre-applies a parallax correction of -0.029mas, and it uses different uncertainty measures and in some cases (bimodal PDF) a different best distance estimator.

Epoch Propagation

The Gaia source catalogue provides, for at least some sources, the six-parameter astrometric solution (Right Ascension, Declination, Parallax, Proper motion in RA and Dec, and Radial Velocity), along with errors on these values and correlations between these errors. While a crude estimate of the position at an earlier or later epoch than that of the measurement can be made by multiplying the proper motion components by epoch difference and adding to the measured position, a more careful treatment is required for accurate propagation between epochs of the astrometric parameters, and if required their errors and correlations. The expressions for this are set out in section 1.5.5 (Volume 1) of The Hipparcos and Tycho Catalogues, ESA SP-1200 (1997) (but see below), and the code is based on an implementation by Alexey Butkevich and Daniel Michalik (DPAC). A correction is applied to the SP-1200 treatment of radial velocity uncertainty following Michalik et al. 2014 2014A&A...571A..85M because of their better handling of small radial velocities or parallaxes.

The calculations give the same results, though not exactly in the same form, as the epoch propagation functions available in the Gaia archive service.

Since:
  2 Mar 2018
Author:
  Mark Taylor

## Field Summary

Fields
Modifier and Type
Field
Description
`static final double`
`AU_YRKMS`
This quantity is A_v, the Astronomical Unit expressed in km.yr/sec.
`static final double`
`C_KMS`
The speed of light in km/s (exact).
`static final double`
`PC_AU`
Parsec in Astronomical Units, equal to 648000/PI.
`static final double`
`PC_YRKMS`
Parsec in units of km.yr/sec.

## Method Summary

Modifier and Type
Method
Description
`static double[]`
`astromUVW(double[] astrom6)`
Calculates Cartesian components of velocity from quantities available in the Gaia source catalogue.
`static double[]`
`astromUVW(double ra, double dec, double pmra, double pmdec, double radial_velocity, double r_parsec, boolean useDoppler)`
Calculates Cartesian components of velocity from the observed position and proper motion, radial velocity and radial distance, with optional light-time correction.
`static double[]`
`astromXYZ(double ra, double dec, double parallax)`
Calculates Cartesian components of position from RA, Declination and parallax.
`static double[]`
`distanceBoundsEdsd(double plxMas, double plxErrorMas, double lPc)`
Calculates the 5th and 95th percentile confidence intervals on the distance estimate using the Exponentially Decreasing Space Density prior.
`static double`
`distanceEstimateEdsd(double plxMas, double plxErrorMas, double lPc)`
Best estimate of distance using the Exponentially Decreasing Space Density prior.
`static double[]`
`distanceQuantilesEdsd(double plxMas, double plxErrorMas, double lPc, double... qpoints)`
Calculates arbitrary quantiles for the distance estimate using the Exponentially Decreasing Space Density prior.
`static double`
`distanceToModulus(double distPc)`
Converts a distance in parsec to a distance modulus.
`static double[]`
`eclToIcrs(double[] xyz)`
Converts a 3-element vector representing ecliptic coordinates to ICRS (equatorial) coordinates.
`static double[]`
`epochProp(double tYr, double[] astrom6)`
Propagates the astrometry parameters, supplied as a 6-element array, to a different epoch.
`static double[]`
`epochPropErr(double tYr, double[] astrom22)`
Propagates the astrometry parameters and their associated errors and correlations, supplied as a 22-element array, to a different epoch.
`static double[]`
`galToIcrs(double[] xyz)`
Converts a 3-element vector representing galactic coordinates to ICRS (equatorial) coordinates.
`static double[]`
`icrsToEcl(double[] xyz)`
Converts a 3-element vector representing ICRS (equatorial) coordinates to ecliptic coordinates.
`static double[]`
`icrsToGal(double[] xyz)`
Converts a 3-element vector representing ICRS (equatorial) coordinates to galactic coordinates.
`static double`
`modulusToDistance(double distmod)`
Converts a distance modulus to a distance in parsec.
`static double[]`
`polarXYZ(double phi, double theta, double r)`
Converts from spherical polar to Cartesian coordinates.
`static double`
`rvKmsToMasyr(double rvKms, double plxMas)`
Converts from unnormalised radial velocity in km/s to normalised radial velocity in mas/year.
`static double`
`rvMasyrToKms(double rvMasyr, double plxMas)`
Converts from normalised radial velocity in mas/year to unnormalised radial velocity in km/s.

### Methods inherited from class Object

`clone, equals, finalize, getClass, hashCode, notify, notifyAll, toString, wait, wait, wait`

## Field Details

### AU_YRKMS

public static final double AU_YRKMS
This quantity is A_v, the Astronomical Unit expressed in km.yr/sec. See the Hipparcos catalogue (ESA SP-1200) table 1.2.2 and Eq. 1.5.24.

See Also:
  - Constant Field Values

### PC_AU

public static final double PC_AU
Parsec in Astronomical Units, equal to 648000/PI.

See Also:
  - Constant Field Values

### PC_YRKMS

public static final double PC_YRKMS
Parsec in units of km.yr/sec.

See Also:
  - Constant Field Values

### C_KMS

public static final double C_KMS
The speed of light in km/s (exact).

See Also:
  - Constant Field Values

## Method Details

### polarXYZ

public static double[] polarXYZ(double phi, double theta, double r)
Converts from spherical polar to Cartesian coordinates.

Parameters:
  `phi` - longitude in degrees
  `theta` - latitude in degrees
  `r` - radial distance
Returns:
  3-element vector giving Cartesian coordinates
Examples:
  `polarXYZ(ra, dec, distance_estimate)`, `polarXYZ(l, b, 3262./parallax)` - calculates vector components in units of light year in the galactic system, on the assumption that distance is the inverse of parallax

### astromXYZ

public static double[] astromXYZ(double ra, double dec, double parallax)
Calculates Cartesian components of position from RA, Declination and parallax. This is a convenience function, equivalent to:

       polarXYZ(ra, dec, 1000./parallax)

Note that this performs distance scaling using a simple inversion of parallax, which is not in general reliable for parallaxes with non-negligable errors. Use at your own risk.

Parameters:
  `ra` - Right Ascension in degrees
  `dec` - Declination in degrees
  `parallax` - parallax in mas
Returns:
  3-element vector giving equatorial space coordinates in parsec
Examples:
  `astromXYZ(ra, dec, parallax)`, `icrsToGal(astromXYZ(ra, dec, parallax))`

### icrsToGal

public static double[] icrsToGal(double[] xyz)
Converts a 3-element vector representing ICRS (equatorial) coordinates to galactic coordinates. This can be used with position or velocity vectors.

The input vector is multiplied by the matrix AG', given in Eq. 3.61 of the Gaia DR2 documentation, following Eq. 1.5.13 of the Hipparcos catalogue.

The output coordinate system is right-handed, with the three components positive in the directions of the Galactic center, Galactic rotation, and the North Galactic Pole respectively.

Parameters:
  `xyz` - 3-element vector giving ICRS Cartesian components
Returns:
  3-element vector giving Galactic Cartesian components
Examples:
  `icrsToGal(polarXYZ(ra, dec, distance))`

### galToIcrs

public static double[] galToIcrs(double[] xyz)
Converts a 3-element vector representing galactic coordinates to ICRS (equatorial) coordinates. This can be used with position or velocity vectors.

The input vector is multiplied by the matrix AG, given in Eq. 3.61 of the Gaia DR2 documentation, following Eq. 1.5.13 of the Hipparcos catalogue.

The input coordinate system is right-handed, with the three components positive in the directions of the Galactic center, Galactic rotation, and the North Galactic Pole respectively.

Parameters:
  `xyz` - 3-element vector giving Galactic Cartesian components
Returns:
  3-element vector giving ICRS Cartesian components
Examples:
  `galToIcrs(polarXYZ(l, b, distance))`

### icrsToEcl

public static double[] icrsToEcl(double[] xyz)
Converts a 3-element vector representing ICRS (equatorial) coordinates to ecliptic coordinates. This can be used with position or velocity vectors.

The transformation corresponds to that between the coordinates `(ra,dec)` and `(ecl_lon,ecl_lat)` in the Gaia source catalogue (DR2).

Parameters:
  `xyz` - 3-element vector giving ICRS Cartesian components
Returns:
  3-element vector giving ecliptic Cartesian components
Examples:
  `icrsToEcl(polarXYZ(ra, dec, distance))`

### eclToIcrs

public static double[] eclToIcrs(double[] xyz)
Converts a 3-element vector representing ecliptic coordinates to ICRS (equatorial) coordinates. This can be used with position or velocity vectors.

The transformation corresponds to that between the coordinates `(ecl_lon,ecl_lat)` and `(ra,dec)` in the Gaia source catalogue (DR2).

Parameters:
  `xyz` - 3-element vector giving ecliptic Cartesian coordinates
Returns:
  3-element vector giving ICRS Cartesian coordinates
Examples:
  `eclToIcrs(polarXYZ(ecl_lon, ecl_lat, distance))`

### astromUVW

public static double[] astromUVW(double[] astrom6)
Calculates Cartesian components of velocity from quantities available in the Gaia source catalogue. The output is in the same coordinate system as the inputs, that is ICRS for the correspondingly-named Gaia quantities.

The input astrometry parameters are represented by a 6-element array, with the following elements:

    index  gaia_source name  unit    description
    -----  ----------------  ----    -----------
      0:   ra                deg     right ascension
      1:   dec               deg     declination
      2:   parallax          mas     parallax
      3:   pmra              mas/yr  proper motion in ra * cos(dec)
      4:   pmdec             mas/yr  proper motion in dec
      5:   radial_velocity   km/s    barycentric radial velocity

The units used by this function are the units used in the `gaia_source` table.

This convenience function just invokes the 7-argument `astromUVW` function using the inverted parallax for the radial distance, and without invoking the Doppler correction. It is exactly equivalent to:

       astromUVW(a[0], a[1], a[3], a[4], a[5], 1000./a[2], false)

Note this naive inversion of parallax to estimate distance is not in general reliable for parallaxes with non-negligable errors.

Parameters:
  `astrom6` - vector of 6 astrometric parameters as provided by the Gaia source catalogue
Returns:
  3-element vector giving equatorial velocity components in km/s
Examples:
  `astromUVW(array(ra, dec, parallax, pmra, pmdec, radial_velocity))`, `icrsToGal(astromUVW(array(ra, dec, parallax, pmra, pmdec, radial_velocity)))`

### astromUVW

public static double[] astromUVW(double ra, double dec, double pmra, double pmdec, double radial_velocity, double r_parsec, boolean useDoppler)
Calculates Cartesian components of velocity from the observed position and proper motion, radial velocity and radial distance, with optional light-time correction. The output is in the same coordinate system as the inputs, that is ICRS for the correspondingly-named Gaia quantities.

The radial distance must be supplied using the `r_parsec` parameter. A naive estimate from quantities in the Gaia source catalogue may be made with the expression `1000./parallax`, though note that this simple inversion of parallax is not in general reliable for parallaxes with non-negligable errors.

The calculations are fairly straightforward, following Eq. 1.5.74 from the Hipparcos catalogue. A (usually small) Doppler factor accounting for light-time effects can also optionally be applied. The effect of this is to multiply the returned vector by a factor of `1/(1-radial_velocity/c)`, as discussed in Eq. 1.2.21 of the Hipparcos catalogue.

Note that no attempt is made to adjust for solar motion.

Parameters:
  `ra` - Right Ascension in degrees
  `dec` - Declination in degrees
  `pmra` - proper motion in RA * cos(dec) in mas/yr
  `pmdec` - proper motion in declination in mas/yr
  `radial_velocity` - radial velocity in km/s
  `r_parsec` - radial distance in parsec
  `useDoppler` - whether to apply the Doppler factor to account for light-time effects
Returns:
  3-element vector giving equatorial velocity components in km/s
Examples:
  `astromUVW(ra, dec, pmra, pmdec, radial_velocity, dist, true)`, `icrsToGal(astromUVW(ra, dec, pmra, pmdec, radial_velocity, 1000./parallax, false))`

### epochProp

public static double[] epochProp(double tYr, double[] astrom6)
Propagates the astrometry parameters, supplied as a 6-element array, to a different epoch.

The input and output astrometry parameters are each represented by a 6-element array, with the following elements:

    index  gaia_source name  unit    description
    -----  ----------------  ----    -----------
      0:   ra                deg     right ascension
      1:   dec               deg     declination
      2:   parallax          mas     parallax
      3:   pmra              mas/yr  proper motion in ra * cos(dec)
      4:   pmdec             mas/yr  proper motion in dec
      5:   radial_velocity   km/s    barycentric radial velocity

The units used by this function are the units used in the `gaia_source` table.

Null values for `parallax`, `pmra`, `pmdec` and `radial_velocity` are treated as if zero for the purposes of propagation. The documentation of the equivalent function in the Gaia archive comments "This is a reasonable choice for most stars because those quantities would be either small (parallax and proper motion) or irrelevant (radial velocity). However, this is not true for stars very close to the Sun, where appropriate values need to be taken from the literature (e.g. average velocity field in the solar neighbourhood)."

The effect is that the output represents the best estimates available for propagated astrometry; proper motions, parallax and RV are applied if present, but if not the output values are calculated or simply copied across as if those quantities were zero.

Parameters:
  `tYr` - epoch difference in years
  `astrom6` - astrometry at time t0, represented by a 6-element array as above (a 5-element array is also permitted where radial velocity is zero or unknown)
Returns:
  astrometry at time `t0+tYr`, represented by a 6-element array as above
Examples:
  `epochProp(-15.5, array(ra,dec,parallax,pmra,pmdec,radial_velocity))` - calculates the astrometry at 2000.0 of gaia_source values that were observed at 2015.5

### epochPropErr

public static double[] epochPropErr(double tYr, double[] astrom22)
Propagates the astrometry parameters and their associated errors and correlations, supplied as a 22-element array, to a different epoch.

The input and output astrometry parameters with associated error and correlation information are each represented by a 22-element array, with the following elements:

    index  gaia_source name      unit    description
    -----  ----------------      ----    -----------
      0:   ra                    deg     right ascension
      1:   dec                   deg     declination
      2:   parallax              mas     parallax
      3:   pmra                  mas/yr  proper motion in RA * cos(dec)
      4:   pmdec                 mas/yr  proper motion in Declination
      5:   radial_velocity       km/s    barycentric radial velocity
      6:   ra_error              mas     error in right ascension
      7:   dec_error             mas     error in declination
      8:   parallax_error        mas     error in parallax
      9:   pmra_error            mas/yr  error in RA proper motion * cos(dec)
     10:   pmdec_error           mas/yr  error in Declination proper motion
     11:   radial_velocity_error km/s    error in barycentric radial velocity
     12:   ra_dec_corr                   correlation between ra and dec
     13:   ra_parallax_corr              correlation between ra and parallax
     14:   ra_pmra_corr                  correlation between ra and pmra
     15:   ra_pmdec_corr                 correlation between ra and pmdec
     16:   dec_parallax_corr             correlation between dec and parallax
     17:   dec_pmra_corr                 correlation between dec and pmra
     18:   dec_pmdec_corr                correlation between dec and pmdec
     19:   parallax_pmra_corr            correlation between parallax and pmra
     20:   parallax_pmdec_corr           correlation between parallax and pmdec
     21:   pmra_pmdec_corr               correlation between pmra and pmdec

Note the correlation coefficients, always in the range -1..1, are dimensionless.

Null values for `parallax`, `pmra`, `pmdec` and `radial_velocity` are treated as if zero for the purposes of propagation. The documentation of the equivalent function in the Gaia archive comments "This is a reasonable choice for most stars because those quantities would be either small (parallax and proper motion) or irrelevant (radial velocity). However, this is not true for stars very close to the Sun, where appropriate values need to be taken from the literature (e.g. average velocity field in the solar neighbourhood)."

The effect is that the output represents the best estimates available for propagated astrometry; proper motions, parallax and RV are applied if present, but if not the output values are calculated or simply copied across as if those quantities were zero.

This transformation is only applicable for radial velocities determined independently of the astrometry, such as those obtained with a spectrometer. It is not applicable for the back-transformation of data already propagated to another epoch.

This is clearly an unwieldy function to invoke, but if you are using it with the gaia_source catalogue itself, or other similar catalogues with the same column names and units, you can invoke it by just copying and pasting the example shown in this documentation.

Parameters:
  `tYr` - epoch difference in years
  `astrom22` - astrometry at time t0, represented by a 22-element array as above
Returns:
  astrometry at time t0+tYr, represented by a 22-element array as above
Examples:
  `epochPropErr(-15.5, array( ra,dec,parallax,pmra,pmdec,radial_velocity, ra_error,dec_error,parallax_error,pmra_error,pmdec_error,radial_velocity_error, ra_dec_corr,ra_parallax_corr,ra_pmra_corr,ra_pmdec_corr, dec_parallax_corr,dec_pmra_corr,dec_pmdec_corr, parallax_pmra_corr,parallax_pmdec_corr, pmra_pmdec_corr))` - calculates the astrometry with all errors and correlations at 2000.0 for gaia_source values that were observed at 2015.5.

### rvMasyrToKms

public static double rvMasyrToKms(double rvMasyr, double plxMas)
Converts from normalised radial velocity in mas/year to unnormalised radial velocity in km/s.

The output is calculated as `AU_YRKMS * rvMasyr / plxMas`, where `AU_YRKMS=4.740470446` is one Astronomical Unit in km.yr/sec.

Parameters:
  `rvMasyr` - normalised radial velocity, in mas/year
  `plxMas` - parallax in mas
Returns:
  radial velocity in km/s

### rvKmsToMasyr

public static double rvKmsToMasyr(double rvKms, double plxMas)
Converts from unnormalised radial velocity in km/s to normalised radial velocity in mas/year.

The output is calculated as `rvKms * plxMas / AU_YRKMS`, where `AU_YRKMS=4.740470446` is one Astronomical Unit in km.yr/sec.

Parameters:
  `rvKms` - unnormalised radial velocity, in mas/year
  `plxMas` - parallax in mas
Returns:
  radial velocity in mas/year

### distanceEstimateEdsd

public static double distanceEstimateEdsd(double plxMas, double plxErrorMas, double lPc)
Best estimate of distance using the Exponentially Decreasing Space Density prior. This estimate is provided by the mode of the PDF.

Parameters:
  `plxMas` - parallax in mas
  `plxErrorMas` - parallax error in mas
  `lPc` - length scale in parsec
Returns:
  best distance estimate in parsec

### distanceBoundsEdsd

public static double[] distanceBoundsEdsd(double plxMas, double plxErrorMas, double lPc)
Calculates the 5th and 95th percentile confidence intervals on the distance estimate using the Exponentially Decreasing Space Density prior.

Note this function has to numerically integrate the PDF to determine quantile values, so it is relatively slow.

Parameters:
  `plxMas` - parallax in mas
  `plxErrorMas` - parallax error in mas
  `lPc` - length scale in parsec
Returns:
  2-element array giving the 5th and 95th percentiles in parsec of the EDSD distance PDF

### distanceQuantilesEdsd

public static double[] distanceQuantilesEdsd(double plxMas, double plxErrorMas, double lPc, double... qpoints)
Calculates arbitrary quantiles for the distance estimate using the Exponentially Decreasing Space Density prior.

Note this function has to numerically integrate the PDF to determine quantile values, so it is relatively slow.

Parameters:
  `plxMas` - parallax in mas
  `plxErrorMas` - parallax error in mas
  `lPc` - length scale in parsec
  `qpoints` - one or more required quantile cut points, each in the range 0..1
Returns:
  array with one element for each of the supplied `qpoints` giving the corresponding distance in parsec
Examples:
  `distanceQuantilesEdsd(parallax, parallax_error, 1350, 0.5)[0]` calculates the median of the EDSD distance PDF using a length scale of 1.35kpc, `distanceQuantilesEdsd(parallax, parallax_error, 3000, 0.01, 0.99)` returns a 2-element array giving the 1st and 99th percentile of the distance estimate using a length scale of 3kpc

### distanceToModulus

public static double distanceToModulus(double distPc)
Converts a distance in parsec to a distance modulus. The formula is `5*log10(distPc)-5`.

Parameters:
  `distPc` - distance in parsec
Returns:
  distance modulus in magnitudes

### modulusToDistance

public static double modulusToDistance(double distmod)
Converts a distance modulus to a distance in parsec. The formula is `10^(1+distmod/5)`.

Parameters:
  `distmod` - distance modulus in magnitudes
Returns:
  distance in parsec