Class Pal

java.lang.Object
uk.ac.starlink.pal.Pal

public class Pal extends Object
Positional Astronomy Library.
Version:
1.0 [Latest Revision: January 2003]

Based on the C version of slalib written by P T Wallace.

Author:
R T Platon (Starlink)
  • Field Summary

    Fields
    Modifier and Type
    Field
    Description
    static final double
    Astronomical unit to kilometers
    static final double
    Light time for 1 AU (sec)
    static final double
    Speed of light (AU per day)
    static final double
    Seconds in a day
    static final double
    Nominal mean sidereal speed of Earth equator in km/s (the actual value is about 0.4651)
    int
    Flag for additional status information
    static final double
    Gravitational radius of the Sun x 2: (2*mu/c**2, au)
    static final double
    Degrees to radians
    static final double
    Ratio between solar and sidereal time
    int
    Current Status flag
    static final double
    Km/s to AU/year
  • Constructor Summary

    Constructors
    Constructor
    Description
    Pal()
     
  • Method Summary

    Modifier and Type
    Method
    Description
    Addet(AngleDR m, double eq)
    Add the e-terms (elliptic component of annual aberration) to a pre IAU 1976 mean place to conform to the old catalogue convention.
    Amp(AngleDR ap, double date, double eq)
    Convert star RA,Dec from geocentric apparent to mean place.
    Ampqk(AngleDR ap, AMParams amprms)
    Convert star RA,Dec from geocentric apparent to mean place.
    Aoppa(UTCdate date, ObsPosition obs, Cartesian pm, double tdk, double pmb, double rh, double wl, double tlr)
    Precompute apparent to observed place parameters required by Aopqk and Oapqk.
    void
    Aoppat(double date, AOParams aoprms)
    Recompute the sidereal time in the apparent to observed place star-independent parameter block.
    double
    Caldj(int iy, int im, int id)
    Gregorian calendar to Modified Julian Date.
    double
    Cldj(int iy, int im, int id)
    Gregorian calendar to Modified Julian Date.
    double
    Daf2r(int ideg, int iamin, double asec)
    Convert degrees, arcminutes, arcseconds to radians.
    double
    Dat(double utc)
    Increment to be applied to Coordinated Universal Time UTC to give International Atomic Time TAI.
    double[][]
    Dav2m(double[] axvec)
    Form the rotation matrix corresponding to a given axial vector.
    double
    Dbjin(palString string, double dreslt)
    Convert free-format input into double precision floating point, using Dfltin but with special syntax extensions.
    Conversion of position & velocity in Cartesian coordinates to spherical coordinates.
    Dcc2s(double[] v)
    Direction cosines to spherical coordinates.
    double[]
    Spherical coordinates to direction cosines.
    Dd2tf(double days)
    Convert an interval in days into hours, minutes, seconds.
    double[][]
    Deuler(String order, double phi, double theta, double psi)
    Form a rotation matrix from the Euler angles - three successive rotations about specified Cartesian axes.
    double
    Dfltin(palString string, double dreslt)
    Convert free-format input into double precision floating point.
    double[]
    Dimxv(double[][] dm, double[] va)
    Performs the 3-d backward unitary transformation.
    Djcal(double djm)
    Modified Julian Date to Gregorian calendar, expressed in a form convenient for formatting messages (namely rounded to a specified precision, and with the fields stored in a single array).
    Djcl(double djm)
    Modified Julian Date to Gregorian year, month, day, and fraction of a day.
    double[]
    Dm2av(double[][] rmat)
    From a rotation matrix, determine the corresponding axial vector.
    double
    Dmat(double[][] a, double[] y)
    Matrix inversion & solution of simultaneous equations.
    double[][]
    Dmxm(double[][] a, double[][] b)
    Product of two 3x3 matrices.
    double[]
    Dmxv(double[][] dm, double[] va)
    Performs the 3-d forward unitary transformation.
    Dr2af(double angle)
    Convert an angle in radians into degrees, arcminutes, arcseconds.
    Dr2tf(double angle)
    Convert an angle in radians to hours, minutes, seconds.
    double
    Drange(double angle)
    Normalize angle into range +/- pi.
    double
    Dranrm(double angle)
    Normalize angle into range 0-2 π.
    double
    Drverot(double phi, AngleDR r, double st)
    Velocity component in a given direction due to Earth rotation.
    double
    Velocity component in a given direction due to the rotation of the Galaxy.
    double
    Drvlg(AngleDR r2000)
    Velocity component in a given direction due to the combination of the rotation of the Galaxy and the motion of the Galaxy relative to the mean motion of the local group.
    double
    Velocity component in a given direction due to the Sun's motion with respect to the dynamical Local Standard of Rest.
    double
    Velocity component in a given direction due to the Sun's motion with respect to an adopted kinematic Local Standard of Rest.
    Conversion of position & velocity in spherical coordinates to Cartesian coordinates.
    Projection of spherical coordinates onto tangent plane ('gnomonic' projection - 'standard coordinates').
    double
    Dt(double epoch)
    Estimate the offset between dynamical time and Universal Time for a given historical epoch.
    double
    Dtf2d(int ihour, int imin, double sec)
    Convert hours, minutes, seconds to days.
    double
    Dtf2r(int ihour, int imin, double sec)
    Convert hours, minutes, seconds to radians.
    Transform tangent plane coordinates into spherical.
    double
    Dtt(double utc)
    Increment to be applied to Coordinated Universal Time UTC to give Terrestrial Time TT (formerly Ephemeris Time ET).
    double
    Dvdv(double[] va, double[] vb)
    Scalar product of two 3-vectors.
    double
    Dvn(double[] v, double[] uv)
    Normalizes a 3-vector also giving the modulus.
    double[]
    Dvxv(double[] va, double[] vb)
    Vector product of two 3-vectors.
    Ecleq(AngleDR dl, double date)
    Transformation from ecliptic coordinates to J2000.0 equatorial coordinates.
    double[][]
    Ecmat(double date)
    Form the equatorial to ecliptic rotation matrix (IAU 1980 theory).
    double
    Epb(double date)
    Conversion of Modified Julian Date to Besselian epoch.
    double
    Epb2d(double epb)
    Conversion of Besselian epoch to Modified Julian Date.
    double
    Epco(char k0, char k, double e)
    Convert an epoch into the appropriate form - 'B' or 'J'.
    double
    Epj(double date)
    Conversion of Modified Julian Date to Julian epoch.
    double
    Epj2d(double epj)
    Conversion of Julian epoch to Modified Julian Date.
    Eqecl(AngleDR d, double date)
    Transformation from J2000.0 equatorial coordinates to ecliptic coordinates.
    double
    Eqeqx(double date)
    Equation of the equinoxes (IAU 1994, double precision).
    Transformation from J2000.0 equatorial coordinates to IAU 1958 Galactic coordinates.
    double[]
    Etrms(double ep)
    Compute the e-terms (elliptic component of annual aberration) vector.
    void
    Evp(double date, double deqx, double[] dvb, double[] dpb, double[] dvh, double[] dph)
    Barycentric and heliocentric velocity and position of the Earth.
    Fk425(Stardata s1950)
    Convert B1950.0 FK4 star data to J2000.0 FK5.
    Fk45z(AngleDR r1950, double bepoch)
    Convert B1950.0 FK4 star data to J2000.0 FK5 assuming zero proper motion in the FK5 frame (double precision).
    Fk524(Stardata j2000)
    Convert J2000.0 FK5 star data to B1950.0 FK4.
    Fk54z(AngleDR r2000, double bepoch)
    Convert a J2000.0 FK5 star position to B1950.0 FK4 assuming zero proper motion and parallax.
    Transformation from IAU 1958 Galactic coordinates to J2000.0 equatorial coordinates.
    Transformation from IAU 1958 Galactic coordinates to De Vaucouleurs supergalactic coordinates.
    double[]
    Geoc(double p, double h)
    Convert geodetic position to geocentric.
    double
    Gmst(double ut1)
    Conversion from Universal Time to Sidereal Time.
    char
    Kbj(int jb, double e)
    Select epoch prefix 'B' or 'J'.
    Map(Stardata sd, double epq, double date)
    Transform star RA,Dec from mean place to geocentric apparent.
    Mappa(double eq, double date)
    Compute star-independent parameters in preparation for conversions between mean place and geocentric apparent place.
    Mapqk(Stardata s, AMParams amprms)
    Quick mean to apparent place: transform a star RA,Dec from mean place to geocentric apparent place, given the star-independent parameters.
    Mapqkz(AngleDR rm, AMParams amprms)
    Quick mean to apparent place: transform a star RA,dec from mean place to geocentric apparent place, given the star-independent parameters, and assuming zero parallax and proper motion.
    double[][]
    Nut(double date)
    Form the matrix of nutation for a given date (IAU 1980 theory).
    double[]
    Nutc(double date)
    Nutation: longitude & obliquity components and mean obliquity (IAU 1980 theory).
    Obs(int n)
    Parameters of selected groundbased observing stations.
    Obs(String id)
    Parameters of selected groundbased observing stations.
    Pm(AngleDR r0, double[] pm, double px, double rv, double ep0, double ep1)
    Apply corrections for proper motion to a star RA,Dec.
    double[][]
    Prebn(double bep0, double bep1)
    Generate the matrix of precession between two epochs, using the old, pre-IAU1976, Bessel-Newcomb model, using Kinoshita's formulation (double precision).
    double[][]
    Prec(double ep0, double ep1)
    Form the matrix of precession between two epochs (IAU 1976, FK5).
    Preces(String sys, double ep0, double ep1, AngleDR d)
    Precession - either FK4 (Bessel-Newcomb, pre-IAU1976) or FK5 (Fricke, post-IAU1976) as required.
    double[][]
    Precl(double ep0, double ep1)
    Form the matrix of precession between two epochs, using the model of Simon et al (1994), which is suitable for long periods of time.
    double[][]
    Prenut(double epoch, double date)
    Form the matrix of precession and nutation (IAU 1976/1980/FK5).
    double[]
    Refco(double hm, double tdk, double pmb, double rh, double wl, double phi, double tlr, double eps)
    Determine constants A and B in atmospheric refraction model dz = A tan z + B tan^3 z.
    double
    Refro(double zobs, double hm, double tdk, double pmb, double rh, double wl, double phi, double tlr, double eps)
    Atmospheric refraction for radio and optical/IR wavelengths.
    Subet(AngleDR rc, double eq)
    Remove the e-terms (elliptic component of annual aberration) from a pre IAU 1976 catalogue RA,Dec to give a mean place.
    Transformation from De Vaucouleurs supergalactic coordinates to IAU 1958 Galactic coordinates.
    double
    Zd(double ha, double dec, double phi)
    HA, Dec to Zenith Distance.

    Methods inherited from class java.lang.Object

    clone, equals, finalize, getClass, hashCode, notify, notifyAll, toString, wait, wait, wait
  • Field Details

    • GR2

      public static final double GR2
      Gravitational radius of the Sun x 2: (2*mu/c**2, au)
      See Also:
    • VF

      public static final double VF
      Km/s to AU/year
      See Also:
    • C

      public static final double C
      Speed of light (AU per day)
      See Also:
    • SOLSID

      public static final double SOLSID
      Ratio between solar and sidereal time
      See Also:
    • ESPEED

      public static final double ESPEED
      Nominal mean sidereal speed of Earth equator in km/s (the actual value is about 0.4651)
      See Also:
    • D2S

      public static final double D2S
      Seconds in a day
      See Also:
    • AUKM

      public static final double AUKM
      Astronomical unit to kilometers
      See Also:
    • AUSEC

      public static final double AUSEC
      Light time for 1 AU (sec)
      See Also:
    • R2D

      public static final double R2D
      Degrees to radians
      See Also:
    • Status

      public int Status
      Current Status flag
    • Flag

      public int Flag
      Flag for additional status information
  • Constructor Details

    • Pal

      public Pal()
  • Method Details

    • Addet

      public AngleDR Addet(AngleDR m, double eq)
      Add the e-terms (elliptic component of annual aberration) to a pre IAU 1976 mean place to conform to the old catalogue convention.
      Explanation:
      Most star positions from pre-1984 optical catalogues (or derived from astrometry using such stars) embody the e-terms. If it is necessary to convert a formal mean place (for example a pulsar timing position) to one consistent with such a star catalogue, then the RA,Dec should be adjusted using this routine.
      Reference:
      Explanatory Supplement to the Astronomical Almanac, ed P.K.Seidelmann (1992), page 169.
      Parameters:
      m - RA,Dec (radians) without e-terms
      eq - Besselian epoch of mean equator and equinox
      Returns:
      RA,dec (radians) with e-terms included
    • Amp

      public AngleDR Amp(AngleDR ap, double date, double eq)
      Convert star RA,Dec from geocentric apparent to mean place.

      The mean coordinate system is the post IAU 1976 system, loosely called FK5.

      Notes:
      1. The distinction between the required TDB and TT is always negligible. Moreover, for all but the most critical applications UTC is adequate.
      2. Iterative techniques are used for the aberration and light deflection corrections so that the routines Amp (or Ampqk) and Map (or Mapqk) are accurate inverses; even at the edge of the Sun's disc the discrepancy is only about 1 nanoarcsecond.
      3. Where multiple apparent places are to be converted to mean places, for a fixed date and equinox, it is more efficient to use the Mappa routine to compute the required parameters once, followed by one call to Ampqk per star.
      4. The accuracy is sub-milliarcsecond, limited by the precession-nutation model (IAU 1976 precession, Shirai & Fukushima 2001 forced nutation and precession corrections).
      5. The accuracy is further limited by the routine Evp, called by Mappa, which computes the Earth position and velocity using the methods of Stumpff. The maximum error is about 0.3 mas.
      References:
      1984 Astronomical Almanac, pp B39-B41. (also Lederle & Schwan, Astron. Astrophys. 134, 1-6, 1984)
      Parameters:
      ap - apparent RA & Dec (radians)
      date - TDB for apparent place (JD-2400000.5)
      eq - equinox: Julian epoch of mean place
      Returns:
      mean RA & Dec (Radians)
    • Ampqk

      public AngleDR Ampqk(AngleDR ap, AMParams amprms)
      Convert star RA,Dec from geocentric apparent to mean place.

      The mean coordinate system is the post IAU 1976 system, loosely called FK5.

      Use of this routine is appropriate when efficiency is important and where many star positions are all to be transformed for one epoch and equinox. The star-independent parameters can be obtained by calling the Mappa routine.

      References:
      1984 Astronomical Almanac, pp B39-B41. (also Lederle & Schwan, Astron. Astrophys. 134, 1-6, 1984)
      Note:
      Iterative techniques are used for the aberration and light deflection corrections so that the routines Amp (or Ampqk) and Map (or Mapqk) are accurate inverses; even at the edge of the Sun's disc the discrepancy is only about 1 nanoarcsecond.
      Parameters:
      ap - apparent RA & Dec (radians)
      amprms - star-independent mean-to-apparent parameters
      Returns:
      mean RA & Dec (radians)
    • Aoppa

      public AOParams Aoppa(UTCdate date, ObsPosition obs, Cartesian pm, double tdk, double pmb, double rh, double wl, double tlr)

      Precompute apparent to observed place parameters required by Aopqk and Oapqk.

      Notes:
      1. It is advisable to take great care with units, as even unlikely values of the input parameters are accepted and processed in accordance with the models used.
      2. The date argument is UTC expressed as an MJD. This is, strictly speaking, improper, because of leap seconds. However, as long as the delta UT and the UTC are consistent there are no difficulties, except during a leap second. In this case, the start of the 61st second of the final minute should begin a new MJD day and the old pre-leap delta UT should continue to be used. As the 61st second completes, the MJD should revert to the start of the day as, simultaneously, the delta UTC changes by one second to its post-leap new value.
      3. The delta UT (UT1-UTC) is tabulated in IERS circulars and elsewhere. It increases by exactly one second at the end of each UTC leap second, introduced in order to keep delta UT within +/- 0.9 seconds.
      4. IMPORTANT -- TAKE CARE WITH THE LONGITUDE SIGN CONVENTION. The longitude required by the present routine is east-positive, in accordance with geographical convention (and right-handed). In particular, note that the longitudes returned by the Obs routine are west-positive, following astronomical usage, and must be reversed in sign before use in the present routine.
      5. The polar coordinates xp,yp can be obtained from IERS circulars and equivalent publications. The maximum amplitude is about 0.3 arcseconds. If xp,yp values are unavailable, use xp=yp=0.0. See page B60 of the 1988 Astronomical Almanac for a definition of the two angles.
      6. The height above sea level of the observing station, hm, can be obtained from the Astronomical Almanac (Section J in the 1988 edition), or via the routine Obs. If p, the pressure in millibars, is available, an adequate estimate of hm can be obtained from the expression

        hm = -29.3 * tsl * log ( p / 1013.25 );

        where tsl is the approximate sea-level air temperature in deg K (See Astrophysical Quantities, C.W.Allen, 3rd edition, section 52). Similarly, if the pressure p is not known, it can be estimated from the height of the observing station, hm as follows:

        p = 1013.25 * exp ( -hm / ( 29.3 * tsl ) );

        Note, however, that the refraction is proportional to the pressure and that an accurate p value is important for precise work.
      7. Repeated, computationally-expensive, calls to Aoppa for times that are very close together can be avoided by calling Aoppa just once and then using Aoppat for the subsequent times. Fresh calls to Aoppa will be needed only when changes in the precession have grown to unacceptable levels or when anything affecting the refraction has changed.
      Parameters:
      date - UTC date/time (Modified Julian Date, JD-2400000.5) & delta UT: UT1-UTC (UTC seconds)
      pm - mean longitude of the observer (radians, east +ve), mean geodetic latitude of the observer (radians), observer's height above sea level (metres) & polar motion x-coordinate (radians)
      tdk - local ambient temperature (DegK; std=273.155)
      pmb - local atmospheric pressure (mB; std=1013.25)
      rh - local relative humidity (in the range 0.0-1.0)
      wl - effective wavelength (micron, e.g. 0.55)
      tlr - tropospheric lapse rate (DegK/metre, e.g. 0.0065)
      Returns:
      aoprms star-independent apparent-to-observed parameters
    • Aoppat

      public void Aoppat(double date, AOParams aoprms)
      Recompute the sidereal time in the apparent to observed place star-independent parameter block.

      For more information, see Aoppa.

      Parameters:
      date - UTC date/time (Modified Julian Date, JD-2400000.5) (see Aoppa source for comments on leap seconds)
      aoprms - star-independent apparent-to-observed parameters
    • Caldj

      public double Caldj(int iy, int im, int id) throws palError
      Gregorian calendar to Modified Julian Date.

      (Includes century default feature: use Cldj for years before 100AD.)

      Parameters:
      iy - Year in Gregorian calendar
      im - Month in Gregorian calendar
      id - Day in Gregorian calendar
      Returns:
      Modified Julian Date (JD-2400000.5) for 0 hrs
      Throws:
      palError - if bad day, month or year
        0 = ok
        1 = bad year   (MJD not computed)
        2 = bad month  (MJD not computed)
        3 = bad day    (MJD computed)
        
        Acceptable years are 00-49, interpreted as 2000-2049,
                             50-99,     "       "  1950-1999,
                             100 upwards, interpreted literally.
        
    • Cldj

      public double Cldj(int iy, int im, int id) throws palError
      Gregorian calendar to Modified Julian Date.

      The year must be -4699 (i.e. 4700BC) or later.

      The algorithm is derived from that of Hatcher 1984 (QJRAS 25, 53-55).

      Parameters:
      iy - Year in Gregorian calendar
      im - Month in Gregorian calendar
      id - Day in Gregorian calendar
      Returns:
      Modified Julian Date (JD-2400000.5) for 0 hrs
      Throws:
      palError - if bad day, month or year
        0 = OK
        1 = bad year   (MJD not computed)
        2 = bad month  (MJD not computed)
        3 = bad day    (MJD computed)
        
    • Daf2r

      public double Daf2r(int ideg, int iamin, double asec) throws palError
      Convert degrees, arcminutes, arcseconds to radians.
      Notes:
      1. The result is computed even if any of the range checks fail.
      2. The sign must be dealt with outside this routine.
      Parameters:
      ideg - Degrees
      iamin - Arcminutes
      asec - Arcseconds
      Returns:
      Angle in radians
      Throws:
      palError - degrees, arcmins or arcsecs out of range
        Status returned:
            1 = ideg outside range 0-359
            2 = iamin outside range 0-59
            3 = asec outside range 0-59.999...
        
    • Dat

      public double Dat(double utc)
      Increment to be applied to Coordinated Universal Time UTC to give International Atomic Time TAI.
      Notes:
      1. The UTC is specified to be a date rather than a time to indicate that care needs to be taken not to specify an instant which lies within a leap second. Though in most cases the utc argument can include the fractional part, correct behaviour on the day of a leap second can only be guaranteed up to the end of the second 23:59:59.
      2. For epochs from 1961 January 1 onwards, the expressions from the file ftp://maia.usno.navy.mil/ser7/tai-utc.dat are used.
      3. The 5ms timestep at 1961 January 1 is taken from 2.58.1 (p87) of the 1992 Explanatory Supplement.
      4. UTC began at 1960 January 1.0 (JD 2436934.5) and it is improper to call the routine with an earlier epoch. However, if this is attempted, the TAI-UTC expression for the year 1960 is used.
      Latest leap second:
      2017 January 1
      Parameters:
      utc - UTC date as a modified JD (JD-2400000.5)
      Returns:
      TAI-UTC in seconds
    • Dav2m

      public double[][] Dav2m(double[] axvec)
      Form the rotation matrix corresponding to a given axial vector.

      A rotation matrix describes a rotation about some arbitrary axis. The axis is called the Euler axis, and the angle through which the reference frame rotates is called the Euler angle. The axial vector supplied to this routine has the same direction as the Euler axis, and its magnitude is the Euler angle in radians.

      If axvec is null, the unit matrix is returned.

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

      Parameters:
      axvec - Axial vector (radians)
      Returns:
      Rotation matrix
    • Dbjin

      public double Dbjin(palString string, double dreslt)
      Convert free-format input into double precision floating point, using Dfltin but with special syntax extensions.

      The purpose of the syntax extensions is to help cope with mixed FK4 and FK5 data. In addition to the syntax accepted by Dfltin, the following two extensions are recognized by dbjin:

      1. A valid non-null field preceded by the character 'B' (or 'b') is accepted.
      2. A valid non-null field preceded by the character 'J' (or 'j') is accepted.

      The calling program is notified of the incidence of either of these extensions through an supplementary status argument. The rest of the arguments are as for Dfltin.

      The Status returned is one of the following:

      -1 or 0 = OK
      1 = null field
      2 = error

      And the additional Flag is one of the following:

      0 = normal Dfltin syntax
      1 = 'B' or 'b'
      2 = 'J' or 'j'

      For details of the basic syntax, see Dfltin.

      Parameters:
      string - String containing field to be decoded
      dreslt - Previous Result
      Returns:
      Result
    • Dc62s

      public Spherical Dc62s(Cartesian v)
      Conversion of position & velocity in Cartesian coordinates to spherical coordinates.
      Parameters:
      v - Cartesian position & velocity vector
      Returns:
      Spherical Coordinates (Radians) - Longitude, Latitude, Radial plus derivitives
    • Dcc2s

      public AngleDR Dcc2s(double[] v)
      Direction cosines to spherical coordinates.

      The spherical coordinates are longitude (+ve anticlockwise looking from the +ve latitude pole) and latitude. The Cartesian coordinates are right handed, with the x axis at zero longitude and latitude, and the z axis at the +ve latitude pole.

      If v is null, zero a and b are returned. At either pole, zero a is returned.

      Parameters:
      v - x, y, z vector
      Returns:
      spherical coordinates in radians (RA, Dec)
    • Dcs2c

      public double[] Dcs2c(AngleDR a)
      Spherical coordinates to direction cosines.

      The spherical coordinates are longitude (+ve anticlockwise looking from the +ve latitude pole) and latitude. The Cartesian coordinates are right handed, with the x axis at zero longitude and latitude, and the z axis at the +ve latitude pole.

      Parameters:
      a - spherical coordinates in radians (RA,Dec)
      Returns:
      x, y, z unit vector
    • Dd2tf

      public palTime Dd2tf(double days)
      Convert an interval in days into hours, minutes, seconds.
      Parameters:
      days - interval in days
      Returns:
      hours, minutes, seconds, fraction
    • Deuler

      public double[][] Deuler(String order, double phi, double theta, double psi)
      Form a rotation matrix from the Euler angles - three successive rotations about specified Cartesian axes.

      A rotation is positive when the reference frame rotates anticlockwise as seen looking towards the origin from the positive region of the specified axis.

      The characters of order define which axes the three successive rotations are about. A typical value is 'zxz', indicating that rmat is to become the direction cosine matrix corresponding to rotations of the reference frame through phi radians about the old z-axis, followed by theta radians about the resulting x-axis, then psi radians about the resulting z-axis.

      The axis names can be any of the following, in any order or combination: x, y, z, uppercase or lowercase, 1, 2, 3. Normal axis labelling/numbering conventions apply; the xyz (=123) triad is right-handed. Thus, the 'zxz' example given above could be written 'zxz' or '313' (or even 'zxz' or '3xz'). Order is terminated by length or by the first unrecognized character.

      Fewer than three rotations are acceptable, in which case the later angle arguments are ignored. Zero rotations leaves rmat set to the identity matrix.

      Parameters:
      order - specifies about which axes the rotations occur
      phi - 1st rotation (radians)
      theta - 2nd rotation ( " )
      psi - 3rd rotation ( " )
      Returns:
      rotation matrix
    • Dfltin

      public double Dfltin(palString string, double dreslt)
      Convert free-format input into double precision floating point.
      Notes:
      1. A tab character is interpreted as a space, and lower case d,e are interpreted as upper case.
      2. The basic format is #^.^@#^ where # means + or -, ^ means a decimal subfield and @ means D or E.
      3. Spaces: Leading spaces are ignored. Embedded spaces are allowed only after # and D or E, and after . where the first ^ is absent. Trailing spaces are ignored; the first signifies end of decoding and subsequent ones are skipped.
      4. Field separators: Any character other than +,-,0-9,.,D,E or space may be used to end a field. Comma is recognized by Dfltin as a special case; it is skipped, leaving the pointer on the next character. See 12, below.
      5. Both signs are optional. The default is +.
      6. The mantissa defaults to 1.
      7. The exponent defaults to e0.
      8. The decimal subfields may be of any length.
      9. The decimal point is optional for whole numbers.
      10. A null field is one that does not begin with +,-,0-9,.,D or E, or consists entirely of spaces. If the field is null, jflag is set to 1 and dreslt is left untouched.
      11. nstrt = 1 for the first character in the string.
      12. On return from Dfltin, nstrt is set ready for the next decode - following trailing blanks and (if used) the comma separator. If a separator other than comma is being used, nstrt must be incremented before the next call to Dfltin.
      13. Errors (jflag=2) occur when: a) A +, -, D or E is left unsatisfied. b) The decimal point is present without at least one decimal subfield. c) An exponent more than 100 has been presented.
      14. When an error has been detected, nstrt is left pointing to the character following the last one used before the error came to light. This may be after the point at which a more sophisticated program could have detected the error. For example, Dfltin does not detect that '1e999' is unacceptable until the whole field has been read.
      15. Certain highly unlikely combinations of mantissa & exponent can cause arithmetic faults during the decode, in some cases despite the fact that they together could be construed as a valid number.
      16. Decoding is left to right, one pass.
      17. End of field may occur in either of two ways: a) As dictated by the string length. b) Detected during the decode. (b overrides a.)
      18. See also Flotin and Intin.
      Parameters:
      string - String containing field to be decoded
      dreslt - Previous result
      Returns:
      Result
    • Dimxv

      public double[] Dimxv(double[][] dm, double[] va)
      Performs the 3-d backward unitary transformation.
      vector vb = (inverse of matrix dm) * vector va

      (n.b. The matrix must be unitary, as this routine assumes that the inverse and transpose are identical).

      Parameters:
      dm - n x n matrix
      va - vector
      Returns:
      vector
    • Djcal

      public mjDate Djcal(double djm) throws palError
      Modified Julian Date to Gregorian calendar, expressed in a form convenient for formatting messages (namely rounded to a specified precision, and with the fields stored in a single array).

      Any date after 4701BC March 1 is accepted.

      Large ndp values risk internal overflows. It is typically safe to use up to ndp=4.

      The algorithm is derived from that of Hatcher 1984 (QJRAS 25, 53-55).

      Parameters:
      djm - Modified Julian Date (JD-2400000.5)
      Returns:
      year, month, day, fraction in Gregorian calendar
      Throws:
      palError
    • Djcl

      public mjDate Djcl(double djm) throws palError
      Modified Julian Date to Gregorian year, month, day, and fraction of a day.

      The algorithm is derived from that of Hatcher 1984 (QJRAS 25, 53-55).

      Parameters:
      djm - Modified Julian Date (JD-2400000.5)
      Returns:
      Year, month, day and fraction of day
      Throws:
      palError - unacceptable date (before 4701BC March 1)
    • Dm2av

      public double[] Dm2av(double[][] rmat)
      From a rotation matrix, determine the corresponding axial vector.

      A rotation matrix describes a rotation about some arbitrary axis. The axis is called the Euler axis, and the angle through which the reference frame rotates is called the Euler angle. The axial vector returned by this routine has the same direction as the Euler axis, and its magnitude is the Euler angle in radians. (The magnitude and direction can be separated by means of the routine Dvn.)

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

      If rmat is null, so is the result.

      Parameters:
      rmat - Rotation matrix
      Returns:
      Axial vector (radians)
    • Dmat

      public double Dmat(double[][] a, double[] y)
      Matrix inversion & solution of simultaneous equations.
      For the set of n simultaneous equations in n unknowns:
      a.y = x
      where:
      a is a non-singular n x n matrix
      y is the vector of n unknowns
      x is the known vector
      Dmat computes:
      the inverse of matrix a
      the determinant of matrix a
      the vector of n unknowns
      Arguments:
      symbol type dimension before after
      n int no. of unknowns unchanged
      *a double [n][n] matrix inverse
      *y double [n] vector solution
      *d double - determinant
      *jf int - # singularity flag
      *iw int [n] - workspace
      # jf is the singularity flag. If the matrix is non-singular, jf=0 is returned. If the matrix is singular, jf=-1 & d=0.0 are returned. In the latter case, the contents of array a on return are undefined.
      Algorithm:
      Gaussian elimination with partial pivoting.
      Speed:
      Very fast.
      Accuracy:
      Fairly accurate - errors 1 to 4 times those of routines optimized for accuracy.
      Parameters:
      a - Matrix
      y - Vector
      Returns:
      Determinant
    • Dmxm

      public double[][] Dmxm(double[][] a, double[][] b)
      Product of two 3x3 matrices.
      matrix c = matrix a x matrix b
      Note:
      the same array may be nominated more than once.
      Parameters:
      a - Matrix
      b - Matrix
      Returns:
      Matrix result
    • Dmxv

      public double[] Dmxv(double[][] dm, double[] va)
      Performs the 3-d forward unitary transformation.
      vector vb = matrix dm * vector va
      Parameters:
      dm - Matrix
      va - Vector
      Returns:
      Result vector
    • Dr2af

      public palTime Dr2af(double angle)
      Convert an angle in radians into degrees, arcminutes, arcseconds. WARNING: broken doesn't preserve sign in "palTime.toString()".
      Parameters:
      angle - angle in radians
      Returns:
      Time as degrees, arcminutes, arcseconds, fraction
    • Dr2tf

      public palTime Dr2tf(double angle)
      Convert an angle in radians to hours, minutes, seconds.
      Parameters:
      angle - Angle in radians
      Returns:
      Time as hours, minutes, seconds, fraction(integer)
    • Drange

      public double Drange(double angle)
      Normalize angle into range +/- pi.
      Parameters:
      angle - Angle in radians
      Returns:
      Angle expressed in the range +/- π
    • Dranrm

      public double Dranrm(double angle)
      Normalize angle into range 0-2 π.
      Parameters:
      angle - Angle in radians
      Returns:
      Angle expressed in the range 0-2 pi
    • Drverot

      public double Drverot(double phi, AngleDR r, double st)
      Velocity component in a given direction due to Earth rotation.
      Sign convention:
      The result is +ve when the observer is receding from the given point on the sky.
      Accuracy:
      The simple algorithm used assumes a spherical Earth, of a radius chosen to give results accurate to about 0.0005 km/s for observing stations at typical latitudes and heights. For applications requiring greater precision, use the routine Pvobs.
      Parameters:
      phi - Latitude of observing station (geodetic)
      r - Apparent RA,Dec (radians)
      st - local apparent sidereal time
      Returns:
      Component of Earth rotation in direction ra,da (km/s)
    • Drvgalc

      public double Drvgalc(AngleDR r2000)
      Velocity component in a given direction due to the rotation of the Galaxy.
      Sign convention:
      The result is +ve when the dynamical LSR is receding from the given point on the sky.
      Note:
      The Local Standard of Rest used here is a point in the vicinity of the Sun which is in a circular orbit around the Galactic centre. Sometimes called the "dynamical" LSR, it is not to be confused with a "kinematical" LSR, which is the mean standard of rest of star catalogues or stellar populations.
      Reference:
      The orbital speed of 220 km/s used here comes from Kerr & Lynden-Bell (1986), MNRAS, 221, p1023.
      Parameters:
      r2000 - J2000.0 mean RA,Dec (radians)
      Returns:
      Component of dynamical LSR motion in direction r2000,d2000 (km/s)
    • Drvlg

      public double Drvlg(AngleDR r2000)
      Velocity component in a given direction due to the combination of the rotation of the Galaxy and the motion of the Galaxy relative to the mean motion of the local group.
      Sign convention:
      The result is +ve when the Sun is receding from the given point on the sky.
      Reference:
      IAU trans 1976, 168, p201.
      Parameters:
      r2000 - J2000.0 mean RA,Dec (radians)
      Returns:
      Component of solar motion in direction r2000,d2000 (km/s)
    • Drvlsrd

      public double Drvlsrd(AngleDR r2000)
      Velocity component in a given direction due to the Sun's motion with respect to the dynamical Local Standard of Rest.
      Sign convention:
      The result is +ve when the Sun is receding from the given point on the sky.
      Note:
      The Local Standard of Rest used here is the "dynamical" LSR, a point in the vicinity of the Sun which is in a circular orbit around the Galactic centre. The Sun's motion with respect to the dynamical LSR is called the "peculiar" solar motion.

      There is another type of LSR, called a "kinematical" LSR. A kinematical LSR is the mean standard of rest of specified star catalogues or stellar populations, and several slightly different kinematical LSRs are in use. The Sun's motion with respect to an agreed kinematical LSR is known as the "standard" solar motion. To obtain a radial velocity correction with respect to an adopted kinematical LSR use the routine Rvlsrk.

      Reference:
      Delhaye (1965), in "Stars and Stellar Systems", vol 5, p73.
      Parameters:
      r2000 - J2000.0 mean RA,Dec (radians)
      Returns:
      Component of "peculiar" solar motion in direction R2000,D2000 (km/s)
    • Drvlsrk

      public double Drvlsrk(AngleDR r2000)
      Velocity component in a given direction due to the Sun's motion with respect to an adopted kinematic Local Standard of Rest.
      Sign convention:
      The result is +ve when the Sun is receding from the given point on the sky.
      Note:
      The Local Standard of Rest used here is one of several "kinematical" LSRs in common use. A kinematical LSR is the mean standard of rest of specified star catalogues or stellar populations. The Sun's motion with respect to a kinematical LSR is known as the "standard" solar motion.

      There is another sort of LSR, the "dynamical" LSR, which is a point in the vicinity of the Sun which is in a circular orbit around the Galactic centre. The Sun's motion with respect to the dynamical LSR is called the "peculiar" solar motion. To obtain a radial velocity correction with respect to the dynamical LSR use the routine Rvlsrd.

      Reference:
      Delhaye (1965), in "Stars and Stellar Systems", vol 5, p73.
      Parameters:
      r2000 - J2000.0 mean RA,Dec (radians)
      Returns:
      Component of "standard" solar motion in direction R2000,D2000 (km/s)
    • Ds2c6

      public Cartesian Ds2c6(Spherical s)
      Conversion of position & velocity in spherical coordinates to Cartesian coordinates.
      Parameters:
      s - Spherical coordinates (longitude, latitude, radial)
      Returns:
      Cartesian position & velocity vector
    • Ds2tp

      public AngleDR Ds2tp(AngleDR r, AngleDR rz) throws palError
      Projection of spherical coordinates onto tangent plane ('gnomonic' projection - 'standard coordinates').
      Parameters:
      r - spherical coordinates of point to be projected
      rz - spherical coordinates of tangent point
      Returns:
      rectangular coordinates on tangent plane (xi, eta)
      Throws:
      palError
    • Dt

      public double Dt(double epoch)
      Estimate the offset between dynamical time and Universal Time for a given historical epoch.
      Notes:
      1. Depending on the epoch, one of three parabolic approximations is used:
        before 979 Stephenson & Morrison's 390 BC to AD 948 model
        979 to 1708 Stephenson & Morrison's 948 to 1600 model
        after 1708 McCarthy & Babcock's post-1650 model
        The breakpoints are chosen to ensure continuity: they occur at places where the adjacent models give the same answer as each other.
      2. The accuracy is modest, with errors of up to 20 sec during the interval since 1650, rising to perhaps 30 min by 1000 BC. Comparatively accurate values from AD 1600 are tabulated in the Astronomical Almanac (see section K8 of the 1995 AA).
      3. The use of double-precision for both argument and result is purely for compatibility with other LIB time routines.
      4. The models used are based on a lunar tidal acceleration value of -26.00 arcsec per century.
      Reference:
      Explanatory Supplement to the Astronomical Almanac, ed P.K.Seidelmann, University Science Books (1992), section 2.553, p83. This contains references to the Stephenson & Morrison and McCarthy & Babcock papers.
      Parameters:
      epoch - (Julian) epoch (e.g. 1850.0)
      Returns:
      Estimate of ET-UT (after 1984, TT-UT1) at the given epoch, in seconds
    • Dtf2d

      public double Dtf2d(int ihour, int imin, double sec) throws palError
      Convert hours, minutes, seconds to days.
      Notes:
      1. The result is computed even if any of the range checks fail.
      2. The sign must be dealt with outside this routine.
      Parameters:
      ihour - Hours
      imin - Minutes
      sec - Seconds
      Returns:
      Interval in days
      Throws:
      palError - Hour, Min or Sec out of range
    • Dtf2r

      public double Dtf2r(int ihour, int imin, double sec) throws palError
      Convert hours, minutes, seconds to radians.
      Notes:
      1. The result is computed even if any of the range checks fail.
      2. The sign must be dealt with outside this routine.
      Parameters:
      ihour - Hours
      imin - Minutes
      sec - Seconds
      Returns:
      Angle in radians
      Throws:
      palError - Hour, Min or Sec out of range
    • Dtp2s

      public AngleDR Dtp2s(AngleDR x, AngleDR rz)
      Transform tangent plane coordinates into spherical.
      Parameters:
      x - Tangent plane rectangular coordinates (xi, eta)
      rz - Spherical coordinates of tangent point (ra, dec)
      Returns:
      Spherical coordinates (0-2pi,+/-pi/2)
    • Dtt

      public double Dtt(double utc)
      Increment to be applied to Coordinated Universal Time UTC to give Terrestrial Time TT (formerly Ephemeris Time ET).
      Notes:
      1. The UTC is specified to be a date rather than a time to indicate that care needs to be taken not to specify an instant which lies within a leap second. Though in most cases UTC can include the fractional part, correct behaviour on the day of a leap second can only be guaranteed up to the end of the second 23:59:59.
      2. Pre 1972 January 1 a fixed value of 10 + ET-TAI is returned.
      3. See also the routine Dt, which roughly estimates ET-UT for historical epochs.
      Parameters:
      utc - UTC date as a modified JD (JD-2400000.5)
      Returns:
      TT-UTC in seconds
    • Dvdv

      public double Dvdv(double[] va, double[] vb)
      Scalar product of two 3-vectors.
      Parameters:
      va - First vector
      vb - Second vector
      Returns:
      The scalar product va.vb
    • Dvn

      public double Dvn(double[] v, double[] uv)
      Normalizes a 3-vector also giving the modulus.
      Note:
      v and uv may be the same array.

      If the modulus of v is zero, uv is set to zero as well.

      Parameters:
      v - Vector
      uv - (Returned) Unit vector in direction of v
      Returns:
      modulus of v
    • Dvxv

      public double[] Dvxv(double[] va, double[] vb)
      Vector product of two 3-vectors.
      Note:
      the same vector may be specified more than once.
      Parameters:
      va - First vector
      vb - Second vector
      Returns:
      Vector result
    • Ecleq

      public AngleDR Ecleq(AngleDR dl, double date)
      Transformation from ecliptic coordinates to J2000.0 equatorial coordinates.
      Parameters:
      dl - ecliptic longitude and latitude (mean of date, IAU 1980 theory, radians)
      date - TDB (loosely ET) as Modified Julian Date (JD-2400000.5)
      Returns:
      J2000.0 mean (RA, Dec) (radians)
    • Ecmat

      public double[][] Ecmat(double date)
      Form the equatorial to ecliptic rotation matrix (IAU 1980 theory).
      References:
      Murray, C.A., Vectorial Astrometry, section 4.3.
      Note:
      The matrix is in the sense v[ecl] = rmat * v[equ]; the equator, equinox and ecliptic are mean of date.
      Parameters:
      date - TDB (loosely ET) as Modified Julian Date (JD-2400000.5)
      Returns:
      Rotation matrix
    • Epb

      public double Epb(double date)
      Conversion of Modified Julian Date to Besselian epoch.
      Reference:
      Lieske,J.H., 1979. Astron. Astrophys.,73,282.
      Parameters:
      date - Modified Julian Date (JD - 2400000.5)
      Returns:
      The Besselian epoch
    • Epb2d

      public double Epb2d(double epb)
      Conversion of Besselian epoch to Modified Julian Date.
      Reference:
      Lieske,J.H., 1979. Astron. Astrophys.,73,282.
      Parameters:
      epb - Besselian epoch
      Returns:
      Modified Julian Date (JD - 2400000.5)
    • Epco

      public double Epco(char k0, char k, double e)
      Convert an epoch into the appropriate form - 'B' or 'J'.
      Notes:
      1. The result is always either equal to or very close to the given epoch e. The routine is required only in applications where punctilious treatment of heterogeneous mixtures of star positions is necessary.
      2. k0 and k are not validated, and only their first characters are used, interpreted as follows:
        • If k0 and k are the same the result is e.
        • If k0 is 'B' or 'b' and k isn't, the conversion is J to B.
        • In all other cases, the conversion is B to J.
      Parameters:
      k0 - Form of result: 'B'=Besselian, 'J'=Julian
      k - Form of given epoch: 'B' or 'J'
      e - Epoch
      Returns:
      Epoch in appropriate form
    • Epj

      public double Epj(double date)
      Conversion of Modified Julian Date to Julian epoch.
      Reference:
      Lieske,J.H., 1979. Astron. Astrophys.,73,282.
      Parameters:
      date - Modified Julian Date (JD - 2400000.5)
      Returns:
      Julian epoch
    • Epj2d

      public double Epj2d(double epj)
      Conversion of Julian epoch to Modified Julian Date.
      Reference:
      Lieske,J.H., 1979. Astron. Astrophys.,73,282.
      Parameters:
      epj - Julian epoch
      Returns:
      Modified Julian Date (JD - 2400000.5)
    • Eqecl

      public AngleDR Eqecl(AngleDR d, double date)
      Transformation from J2000.0 equatorial coordinates to ecliptic coordinates.
      Parameters:
      d - J2000.0 mean RA,Dec (radians)
      date - TDB (loosely ET) as Modified Julian Date (JD-2400000.5)
      Returns:
      ecliptic longitude and latitude (mean of date, IAU 1980 theory, radians)
    • Eqeqx

      public double Eqeqx(double date)
      Equation of the equinoxes (IAU 1994, double precision).

      Greenwich apparent ST = Greenwich mean ST + equation of the equinoxes

      References:
      IAU Resolution C7, Recommendation 3 (1994) Capitaine, N. & Gontier, A.-M., Astron. Astrophys., 275, 645-650 (1993)
      Parameters:
      date - TDB (loosely ET) as Modified Julian Date (JD-2400000.5)
      Returns:
      Equation of the equinoxes (in radians)
    • Eqgal

      public Galactic Eqgal(AngleDR dr)
      Transformation from J2000.0 equatorial coordinates to IAU 1958 Galactic coordinates.
      Note:
      The equatorial coordinates are J2000.0. Use the routine slaEg50 if conversion from B1950.0 'FK4' coordinates is required.
      Reference:
      Blaauw et al, Mon.Not.R.astron.Soc.,121,123 (1960)
      Parameters:
      dr - J2000.0 (RA, Dec) (in radians)
      Returns:
      Galactic longitude and latitude (l2, b2) (in radians)
    • Etrms

      public double[] Etrms(double ep)
      Compute the e-terms (elliptic component of annual aberration) vector.
      References:
      1. Smith, C.A. et al, 1989. "The transformation of astrometric catalog systems to the equinox J2000.0". Astron.J. 97, 265.
      2. Yallop, B.D. et al, 1989. "Transformation of mean star places from FK4 B1950.0 to FK5 J2000.0 using matrices in 6-space". Astron.J. 97, 274.
      Note the use of the J2000 aberration constant (20.49552 arcsec). This is a reflection of the fact that the e-terms embodied in existing star catalogues were computed from a variety of aberration constants. Rather than adopting one of the old constants the latest value is used here.
      Parameters:
      ep - Besselian epoch
      Returns:
      E-terms as ( dx, dy, dz )
    • Evp

      public void Evp(double date, double deqx, double[] dvb, double[] dpb, double[] dvh, double[] dph)
      Barycentric and heliocentric velocity and position of the Earth.
      Accuracy:
      The maximum deviations from the JPL DE96 ephemeris are as follows:
      barycentric velocity 42 cm/s
      barycentric position 6900 km
      heliocentric velocity 42 cm/s
      heliocentric position 1600 km

      This routine is adapted from the BARVEL and BARCOR Fortran subroutines of P.Stumpff, which are described in Astron. Astrophys. Suppl. Ser. 41, 1-8 (1980). The present routine uses double precision throughout; most of the other changes are essentially cosmetic and do not affect the results. However, some adjustments have been made so as to give results that refer to the new (IAU 1976 "FK5") equinox and precession, although the differences these changes make relative to the results from Stumpff's original "FK4" version are smaller than the inherent accuracy of the algorithm. One minor shortcoming in the original routines that has not been corrected is that better numerical accuracy could be achieved if the various polynomial evaluations were nested. Note also that one of Stumpff's precession constants differs by 0.001 arcsec from the value given in the Explanatory Supplement to the A.E.

      (Units are AU/s for velocity and AU for position)

      Parameters:
      date - TDB (loosely ET) as a Modified Julian Date (JD-2400000.5)
      deqx - Julian epoch (e.g. 2000.0) of mean equator and equinox of the vectors returned. If deqx <= 0.0, all vectors are referred to the mean equator and equinox (FK5) of epoch date
      dvb - (Returned) barycentric velocity
      dpb - (Returned) barycentric position
      dvh - (Returned) heliocentric velocity
      dph - (Returned) heliocentric position
    • Fk425

      public Stardata Fk425(Stardata s1950)
      Convert B1950.0 FK4 star data to J2000.0 FK5.

      This routine converts stars from the old, Bessel-Newcomb, FK4 system to the new, IAU 1976, FK5, Fricke system. The precepts of Smith et al (Ref 1) are followed, using the implementation by Yallop et al (Ref 2) of a matrix method due to Standish. Kinoshita's development of Andoyer's post-Newcomb precession is used. The numerical constants from Seidelmann et al (Ref 3) are used canonically.

      Notes:
      1. The proper motions in RA are dRA/dt rather than cos(Dec)*dRA/dt, and are per year rather than per century.
      2. Conversion from Besselian epoch 1950.0 to Julian epoch 2000.0 only is provided for. Conversions involving other epochs will require use of the appropriate precession, proper motion, and E-terms routines before and/or after FK425 is called.
      3. In the FK4 catalogue the proper motions of stars within 10 degrees of the poles do not embody the differential E-term effect and should, strictly speaking, be handled in a different manner from stars outside these regions. However, given the general lack of homogeneity of the star data available for routine astrometry, the difficulties of handling positions that may have been determined from astrometric fields spanning the polar and non-polar regions, the likelihood that the differential E-terms effect was not taken into account when allowing for proper motion in past astrometry, and the undesirability of a discontinuity in the algorithm, the decision has been made in this routine to include the effect of differential E-terms on the proper motions for all stars, whether polar or not. At epoch 2000, and measuring on the sky rather than in terms of dRA, the errors resulting from this simplification are less than 1 milliarcsecond in position and 1 milliarcsecond per century in proper motion.
      References:
      1. Smith, C.A. et al, 1989. "The transformation of astrometric catalog systems to the equinox J2000.0". Astron.J. 97, 265.
      2. Yallop, B.D. et al, 1989. "Transformation of mean star places from FK4 B1950.0 to FK5 J2000.0 using matrices in 6-space". Astron.J. 97, 274.
      3. Seidelmann, P.K. (ed), 1992. "Explanatory Supplement to the Astronomical Almanac", ISBN 0-935702-68-7.
      Parameters:
      s1950 - B1950.0 RA,dec (rad), proper motions (rad/trop.yr), parallax (arcsec), radial velocity (km/s, +ve = moving away)
      Returns:
      J2000.0 RA,dec (rad), proper motions (rad/trop.yr), parallax (arcsec), radial velocity (km/s, +ve = moving away)
    • Fk45z

      public AngleDR Fk45z(AngleDR r1950, double bepoch)
      Convert B1950.0 FK4 star data to J2000.0 FK5 assuming zero proper motion in the FK5 frame (double precision).

      This routine converts stars from the old, Bessel-Newcomb, FK4 system to the new, IAU 1976, FK5, Fricke system, in such a way that the FK5 proper motion is zero. Because such a star has, in general, a non-zero proper motion in the FK4 system, the routine requires the epoch at which the position in the FK4 system was determined.

      The method is from Appendix 2 of Ref 1, but using the constants of Ref 4.

      Notes:
      1. The epoch BEPOCH is strictly speaking Besselian, but if a Julian epoch is supplied the result will be affected only to a negligible extent.
      2. Conversion from Besselian epoch 1950.0 to Julian epoch 2000.0 only is provided for. Conversions involving other epochs will require use of the appropriate precession, proper motion, and E-terms routines before and/or after FK45Z is called.
      3. In the FK4 catalogue the proper motions of stars within 10 degrees of the poles do not embody the differential E-term effect and should, strictly speaking, be handled in a different manner from stars outside these regions. However, given the general lack of homogeneity of the star data available for routine astrometry, the difficulties of handling positions that may have been determined from astrometric fields spanning the polar and non-polar regions, the likelihood that the differential E-terms effect was not taken into account when allowing for proper motion in past astrometry, and the undesirability of a discontinuity in the algorithm, the decision has been made in this routine to include the effect of differential E-terms on the proper motions for all stars, whether polar or not. At epoch 2000, and measuring on the sky rather than in terms of dRA, the errors resulting from this simplification are less than 1 milliarcsecond in position and 1 milliarcsecond per century in proper motion.
      References:
      1. Aoki,S., et al, 1983. Astron. Astrophys., 128, 263.
      2. Smith, C.A. et al, 1989. "The transformation of astrometric catalog systems to the equinox J2000.0". Astron.J. 97, 265.
      3. Yallop, B.D. et al, 1989. "Transformation of mean star places from FK4 B1950.0 to FK5 J2000.0 using matrices in 6-space". Astron.J. 97, 274.
      4. Seidelmann, P.K. (ed), 1992. "Explanatory Supplement to the Astronomical Almanac", ISBN 0-935702-68-7.
      Parameters:
      r1950 - B1950.0 FK4 RA,Dec at epoch (rad)
      bepoch - Besselian epoch (e.g. 1979.3)
      Returns:
      J2000.0 FK5 RA,Dec (rad)
    • Fk524

      public Stardata Fk524(Stardata j2000)
      Convert J2000.0 FK5 star data to B1950.0 FK4.

      This routine converts stars from the new, IAU 1976, FK5, Fricke system, to the old, Bessel-Newcomb, FK4 system. The precepts of Smith et al (Ref 1) are followed, using the implementation by Yallop et al (Ref 2) of a matrix method due to Standish. Kinoshita's development of Andoyer's post-Newcomb precession is used. The numerical constants from Seidelmann et al (Ref 3) are used canonically.

      Notes:
      1. The proper motions in RA are dRA/dt rather than cos(Dec)*dRA/dt, and are per year rather than per century.
      2. Note that conversion from Julian epoch 2000.0 to Besselian epoch 1950.0 only is provided for. Conversions involving other epochs will require use of the appropriate precession, proper motion, and E-terms routines before and/or after FK524 is called.
      3. In the FK4 catalogue the proper motions of stars within 10 degrees of the poles do not embody the differential E-term effect and should, strictly speaking, be handled in a different manner from stars outside these regions. However, given the general lack of homogeneity of the star data available for routine astrometry, the difficulties of handling positions that may have been determined from astrometric fields spanning the polar and non-polar regions, the likelihood that the differential E-terms effect was not taken into account when allowing for proper motion in past astrometry, and the undesirability of a discontinuity in the algorithm, the decision has been made in this routine to include the effect of differential E-terms on the proper motions for all stars, whether polar or not. At epoch 2000, and measuring on the sky rather than in terms of dRA, the errors resulting from this simplification are less than 1 milliarcsecond in position and 1 milliarcsecond per century in proper motion.
      References:
      1. Smith, C.A. et al, 1989. "The transformation of astrometric catalog systems to the equinox J2000.0". Astron.J. 97, 265.
      2. Yallop, B.D. et al, 1989. "Transformation of mean star places from FK4 B1950.0 to FK5 J2000.0 using matrices in 6-space". Astron.J. 97, 274.
      3. Seidelmann, P.K. (ed), 1992. "Explanatory Supplement to the Astronomical Almanac", ISBN 0-935702-68-7.
      Parameters:
      j2000 - J2000.0 RA,Dec (rad), J2000.0 proper motions (rad/Jul.yr), parallax (arcsec), radial velocity (km/s, +ve = moving away)
      Returns:
      B1950.0 RA,Dec (rad), proper motions (rad/trop.yr), parallax (arcsec), radial velocity (km/s, +ve = moving away)
    • Fk54z

      public Stardata Fk54z(AngleDR r2000, double bepoch)
      Convert a J2000.0 FK5 star position to B1950.0 FK4 assuming zero proper motion and parallax.

      This routine converts star positions from the new, IAU 1976, FK5, Fricke system to the old, Bessel-Newcomb, FK4 system.

      Notes:
      1. The proper motion in RA is dRA/dt rather than cos(Dec)*dRA/dt.
      2. Conversion from Julian epoch 2000.0 to Besselian epoch 1950.0 only is provided for. Conversions involving other epochs will require use of the appropriate precession routines before and after this routine is called.
      3. Unlike in the FK524 routine, the FK5 proper motions, the parallax and the radial velocity are presumed zero.
      4. It is the intention that FK5 should be a close approximation to an inertial frame, so that distant objects have zero proper motion; such objects have (in general) non-zero proper motion in FK4, and this routine returns those fictitious proper motions.
      5. The position returned by this routine is in the B1950 reference frame but at Besselian epoch bepoch. For comparison with catalogues the bepoch argument will frequently be 1950.0.
      Parameters:
      r2000 - J2000.0 FK5 RA,Dec (rad)
      bepoch - Besselian epoch (e.g. 1950)
      Returns:
      B1950.0 FK4 RA,Dec (rad) at epoch BEPOCH, B1950.0 FK4 proper motions (rad/trop.yr)
    • Galeq

      public AngleDR Galeq(Galactic gl)
      Transformation from IAU 1958 Galactic coordinates to J2000.0 equatorial coordinates.
      Note:
      The equatorial coordinates are J2000.0. Use the routine Ge50 if conversion to B1950.0 'FK4' coordinates is required.
      Reference:
      Blaauw et al, Mon.Not.R.astron.Soc.,121,123 (1960)
      Parameters:
      gl - Galactic longitude and latitude l2, b2
      Returns:
      J2000.0 RA, dec
    • Galsup

      public Galactic Galsup(Galactic gl)
      Transformation from IAU 1958 Galactic coordinates to De Vaucouleurs supergalactic coordinates.
      References:
      1. De Vaucouleurs, De Vaucouleurs, & Corwin, Second reference catalogue of bright galaxies, U. Texas, page 8.
      2. Systems & Applied Sciences Corp., Documentation for the machine-readable version of the above catalogue, contract NAS 5-26490.
      (These two references give different values for the Galactic longitude of the Supergalactic origin. Both are wrong; the correct value is l2 = 137.37.)
      Parameters:
      gl - Galactic longitude and latitude l2,b2
      Returns:
      Supergalactic longitude and latitude
    • Geoc

      public double[] Geoc(double p, double h)
      Convert geodetic position to geocentric.
      Notes:
      1. Geocentric latitude can be obtained by evaluating atan2(z,r).
      2. IAU 1976 constants are used.
      Reference:
      Green,R.M., Spherical Astronomy, CUP 1985, p98.
      Parameters:
      p - latitude (geodetic, radians)
      h - height above reference spheroid (geodetic, metres)
      Returns:
      distance from Earth axis (AU), distance from plane of Earth equator (AU)
    • Gmst

      public double Gmst(double ut1)
      Conversion from Universal Time to Sidereal Time.

      The IAU 1982 expression (see page S15 of the 1984 Astronomical Almanac) is used, but rearranged to reduce rounding errors. This expression is always described as giving the GMST at 0 hours UT. In fact, it gives the difference between the GMST and the UT, which happens to equal the GMST (modulo 24 hours) at 0 hours UT each day. In this routine, the entire UT is used directly as the argument for the standard formula, and the fractional part of the UT is added separately; note that the factor 1.0027379... does not appear.

      See also the routine Gmsta, which delivers better numerical precision by accepting the UT date and time as separate arguments.

      Parameters:
      ut1 - Universal Time (strictly UT1) expressed as Modified Julian Date (JD-2400000.5)
      Returns:
      Greenwich Mean Sidereal Time (radians)
    • Kbj

      public char Kbj(int jb, double e) throws palError
      Select epoch prefix 'B' or 'J'.
      Parameters:
      jb - Dbjin prefix status: 0=none, 1='B', 2='J'
      e - epoch - Besselian or Julian
      Returns:
      'B' or 'J'
      Throws:
      palError - Illegal prefix

      If jb=0, B is assumed for e < 1984.0, otherwise J.

    • Map

      public AngleDR Map(Stardata sd, double epq, double date)
      Transform star RA,Dec from mean place to geocentric apparent.

      The reference frames and timescales used are post IAU 1976.

      Notes:
      1. eq is the Julian epoch specifying both the reference frame and the epoch of the position - usually 2000. For positions where the epoch and equinox are different, use the routine Pm to apply proper motion corrections before using this routine.
      2. The distinction between the required TDB and TT is always negligible. Moreover, for all but the most critical applications UTC is adequate.
      3. The proper motions in RA are dRA/dt rather than cos(Dec)*dRA/dt.
      4. This routine may be wasteful for some applications because it recomputes the Earth position/velocity and the precession- nutation matrix each time, and because it allows for parallax and proper motion. Where multiple transformations are to be carried out for one epoch, a faster method is to call the Mappa routine once and then either the Mapqk routine (which includes parallax and proper motion) or Mapqkz (which assumes zero parallax and proper motion).
      5. The accuracy is sub-milliarcsecond, limited by the precession-nutation model (IAU 1976 precession, Shirai & Fukushima 2001 forced nutation and precession corrections).
      6. The accuracy is further limited by the routine Evp, called by Mappa, which computes the Earth position and velocity using the methods of Stumpff. The maximum error is about 0.3 mas.
      References:
      1984 Astronomical Almanac, pp B39-B41. (also Lederle & Schwan, Astron. Astrophys. 134, 1-6, 1984)
      Parameters:
      sd - mean RA,Dec (rad) proper motions (RA,Dec changes per Julian year), parallax (arcsec), radial velocity (km/sec, +ve if receding)
      epq - Epoch and equinox of star data (Julian)
      date - TDB for apparent place (JD-2400000.5)
      Returns:
      Apparent RA,Dec (rad)
    • Mappa

      public AMParams Mappa(double eq, double date)
      Compute star-independent parameters in preparation for conversions between mean place and geocentric apparent place.

      The parameters produced by this routine are required in the parallax, light deflection, aberration, and precession/nutation parts of the mean/apparent transformations.

      The reference frames and timescales used are post IAU 1976.

      Notes:
      1. For date, the distinction between the required TDB and TT is always negligible. Moreover, for all but the most critical applications UTC is adequate.
      2. The vectors amprms(1-3) and amprms(4-6) are referred to the mean equinox and equator of epoch eq.
      3. The parameters AMPRMS produced by this routine are used by Ampqk, Mapqk and Mapqkz.
      4. The accuracy is sub-milliarcsecond, limited by the precession-nutation model (IAU 1976 precession, Shirai & Fukushima 2001 forced nutation and precession corrections).
      5. A further limit to the accuracy of routines using the parameter array AMPRMS is imposed by the routine Evp, used here to compute the Earth position and velocity by the methods of Stumpff. The maximum error in the resulting aberration corrections is about 0.3 milliarcsecond.
      References:
      1984 Astronomical Almanac, pp B39-B41. (also Lederle & Schwan, Astron. Astrophys. 134, 1-6, 1984)
      Parameters:
      eq - epoch of mean equinox to be used (Julian)
      date - TDB (JD-2400000.5)
      Returns:
      star-independent mean-to-apparent parameters
    • Mapqk

      public AngleDR Mapqk(Stardata s, AMParams amprms)
      Quick mean to apparent place: transform a star RA,Dec from mean place to geocentric apparent place, given the star-independent parameters.

      Use of this routine is appropriate when efficiency is important and where many star positions, all referred to the same equator and equinox, are to be transformed for one epoch. The star-independent parameters can be obtained by calling the Mappa routine.

      If the parallax and proper motions are zero the Mapqkz routine can be used instead.

      The reference frames and timescales used are post IAU 1976.

      Notes:
      1. The vectors amprms(1-3) and amprms(4-6) are referred to the mean equinox and equator of epoch eq.
      2. Strictly speaking, the routine is not valid for solar-system sources, though the error will usually be extremely small. However, to prevent gross errors in the case where the position of the Sun is specified, the gravitational deflection term is restrained within about 920 arcsec of the centre of the Sun's disc. The term has a maximum value of about 1.85 arcsec at this radius, and decreases to zero as the centre of the disc is approached.
      References:
      1984 Astronomical Almanac, pp B39-B41. (also Lederle & Schwan, Astron. Astrophys. 134, 1-6, 1984)
      Parameters:
      s - Mean RA,Dec (rad), proper motions (RA,Dec changes per Julian year), parallax (arcsec), radial velocity (km/sec, +ve if receding)
      amprms - star-independent mean-to-apparent parameters
      Returns:
      Apparent RA,Dec (rad)
    • Mapqkz

      public AngleDR Mapqkz(AngleDR rm, AMParams amprms)
      Quick mean to apparent place: transform a star RA,dec from mean place to geocentric apparent place, given the star-independent parameters, and assuming zero parallax and proper motion.

      Use of this routine is appropriate when efficiency is important and where many star positions, all with parallax and proper motion either zero or already allowed for, and all referred to the same equator and equinox, are to be transformed for one epoch. The star-independent parameters can be obtained by calling the Mappa routine.

      The corresponding routine for the case of non-zero parallax and proper motion is Mapqk.

      The reference frames and timescales used are post IAU 1976.

      Notes:
      1. The vectors amprms(1-3) and amprms(4-6) are referred to the mean equinox and equator of epoch eq.
      2. Strictly speaking, the routine is not valid for solar-system sources, though the error will usually be extremely small. However, to prevent gross errors in the case where the position of the Sun is specified, the gravitational deflection term is restrained within about 920 arcsec of the centre of the Sun's disc. The term has a maximum value of about 1.85 arcsec at this radius, and decreases to zero as the centre of the disc is approached.
      References:
      1984 Astronomical Almanac, pp B39-B41. (also Lederle & Schwan, Astron. Astrophys. 134, 1-6, 1984)
      Parameters:
      rm - Mean RA,dec (rad)
      amprms - Star-independent mean-to-apparent parameters
      Returns:
      Apparent RA,dec (rad)
    • Nut

      public double[][] Nut(double date)
      Form the matrix of nutation for a given date (IAU 1980 theory).
      Parameters:
      date - TDB (loosely ET) as Modified Julian Date (=JD-2400000.5)
      Returns:
      Nutation matrix

      The matrix is in the sense v(true) = rmatn * v(mean) .

    • Nutc

      public double[] Nutc(double date)
      Nutation: longitude & obliquity components and mean obliquity (IAU 1980 theory).
      Notes:
      1. The routine predicts forced nutation (but not free core nutation) plus corrections to the IAU 1976 precession model.
      2. Earth attitude predictions made by combining the present nutation model with IAU 1976 precession are accurate to 1 mas (with respect to the ICRF) for a few decades around 2000.
      3. The slaNutc80 routine is the equivalent of the present routine but using the IAU 1980 nutation theory. The older theory is less accurate, leading to errors as large as 350 mas over the interval 1900-2100, mainly because of the error in the IAU 1976 precession.
      References:
      1. Shirai, T. & Fukushima, T., Astron.J. 121, 3270-3283 (2001).
      2. Fukushima, T., 1991, Astron.Astrophys. 244, L11 (1991).
      3. Simon, J. L., Bretagnon, P., Chapront, J., Chapront-Touze, M., Francou, G. & Laskar, J., Astron.Astrophys. 282, 663 (1994).
      Parameters:
      date - TDB (loosely ET) as Modified Julian Date (JD-2400000.5)
      Returns:
      Nutation in longitude, obliquity, and mean obliquity
    • Obs

      public Observatory Obs(int n)
      Parameters of selected groundbased observing stations.
      Parameters:
      n - Number specifying observing station
      Returns:
      Name of specified observing station, longitude (radians, West +ve), geodetic latitude (radians, North +ve), height above sea level (metres)
    • Obs

      public Observatory Obs(String id)
      Parameters of selected groundbased observing stations.
      Parameters:
      id - Identifier specifying observing station
      Returns:
      Name of specified observing station, longitude (radians, West +ve), geodetic latitude (radians, North +ve), height above sea level (metres)
    • Pm

      public AngleDR Pm(AngleDR r0, double[] pm, double px, double rv, double ep0, double ep1)
      Apply corrections for proper motion to a star RA,Dec.
      Notes:
      1. The proper motions in RA are dRA/dt rather than cos(Dec)*dRA/dt, and are in the same coordinate system as R0,D0.
      2. If the available proper motions are pre-FK5 they will be per tropical year rather than per Julian year, and so the epochs must both be Besselian rather than Julian. In such cases, a scaling factor of 365.2422D0/365.25D0 should be applied to the radial velocity before use.
      Parameters:
      r0 - RA,Dec at epoch ep0 (rad)
      pm - proper motions: RA,Dec changes per year of epoch
      px - parallax (arcsec)
      rv - radial velocity (km/sec, +ve if receding)
      ep0 - start epoch in years (e.g Julian epoch)
      ep1 - end epoch in years (same system as ep0)
      Returns:
      RA,Dec at epoch ep1 (rad)
    • Prebn

      public double[][] Prebn(double bep0, double bep1)
      Generate the matrix of precession between two epochs, using the old, pre-IAU1976, Bessel-Newcomb model, using Kinoshita's formulation (double precision).

      The matrix is in the sense v(bep1) = rmatp * v(bep0)

      Reference:
      Kinoshita, H. (1975) 'Formulas for precession', SAO Special Report No. 364, Smithsonian Institution Astrophysical Observatory, Cambridge, Massachusetts.
      Parameters:
      bep0 - Beginning Besselian epoch
      bep1 - Ending Besselian epoch
      Returns:
      Precession matrix
    • Prec

      public double[][] Prec(double ep0, double ep1)
      Form the matrix of precession between two epochs (IAU 1976, FK5).
      Notes:
      1. The epochs are TDB (loosely ET) Julian epochs.
      2. The matrix is in the sense v(ep1) = rmatp * v(ep0).
      3. Though the matrix method itself is rigorous, the precession angles are expressed through canonical polynomials which are valid only for a limited time span. There are also known errors in the IAU precession rate. The absolute accuracy of the present formulation is better than 0.1 arcsec from 1960AD to 2040AD, better than 1 arcsec from 1640AD to 2360AD, and remains below 3 arcsec for the whole of the period 500BC to 3000AD. The errors exceed 10 arcsec outside the range 1200BC to 3900AD, exceed 100 arcsec outside 4200BC to 5600AD and exceed 1000 arcsec outside 6800BC to 8200AD. The LIB routine Precl implements a more elaborate model which is suitable for problems spanning several thousand years.
      References:
      1. Lieske,J.H., 1979. Astron. Astrophys.,73,282. equations (6) & (7), p283.
      2. Kaplan,G.H., 1981. USNO circular no. 163, pa2.
      Parameters:
      ep0 - Beginning epoch
      ep1 - Ending epoch
      Returns:
      Precession matrix
    • Preces

      public AngleDR Preces(String sys, double ep0, double ep1, AngleDR d)
      Precession - either FK4 (Bessel-Newcomb, pre-IAU1976) or FK5 (Fricke, post-IAU1976) as required.
      Notes:
      1. The epochs are Besselian if sys='FK4' and Julian if 'FK5'. For example, to precess coordinates in the old system from equinox 1900.0 to 1950.0 the call would be: Preces ( "FK4", 1900.0, 1950.0, &ra, &dc )
      2. This routine will not correctly convert between the old and the new systems - for example conversion from B1950 to J2000. For these purposes see Fk425, Fk524, Fk45z and Fk54z.
      3. If an invalid sys is supplied, values of -99.0,-99.0 will be returned for both ra and dc.
      Parameters:
      sys - Precession to be applied: "FK4" or "FK5"
      ep0 - Starting epoch
      ep1 - Ending epoch
      d - RA,Dec, mean equator & equinox of epoch ep0
      Returns:
      RA,Dec, mean equator & equinox of epoch ep1
    • Precl

      public double[][] Precl(double ep0, double ep1)
      Form the matrix of precession between two epochs, using the model of Simon et al (1994), which is suitable for long periods of time.
      Notes:
      1. The epochs are TDB (loosely ET) Julian epochs.
      2. The matrix is in the sense v(ep1) = rmatp * v(ep0).
      3. The absolute accuracy of the model is limited by the uncertainty in the general precession, about 0.3 arcsec per 1000 years. The remainder of the formulation provides a precision of 1 mas over the interval from 1000AD to 3000AD, 0.1 arcsec from 1000BC to 5000AD and 1 arcsec from 4000BC to 8000AD.
      Reference:
      Simon, J.L., et al., 1994. Astron.Astrophys., 282, 663-683.
      Parameters:
      ep0 - Beginning epoch
      ep1 - Ending epoch
      Returns:
      Precession matrix
    • Prenut

      public double[][] Prenut(double epoch, double date)
      Form the matrix of precession and nutation (IAU 1976/1980/FK5).
      Notes:
      1. The epoch and date are TDB (loosely ET).
      2. The matrix is in the sense v(true) = rmatpn * v(mean).
      Parameters:
      epoch - Julian epoch for mean coordinates
      date - Modified Julian Date (JD-2400000.5) for true coordinates
      Returns:
      Combined precession/nutation matrix
    • Refco

      public double[] Refco(double hm, double tdk, double pmb, double rh, double wl, double phi, double tlr, double eps)
      Determine constants A and B in atmospheric refraction model dz = A tan z + B tan^3 z.

      z is the "observed" zenith distance (i.e. affected by refraction) and dz is what to add to z to give the "topocentric" (i.e. in vacuo) zenith distance.

      Notes:
      1. Typical values for the tlr and eps arguments might be 0.0065 and 1e-10 respectively.
      2. The radio refraction is chosen by specifying wl > 100 micrometres.
      3. The routine is a slower but more accurate alternative to the Refcoq routine. The constants it produces give perfect agreement with Refro at zenith distances arctan(1) (45 deg) and arctan(4) (about 76 deg). It achieves 0.5 arcsec accuracy for ZD < 80 deg, 0.01 arcsec accuracy for ZD < 60 deg, and 0.001 arcsec accuracy for ZD < 45 deg.
      Parameters:
      hm - Height of the observer above sea level (metre)
      tdk - Ambient temperature at the observer (deg k)
      pmb - Pressure at the observer (millibar)
      rh - Relative humidity at the observer (range 0-1)
      wl - Effective wavelength of the source (micrometre)
      phi - Latitude of the observer (radian, astronomical)
      tlr - Temperature lapse rate in the troposphere (degk/metre)
      eps - Precision required to terminate iteration (radian)
      Returns:
      tan z coefficient (radian), tan^3 z coefficient (radian)
    • Refro

      public double Refro(double zobs, double hm, double tdk, double pmb, double rh, double wl, double phi, double tlr, double eps)
      Atmospheric refraction for radio and optical/IR wavelengths.
      Notes:
      1. A suggested value for the tlr argument is 0.0065. The refraction is significantly affected by tlr, and if studies of the local atmosphere have been carried out a better tlr value may be available.
      2. A suggested value for the eps argument is 1e-8. The result is usually at least two orders of magnitude more computationally precise than the supplied eps value.
      3. The routine computes the refraction for zenith distances up to and a little beyond 90 deg using the method of Hohenkerk and Sinclair (NAO Technical Notes 59 and 63, subsequently adopted in the Explanatory Supplement, 1992 edition - see section 3.281).
      4. The C code is an adaptation of the Fortran optical/IR refraction subroutine AREF of C.Hohenkerk (HMNAO, September 1984), with extensions to support the radio case. The following modifications to the original HMNAO optical/IR refraction algorithm have been made:
        • The angle arguments have been changed to radians.
        • Any value of zobs is allowed (see note 6, below).
        • Other argument values have been limited to safe values.
        • Murray's values for the gas constants have been used (Vectorial Astrometry, Adam Hilger, 1983).
        • The numerical integration phase has been rearranged for extra clarity.
        • A better model for Ps(T) has been adopted (taken from Gill, Atmosphere-Ocean Dynamics, Academic Press, 1982).
        • More accurate expressions for Pwo have been adopted (again from Gill 1982).
        • Provision for radio wavelengths has been added using expressions devised by A.T.Sinclair, RGO (private communication 1989), based on the Essen & Froome refractivity formula adopted in Resolution 1 of the 13th International Geodesy Association General Assembly (Bulletin Geodesique 70 p390, 1963).
        • Various small changes have been made to gain speed.
        None of the changes significantly affects the optical/IR results with respect to the algorithm given in the 1992 Explanatory Supplement. For example, at 70 deg zenith distance the present routine agrees with the ES algorithm to better than 0.05 arcsec for any reasonable combination of parameters. However, the improved water-vapour expressions do make a significant difference in the radio band, at 70 deg zenith distance reaching almost 4 arcsec for a hot, humid, low-altitude site during a period of low pressure.
      5. The radio refraction is chosen by specifying wl > 100 micrometres. Because the algorithm takes no account of the ionosphere, the accuracy deteriorates at low frequencies, below about 30 MHz.
      6. Before use, the value of zobs is expressed in the range +/- pi. If this ranged zobs is -ve, the result ref is computed from its absolute value before being made -ve to match. In addition, if it has an absolute value greater than 93 deg, a fixed ref value equal to the result for zobs = 93 deg is returned, appropriately signed.
      7. As in the original Hohenkerk and Sinclair algorithm, fixed values of the water vapour polytrope exponent, the height of the tropopause, and the height at which refraction is negligible are used.
      8. The radio refraction has been tested against work done by Iain Coulson, JACH, (private communication 1995) for the James Clerk Maxwell Telescope, Mauna Kea. For typical conditions, agreement at the 0.1 arcsec level is achieved for moderate ZD, worsening to perhaps 0.5-1.0 arcsec at ZD 80 deg. At hot and humid sea-level sites the accuracy will not be as good.
      9. It should be noted that the relative humidity rh is formally defined in terms of "mixing ratio" rather than pressures or densities as is often stated. It is the mass of water per unit mass of dry air divided by that for saturated air at the same temperature and pressure (see Gill 1982).
      Parameters:
      zobs - Observed zenith distance of the source (radian)
      hm - Height of the observer above sea level (metre)
      tdk - Ambient temperature at the observer (deg K)
      pmb - Pressure at the observer (millibar)
      rh - Relative humidity at the observer (range 0-1)
      wl - Effective wavelength of the source (micrometre)
      phi - Latitude of the observer (radian, astronomical)
      tlr - Tropospheric lapse rate (degK/metre)
      eps - Precision required to terminate iteration (radian)
      Returns:
      Refraction: in vacuo ZD minus observed ZD (radian)
    • Subet

      public AngleDR Subet(AngleDR rc, double eq)
      Remove the e-terms (elliptic component of annual aberration) from a pre IAU 1976 catalogue RA,Dec to give a mean place.
      Explanation:
      Most star positions from pre-1984 optical catalogues (or derived from astrometry using such stars) embody the e-terms. This routine converts such a position to a formal mean place (allowing, for example, comparison with a pulsar timing position).
      Reference:
      Explanatory Supplement to the Astronomical Ephemeris, section 2D, page 48.
      Parameters:
      rc - RA,Dec (radians) with e-terms included
      eq - Besselian epoch of mean equator and equinox
      Returns:
      RA,Dec (radians) without e-terms
    • Supgal

      public Galactic Supgal(Galactic ds)
      Transformation from De Vaucouleurs supergalactic coordinates to IAU 1958 Galactic coordinates.
      References:
      1. De Vaucouleurs, De Vaucouleurs, & Corwin, Second Reference Catalogue of Bright Galaxies, U. Texas, page 8.
      2. Systems & Applied Sciences Corp., Documentation for the machine-readable version of the above catalogue, contract NAS 5-26490.

      (These two references give different values for the Galactic longitude of the supergalactic origin. Both are wrong; the correct value is l2=137.37.)

      Parameters:
      ds - Supergalactic longitude and latitude
      Returns:
      Galactic longitude and latitude l2,b2
    • Zd

      public double Zd(double ha, double dec, double phi)
      HA, Dec to Zenith Distance.
      Notes:
      1. The latitude must be geodetic. In critical applications, corrections for polar motion should be applied.
      2. In some applications it will be important to specify the correct type of hour angle and declination in order to produce the required type of zenith distance. In particular, it may be important to distinguish between the zenith distance as affected by refraction, which would require the "observed" HA,Dec, and the zenith distance in vacuo, which would require the "topocentric" HA,Dec. If the effects of diurnal aberration can be neglected, the "apparent" HA,Dec may be used instead of the topocentric HA,Dec.
      3. No range checking of arguments is done.
      4. In applications which involve many zenith distance calculations, rather than calling the present routine it will be more efficient to use inline code, having previously computed fixed terms such as sine and cosine of latitude, and perhaps sine and cosine of declination.
      Parameters:
      ha - Hour Angle in radians
      dec - Declination in radians
      phi - Observatory latitude in radians
      Returns:
      Zenith distance in the range 0 to pi