Source code for pydl.goddard.astro

# Licensed under a 3-clause BSD style license - see LICENSE.rst
# -*- coding: utf-8 -*-
"""This module corresponds to the goddard/astro directory in idlutils.
"""
from time import time
import numpy as np
from astropy.units import Angstrom


def airtovac(air):
    """Convert air wavelengths to wavelengths in vacuum.

    Parameters
    ----------
    air : array-like
        Values of wavelength in air in Angstroms.
        :class:`~astropy.units.Quantity` objects with valid length
        dimensions will be internally converted to Angstrom.

    Returns
    -------
    array-like
        Values of wavelength in vacuum in Angstroms.  If a
        :class:`~astropy.units.Quantity` object was passed in, the output
        will be converted to the same units as the input.

    Notes
    -----
    * Formula from `P. E. Ciddor, Applied Optics, 35, 1566 (1996)
      <https://ui.adsabs.harvard.edu/abs/1996ApOpt..35.1566C/abstract>`_.
    * Values of wavelength below 2000 Å are not converted.
    """
    try:
        u = air.unit
    except AttributeError:
        u = None
    try:
        t = air.dtype
    except AttributeError:
        # Most likely, air is simply a float.
        t = None
    if t is None:
        if air < 2000.0:
            return air
        vacuum = air
        a = air
        g = None
    else:
        try:
            a = air.to(Angstrom).value
        except AttributeError:
            a = air
        g = a < 2000.0
        if g.all():
            return air
        vacuum = np.zeros(air.shape, dtype=t) + a
    for k in range(2):
        sigma2 = (1.0e4/vacuum)**2
        fact = (1.0 + 5.792105e-2/(238.0185 - sigma2) +
                1.67917e-3/(57.362 - sigma2))
        vacuum = a * fact
    if g is not None:
        vacuum[g] = a[g]
    if u is not None:
        vacuum = (vacuum * Angstrom).to(u)
    return vacuum


def gcirc(ra1, dec1, ra2, dec2, units=2):
    """Computes rigorous great circle arc distances.

    Parameters
    ----------
    ra1, dec1, ra2, dec2 : :class:`float` or array-like
        RA and Dec of two points.
    units : { 0, 1, 2 }, optional
        * units = 0: everything is already in radians
        * units = 1: RA in hours, dec in degrees, distance in arcsec.
        * units = 2: RA, dec in degrees, distance in arcsec (default)

    Returns
    -------
    :class:`float` or array-like
        The angular distance.  Units of the value returned depend on the
        input value of `units`.

    Notes
    -----
    The formula below is the one best suited to handling small angular
    separations.  See:
    https://en.wikipedia.org/wiki/Great-circle_distance
    """
    if units == 0:
        rarad1 = ra1
        dcrad1 = dec1
        rarad2 = ra2
        dcrad2 = dec2
    elif units == 1:
        rarad1 = np.deg2rad(15.0*ra1)
        dcrad1 = np.deg2rad(dec1)
        rarad2 = np.deg2rad(15.0*ra2)
        dcrad2 = np.deg2rad(dec2)
    elif units == 2:
        rarad1 = np.deg2rad(ra1)
        dcrad1 = np.deg2rad(dec1)
        rarad2 = np.deg2rad(ra2)
        dcrad2 = np.deg2rad(dec2)
    else:
        raise ValueError('units must be 0, 1 or 2!')
    deldec2 = (dcrad2-dcrad1)/2.0
    delra2 = (rarad2-rarad1)/2.0
    sindis = np.sqrt(np.sin(deldec2)*np.sin(deldec2) +
                  np.cos(dcrad1)*np.cos(dcrad2)*np.sin(delra2)*np.sin(delra2))
    dis = 2.0*np.arcsin(sindis)
    if units == 0:
        return dis
    else:
        return np.rad2deg(dis)*3600.0


def get_juldate(seconds=None):
    """Returns the current Julian date.

    Uses the MJD trick & adds the offset to get JD.

    Parameters
    ----------
    seconds : :class:`int` or :class:`float`, optional
        Time in seconds since the UNIX epoch.  This should only be used
        for testing.

    Returns
    -------
    :class:`float`
        The Julian Day number as a floating point number.

    Notes
    -----
    Do not use this function if high precision is required, or if you are
    concerned about the distinction between UTC & TAI.
    """
    if seconds is None:
        t = time()
    else:
        t = seconds
    mjd = t/86400.0 + 40587.0
    return mjd + 2400000.5


def get_juldate_main():  # pragma: no cover
    """Entry point for the get_juldate command-line script.
    """
    print(get_juldate())
    return 0


def vactoair(vacuum):
    """Convert vacuum wavelengths to wavelengths in air.

    Parameters
    ----------
    vacuum : array-like
        Values of wavelength in vacuum in Angstroms.
        :class:`~astropy.units.Quantity` objects with valid length
        dimensions will be internally converted to Angstrom.

    Returns
    -------
    array-like
        Values of wavelength in air in Angstroms.
        :class:`~astropy.units.Quantity` object was passed in, the output
        will be converted to the same units as the input.

    Notes
    -----
    * Formula from `P. E. Ciddor, Applied Optics, 35, 1566 (1996)
      <https://ui.adsabs.harvard.edu/abs/1996ApOpt..35.1566C/abstract>`_.
    * Values of wavelength below 2000 Å are not converted.
    """
    try:
        u = vacuum.unit
    except AttributeError:
        u = None
    try:
        t = vacuum.dtype
    except AttributeError:
        # Most likely, vacuum is simply a float.
        t = None
    if t is None:
        if vacuum < 2000.0:
            return vacuum
        air = vacuum
        v = vacuum
        g = None
    else:
        try:
            v = vacuum.to(Angstrom).value
        except AttributeError:
            v = vacuum
        g = v < 2000.0
        if g.all():
            return vacuum
        air = np.zeros(vacuum.shape, dtype=t) + v
    sigma2 = (1.0e4/v)**2
    fact = (1.0 + 5.792105e-2/(238.0185 - sigma2) +
            1.67917e-3/(57.362 - sigma2))
    air = v / fact
    if g is not None:
        air[g] = v[g]
    if u is not None:
        air = (air * Angstrom).to(u)
    return air