Module barycentric_correction

Calculates heliocentric and barycentric velocity components of Earth.

This function takes into account the Earth-Moon motion, and is useful for
radial velocity work to an accuracy of  ~1 m/s.

dvel_hel, dvel_bary = baryvel(dje, deq)
 
Inputs:
@param dje: Julian ephemeris date.
@type dje: array
@param deq: epoch of mean equinox of dvelh and dvelb. If deq = 0 then deq is assumed to be equal to dje
@type deq: float
@return: heliocentric and barycentric velocity components in km/s (The 3-vectors dvelh and dvelb are given in a right-handed coordinate system with the +X axis toward the Vernal Equinox, and +Z axis toward the celestial pole.)
@rtype: array (2X3)          

Functions called:
  premat() -- computes precession matrix
  
NOTES:
  Algorithm taken from FORTRAN program of Stumpff (1980, A&A Suppl, 41,1)
  Stumpf claimed an accuracy of 42 cm/s for the velocity. A comparison with
  the JPL FORTRAN planetary ephemeris program PLEPH found agreement to within
  about 65 cm/s between 1986 and 1994
 
REVISION HISTORY:
  Jeff Valenti,  U.C. Berkeley    Translated BARVEL.FOR to IDL.
  W. Landsman, Cleaned up program sent by Chris McCarthy (SfSU) June 1994
  Converted to IDL V5.0   W. Landsman   September 1997
  Added /JPL keyword  W. Landsman   July 2001
  Documentation update W. Landsman Dec 2005
  Converted to Python S. Koposov 2009-2010
  Adapted for IvS library

Precess positions from J2000.0 (FK5) to B1950.0 (FK4)

Calculates the mean place of a star at B1950.0 on the FK4 system from the mean
place at J2000.0 on the FK5 system.

Input:
  @param ra0: J2000 right ascension in degrees
  @type ra0: N array
  @param dec0: J2000 declination in degrees
  @type dec0: N array
  @keyword mu_radec: array containing the proper motion in seconds of arc per tropical *century* in right ascension and declination.
  @type mu_radec: 2xN array
  @keyword parallax: array giving stellar parallax (seconds of arc)
  @type parallax: N array
  @keyword rad_vel: array giving radial velocity in km/s
  @type rad_vel: N array
  @keyword epoch: scalar giving epoch of original observations, default 2000. This keyword value is only used if mu_radec is not given
  @type epoch:
  @return: The corresponding B1950 right ascension and declination in *degrees*.
  @rtype: N array, N array
  
NOTES:
  The algorithm is taken from the Explanatory Supplement to the Astronomical Almanac 1992, page 186. Also see Aoki et al (1983), A&A, 128,263
  
  BPRECESS distinguishes between the following two cases:
  (1) The proper motion is known and non-zero
  (2) the proper motion is unknown or known to be exactly zero (i.e.
      extragalactic radio sources).   In this case, the reverse of the
      algorithm in Appendix 2 of Aoki et al. (1983) is used to ensure that the
      output proper motion is exactly zero. Better precision can be achieved
      in this case by inputting the epoch of the original observations.
      
REVISION HISTORY:
  Written (W. Landsman, october 1992)
  Vectorized (W. Landsman, february, 1994)
  Treat case where proper motion not known or exactly zero  (november 1994)
  Handling of arrays larger than 32767   (Lars L. Christensen, march 1995)
  Converted to IDL V5.0   (W. Landsman, september 1997)
  Fixed bug where A term not initialized for vector input (W. Landsman, february 2000)
  Converted to python (Sergey Koposov, july 2010)

Convolution of an image with a Point Spread Function (PSF)
The default is to compute the convolution using a product of Fourier transforms (for speed).

Inputs:
  @param image: 2-D array (matrix) to be convolved with psf
  @type image: array
  @param psf: the Point Spread Function, (size < or = to size of image).
  @type psf array
  @keyword ft_psf: passes out/in the Fourier transform of the PSF, (so that it can be re-used the next time function is called).
  @type ft_psc: bool
  @param ft_image: passes out/in the Fourier transform of image.
  @type ft_image: bool
  @param correlate: uses the conjugate of the Fourier transform of PSF to compute the cross-correlation of image and PSF (equivalent to IDL function convol() with NO rotation of PSF)
  @type correlate: bool
  @param auto_corr: computes the auto-correlation function of image using FFT.
  @type auto_corr: bool

METHOD:
  When using FFT, PSF is centered & expanded to size of image.

