Class Pal
- 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
FieldsModifier and TypeFieldDescriptionstatic final double
Astronomical unit to kilometersstatic final double
Light time for 1 AU (sec)static final double
Speed of light (AU per day)static final double
Seconds in a daystatic final double
Nominal mean sidereal speed of Earth equator in km/s (the actual value is about 0.4651)int
Flag for additional status informationstatic final double
Gravitational radius of the Sun x 2: (2*mu/c**2, au)static final double
Degrees to radiansstatic final double
Ratio between solar and sidereal timeint
Current Status flagstatic final double
Km/s to AU/year -
Constructor Summary
Constructors -
Method Summary
Modifier and TypeMethodDescriptionAdd the e-terms (elliptic component of annual aberration) to a pre IAU 1976 mean place to conform to the old catalogue convention.Convert star RA,Dec from geocentric apparent to mean place.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
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
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[][]
Form a rotation matrix from the Euler angles - three successive rotations about specified Cartesian axes.double
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
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
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.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.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.Convert B1950.0 FK4 star data to J2000.0 FK5.Convert B1950.0 FK4 star data to J2000.0 FK5 assuming zero proper motion in the FK5 frame (double precision).Convert J2000.0 FK5 star data to B1950.0 FK4.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'.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.Quick mean to apparent place: transform a star RA,Dec from mean place to geocentric apparent place, given the star-independent parameters.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.Parameters of selected groundbased observing stations.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).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.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.
-
Field Details
-
GR2
public static final double GR2Gravitational radius of the Sun x 2: (2*mu/c**2, au)- See Also:
-
VF
public static final double VFKm/s to AU/year- See Also:
-
C
public static final double CSpeed of light (AU per day)- See Also:
-
SOLSID
public static final double SOLSIDRatio between solar and sidereal time- See Also:
-
ESPEED
public static final double ESPEEDNominal mean sidereal speed of Earth equator in km/s (the actual value is about 0.4651)- See Also:
-
D2S
public static final double D2SSeconds in a day- See Also:
-
AUKM
public static final double AUKMAstronomical unit to kilometers- See Also:
-
AUSEC
public static final double AUSECLight time for 1 AU (sec)- See Also:
-
R2D
public static final double R2DDegrees to radians- See Also:
-
Status
public int StatusCurrent Status flag -
Flag
public int FlagFlag for additional status information
-
-
Constructor Details
-
Pal
public Pal()
-
-
Method Details
-
Addet
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-termseq
- Besselian epoch of mean equator and equinox- Returns:
- RA,dec (radians) with e-terms included
-
Amp
Convert star RA,Dec from geocentric apparent to mean place.The mean coordinate system is the post IAU 1976 system, loosely called FK5.
- Notes:
-
- The distinction between the required TDB and TT is always negligible. Moreover, for all but the most critical applications UTC is adequate.
- 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.
- 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.
- The accuracy is sub-milliarcsecond, limited by the precession-nutation model (IAU 1976 precession, Shirai & Fukushima 2001 forced nutation and precession corrections).
- 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
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:
- 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.
- 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.
- 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.
- 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.
- 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.
- 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. - 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
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
Gregorian calendar to Modified Julian Date.(Includes century default feature: use Cldj for years before 100AD.)
- Parameters:
iy
- Year in Gregorian calendarim
- Month in Gregorian calendarid
- Day in Gregorian calendar- Returns:
- Modified Julian Date (JD-2400000.5) for 0 hrs
- Throws:
palError
- if bad day, month or year0 = 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
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 calendarim
- Month in Gregorian calendarid
- Day in Gregorian calendar- Returns:
- Modified Julian Date (JD-2400000.5) for 0 hrs
- Throws:
palError
- if bad day, month or year0 = OK 1 = bad year (MJD not computed) 2 = bad month (MJD not computed) 3 = bad day (MJD computed)
-
Daf2r
Convert degrees, arcminutes, arcseconds to radians.- Notes:
- The result is computed even if any of the range checks fail.
- The sign must be dealt with outside this routine.
- Parameters:
ideg
- Degreesiamin
- Arcminutesasec
- Arcseconds- Returns:
- Angle in radians
- Throws:
palError
- degrees, arcmins or arcsecs out of rangeStatus 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:
-
- 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.
- For epochs from 1961 January 1 onwards, the expressions from the file ftp://maia.usno.navy.mil/ser7/tai-utc.dat are used.
- The 5ms timestep at 1961 January 1 is taken from 2.58.1 (p87) of the 1992 Explanatory Supplement.
- 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
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:
- A valid non-null field preceded by the character 'B' (or 'b') is accepted.
- 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 decodeddreslt
- Previous Result- Returns:
- Result
-
Dc62s
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
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
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
Convert an interval in days into hours, minutes, seconds.- Parameters:
days
- interval in days- Returns:
- hours, minutes, seconds, fraction
-
Deuler
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 occurphi
- 1st rotation (radians)theta
- 2nd rotation ( " )psi
- 3rd rotation ( " )- Returns:
- rotation matrix
-
Dfltin
Convert free-format input into double precision floating point.- Notes:
-
- A tab character is interpreted as a space, and lower case d,e are interpreted as upper case.
- The basic format is #^.^@#^ where # means + or -, ^ means a decimal subfield and @ means D or E.
- 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.
- 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.
- Both signs are optional. The default is +.
- The mantissa defaults to 1.
- The exponent defaults to e0.
- The decimal subfields may be of any length.
- The decimal point is optional for whole numbers.
- 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.
- nstrt = 1 for the first character in the string.
- 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.
- 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.
- 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.
- 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.
- Decoding is left to right, one pass.
- 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.)
- See also Flotin and Intin.
- Parameters:
string
- String containing field to be decodeddreslt
- 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 matrixva
- vector- Returns:
- vector
-
Djcal
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
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
- Matrixy
- 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
- Matrixb
- 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
- Matrixva
- Vector- Returns:
- Result vector
-
Dr2af
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
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
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
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
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
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
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
Conversion of position & velocity in spherical coordinates to Cartesian coordinates.- Parameters:
s
- Spherical coordinates (longitude, latitude, radial)- Returns:
- Cartesian position & velocity vector
-
Ds2tp
Projection of spherical coordinates onto tangent plane ('gnomonic' projection - 'standard coordinates').- Parameters:
r
- spherical coordinates of point to be projectedrz
- 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:
-
- 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 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).
- The use of double-precision for both argument and result is purely for compatibility with other LIB time routines.
- The models used are based on a lunar tidal acceleration value of -26.00 arcsec per century.
- Depending on the epoch, one of three parabolic approximations
is used:
- 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
Convert hours, minutes, seconds to days.- Notes:
-
- The result is computed even if any of the range checks fail.
- The sign must be dealt with outside this routine.
- Parameters:
ihour
- Hoursimin
- Minutessec
- Seconds- Returns:
- Interval in days
- Throws:
palError
- Hour, Min or Sec out of range
-
Dtf2r
Convert hours, minutes, seconds to radians.- Notes:
-
- The result is computed even if any of the range checks fail.
- The sign must be dealt with outside this routine.
- Parameters:
ihour
- Hoursimin
- Minutessec
- Seconds- Returns:
- Angle in radians
- Throws:
palError
- Hour, Min or Sec out of range
-
Dtp2s
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:
-
- 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.
- Pre 1972 January 1 a fixed value of 10 + ET-TAI is returned.
- 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 vectorvb
- 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
- Vectoruv
- (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 vectorvb
- Second vector- Returns:
- Vector result
-
Ecleq
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:
-
- 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.
- 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'=Juliank
- 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
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
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:
-
- Smith, C.A. et al, 1989. "The transformation of astrometric catalog systems to the equinox J2000.0". Astron.J. 97, 265.
- 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 datedvb
- (Returned) barycentric velocitydpb
- (Returned) barycentric positiondvh
- (Returned) heliocentric velocitydph
- (Returned) heliocentric position
-
Fk425
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:
-
- The proper motions in RA are dRA/dt rather than cos(Dec)*dRA/dt, and are per year rather than per century.
- 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.
- 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:
-
- Smith, C.A. et al, 1989. "The transformation of astrometric catalog systems to the equinox J2000.0". Astron.J. 97, 265.
- 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.
- 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
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:
-
- The epoch BEPOCH is strictly speaking Besselian, but if a Julian epoch is supplied the result will be affected only to a negligible extent.
- 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.
- 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:
-
- Aoki,S., et al, 1983. Astron. Astrophys., 128, 263.
- Smith, C.A. et al, 1989. "The transformation of astrometric catalog systems to the equinox J2000.0". Astron.J. 97, 265.
- 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.
- 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
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:
-
- The proper motions in RA are dRA/dt rather than cos(Dec)*dRA/dt, and are per year rather than per century.
- 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.
- 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:
-
- Smith, C.A. et al, 1989. "The transformation of astrometric catalog systems to the equinox J2000.0". Astron.J. 97, 265.
- 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.
- 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
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:
-
- The proper motion in RA is dRA/dt rather than cos(Dec)*dRA/dt.
- 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.
- Unlike in the FK524 routine, the FK5 proper motions, the parallax and the radial velocity are presumed zero.
- 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.
- 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
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
Transformation from IAU 1958 Galactic coordinates to De Vaucouleurs supergalactic coordinates.- References:
-
- De Vaucouleurs, De Vaucouleurs, & Corwin, Second reference catalogue of bright galaxies, U. Texas, page 8.
- 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:
-
- Geocentric latitude can be obtained by evaluating atan2(z,r).
- 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
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 prefixIf jb=0, B is assumed for e < 1984.0, otherwise J.
-
Map
Transform star RA,Dec from mean place to geocentric apparent.The reference frames and timescales used are post IAU 1976.
- Notes:
-
- 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.
- The distinction between the required TDB and TT is always negligible. Moreover, for all but the most critical applications UTC is adequate.
- The proper motions in RA are dRA/dt rather than cos(Dec)*dRA/dt.
- 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).
- The accuracy is sub-milliarcsecond, limited by the precession-nutation model (IAU 1976 precession, Shirai & Fukushima 2001 forced nutation and precession corrections).
- 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
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:
-
- For date, the distinction between the required TDB and TT is always negligible. Moreover, for all but the most critical applications UTC is adequate.
- The vectors amprms(1-3) and amprms(4-6) are referred to the mean equinox and equator of epoch eq.
- The parameters AMPRMS produced by this routine are used by Ampqk, Mapqk and Mapqkz.
- The accuracy is sub-milliarcsecond, limited by the precession-nutation model (IAU 1976 precession, Shirai & Fukushima 2001 forced nutation and precession corrections).
- 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
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:
-
- The vectors amprms(1-3) and amprms(4-6) are referred to the mean equinox and equator of epoch eq.
- 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
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:
-
- The vectors amprms(1-3) and amprms(4-6) are referred to the mean equinox and equator of epoch eq.
- 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:
-
- The routine predicts forced nutation (but not free core nutation) plus corrections to the IAU 1976 precession model.
- 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.
- 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:
-
- Shirai, T. & Fukushima, T., Astron.J. 121, 3270-3283 (2001).
- Fukushima, T., 1991, Astron.Astrophys. 244, L11 (1991).
- 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
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
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
Apply corrections for proper motion to a star RA,Dec.- Notes:
-
- The proper motions in RA are dRA/dt rather than cos(Dec)*dRA/dt, and are in the same coordinate system as R0,D0.
- 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 epochpx
- 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 epochbep1
- 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:
-
- The epochs are TDB (loosely ET) Julian epochs.
- The matrix is in the sense v(ep1) = rmatp * v(ep0).
- 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:
-
- Lieske,J.H., 1979. Astron. Astrophys.,73,282. equations (6) & (7), p283.
- Kaplan,G.H., 1981. USNO circular no. 163, pa2.
- Parameters:
ep0
- Beginning epochep1
- Ending epoch- Returns:
- Precession matrix
-
Preces
Precession - either FK4 (Bessel-Newcomb, pre-IAU1976) or FK5 (Fricke, post-IAU1976) as required.- Notes:
-
- 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 )
- 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.
- 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 epochep1
- Ending epochd
- 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:
-
- The epochs are TDB (loosely ET) Julian epochs.
- The matrix is in the sense v(ep1) = rmatp * v(ep0).
- 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 epochep1
- 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:
-
- The epoch and date are TDB (loosely ET).
- The matrix is in the sense v(true) = rmatpn * v(mean).
- Parameters:
epoch
- Julian epoch for mean coordinatesdate
- 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:
-
- Typical values for the tlr and eps arguments might be 0.0065 and 1e-10 respectively.
- The radio refraction is chosen by specifying wl > 100 micrometres.
- 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:
-
- 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.
- 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.
- 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).
- 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.
- 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.
- 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.
- 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.
- 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.
- 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
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 includedeq
- Besselian epoch of mean equator and equinox- Returns:
- RA,Dec (radians) without e-terms
-
Supgal
Transformation from De Vaucouleurs supergalactic coordinates to IAU 1958 Galactic coordinates.- References:
-
- De Vaucouleurs, De Vaucouleurs, & Corwin, Second Reference Catalogue of Bright Galaxies, U. Texas, page 8.
- 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:
-
- The latitude must be geodetic. In critical applications, corrections for polar motion should be applied.
- 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.
- No range checking of arguments is done.
- 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 radiansdec
- Declination in radiansphi
- Observatory latitude in radians- Returns:
- Zenith distance in the range 0 to pi
-