HISTORY:
  Written (Frank Varosi, 1992)
  Appropriate precision type for result depending on input image (Markus Hundertmark, february 2006)
  Fix the bug causing the recomputation of FFT(psf) and/or FFT(image) (Sergey Koposov, december 2006)

 Converts Julian dates to Gregorian calendar dates
 
 @param xjd: Julian date, positive double precision scalar or vector
 @type: array
 
OUTPUTS:
  YR: Year (Integer)
  MN: Month (Integer)
  DAY: Day (Integer)
  HR: Hours and fractional hours (Real). If XJD is a vector, then YR,MN,DAY and HR will be vectors of the same length.
  
REVISION HISTORY:
  Converted to IDL from Yeoman's Comet Ephemeris Generator, B. Pfarr, STX, 6/16/88
  Converted to IDL V5.0   W. Landsman   September 1997
 

Transform between Galactic, celestial, and ecliptic coordinates.

Input:
  @param ai: Input longitude in degrees
  @type ai: float
  @param bi: Input latitude in degrees
  @type bi: float
  @keyword select: Integer (1-6) specifying type of coordinate transformation.
  @type select: int 1-6
  @keyword fk4: set the fk4, otherwise use J2000
  @type fk4: bool
  
SELECT   From          To        |   SELECT      From            To
 1     RA-Dec (2000)  Galactic   |     4       Ecliptic      RA-Dec
 2     Galactic       RA-DEC     |     5       Ecliptic      Galactic
 3     RA-Dec         Ecliptic   |     6       Galactic      Ecliptic

Calculate the Galactic space velocity (U,V,W) of star

Calculates the Galactic space velocity U, V, W of star given its
(1) coordinates, (2) proper motion, (3) distance (or parallax), and (4) radial velocity.

OUTPUT PARAMETERS:
  U - Velocity (km/s) positive toward the Galactic *anti*center
  V - Velocity (km/s) positive in the direction of Galactic rotation
  W - Velocity (km/s) positive toward the North Galactic Pole
INPUT KEYWORDS:
  User must supply a position, proper motion,radial velocity and distance
  (or parallax).    Either scalars or vectors can be supplied.
 (1) Position:
  RA - Right Ascension in *Degrees*
  Dec - Declination in *Degrees*
 (2) Proper Motion
  PMRA = Proper motion in RA in arc units (typically milli-arcseconds/yr)
  PMDEC = Proper motion in Declination (typically mas/yr)
 (3) Radial Velocity
  VRAD = radial velocity in km/s
 (4) Distance or Parallax
  DISTANCE - distance in parsecs
             or
  PLX - parallax with same distance units as proper motion measurements
        typically milliarcseconds (mas)

REVISION HISTORY:
  Written, W. Landsman                       December   2000
  fix the bug occuring if the input arrays are longer than 32767
    and update the Sun velocity           Sergey Koposov June 2008
    vectorization of the loop -- performance on large arrays
    is now 10 times higher                Sergey Koposov December 2008

calculates heliocentric Julian date, baricentric and heliocentric radial
velocity corrections from:

INPUT:
@param ra2000: Right ascension of object for epoch 2000.0 (hours)
@type ra2000: float
@param dec2000: declination of object for epoch 2000.0 (degrees)
@type dec2000: float
@param jd: julian date for the middle of exposure
@type jd: float
@keyword obs_long: longitude of observatory (degrees, western direction is positive)
@type obs_long: float
@keyword obs_lat: latitude of observatory (degrees)
@type obs_lat: float
@keyword obs_alt: altitude of observatory (meters)
@type obs_alt: float

Algorithms used are taken from the IRAF task noao.astutils.rvcorrect
and some procedures of the IDL Astrolib are used as well.
Accuracy is about 0.5 seconds in time and about 1 m/s in velocity.

History:
written by Peter Mittermayer, Nov 8,2003
2005-January-13   Kudryavtsev   Made more accurate calculation of the sideral time.
                                Conformity with MIDAS compute/barycorr is checked.
2005-June-20      Kochukhov Included precession of RA2000 and DEC2000 to current epoch

Convert geocentric (reduced) Julian date to heliocentric Julian date

This procedure correct for the extra light travel time between the Earth and the Sun. An online calculator for this quantity is available at http://www.physics.sfasu.edu/astro/javascript/hjd.html

Input:

Parameters:

• date (float) - Reduced Julian date (= JD - 2400000)
• ra (float) - J2000 right ascension in DEGREES (unless B1950 keyword is set)
• dec (float) - J2000 declination in DEGREES (unless B1950 keyword is set)
• b1950 (bool) - if set, then input coordinates are assumed to be in equinox B1950 coordinates
• time_diff (bool) - if set, then HELIO_JD() returns the time difference (heliocentric JD - geocentric JD ) in seconds

Returns:

heliocentric reduced Julian date. If time_diff is True, then helio_jd() instead returns the time difference in seconds between the geocentric and heliocentric Julian date.

Precess coordinates from EQUINOX1 to EQUINOX2.

For interactive display, one can use the procedure ASTRO which calls PRECESS
or use the /PRINT keyword. The default (RA,DEC) system is FK5 based on epoch
J2000.0 but FK4 based on B1950.0 is available via the /FK4 keyword.

Input:
  @param ra0: Input right ascension in DEGREES
  @type ra0: float
  @param dec0: Input declination in degrees
  @type dec0: float
  @param equinox1: first equinox
  @type equinox1: float
  @param equinox2: second equinox
  @type equinox2: float
  @keyword fk4: If this keyword is set and non-zero, the FK4 (B1950.0) system will be used otherwise FK5 (J2000.0) will be used instead.
  @type fk4: float
  @keyword radian: If this keyword is set, input is in radian instead of degrees
  @type radian
  
Restrictions:
  Accuracy of precession decreases for declination values near 90 degrees.
  PRECESS should not be used more than 2.5 centuries from 2000 on the FK5 system (1950.0 on the FK4 system).

Procedure:
  Algorithm from Computational Spherical Astronomy by Taff (1983), p24. (FK4).
  FK5 constants from "Astronomical Almanac Explanatory Supplement 1992, page 104 Table 3.211.1.

Revision history
  Written, Wayne Landsman, STI Corporation  August 1986
  Correct negative output RA values   February 1989
  Provided FK5 (J2000.0)  I. Freedman   January 1994
  Precession Matrix computation now in PREMAT   W. Landsman June 1994
  Work for arrays, not just vectors (W. Landsman, september 2003)
  Convert to Python (S. Koposov, july 2010)
  Converted for use at IvS (K. Smolders)

Precess equatorial geocentric rectangular coordinates.

Input:
@param x: heliocentric rectangular x-coordinates
@type x: float
@param y: heliocentric rectangular y-coordinates
@type y: float
@param z: heliocentric rectangular z-coordinates
@type z: float
@param equinox1: first equinox
@type equinox1: float
@param equinox2: second equinox
@type equinox2: float

NOTES:
  The equatorial geocentric rectangular coords are converted to RA and Dec,
  precessed in the normal way, then changed back to x, y and z using unit
  vectors.

Return the precession matrix needed to go from EQUINOX1 to EQUINOX2.

This matrix is used by the procedures PRECESS and BARYVEL to precess astronomical coordinates
@param equinox1: Original equinox of coordinates, numeric scalar.
@type equinox1: float
@param equinox2: Equinox of precessed coordinates.
@type equinox2: float
@return: double precision 3 x 3 precession matrix, used to precess equatorial rectangular coordinates
@rtype: 3 by 3 array
@keyword fk4: If this keyword is set, the FK4 (B1950.0) system precession angles are used to compute the precession matrix. The default is to use FK5 (J2000.0) precession angles

Revision history:
  Written, Wayne Landsman, HSTX Corporation, June 1994
  Converted to IDL V5.0   W. Landsman   September 1997

Measures the spherical distance in degrees. The input has to be in degrees too

Calculate geocentric X,Y, and Z  and velocity coordinates of the Sun

Calculates geocentric X,Y, and Z vectors and velocity coordinates
(dx, dy and dz) of the Sun.   (The positive X axis is directed towards the
equinox, the y-axis, towards the point on the equator at right ascension 6h,
and the z axis toward the north pole of the equator). Typical position
accuracy is <1e-4 AU (15000 km).

Input:
@param date: reduced julian date (=JD - 2400000)
@type date: float
@keyword equinox: equinox of output. Default is 1950.
@type equinox: float

Output:
  x,y,z: scalars or vectors giving heliocentric rectangular coordinates
         (in A.U) for each date supplied. Note that sqrt(x^2 + y^2 + z^2)
         gives the Earth-Sun distance for the given date.
  xvel, yvel, zvel: velocity vectors corresponding to X, Y and Z.