pal.pyx.in

# Copyright (C) 2012 Tim Jenness and Science and Technology
# Facilities Council.
# Copyright (C) 2014 Tim Jenness

# This program is free software; you can redistribute it and/or
# modify it under the terms of the GNU General Public License as
# published by the Free Software Foundation; either version 3 of
# the License, or (at your option) any later version.
#
# This program is distributed in the hope that it will be
# useful, but WITHOUT ANY WARRANTY; without even the implied
# warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
# PURPOSE. See the GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program.  If not, see <http://www.gnu.org/licenses/>.

#cython: language_level=3
#distutils: language=c

cimport cpal

from libc.stdlib cimport malloc, free

cimport numpy as np
import numpy as np

D2PI = cpal.PAL__D2PI
DAS2R = cpal.PAL__DAS2R
DD2R = cpal.PAL__DD2R
DH2R = cpal.PAL__DH2R
DPI = cpal.PAL__DPI
DPIBY2 = cpal.PAL__DPIBY2
DR2AS = cpal.PAL__DR2AS
DR2D = cpal.PAL__DR2D
DR2H = cpal.PAL__DR2H
DR2S = cpal.PAL__DR2S
DS2R = cpal.PAL__DS2R

__version__ = @PALPY_VERSION@

def palvers():
    """
    versiondict = palvers()

    Returns the version information from the underlying PAL C
    library used to build palpy. Version information is returned in
    a dict with keys:

    verstring : Version in string form, e.g "1.6.2"
    major : Major version number as integer, e.g. 1
    minor : Minor version number as integer, e.g. 6
    patchlevel : Patch level as integer, e.g. 2
    """
    return dict( { "verstring": @VERSTRING@,
                    "major": @MAJOR_VERSION@,
                    "minor": @MINOR_VERSION@,
                    "patchlevel": @PATCH_VERSION@ } )

def addet( double rm, double dm, double eq ):
     """
     (rc, dc) = addet( rm, dm, eq)

     @palAddet@
     """
     cdef double rc
     cdef double dc
     cpal.palAddet( rm, dm, eq, &rc, &dc )
     return ( rc, dc )

def airmas( double zd ):
     """
     airmass = airmas( zd )

     @palAirmas@
     """
     return cpal.palAirmas( zd )

def airmasVector(np.ndarray[double, ndim=1] zd not None):
     """
     airmass = airmasVector( zd )

     where zd and airmass are numpy arrays

     @palAirmas@
     """
     cdef int length = len(zd)
     cdef np.ndarray airmass_out = np.empty(length, dtype=np.float64)
     for ii in range(length):
         airmass_out[ii] = cpal.palAirmas(zd[ii])
     return airmass_out


def altaz(double ha, double dec, double phi):
     """
     (az, azd, azdd, el, eld, eldd, pa, pad, padd) = altaz( ha, dec, phi )

     @palAltaz@
     """
     cdef double az
     cdef double azd
     cdef double azdd
     cdef double el
     cdef double eld
     cdef double eldd
     cdef double pa
     cdef double pad
     cdef double padd
     cpal.palAltaz (ha, dec, phi, &az, &azd, &azdd, &el, &eld, &eldd, &pa, &pad, &padd )
     return (az, azd, azdd, el, eld, eldd, pa, pad, padd )

def altazVector(np.ndarray[double, ndim=1] ha not None,
                np.ndarray[double, ndim=1] dec not None,
                double phi):
     """
     (az, azd, azdd, el, eld, eldd, pa, pad, padd) = altazVector( ha, dec, phi )

     where ha, dec and all outputs are numpy arrays

     @palAltaz@

     Raises
     ------
     ValueError if input numpy arrays are not of the same length
     """
     cdef int length = len(ha)

     if len(dec)!=length:
          raise ValueError("You did not pass as many Decs as " \
                           + "Hour Angles to altazVector")

     cdef np.ndarray azout = np.zeros(length, dtype=np.float64)
     cdef double az
     cdef np.ndarray azdout = np.zeros(length, dtype=np.float64)
     cdef double azd
     cdef np.ndarray azddout = np.zeros(length, dtype=np.float64)
     cdef double azdd
     cdef np.ndarray elout = np.zeros(length, dtype=np.float64)
     cdef double el
     cdef np.ndarray eldout = np.zeros(length, dtype=np.float64)
     cdef double eld
     cdef np.ndarray elddout = np.zeros(length, dtype=np.float64)
     cdef double eldd
     cdef np.ndarray paout = np.zeros(length, dtype=np.float64)
     cdef double pa
     cdef np.ndarray padout = np.zeros(length, dtype=np.float64)
     cdef double pad
     cdef np.ndarray paddout = np.zeros(length, dtype=np.float64)
     cdef double padd
     cdef int i
     for i in range(length):
          cpal.palAltaz(ha[i], dec[i], phi, &az, &azd, &azdd, &el, &eld, &eldd, &pa, &pad, &padd)
          azout[i] = az
          azdout[i] = azd
          azddout[i] = azdd
          elout[i] = el
          eldout[i] = eld
          elddout[i] = eldd
          paout[i] = pa
          padout[i] = pad
          paddout[i] = padd
     return (azout, azdout, azddout, elout, eldout, elddout, paout, padout, paddout)

def amp( double ra, double da, double date, double eq):
     """
     (rm, dm) = amp( ra, da, date, eq )

     @palAmp@
     """
     cdef double rm
     cdef double dm
     cpal.palAmp( ra, da, date, eq, &rm, &dm )
     return ( rm, dm )

def ampqk( double ra, double da, np.ndarray[double, ndim=1] amprms not None ):
     """
     (rm, dm) = ampqk( ra, da, amprms )

     @palAmpqk@
     """
     cdef double rm
     cdef double dm
     cdef double camprms[21]
     for i in range(21):
          camprms[i] = amprms[i]
     cpal.palAmpqk( ra, da, camprms, &rm, &dm )
     return (rm, dm)


def ampqkVector(np.ndarray[double, ndim=1] ra not None,
                np.ndarray[double, ndim=1] dec not None,
                np.ndarray[double, ndim=1] amprms not None):
     """
     (rm, dm) = ampqkVector( ra, dec, amprms )

     where ra, dec, rm, and dm are all numpy arrays

     @palAmpqk@

     Raises
     ------
     ValueError if input arrays of RA and Dec are not of the same
     length
     """

     cdef int length = len(ra)

     if len(dec)!=length:
         raise ValueError("You did not pass as many Decs as RAs to " \
                          + "ampqkVector")

     cdef np.ndarray rm_out = np.empty(length, dtype=np.float64)
     cdef np.ndarray dm_out = np.empty(length, dtype=np.float64)

     cdef double rm
     cdef double dm
     cdef double camprms[21]

     for ii in range(21):
          camprms[ii] = amprms[ii]

     for ii in range(length):
          cpal.palAmpqk(ra[ii], dec[ii], camprms, &rm, &dm)
          rm_out[ii] = rm
          dm_out[ii] = dm

     return (rm_out, dm_out)

def aop( double rap, double dap, double date, double dut,
         double elongm, double phim, double hm, double xp,
         double yp, double tdk, double pmb, double rh,
         double wl, double tlr ):
     """
     (aob, zob, hob, dob, rob) = aop( rap, dap, date, dut,
                                      elongm, phim, hm, xp,
                                      yp, tdk, pmb, rh,
                                      wl, tlr )

     @palAop@
     """
     cdef double aob
     cdef double zob
     cdef double hob
     cdef double dob
     cdef double rob
     cpal.palAop( rap, dap, date, dut, elongm, phim, hm,
                  xp, yp, tdk, pmb, rh, wl, tlr,
                  &aob, &zob, &hob, &dob, &rob)
     return (aob, zob, hob, dob, rob )

def aoppa(double date, double dut, double elongm, double phim,
          double hm, double xp, double yp, double tdk, double pmb,
          double rh, double wl, double tlr ):
     """
     aoprms = aoppa( date, dut, elongm, phim, hm, xp, yp,
                     tdk, pmb, rh, wl, tlr )

     @palAoppa@
     """
     cdef double caoprms[14]
     cdef np.ndarray aoprms = np.zeros( [14], dtype=np.float64 )
     cpal.palAoppa( date, dut, elongm, phim, hm, xp, yp, tdk, pmb,
                    rh, wl, tlr, caoprms )
     for i in range(14):
          aoprms[i] = caoprms[i]
     return aoprms

def aoppat( double date, np.ndarray[double, ndim=1] aoprms not None ):
     """
     aoprms_updated = aoppat( date, aoprms )

     @palAoppat@
     """
     # We can either copy the array or modify in place.
     # For now we return a new copy
     cdef double caoprms[14]
     for i in range(14):
          caoprms[i] = aoprms[i]
     cpal.palAoppat( date, caoprms )
     cdef np.ndarray result = np.zeros( [14], dtype=np.float64 )
     for i in range(14):
          result[i] = caoprms[i]
     return result

def aopqk(double rap, double dap, np.ndarray[double, ndim=1] aoprms not None):
    """
    (aob, zob, hob, dob, rob) = aopqk( rap, dap, aoprms )

    @palAopqk@
    """
    cdef double aob
    cdef double zob
    cdef double hob
    cdef double dob
    cdef double rob

    cdef double caoprms[14]
    for i in range(14):
        caoprms[i]=aoprms[i]

    cpal.palAopqk(rap,dap,caoprms,&aob,&zob,&hob,&dob,&rob)

    return (aob,zob,hob,dob,rob)

def aopqkVector(np.ndarray[double, ndim=1] rap not None,
                np.ndarray[double, ndim=1] dap not None,
                np.ndarray[double, ndim=1] aoprms not None):
    """
    (aob, zob, hob, dob, rob) = aopqkVector( rap, dap, aoprms )

    where rap, dap and all outputs are numpy arrays

    @palAopqk@

    Raises
    ------
    ValueError if input arrays are not of the same length
    """
    cdef int i
    cdef int length=len(rap)

    if len(dap)!=length:
         raise ValueError("You did not pass the same number of Decs " \
                          + "as RAs to aopqkVector")

    cdef np.ndarray aobout = np.zeros(length,dtype=np.float64)
    cdef np.ndarray zobout = np.zeros(length,dtype=np.float64)
    cdef np.ndarray hobout = np.zeros(length,dtype=np.float64)
    cdef np.ndarray dobout = np.zeros(length,dtype=np.float64)
    cdef np.ndarray robout = np.zeros(length,dtype=np.float64)
    cdef double aob
    cdef double zob
    cdef double hob
    cdef double dob
    cdef double rob

    cdef double caoprms[14]
    for i in range(14):
         caoprms[i]=aoprms[i]
    for i in range(length):
         cpal.palAopqk(rap[i],dap[i],caoprms,&aob,&zob,&hob,&dob,&rob)
         aobout[i]=aob
         zobout[i]=zob
         hobout[i]=hob
         dobout[i]=dob
         robout[i]=rob
    return (aobout,zobout,hobout,dobout,robout)

def atmdsp(double tdk, double pmb, double rh, double wl1,
                 double a1, double b1, double wl2):
    """
    (a2, b2) = atmdsp( tdk, pmb, rh, wl1, a1, b1, wl2 )

    @palAtmdsp@
    """
    cdef double a2
    cdef double b2
    cpal.palAtmdsp(tdk, pmb, rh, wl1, a1, b1, wl2, &a2, &b2)
    return (a2, b2)

def caldj( int iy, int im, int id ):
     """
     djm = caldj( iy, im, id )

     @palCaldj@

     Raises
     ------
     ValueError if arguments are out of range.
     """
     cdef double djm
     cdef int j
     cpal.palCaldj( iy, im, id, &djm, &j )
     if j==0:
          return djm
     else:
          bad = ""
          if j==-1:
               bad = "year"
          elif j==-2:
               bad = "month"
          else:
               bad = "day"
          raise ValueError( "Julian date not computed. Bad {}".format(bad) )


def caldjVector(np.ndarray[long, ndim=1] iy not None,
                np.ndarray[long, ndim=1] im not None,
                np.ndarray[long, ndim=1] id not None):
     """
     djm = caldjVector( iy, im, id )

     where iy, im, id and djm are vectors

     @palCaldj@

     Raises
     ------
     ValueError if input arrays have different lengths

     Bad input dates result in numpy.NaN for djm.

     Note: This algorithm does not work for dates before
     -4799 January 1

     Note: If you pass in a two-digit year, this algorithm will
     interpret it as being between 1949 and 2049
     """

     cdef int length = len(iy)

     if len(im)!=length:
         raise ValueError("You did not pass as many months as years " \
                           + "to caldjVector")

     if len(id)!=length:
         raise ValueError("You did not pass as many days as years " \
                           + "to caldjVector")

     cdef np.ndarray mjd_out = np.empty(length, dtype=np.float64)
     cdef double mjd
     cdef int flag

     for ii in range(length):
         cpal.palCaldj(iy[ii], im[ii], id[ii], &mjd, &flag)
         if flag == 0:
             mjd_out[ii] = mjd
         else:
             mjd_out[ii] = np.nan

     return mjd_out

def cldj( int iy, int im, int id ):
     """
     djm = cldj( iy, im, id )

     @palCldj@

     Raises
     ------
     ValueError if arguments are out of range.
     """
     cdef int j
     cdef double djm
     cpal.palCldj( iy, im, id, &djm, &j )
     if j==0:
          return djm
     else:
          bad = ""
          if j==-1:
               bad = "year"
          elif j==-2:
               bad = "month"
          else:
               bad = "day"
          raise ValueError( "Bad {} argument".format(bad) )

def dafin( string, int ipos ):
     """
     Decode a sexagesimal string into an angle.

       (angle, endpos) = pal.dafin( string, startpos )

     where string is the string to be analyzed, startpos
     is the position in the string to start looking
     (starts at 1), endpos is the position in the string
     after the search and angle is the decoded angle in
     radians.

     A ValueError exception is thrown if no angle can
     be found.

     A ValueError exception is thrown if an angle can
     be located but is numerically out of range.

     The interface for this routine is experimental. In
     particular the startpos and endpos variables need
     some thought. startpos could simply be ignored
     as a python slice would work just as well. endpos
     is useful in that you know where to slice the string
     to continue searching for angles. Zero- versus one-
     based indexing is an issue.
     """
     byte_string = string.encode('ascii')
     cdef char * cstring = byte_string
     cdef int cipos = ipos
     cdef double a
     cdef int j
     cpal.palDafin( cstring, &cipos, &a, &j )
     if j==0:
          return ( a, cipos )
     elif j==1:
          raise ValueError( "No angle could be located in string '{}'".format(string) )
     else:
          bad = "unknown"
          if j==-1:
               bad = "degrees"
          elif j==-2:
               bad = "arcminutes"
          elif j==-3:
               bad = "arcseconds"
          raise ValueError( "Bad {} in input string".format(bad) )

def dat( double dju ):
     """
     deltat = dat( dju )

     @palDat@
     """
     return cpal.palDat( dju )

def dav2m( np.ndarray[double, ndim=1] axvec not None ):
     """
     rmat = dav2m( axvec )

     @palDav2m@
     """
     cdef double c_rmat[3][3]
     cdef double c_axvec[3]
     for i in range(3):
         c_axvec[i] = axvec[i]
     cpal.palDav2m( c_axvec, c_rmat )

     cdef np.ndarray rmat = np.zeros([3,3], dtype=np.float64)
     for i in range(3):
         for j in range(3):
             rmat[i][j] = c_rmat[i][j]
     return rmat

def dbear( double a1, double b1, double a2, double b2 ):
     """
     b = dbear( a1, b1, a2, b2 )

     @palDbear@
     """
     return cpal.palDbear( a1, b1, a2, b2 )

def dbearVector( np.ndarray[ double, ndim=1] a1 not None,
                 np.ndarray[ double, ndim=1] b1 not None,
                 np.ndarray[ double, ndim=1] a2 not None,
                 np.ndarray[ double, ndim=1] b2 not None):
    """
    b_out = dbearVector( a1, b1, a2, b2)

    where a1 and b1 are numpy arrays of longitude.  a2 and b2 are
    also numpy arrays.  a2 and b2 can have either one element, or the
    same number of elements as a1, b1.  If a2 and b2 have only one element,
    b_out will  be a numpy array of the bearings from the points in a1, b1
    to the single point in a2, b2.  If a2 and b2 have the same number
    of elements as a1 and b1, b_out will be a numpy array of the bearings
    between each corresponding pair of points, i.e.

    b_out[0] = bearing(a1[0], b1[0], a2[0], b2[0])
    b_out[1] = bearing(a1[1], b1[1], a2[1], b2[1])
    ...
    b_out[n] = bearing(a1[n], b1[n], a2[n], b2[n])

    @palDbear@

    Raises
    ------
    Value error if input arrays have inconsistent lengths
    """
    cdef int length1 = len(a1)
    cdef int length2 = len(a2)
    if len(b1) != length1:
        raise ValueError("The first two arguments of dbearVector must " \
                         + "be numpy arrays with the same number of elements")

    if len(b2) != length2:
        raise ValueError("The second two arguments of dbearVector must " \
                         + "be numpy arrays with the same number of elements")

    if length2!=1 and length2!=length1:
        raise ValueError("The second two arguments of dbearVector must " \
                         + "either have the same number of elements " \
                         + "as the first two arguments, or they must have " \
                         + "only one element")

    cdef np.ndarray b_out = np.zeros(length1, dtype=np.float64)

    if length1==length2:
        for i in range(length1):
            b_out[i] = cpal.palDbear(a1[i], b1[i], a2[i], b2[i])
    else:
        for i in range(length1):
            b_out[i] = cpal.palDbear(a1[i], b1[i], a2[0], b2[0])

    return b_out

def daf2r( int ideg, int iamin, double asec ):
     """
     ang = daf2r( ideg, iamin, asec )

     @palDaf2r@

     Raises
     ------
     ValueError if the input arguments are out of range.
     """
     cdef double rad
     cdef int j
     cpal.palDaf2r( ideg, iamin, asec, &rad, &j )
     if j==0:
          return rad
     else:
          bad = ""
          if j==1:
               bad = "Degree argument outside range 0-359"
          elif j==2:
               bad = "Minute argument outside range 0-59"
          else:
               bad = "Arcsec argument outside range 0-59.9999..."
          raise ValueError( bad )

def daf2rVector(np.ndarray[long, ndim=1] ideg_in not None,
                np.ndarray[long, ndim=1] iamin_in not None,
                np.ndarray[double, ndim=1] asec_in not None):
    """
    ang = daf2rVector( ideg, iamin, asec)

    where ideg, iamin, and asec are numpy arrays of degrees,
    minutes, and seconds.  ideg and iamin are arrays of ints.  asec
    is an array of doubles.  ang is a numpy array of the same angles
    converted to radians.

    Bad values are mapped to numpy.NaNs in ang.

    Bad values are defined as:

    ideg > 359
    iamin > 59
    asec > 59.99999

    @palDaf2r@

    Raises
    ------
    ValueError if input arrays have inconsistent lengths
    """
    cdef int length = len(ideg_in)

    if len(iamin_in)!=length or len(asec_in)!=length:
        raise ValueError("You need to pass the same number of degrees as " \
                         + "minutes and seconds to daf2rVector")

    cdef np.ndarray ang_out = np.zeros(length, dtype=np.float64)
    cdef double rad
    cdef int flag

    for i in range(length):
       cpal.palDaf2r(ideg_in[i], iamin_in[i], asec_in[i], &rad, &flag)
       if flag==0:
           ang_out[i] = rad
       else:
           ang_out[i] = np.nan

    return ang_out

def dcc2s( np.ndarray[double, ndim=1] v not None ):
     """
     (a, b) = dcc2s( v )

     @palDcc2s@
     """
     cdef double cv[3]
     for i in range(3):
          cv[i] = v[i]
     cdef double a
     cdef double b
     cpal.palDcc2s( cv, &a, &b )
     return (a, b)

def dcc2sVector( np.ndarray[double, ndim=2] v_in not None):
    """
    (a, b) = dcc2sVector( v )

    where v is a 2-D numpy array in which each row is a different
    point in Cartesian coordinates.  a and b will be numpy arrays
    containing the corresponding spherical coordinates

    @palDcc2s@

    Raises
    ------
    ValueError if you do not pass in a 2-D numpy array with
    three columns
    """

    if v_in.shape[1] != 3:
        raise ValueError("The input to dcc2sVector should be " \
                         + "a 2-D numpy array in which each row " \
                         + "is a point in 3-D Cartesian coordinates." \
                         + " The rows of your input" \
                         + " have %d dimensions" % v_in.shape[1])

    cdef int length = v_in.shape[0]
    cdef double cv[3]
    cdef double a
    cdef double b
    cdef np.ndarray a_out = np.zeros(length, dtype=np.float64)
    cdef np.ndarray b_out = np.zeros(length, dtype=np.float64)

    for row in range(length):
        for i in range(3):
            cv[i] = v_in[row][i]
        cpal.palDcc2s( cv, &a, &b)
        a_out[row] = a
        b_out[row]= b

    return (a_out, b_out)

def dcs2c( double a, double b ):
     """
     v = dcs2c( a, b )

     @palDcs2c@
     """
     cdef double cv[3]
     cpal.palDcs2c( a, b, cv )
     cdef np.ndarray v = np.zeros( [3], dtype=np.float64 )
     for i in range(3):
          v[i] = cv[i]
     return v

def dcs2cVector( np.ndarray[double, ndim=1] a_in not None,
                 np.ndarray[double, ndim=1] b_in not None):
    """
    v = dcs2cVector(a, b)

    where a and b are numpy arrays of spherical coordinates.
    v is a 2-D numpy array in which each row contains the
    directional cosines of the corresponding Cartesian coordinates

    @palDcs2c@

    Raises
    ------
    ValueError if input arrays have inconsisten lengths
    """

    cdef int length = len(a_in)
    if len(b_in) != length:
        raise ValueError("You did not pass as many latitudes as longitudes " \
                         + "into dcs2sVector")

    cdef double cv[3]
    cdef np.ndarray v_out = np.zeros((length, 3), dtype=np.float64)
    for row in range(length):
        cpal.palDcs2c(a_in[row], b_in[row], cv)
        for i in range(3):
            v_out[row][i] = cv[i]

    return v_out


def dd2tf( int ndp, double days ):
     """
     (sign, ih, im, is, frac) = dd2tf( ndp, days )

     @palDd2tf@
     """
     cdef char csign = b' '
     cdef int ihmsf[4]
     cpal.palDd2tf( ndp, days, &csign, ihmsf )
     sign = chr(csign)
     return ( sign, ihmsf[0], ihmsf[1], ihmsf[2], ihmsf[3] )

def dd2tfVector( int ndp, np.ndarray[double, ndim=1] days not None):
    """
    (sign, ih, im, is ,frac) = ddtfVector( ndp, dayList)

    where ndp is an integer denoting the resolution to which the
    fraction of seconds should be taken and dayList is a numpy
    array of values to be converted. sign, ih, im, is, and frac
    will also be numpy arrays corresponding to the elements in
    dayList.

    @palDd2tf@
    """

    cdef int length = len(days)

    cdef np.ndarray sign_out = np.empty(length, dtype = (unicode, 1))
    cdef np.ndarray ih_out = np.empty(length, dtype = np.int64)
    cdef np.ndarray im_out = np.empty(length, dtype = np.int64)
    cdef np.ndarray is_out = np.empty(length, dtype = np.int64)
    cdef np.ndarray frac_out = np.empty(length, dtype = np.int64)

    cdef char csign = b' '
    cdef int ihmsf[4]
    for row in range(length):
        cpal.palDd2tf( ndp, days[row], &csign, ihmsf)
        sign_out[row] = chr(csign)
        ih_out[row] = ihmsf[0]
        im_out[row] = ihmsf[1]
        is_out[row] = ihmsf[2]
        frac_out[row] = ihmsf[3]

    return (sign_out, ih_out, im_out, is_out, frac_out)

def de2h( double ha, double dec, double phi ):
     """
     (az, el) = de2h( ha, dec, phi )

     @palDe2h@
     """
     cdef double az
     cdef double el
     cpal.palDe2h( ha, dec, phi, &az, &el )
     return (az, el)

def de2hVector( np.ndarray[double, ndim=1] ha not None,
                np.ndarray[double, ndim=1] dec not None,
                double phi ):
     """
     (az, el) = de2hVector( ha, dec, phi )

     where ha, dec, az, and el are numpy arrays

     @palDe2h@

     Raises
     ------
     ValueError if input arrays are not of the same length
     """
     cdef int i
     cdef int length=len(ha)

     if len(dec)!=length:
         raise ValueError("You did not pass the same number of " \
                          + "Decs as Hour Angles to de2hVector")

     cdef np.ndarray azout = np.zeros(length,dtype=np.float64)
     cdef np.ndarray elout = np.zeros(length,dtype=np.float64)
     cdef double az
     cdef double el
     for i in range(length):
          cpal.palDe2h( ha[i], dec[i], phi, &az, &el )
          azout[i]=az
          elout[i]=el
     return (azout, elout)



def deuler( order, double phi, double theta, double psi ):
    """
    rmat = deuler( order, phi, theta, psi )

    @palDeuler@
    """
    cdef double c_rmat[3][3]
    byte_order = order.encode('ascii')
    cdef char * c_order = byte_order
    cpal.palDeuler( c_order, phi, theta, psi, c_rmat )
    cdef np.ndarray rmat = np.zeros([3,3], dtype=np.float64)
    for i in range(3):
        for j in range(3):
            rmat[i][j] = c_rmat[i][j]
    return rmat

# dfltin() not implemented -- not necessary for python

def dh2e( double az, double el, double phi ):
     """
     (ha, dec) = dh2e( az, el, phi )

     @palDh2e@
     """
     cdef double ha
     cdef double dec
     cpal.palDh2e( az, el, phi, &ha, &dec )
     return (ha, dec)

def dh2eVector(np.ndarray[double, ndim=1] az not None,
               np.ndarray[double, ndim=1] el not None,
               double phi):
     """
     (ha, dec) = dh2eVector( az, el, phi)

     where az, el, ha, and dec are all numpy arrays

     @palDh2e@

     Raises
     ------
     ValueError if the input arrays are of different lengths
     """

     cdef int length = len(az)

     if len(el)!=length:
         raise ValueError("You did not pass as many elevations as " \
                          + "azimuths to dh2eVector")

     cdef np.ndarray ha_out = np.empty(length, dtype=np.float64)
     cdef np.ndarray dec_out = np.empty(length, dtype=np.float64)

     cdef double ha
     cdef double dec

     for ii in range(length):
         cpal.palDh2e(az[ii], el[ii], phi, &ha, &dec)
         ha_out[ii] = ha
         dec_out[ii] = dec

     return (ha_out, dec_out)

def dimxv( np.ndarray[double, ndim=2] dm not None, np.ndarray[double, ndim=1] va not None):
     """
     vb = dimxv( dm, va )

     @palDimxv@
     """
     cdef double c_dm[3][3]
     cdef double c_va[3]
     cdef double c_vb[3]
     for i in range(3):
          for j in range(3):
                c_dm[i][j] = dm[i][j]
          c_va[i] = va[i]
     cpal.palDimxv(c_dm, c_va, c_vb)

     cdef np.ndarray vb = np.zeros([3], dtype=np.float64)
     for i in range(3):
          vb[i] = c_vb[i]
     return vb

def djcal( int ndp, double djm ):
     """
     (iy, im, id, frac) = djcal( ndp, djm )

     @palDjcal@

     Raises
     ------
     ValueError for unacceptable dates.
     """
     cdef int iymdf[4]
     cdef int j
     cpal.palDjcal( ndp, djm, iymdf, &j )
     if j==0:
          return ( iymdf[0], iymdf[1], iymdf[2], iymdf[3] )
     else:
          raise ValueError( "Unacceptable date" )


def djcalVector(int ndp, np.ndarray[double, ndim=1] djm not None):
     """
     (iy, im, id, frac) = djcalVector( ndp, djm)

     where djm and outputs are numpy arrays

     @palDjcal@

     Raises
     ------
     Does not raise any exceptions; unacceptable dates (mjd<-2468570
     or mjd>1e9) result in -1 in iy, im, id and frac.

     Note: -1 is a valid result of iy (not for im, id, or frac), so
     any filters on unacceptable dates should examine the other
     outputs.
     """

     cdef int length = len(djm)
     cdef np.ndarray iy_out = np.empty(length, dtype=np.int64)
     cdef np.ndarray im_out = np.empty(length, dtype=np.int64)
     cdef np.ndarray id_out = np.empty(length, dtype=np.int64)
     cdef np.ndarray frac_out = np.empty(length, dtype=np.int64)
     cdef int output[4]
     cdef int flag

     for ii in range(length):
          cpal.palDjcal(ndp, djm[ii], output, &flag)
          if flag==0:
              iy_out[ii] = output[0]
              im_out[ii] = output[1]
              id_out[ii] = output[2]
              frac_out[ii] = output[3]
          else:
              iy_out[ii] = -1
              im_out[ii] = -1
              id_out[ii] = -1
              frac_out[ii] = -1

     return (iy_out, im_out, id_out, frac_out)

def djcl( double djm ):
     """
     (iy, im, id, frac) = djcl( djm )

     @palDjcl@

     Raises
     ------
     ValueError for unacceptable date.
     """
     cdef int iy
     cdef int im
     cdef int id
     cdef double fd
     cdef int j
     cpal.palDjcl( djm, &iy, &im, &id, &fd, &j )
     if j==0:
          return ( iy, im, id, fd )
     else:
          raise ValueError( "Unacceptable date" )

def dmat( np.ndarray[double, ndim=2] a not None,
          np.ndarray[double, ndim=1] y not None ):
     """
     (na, ny, d) = dmat( a, y )

     @palDmat@

     Raises
     ------
     MemoryError: Unable to get workspace

     ArithmeticError: Singular matrix

     ValueError: Shapes of input matrices incompatible
     """
     # validate the arguments and get the dimension
     ashape = a.shape
     yshape = y.shape
     if ashape[0] != ashape[1]:
          raise ValueError( "Matrix must be square" )
     if yshape[0] != ashape[0]:
          raise ValueError( "Matrix must match number of elements in supplied vector" )

     cdef int n = y.size
     cdef int j
     cdef double d
     cdef int *iw = <int *>malloc( n * sizeof(int) )
     cdef double *ca = <double *>malloc( n * n * sizeof(double))
     cdef double *cy = <double *>malloc( n * sizeof(double))

     if not ca or not iw or not cy:
          if ca:
               free(ca)
          if iw:
               free(iw)
          if cy:
               free(cy)
          raise MemoryError( "Could not get dynamic memory for matrix" )

     # Need to flatten the input 2d matrix
     k = 0;
     for i in range(n):
          cy[i] = y[i]
          for j in range(n):
               ca[k] = a[ i, j ]
               k = k + 1

     cpal.palDmat( n, ca, cy, &d, &j, iw )
     free(iw)

     cdef np.ndarray na = np.zeros( [n,n], dtype=np.float64 )
     cdef np.ndarray ny = np.zeros( [n], dtype=np.float64 )

     if j==0:
          k = 0
          for i in range(n):
               ny[i] = cy[i]
               for j in range(n):
                    na[i,j] = ca[k]
                    k = k + 1
          free(ca)
          free(cy)
          return ( na, ny, d )
     else:
          free(ca)
          free(cy)
          raise ArithmeticError( "Matrix is singular" )

def dmoon( double date ):
     """
     pv = dmoon( date )

     @palDmoon@
     """
     cdef double cpv[6]
     cpal.palDmoon( date, cpv )
     cdef np.ndarray pv = np.zeros( [6], dtype=np.float64 )
     for i in range(6):
          pv[i] = cpv[i]
     return pv


def dmoonVector(np.ndarray[double, ndim=1] date):
     """
     (x, y, z, xd, yd, zd) = dmoonVector(date)

     where date and all of the outputs are numpy arrays

     @palDmoon@
     """

     cdef int length = len(date)
     cdef np.ndarray x_out = np.empty(length, dtype=np.float64)
     cdef np.ndarray y_out = np.empty(length, dtype=np.float64)
     cdef np.ndarray z_out = np.empty(length, dtype=np.float64)
     cdef np.ndarray xd_out = np.empty(length, dtype=np.float64)
     cdef np.ndarray yd_out = np.empty(length, dtype=np.float64)
     cdef np.ndarray zd_out = np.empty(length, dtype=np.float64)
     cdef double cpv [6]

     for ii in range(length):
         cpal.palDmoon(date[ii], cpv)
         x_out[ii] = cpv[0]
         y_out[ii] = cpv[1]
         z_out[ii] = cpv[2]
         xd_out[ii] = cpv[3]
         yd_out[ii] = cpv[4]
         zd_out[ii] = cpv[5]

     return (x_out, y_out, z_out, xd_out, yd_out, zd_out)

def dmxm(np.ndarray[double, ndim=2] a not None, np.ndarray[double, ndim=2] b not None):
    """
    c = dmxm( a, b )

    @palDmxm@
    """
    cdef double c_a[3][3]
    cdef double c_b[3][3]
    cdef double c_c[3][3]
    for i in range(3):
        for j in range(3):
            c_a[i][j] = a[i][j]
            c_b[i][j] = b[i][j]
    cpal.palDmxm(c_a, c_b, c_c)

    cdef np.ndarray c = np.zeros([3,3], dtype=np.float64)
    for i in range(3):
        for j in range(3):
            c[i][j] = c_c[i][j]
    return c

def dmxv(np.ndarray[double, ndim=2] dm not None, np.ndarray[double, ndim=1] va not None):
     """
     vb = dmxv( dm, va )

     @palDmxv@
     """
     cdef double c_dm[3][3]
     cdef double c_va[3]
     cdef double c_vb[3]
     for i in range(3):
          for j in range(3):
                c_dm[i][j] = dm[i][j]
          c_va[i] = va[i]
     cpal.palDmxv(c_dm, c_va, c_vb)

     cdef np.ndarray vb = np.zeros([3], dtype=np.float64)
     for i in range(3):
          vb[i] = c_vb[i]
     return vb

def dm2av(np.ndarray[double, ndim=2] rmat not None):
    """
    axvec = dm2av( rmat )

    @palDm2av@
    """
    cdef double c_rmat[3][3]
    cdef double c_axvec[3]
    for i in range(3):
        for j in range(3):
            c_rmat[i][j] = rmat[i][j]
    cpal.palDm2av( c_rmat, c_axvec )

    cdef np.ndarray axvec = np.zeros([3], dtype=np.float64)
    for i in range(3):
        axvec[i] = c_axvec[i]
    return axvec

def dpav( np.ndarray[double, ndim=1] v1 not None, np.ndarray[double, ndim=1] v2 not None ):
     """
     pa = dpav( v1, v2 )

     @palDpav@
     """
     cdef double cv1[3]
     cdef double cv2[3]
     cdef double result
     for i in range(3):
          cv1[i] = v1[i]
          cv2[i] = v2[i]
     result = cpal.palDpav( cv1, cv2 )
     return result

def dpavVector( np.ndarray[double, ndim=2] v1_in not None,
                np.ndarray v2_in not None):
    """
    pa = dpavVector(v1, v2)

    where v1 is a 2-D numpy array in which each row is a different
    point on the unit sphere.  v2 is either a 1-D numpy array
    corresponding to one point on the unit sphere, in which case
    pa is a numpy array of position angles between the points in
    v1 and the point v2, or v2 is a 2-D numpy array of the same
    shape as v1, in which case pa is a numpy array of position angles
    between the points in v1 and the corresponding points in v2, i.e.

    pa[0] = dpav(v1[0], v2[0])
    pa[1] = dpav(v1[1], v2[1])
    ...
    pa[n] = dpav(v1[n], v2[n])

    @palDpav@

    Raises
    ------
    ValueError if the inputs are of the wrong shape
    """

    cdef nPts = v1_in.shape[0]
    cdef nDim2 = len((<object> v2_in).shape)

    if nDim2!=1 and v2_in.shape[0]!=nPts:
        raise ValueError("In dpavVector, v2 must either be a numpy array of " \
                         + "shape (3,) or a numpy array of shape (n, 3) "\
                         + "where n is the number of rows in v1")

    if v1_in.shape[1]!=3:
        raise ValueError("In dpavVector, v1 must be a numpy array in which " \
                         + "each row has 3 elements. " \
                         + "Your rows have %d elements." % v1_in.shape[1])

    if nDim2!=1 and v2_in.shape[1]!=3:
        raise ValueError("In dpavVector, v2 must be a numpy array in which " \
                         + "each row has 3 elements. " \
                         + "Your rows have %d elements." % v2_in.shape[1])

    if nDim2==1 and v2_in.shape[0]!=3:
        raise ValueError("The v2 you passed into dpavVector represents a single " \
                         "point with %d elements; it should have 3 elements" \
                         % v2_in.shape[0])

    cdef double cv1[3]
    cdef double cv2[3]
    cdef np.ndarray pa_out = np.ndarray(nPts, dtype=np.float64)

    if nDim2 == 1:
        for j in range(3):
            cv2[j] = v2_in[j]
        for i in range(nPts):
            for j in range(3):
                cv1[j] = v1_in[i][j]

            pa_out[i] = cpal.palDpav(cv1, cv2)
    else:
        for i in range(nPts):
            for j in range(3):
                cv1[j] = v1_in[i][j]
                cv2[j] = v2_in[i][j]

            pa_out[i] = cpal.palDpav(cv1, cv2)

    return pa_out

def dr2af( int ndp, double angle ):
     """
     (sign, id, im, is, frac) = dr2af( ndp, angle )
     @palDr2af@
     """
     cdef char csign = b' '
     cdef int idmsf[4]
     cpal.palDr2af( ndp, angle, &csign, idmsf )
     sign = chr(csign)
     return ( sign, idmsf[0], idmsf[1], idmsf[2], idmsf[3] )

def dr2afVector( int ndp, np.ndarray[double, ndim=1] angle not None):
    """
    (sign, id, im, is, frac) = dr2afVector( ndp, angle)

    where angle is a numpy array of angles in radians
    and all outputs are numpy arrays of the same size

    @palDr2af@
    """
    cdef int length = len(angle)
    cdef np.ndarray sign_out = np.ndarray(length, dtype=(unicode, 1))
    cdef np.ndarray id_out = np.ndarray(length, dtype=int)
    cdef np.ndarray im_out = np.ndarray(length, dtype=int)
    cdef np.ndarray isec_out = np.ndarray(length, dtype=int)
    cdef np.ndarray frac_out = np.ndarray(length, dtype=int)
    cdef char csign = b' '
    cdef int idmsf[4]

    for ii in range(length):
        cpal.palDr2af(ndp, angle[ii], &csign, idmsf)
        sign_out[ii] = chr(csign)
        id_out[ii] = idmsf[0]
        im_out[ii] = idmsf[1]
        isec_out[ii] = idmsf[2]
        frac_out[ii] = idmsf[3]

    return (sign_out, id_out, im_out, isec_out, frac_out)

def dr2tf( int ndp, double angle ):
     """
     (sign, ih, im, is, frac) = dr2tf( ndp, angle )
     @palDr2tf@
     """
     cdef char csign = b' '
     cdef int ihmsf[4]
     cpal.palDr2tf( ndp, angle, &csign, ihmsf )
     sign = chr(csign)
     return ( sign, ihmsf[0], ihmsf[1], ihmsf[2], ihmsf[3] )

def dr2tfVector( int ndp, np.ndarray[double, ndim=1] angle not None):
    """
    (sign, ih, im, is, frac) = dr2tfVector( ndp, angle)

    where angle is a numpy array of angles in radians
    and all outputs are numpy arrays of the same size

    @palDr2tf@
    """
    cdef int length = len(angle)
    cdef np.ndarray sign_out = np.ndarray(length, dtype=(unicode, 1))
    cdef np.ndarray ih_out = np.ndarray(length, dtype=int)
    cdef np.ndarray im_out = np.ndarray(length, dtype=int)
    cdef np.ndarray isec_out = np.ndarray(length, dtype=int)
    cdef np.ndarray frac_out = np.ndarray(length, dtype=int)
    cdef char csign = b' '
    cdef int idmsf[4]

    for ii in range(length):
        cpal.palDr2tf(ndp, angle[ii], &csign, idmsf)
        sign_out[ii] = chr(csign)
        ih_out[ii] = idmsf[0]
        im_out[ii] = idmsf[1]
        isec_out[ii] = idmsf[2]
        frac_out[ii] = idmsf[3]

    return (sign_out, ih_out, im_out, isec_out, frac_out)

def drange( double angle ):
     """
     a = drange( angle )

     @palDrange@
     """
     return cpal.palDrange( angle )

def drangeVector(np.ndarray[double, ndim=1] angle not None):
     """
     a = drangeVector(angle)

     where angle and a are numpy arrays

     @palDrange@
     """
     cdef int length = len(angle)
     cdef np.ndarray angle_out = np.empty(length, dtype=np.float64)
     for ii in range(length):
         angle_out[ii] = cpal.palDrange(angle[ii])
     return angle_out

def dranrm( double angle ):
     """
     a = dranrm( angle )

     @palDranrm@
     """
     return cpal.palDranrm( angle )

def dranrmVector(np.ndarray[double, ndim=1] angle not None):
    """
    a = dranrmVector( angle )

    where a and angle are numpy arrays

    @palDranrm@
    """
    cdef int length=len(angle)
    cdef np.ndarray angle_out = np.empty(length, dtype=np.float64)
    for ii in range(length):
        angle_out[ii] = cpal.palDranrm(angle[ii])
    return angle_out

def ds2tp( double ra, double dec, double raz, double decz ):
     """
     (xi, eta) = ds2tp( ra, dec, raz, decz )

     @palDs2tp@

     Raises
     ------
     ValueError: Bad input arguments
     """
     cdef double xi
     cdef double eta
     cdef int j
     cpal.palDs2tp( ra, dec, raz, decz, &xi, &eta, &j )
     if j==0:
          return (xi, eta)
     elif j==1:
          raise ValueError( "Star too far from axis" )
     elif j==2:
          raise ValueError( "Antistar on tangent plane" )
     else:
          raise ValueError( "Antistart too far from axis" )

def ds2tpVector( np.ndarray[double, ndim=1] ra not None,
                 np.ndarray[double, ndim=1] dec not None,
                 double raz, double decz ):
     """
     (xi, eta) = ds2tpVector( ra, dec, raz, decz )

     @palDs2tp@

     Raises
     ------
     ValueError in input arrays are not of the same length

     If bad input values are given, NaNs will be placed in
     the output arrays
     """
     cdef int i
     cdef int j
     cdef int length=len(ra)

     if len(dec)!=length:
         raise ValueError("You did not pass as many RAs as Decs " \
                          + "to ds2tpVector")

     cdef double xi
     cdef double eta
     cdef np.ndarray xiout=np.zeros(length,dtype=np.float64)
     cdef np.ndarray etaout=np.zeros(length,dtype=np.float64)

     for i in range(length):
          cpal.palDs2tp( ra[i], dec[i], raz, decz, &xi, &eta, &j )
          if j==0:
               xiout[i] = xi
               etaout[i] = eta
          else:
              xiout[i] = np.nan
              etaout[i] = np.nan

     return xiout, etaout

def dsep( double a1, double b1, double a2, double b2 ):
     """
     s = dsep( a1, b1, a2, b2 )

     @palDsep@
     """
     return cpal.palDsep( a1, b1, a2, b2 )

def dsepVector(np.ndarray a1, np.ndarray b1, np.ndarray a2, np.ndarray b2):
     """
     s = dsepVector( a1, b1, a2, b2 )

     where a1, b1, a2, b2, and s are numpy arrays

     @palDsep@

     Raises
     ------
     ValueError if input arrays have different lengths
     """
     cdef int length=len(a1)

     if len(b1)!=length or len(a2)!=length or len(b2)!=length:
         raise ValueError("The arrays you passed to dsepVector " \
                          + "are not of the same length")

     cdef np.ndarray dd = np.zeros(length, dtype=np.float64)
     cdef double d
     cdef int i
     for i in range(length):
          d = cpal.palDsep(a1[i], b1[i], a2[i], b2[i])
          dd[i] = d
     return dd

def dsepv( np.ndarray[double, ndim=1] v1 not None, np.ndarray[double, ndim=1] v2 not None ):
     """
     s = dsepv( v1, v2 )

     @palDsepv@
     """
     cdef double cv1[3]
     cdef double cv2[3]
     cdef double result
     for i in range(3):
          cv1[i] = v1[i]
          cv2[i] = v2[i]
     result = cpal.palDsepv( cv1, cv2 )
     return result

def dsepvVector(np.ndarray[double, ndim=2] v1 not None, np.ndarray v2 not None):
    """
    s = dsepvVector(v1, v2)

    where v1 is a 2-D numpy array in which each row is a different 3-vector.
    v2 is either a 1-D numpy array representing a 3-vector, in which case
    s is a numpy array of the angular separation between each row of v1 and
    v2, or v2 is a 2-D numpy array with the same shape as v1, in which case
    s is a numpy array ofthe separation between the points in v1 and the
    corresponding points in v2, i.e.

    s[0] = dsepv(v1[0], v2[0])
    s[1] = dsepv(v1[1], v2[1])
    ...
    s[n] = dsepv(v1[n], v2[n])

    @palDsepv@

    Raises
    ------
    ValueError if the inputs are not formatted correctly
    """

    cdef int nPts = v1.shape[0]
    cdef int nDim2 = len((<object> v2).shape)

    if nDim2!=1 and v2.shape[0]!=nPts:
        raise ValueError("In dsepvVector v2 must either have one row or " \
                         + "n rows, where n is the number of rows in v1")

    if v1.shape[1]!=3:
        raise ValueError("In dsepvVector v1 must be a numpy array " \
                         + "in which each row corresponds to a direction " \
                         + "in 3-D space.  Your v1 has %d columns." % v1.shape[1])

    if nDim2!=1 and v2.shape[1]!=3:
        raise ValueError("In dsepvVector v2 must be a numpy array " \
                         + "in which each row corresponds to a direction " \
                         + "in 3-D space.  Your v2 has %d columns." % v2.shape[1])

    if nDim2==1 and v2.shape[0]!=3:
        raise ValueError("The v2 you passed to dsepvVector is a single point " \
                          + "with %d elements; it should have 3 elements" % v2.shape[0])


    cdef double cv1[3]
    cdef double cv2[3]

    cdef np.ndarray sep = np.empty(nPts, dtype=np.float64)

    if nDim2==1:
        for i in range(3):
            cv2[i] = v2[i]
        for i in range(nPts):
            for j in range(3):
                cv1[j] = v1[i][j]
            sep[i] = cpal.palDsepv(cv1, cv2)
    else:
        for i in range(nPts):
            for j in range(3):
                cv1[j] = v1[i][j]
                cv2[j] = v2[i][j]
            sep[i] = cpal.palDsepv(cv1, cv2)

    return sep


def dt( double epoch ):
     """
     d = dt( epoch )

     @palDt@
     """
     return cpal.palDt( epoch )

def dtf2d( int ihour, int imin, double sec ):
     """
     days = dtf2d( ihour, imin, sec )

     @palDtf2d@

     Raises
     ------
     ValueError: Arguments out of range.
     """
     cdef double days
     cdef int j
     cpal.palDtf2d( ihour, imin, sec, &days, &j )
     if j==0:
          return days
     else:
          bad = ""
          if j==1:
               bad = "Hour argument outside range 0-23"
          elif j==2:
               bad = "Minute argument outside range 0-59"
          else:
               bad = "Arcsec argument outside range 0-59.9999..."
          raise ValueError( bad )

def dtf2dVector( np.ndarray[long, ndim=1] ihour not None,
                 np.ndarray[long, ndim=1] imin not None,
                 np.ndarray[double, ndim=1] sec not None):
    """
    days = dtf2dVector(ihour, imin, sec)

    where ihour, imin, sec, and days are all numpy arrays

    @palDtf2d@

    Raises
    ------
    Raises a ValueError if inputs have different sizes.

    Bad values (hour outside of 0-23 range; minute outside of
    0-59 range, or second outside of 0-59.9999 range)
    will result in numpy.NaNs in the output.
    """
    cdef int length = len(ihour)

    if len(imin)!=length:
        raise ValueError("In dtf2dVector, imin does not have the same length " \
                         + "as ihour")

    if len(sec)!=length:
        raise ValueError("In dtf2dVector, sec does not have the same length " \
                         + "as ihour")

    cdef np.ndarray days = np.empty(length, dtype=np.float64)
    cdef int flag
    cdef double value

    for ii in range(length):
        cpal.palDtf2d( ihour[ii], imin[ii], sec[ii], &value, &flag )
        if flag==0:
            days[ii] = value
        else:
            days[ii] = np.nan

    return days


def dtf2r( int ihour, int imin, double sec ):
     """
     rad = dtf2r( ihour, imin, sec )

     @palDtf2r@

     Raises
     ------
     ValueError: Arguments out of range.
     """
     cdef double rad
     cdef int j
     cpal.palDtf2r( ihour, imin, sec, &rad, &j )
     if j==0:
          return rad
     else:
          bad = ""
          if j==1:
               bad = "Hour argument outside range 0-23"
          elif j==2:
               bad = "Minute argument outside range 0-59"
          else:
               bad = "Arcsec argument outside range 0-59.9999..."
          raise ValueError( bad )

def dtf2rVector(np.ndarray[long, ndim=1] ihour, np.ndarray[long, ndim=1] imin,
                np.ndarray[double, ndim=1] sec):
    """
    rad = dtf2rVector(ihour, imin, sec)

    where ihour, imin, sec, and rad are all numpy arrays

    @palDtf2r@

    Raises
    ------
    ValueError if input arrays do not have the same length

    If bad values (ihour outside of 0-23, imin outside of 0-59,
    sec outside of 0-59.999999) are passed in, these result in
    numpy.NaNs in the output.
    """

    cdef int length = len(ihour)

    if len(imin)!=length:
        raise ValueError("In dtf2rVector, imin does not have the same length " \
                         + "as ihour")

    if len(sec)!=length:
        raise ValueError("In dtf2rVector, sec does not have the same length " \
                         + "as ihour")

    cdef np.ndarray radians = np.empty(length, dtype=np.float64)
    cdef int flag
    cdef double value

    for ii in range(length):
        cpal.palDtf2r( ihour[ii], imin[ii], sec[ii], &value, &flag )
        if flag==0:
            radians[ii] = value
        else:
            radians[ii] = np.nan

    return radians



def dtp2s( double xi, double eta, double raz, double decz):
     """
     (ra,dec) = dtp2s( xi, eta, raz, decz )

     @palDtp2s@
     """
     cdef double ra
     cdef double dec
     cpal.palDtp2s( xi, eta, raz, decz, &ra, &dec )
     return (ra,dec)

def dtp2sVector(np.ndarray[double, ndim=1] xi not None,
                np.ndarray[double, ndim=1] eta not None,
                double raz, double decz):
     """
     (ra, dec) = dtp2sVector( xi, eta, raz, decz)

     where xi, eta, ra, and dex (NOT raz, decz) are numpy arrays

     @palDtp2s@

     Raises
     ------
     ValueError if input arrays are of different lengths
     """

     cdef int length = len(xi)

     if len(eta)!=length:
          raise ValueError("The arrays of tangent plane coordinates " \
                           + "you passed to dtp2sVector are of " \
                           + "different lengths")

     cdef np.ndarray ra_out = np.empty(length, dtype=np.float64)
     cdef np.ndarray dec_out = np.empty(length, dtype=np.float64)

     cdef double ra
     cdef double dec

     for ii in range(length):
         cpal.palDtp2s(xi[ii], eta[ii], raz, decz, &ra, &dec)
         ra_out[ii] = ra
         dec_out[ii] = dec

     return (ra_out, dec_out)


def dtps2c( double xi, double eta, double ra, double dec ):
     """
     (raz1, decz1, raz2, decz2) = dtps2c( xi, eta, ra, dec )

     @palDtps2c@

     Python notes
     ------------
     If some solutions are missing None is returned for each value.
     """
     cdef double raz1
     cdef double decz1
     cdef double raz2
     cdef double decz2
     cdef int n
     cpal.palDtps2c( xi, eta, ra, dec, &raz1, &decz1, &raz2, &decz2, &n )
     if n==0:
          return (None, None, None, None)
     elif n==1:
          return (raz1, decz1, None, None)
     else:
          return (raz1, decz1, raz2, decz2 )

def dtt( double dju ):
     """
     d = dtt( dju )

     @palDtt@
     """
     return cpal.palDtt( dju )

def dvn( np.ndarray[double, ndim=1] v not None):
    """
    (uv, vm) = dvn( v )

    @palDvn@
    """
    cdef double c_v[3]
    cdef double c_uv[3]
    cdef double vm
    for i in range(3):
        c_v[i] = v[i]
    cpal.palDvn( c_v, c_uv, &vm )

    cdef np.ndarray uv = np.zeros([3], dtype=np.float64)
    for i in range(3):
        uv[i] = c_uv[i]
    return (uv, vm)

def dvxv( np.ndarray[double, ndim=1] va not None, np.ndarray[double, ndim=1] vb not None):
    """
    vc = dvxv( va, vb )

    @palDvxv@
    """
    cdef double c_va[3]
    cdef double c_vb[3]
    cdef double c_vc[3]
    for i in range(3):
        c_va[i] = va[i]
        c_vb[i] = vb[i]
    cpal.palDvxv( c_va, c_vb, c_vc )

    cdef np.ndarray vc = np.zeros([3], dtype=np.float64)
    for i in range(3):
        vc[i] = c_vc[i]
    return vc

def ecleq( double dl, double db, double date ):
     """
     (dr, dd) = ecleq( dl, db, date )

     @palEcleq@
     """
     cdef double dr
     cdef double dd
     cpal.palEcleq( dl, db, date, &dr, &dd )
     return (dr, dd)

def ecleqVector(np.ndarray[double, ndim=1] dl not None,
                np.ndarray[double, ndim=1] db not None,
                double date):
    """
    (dr, dd) = ecleqVector(dl, db, date)

    where dl, db, dr, and dd are all numpy arrays

    @palEcleq@

    Raises
    ------
    ValueError if input arrays are of different lengths
    """

    cdef int length = len(dl)

    if len(db)!=length:
        raise ValueError("You did not pass the same number of " \
                         + "longitudes as latitudes to ecleqVector")

    cdef np.ndarray rr_out = np.empty(length, dtype=np.float64)
    cdef np.ndarray dd_out = np.empty(length, dtype=np.float64)

    cdef double rr
    cdef double dd

    for ii in range(length):
        cpal.palEcleq(dl[ii], db[ii], date, &rr, &dd)
        rr_out[ii] = rr
        dd_out[ii] = dd

    return (rr_out, dd_out)

def ecmat( double date ):
     """
     rmat = ecmat( date )

     @palEcmat@
     """
     cdef double crmat[3][3]
     cpal.palEcmat( date, crmat )
     cdef np.ndarray rmat = np.zeros( [3,3], dtype=np.float64 )
     for i in range(3):
          for j in range(3):
               rmat[i,j] = crmat[i][j]
     return rmat

def el2ue( double date, int jform, double epoch, double orbinc,
            double anode, double perih, double aorq,  double e,
            double aorl, double dm ):
    """
    u = el2ue( date, jform, epoch, orbinc, anode, perih, aorq, e, aorl, dm )

    @palEl2ue@

    Raises
    ------
    ValueError: Illegal input arguments
    ArithmeticError: Numerical error
    """
    cdef int jstat
    cdef double cu[13]
    cpal.palEl2ue( date, jform, epoch, orbinc, anode, perih, aorq,
                    e, aorl, dm, cu, &jstat )
    if jstat == -1:
        raise ValueError( "Illegal jform" )
    elif jstat == -2:
        raise ValueError( "Illegal e" )
    elif jstat == -3:
        raise ValueError( "Illegal aorq" )
    elif jstat == -4:
        raise ValueError( "Illegal dm" )
    elif jstat == -5:
        raise ArithmeticError( "Numerical error" )

    cdef np.ndarray u = np.zeros( [13], dtype=np.float64 )
    for i in range(13):
        u[i] = cu[i]
    return u

def epb( double date ):
     """
     e = epb( date )

     @palEpb@
     """
     return cpal.palEpb(date)

def epb2d( double epb ):
     """
     d = epb2d( epb )

     @palEpb2d@
     """
     return cpal.palEpb2d(epb)

def epco( k0, k, double e ):
     """
     e = epco( k0, k, e )

     @palEpco@
     """
     k0_bytes = k0.encode('ascii')
     k_bytes = k.encode('ascii')
     cdef char * ck0 = k0_bytes
     cdef char * ck = k_bytes
     return cpal.palEpco( ck0[0], ck[0], e )

def epj( double date ):
     """
     e = epj( date )

     @palEpj@
     """
     return cpal.palEpj(date)

def epjVector( np.ndarray[double, ndim=1] date not None):
    """
    e = epjVector( date )

    where e and date are both numpy arrays

    @palEpj@
    """
    cdef int length = len(date)
    cdef np.ndarray epj = np.empty(length, dtype=np.float64)
    for ii in range(length):
        epj[ii] = cpal.palEpj(date[ii])
    return epj


def epj2d( double epj ):
     """
     d = epj2d( epj )

     @palEpj2d@
     """
     return cpal.palEpj2d(epj)

def epj2dVector(np.ndarray[double, ndim=1] epj not None):
    """
    d = epj2dVector(epj)

    where both d and epj are numpy arrays

    @palEpj2d@
    """
    cdef int length = len(epj)
    cdef np.ndarray mjd = np.empty(length, dtype=np.float64)
    for ii in range(length):
        mjd[ii] = cpal.palEpj2d(epj[ii])
    return mjd

# epv goes here
def epv( double date ):
    """
    (ph, vh, pb, vb) = epv( date )

    @palEpv@
    """
    cdef double cph[3]
    cdef double cvh[3]
    cdef double cpb[3]
    cdef double cvb[3]
    cpal.palEpv( date, cph, cvh, cpb, cvb )

    cdef np.ndarray ph = np.zeros( [3], dtype=np.float64 )
    cdef np.ndarray vh = np.zeros( [3], dtype=np.float64 )
    cdef np.ndarray pb = np.zeros( [3], dtype=np.float64 )
    cdef np.ndarray vb = np.zeros( [3], dtype=np.float64 )
    for i in range(3):
        ph[i] = cph[i]
        vh[i] = cvh[i]
        pb[i] = cpb[i]
        vb[i] = cvb[i]
    return (ph, vh, pb, vb)

def eqecl( double dr, double dd, double date ):
     """
     (dl, db) = eqecl( dr, dd, date )

     @palEqecl@
     """
     cdef double dl
     cdef double db
     cpal.palEqecl( dr, dd, date, &dl, &db )
     return (dl, db)

def eqeclVector(np.ndarray[double, ndim=1] dr not None,
                np.ndarray[double, ndim=1] dd not None,
                double date):
    """
    (dl, db) = eqeclVector(dr, dd, date)

    where dr, dd, dl, and db are all numpy arrays

    @palEqEcl@

    Raises
    ------
    ValueError if input arrays have different lengths
    """

    cdef int length = len(dr)

    if len(dd)!=length:
        raise ValueError("You did not pass the same number of Decs " \
                         + "as RAs to eqeclVector")

    cdef np.ndarray dl_out = np.empty(length, dtype=np.float64)
    cdef np.ndarray db_out = np.empty(length, dtype=np.float64)

    cdef double dl
    cdef double db

    for ii in range(length):
        cpal.palEqecl(dr[ii], dd[ii], date, &dl, &db)
        dl_out[ii] = dl
        db_out[ii] = db

    return (dl_out, db_out)

def eqeqx( double date ):
     """
     eq = eqeqx( date )

     @palEqeqx@
     """
     return cpal.palEqeqx( date )

def eqeqxVector(np.ndarray date):
     """
     eq = eqeqxVector( date )

     where date and eq are numpy arrays

     @palEqeqx@
     """
     cdef int length=len(date)
     cdef np.ndarray eqout = np.zeros(length, dtype=np.float64)
     cdef double eq
     cdef int i
     for i in range(length):
          eq = cpal.palEqeqx(date[i])
          eqout[i] = eq
     return eqout

def eqgal( double dr, double dd ):
     """
     (dl, db) = eqgal( dr, dd )

     @palEqgal@
     """
     cdef double dl
     cdef double db
     cpal.palEqgal( dr, dd, &dl, &db )
     return (dl, db)

def eqgalVector(np.ndarray[double, ndim=1] dr not None,
                np.ndarray[double, ndim=1] dd not None):
     """
     (dl, db) = eqgalVector( dr, dd )

     where dr, dd, dl, db are all numpy arrays

     @palEqgal@

     Raises
     ------
     ValueError if input arrays are not of the same length
     """

     cdef int length=len(dr)

     if len(dd) != length:
         raise ValueError("You did not pass the same number of RAs " \
                          + "as Decs to eqgalVector")

     cdef np.ndarray dlout = np.zeros(length, dtype=np.float64)
     cdef np.ndarray dbout = np.zeros(length, dtype=np.float64)
     cdef double dl
     cdef double db

     for i in range(length):
         cpal.palEqgal(dr[i], dd[i], &dl, &db)
         dlout[i] = dl
         dbout[i] = db

     return (dlout, dbout)

def etrms( double ep ):
     """
     ev = etrms( ep )

     @palEtrms@
     """
     cdef double cev[3]
     cpal.palEtrms( ep, cev )
     cdef np.ndarray ev = np.zeros( [3], dtype=np.float64 )
     for i in range(3):
          ev[i] = cev[i]
     return ev

def evp(double date, double deqx):
    """
    (dvb, dpb, dvh, dph) = evp( date, deqx )

    @palEvp@
    """
    cdef double cdvb[3]
    cdef double cdpb[3]
    cdef double cdvh[3]
    cdef double cdph[3]
    cpal.palEvp( date, deqx, cdvb, cdpb, cdvh, cdph )

    cdef np.ndarray dvb = np.zeros( [3], dtype=np.float64 )
    cdef np.ndarray dpb = np.zeros( [3], dtype=np.float64 )
    cdef np.ndarray dvh = np.zeros( [3], dtype=np.float64 )
    cdef np.ndarray dph = np.zeros( [3], dtype=np.float64 )
    for i in range(3):
        dvb[i] = cdvb[i]
        dpb[i] = cdpb[i]
        dvh[i] = cdvh[i]
        dph[i] = cdph[i]
    return (dvb, dpb, dvh, dph)

def fk45z( double r1950, double d1950, double bepoch ):
     """
     (r2000, d2000) = fk45z( r1950, d1950, bepoch )

     @palFk45z@
     """
     cdef double r2000
     cdef double d2000
     cpal.palFk45z( r1950, d1950, bepoch, &r2000, &d2000 )
     return (r2000, d2000)

def fk45zVector(np.ndarray[double, ndim=1] r1950 not None,
                np.ndarray[double, ndim=1] d1950 not None,
                double bepoch):
    """
    (r2000, d2000) = fk45zVector(r1950, d1950, bepoch)

    where r1950, d1950, r2000, d2000 are numpy arrays

    @palFk45z@

    Raises
    ------
    ValueError if input arrays are not of the same length
    """

    cdef int length = len(r1950)

    if len(d1950)!=length:
        raise ValueError("You did not pass the same number of Decs as " \
                         + "RAs to fk45zVector")

    cdef np.ndarray r2000_out = np.empty(length, dtype=np.float64)
    cdef np.ndarray d2000_out = np.empty(length, dtype=np.float64)

    cdef double r2000
    cdef double d2000

    for ii in range(length):
        cpal.palFk45z(r1950[ii], d1950[ii], bepoch, &r2000, &d2000)
        r2000_out[ii] = r2000
        d2000_out[ii] = d2000

    return (r2000_out, d2000_out)

def fk524( double r2000, double d2000, double dr2000,
           double dd2000, double p2000, double v2000 ):
     """
     (r1950, d1950, dr1950, dd1950, p1950, v1950) = fk524( r2000, d2000, dr2000, dd2000, p2000, v2000 )

     @palFk524@
     """
     cdef double r1950
     cdef double d1950
     cdef double dr1950
     cdef double dd1950
     cdef double p1950
     cdef double v1950
     cpal.palFk524( r2000, d2000, dr2000, dd2000, p2000, v2000,
                    &r1950, &d1950, &dr1950, &dd1950,
                    &p1950, &v1950 )
     return (r1950, d1950, dr1950, dd1950, p1950, v1950 )

def fk524Vector(np.ndarray[double, ndim=1] r2000 not None,
                np.ndarray[double, ndim=1] d2000 not None,
                np.ndarray[double, ndim=1] dr2000 not None,
                np.ndarray[double, ndim=1] dd2000 not None,
                np.ndarray[double, ndim=1] p2000 not None,
                np.ndarray[double, ndim=1] v2000 not None):
    """
    (r1950, d1950, dr1950, dd1950, p1950, v1950) = fk524Vector( r2000, d2000, dr2000, dd2000, p2000, v2000 )

    where r2000, d2000, dr2000, dd2000, p2000, v2000,
    r1950, d1950, dr1950, dd1950, p1950, and v1950 are
    all numpy arrays

    @palFk524@

    Raises
    ------
    ValueError if input arrays do not all of the same length
    """
    cdef int length = len(r2000)

    if len(d2000)!= length:
        raise ValueError("You did not pass as many Decs as RAs into " \
                         + "fk524Vector")

    if len(dr2000)!= length:
        raise ValueError("You did not pass as many RA proper motions as RAs into " \
                         + "fk524Vector")

    if len(dd2000)!= length:
        raise ValueError("You did not pass as many Dec proper motions as RAs into " \
                         + "fk524Vector")

    if len(p2000)!= length:
        raise ValueError("You did not pass as many parallaxes as RAs into " \
                         + "fk524Vector")

    if len(v2000)!= length:
        raise ValueError("You did not pass as many radial velocities as RAs into " \
                         + "fk524Vector")


    cdef np.ndarray r1950_out = np.empty(length, dtype=np.float64)
    cdef np.ndarray d1950_out = np.empty(length, dtype=np.float64)
    cdef np.ndarray dr1950_out = np.empty(length, dtype=np.float64)
    cdef np.ndarray dd1950_out = np.empty(length, dtype=np.float64)
    cdef np.ndarray p1950_out = np.empty(length, dtype=np.float64)
    cdef np.ndarray v1950_out = np.empty(length, dtype=np.float64)

    cdef double r1950
    cdef double d1950
    cdef double dr1950
    cdef double dd1950
    cdef double p1950
    cdef double v1950

    for ii in range(length):
        cpal.palFk524(r2000[ii], d2000[ii], dr2000[ii], dd2000[ii],
                      p2000[ii], v2000[ii], &r1950, &d1950,
                      &dr1950, &dd1950, &p1950, &v1950)

        r1950_out[ii] = r1950
        d1950_out[ii] = d1950
        dr1950_out[ii] = dr1950
        dd1950_out[ii] = dd1950
        p1950_out[ii] = p1950
        v1950_out[ii] = v1950

    return (r1950_out, d1950_out, dr1950_out, dd1950_out,
            p1950_out, v1950_out)

def fk54z(double r2000, double d2000, double bepoch):
     """
     (r1950, d1950, dr1950, dd1950) = fk54z( r2000, d2000, bepoch )

     @palFk54z@
     """
     cdef double r1950
     cdef double d1950
     cdef double dr1950
     cdef double dd1950
     cpal.palFk54z( r2000, d2000, bepoch, &r1950, &d1950,
                    &dr1950, &dd1950 )
     return (r1950, d1950, dr1950, dd1950 )

def fk54zVector(np.ndarray[double, ndim=1] r2000 not None,
                np.ndarray[double, ndim=1] d2000 not None,
                double bepoch):
    """
    (r1950, d1950, dr1950, dd1950) = fk45zVector( r2000, d2000, bepoch )

    where r2000, d2000, and the outputs are numpy arrays

    @palFk54z@

    Raises
    ------
    Value error if input arrays have different lengths
    """

    cdef int length = len(r2000)

    if len(d2000)!=length:
        raise ValueError("You did not pass the same number of Decs as " \
                         + "RAs to fk54zVector")

    cdef np.ndarray r1950_out = np.empty(length, dtype=np.float64)
    cdef np.ndarray d1950_out = np.empty(length, dtype=np.float64)
    cdef np.ndarray dr1950_out = np.empty(length, dtype=np.float64)
    cdef np.ndarray dd1950_out = np.empty(length, dtype=np.float64)

    cdef double r1950
    cdef double d1950
    cdef double dr1950
    cdef double dd1950

    for ii in range(length):
        cpal.palFk54z(r2000[ii], d2000[ii], bepoch, \
                      &r1950, &d1950, &dr1950, &dd1950)

        r1950_out[ii] = r1950
        d1950_out[ii] = d1950
        dr1950_out[ii] = dr1950
        dd1950_out[ii] = dd1950

    return (r1950_out, d1950_out, dr1950_out, dd1950_out)

def fk5hz( double r5, double d5, double epoch):
     """
     (rh, dh) = fk5hz( r5, d5, epoch )

     @palFk5hz@
     """
     cdef double rh
     cdef double dh
     cpal.palFk5hz( r5, d5, epoch, &rh, &dh )
     return (rh, dh)

def fk5hzVector(np.ndarray[double, ndim=1] r5 not None,
                np.ndarray[double, ndim=1] d5 not None,
                double epoch):
    """
    (rh, dh) = fk5hzVector(r5, d5, epoch)

    where rh, dh, r5, and d5 are numpy arrays

    @palFk5hz@

    Raises
    ------
    ValueError if r5 and d5 have different lengths
    """

    cdef int length = len(r5)

    if len(d5)!=length:
        raise ValueError("You did not pass the same number of " \
                         + "RAs as Decs to fk5hzVector")

    cdef np.ndarray ra_out = np.empty(length, dtype=np.float64)
    cdef np.ndarray dec_out = np.empty(length, dtype=np.float64)
    cdef double rh
    cdef double dh

    for ii in range(length):
        cpal.palFk5hz( r5[ii], d5[ii], epoch, &rh, &dh )
        ra_out[ii] = rh
        dec_out[ii] = dh

    return (ra_out, dec_out)

def galeq( double dl, double db ):
     """
     (dr, dd) = galeq( dl, db )

     @palGaleq@
     """
     cdef double dr
     cdef double dd
     cpal.palGaleq( dl, db, &dr, &dd )
     return (dr, dd)

def galeqVector(np.ndarray[double, ndim=1] dl not None,
                np.ndarray [double, ndim=1] db not None):
     """
     (dr, dd) = galeqVector( dl, db )

     where dl, db, dr, and dd are numpy arrays

     @palGaleq@

     Raises
     ------
     ValueError if input arrays are of different lengths
     """
     cdef int length=len(dl)

     if len(db)!=length:
         raise ValueError("You did not pass the same number of " \
                          + "longitudes as latitudes to galeqVector")

     cdef np.ndarray drout = np.zeros(length, dtype=np.float64)
     cdef np.ndarray ddout = np.zeros(length, dtype=np.float64)
     cdef double dr
     cdef double dd
     cdef int i
     for i in range(length):
          cpal.palGaleq(dl[i], db[i], &dr, &dd)
          drout[i] = dr
          ddout[i] = dd
     return (drout, ddout)

def galsup( double dl, double db ):
     """
     (dsl, dsb) = galsup( dl, db )

     @palGalsup@
     """
     cdef double dsl
     cdef double dsb
     cpal.palGalsup( dl, db, &dsl, &dsb )
     return (dsl, dsb)

def galsupVector(np.ndarray[double, ndim=1] dl not None,
                 np.ndarray[double, ndim=1] db not None):
    """
    (dsl, dsb) = galsupVector( dl, db )

    where dl, db, dsl, and dsb are all numpy arrays

    @palGalsup@

    Raises
    ------
    ValueError if input arrays have different lengths
    """

    cdef int length=len(dl)

    if len(db)!=length:
        raise ValueError("You did not pass the same number of " \
                         + "Galactic longitudes as Galactic latitudes " \
                         + "into galsupVector")

    cdef np.ndarray dsl_out = np.empty(length, dtype=np.float64)
    cdef np.ndarray dsb_out = np.empty(length, dtype=np.float64)
    cdef double dsl
    cdef double dsb

    for ii in range(length):
        cpal.palGalsup(dl[ii], db[ii], &dsl, &dsb)
        dsl_out[ii] = dsl
        dsb_out[ii] = dsb

    return (dsl_out, dsb_out)

def ge50( double dl, double db ):
     """
     (rd, dd) = ge50( dl, db )

     @palGe50@
     """
     cdef double dr
     cdef double dd
     cpal.palGe50( dl, db, &dr, &dd )
     return (dr, dd)

def ge50Vector( np.ndarray[double, ndim=1] dl not None,
                np.ndarray[double, ndim=1] db not None):
    """
    (rd, dd) = ge50( dl, db )

    where dl, db, rd, and dd are all numpy arrays

    @palGe50@

    Raises
    ------
    ValueError if input arrays have different lengths
    """

    cdef int length = len(dl)

    if len(db)!=length:
        raise ValueError("You did not pass the same number of " \
                         + "Galactic longitudes as Galactic latitudes " \
                         + "into ge50Vector")

    cdef np.ndarray dr_out = np.empty(length, dtype=np.float64)
    cdef np.ndarray dd_out = np.empty(length, dtype=np.float64)
    cdef double dr
    cdef double dd

    for ii in range(length):
        cpal.palGe50(dl[ii], db[ii], &dr, &dd)
        dr_out[ii] = dr
        dd_out[ii] = dd

    return (dr_out, dd_out)

def geoc( double p, double h ):
     """
     (r, z) = geoc( p, h )

     @palGeoc@
     """
     cdef double r
     cdef double z
     cpal.palGeoc( p, h, &r, &z )
     return (r, z)

def gmst( double ut1 ):
     """
     t = gmst( ut1 )

     @palGmst@
     """
     return cpal.palGmst( ut1 )

def gmstVector(np.ndarray ut1):
     """
     t = gmstVector( ut1 )

     where ut1 and t are both numpy arrays

     @palGmst@
     """
     cdef int length = len(ut1)
     cdef np.ndarray gmout = np.zeros(length, dtype=np.float64)
     cdef double gm
     cdef int i
     for i in range(length):
          gm = cpal.palGmst(ut1[i])
          gmout[i] = gm
     return gmout

def gmsta( double date, double ut1 ):
     """
     t = gmsta( date, ut1 )

     @palGmsta@
     """
     return cpal.palGmsta( date, ut1 )

def gmstaVector(np.ndarray date, np.ndarray ut1):
     """
     t = gmsta( date, ut1 )

     where date, ut1, and t are all numpy arrays

     @palGmsta@

     Raises
     ------
     ValueError if input arrays are of different lengths
     """
     cdef int length = len(date)

     if len(ut1)!=length:
         raise ValueError("The arrays you passed into gmstaVector " \
                          + "are of different lengths")

     cdef np.ndarray gmout = np.zeros(length, dtype=np.float64)
     cdef int i
     for i in range(length):
          gmout[i] = cpal.palGmsta(date[i], ut1[i])
     return gmout

def hfk5z( double rh, double dh, double epoch ):
     """
     (r5, d5, dr5, dd5) = hfk5z( rh, dh, epoch )

     @palHfk5z@
     """
     cdef double r5
     cdef double d5
     cdef double dr5
     cdef double dd5
     cpal.palHfk5z( rh, dh, epoch, &r5, &d5, &dr5, &dd5 )
     return (r5, d5, dr5, dd5)

def hfk5zVector(np.ndarray[double, ndim=1] rh not None,
                np.ndarray[double, ndim=1] dh not None,
                double epoch):
    """
    (r5, d5, dr5, dd5) = hfk5zVector( rh, dh, epoch)

    where rh, dh, r5, d5, dr5 and dd5 are all numpy arrays

    @palHfk5z@

    Raises
    ------
    ValueError if input RA and Dec arrays do not have the same length
    """

    cdef int length = len(rh)

    if len(dh)!=length:
        raise ValueError("You did not pass the same number of RAs as " \
                         + "Decs to hfk5zVector")

    cdef np.ndarray ra_out = np.empty(length, dtype=np.float64)
    cdef np.ndarray dec_out = np.empty(length, dtype=np.float64)
    cdef np.ndarray dr_out = np.empty(length, dtype=np.float64)
    cdef np.ndarray dd_out = np.empty(length, dtype=np.float64)

    cdef double r5
    cdef double d5
    cdef double dr5
    cdef double dd5

    for ii in range(length):
        cpal.palHfk5z( rh[ii], dh[ii], epoch, &r5, &d5, &dr5, &dd5 )
        ra_out[ii] = r5
        dec_out[ii] = d5
        dr_out[ii] = dr5
        dd_out[ii] = dd5

    return (ra_out, dec_out, dr_out, dd_out)

# We need to return the sign in order to work out whether -0
# is a negative integer (important when parsing "-0 22 33.0"
# sexagesimal format. I'm assuming that no-one is going to really
# use the intin() function.
# We also need to handle overflow. If we raise an exception on
# overflow we can't continue to traverse the string so we have
# to return the position to continue but then document the
# magic value of LONG_MAX and LONG_MIN somehow but where would
# that go? For now raise an exception.
def intin( string, int nstrt ):
     """
     (ireslt, nstrt, sign) = intin( string, nstrt )

     @palIntin@

     Raises
     ------
     OverflowError: Integer too large to fit in 32-bit integer
     """
     cdef long ireslt
     cdef int j
     string_bytes = string.encode('ascii')
     cdef char * cstring = string_bytes
     cpal.palIntin( cstring, &nstrt, &ireslt, &j)
     sign = 0
     if j==0:
          sign = 1
     elif j==-1:
          sign = -1
     elif j==1:
          return (None, nstrt, None)
     else:
          raise OverflowError( "Integer too large for pal.intin() function" )
     return ( ireslt, nstrt, sign )

def map(double rm, double dm, double pr, double pd, double px, double rv, double eq, double date ):
     """
     (ra, da) = map( rm, dm, pr, pd, px, rv, eq, date )

     @palMap@
     """
     cdef double ra
     cdef double da
     cpal.palMap( rm, dm, pr, pd, px, rv, eq, date, &ra, &da )
     return (ra, da)

def mappa( double eq, double date ):
     """
     amprms = mappa( eq, date )

     @palMappa@
     """
     cdef double camprms[21]
     cdef np.ndarray amprms = np.zeros( [21], dtype=np.float64 )
     cpal.palMappa( eq, date, camprms )
     for i in range(21):
          amprms[i] = camprms[i]
     return amprms

def mapqk( double rm, double dm, double pr, double pd, double px, double rv, np.ndarray[double, ndim=1] amprms not None):
     """
     (ra, da) = mapqk( rm, dm, pr, pd, px, rv, amprms )

     @palMapqk@
     """
     cdef double ra
     cdef double da
     cdef double camprms[21]
     for i in range(21):
          camprms[i] = amprms[i]
     cpal.palMapqk( rm, dm, pr, pd, px, rv, camprms, &ra, &da )
     return (ra, da)

def mapqkVector( np.ndarray[double, ndim=1] rm not None,
                 np.ndarray[double, ndim=1] dm not None,
                 np.ndarray[double, ndim=1] pr not None,
                 np.ndarray[double, ndim=1] pd not None,
                 np.ndarray[double, ndim=1] px not None,
                 np.ndarray[double, ndim=1] rv not None,
                 np.ndarray[double, ndim=1] amprms not None):
     """
     (ra, da) = mapqkVector( rm, dm, pr, pd, px, rv, amprms )

     where rm, dm, pr, pd, px, rv, ra, and da are all numpy arrays

     @palMapqk@

     Raises
     ------
     ValueError if input arrays are not of the same length
     """
     cdef int i
     cdef int length=len(rm)

     if len(dm)!=length or len(pr)!=length or len(pd)!=length \
     or len(px)!=length or len(rv)!=length:
         raise ValueError("The arrays you passed into mapqkVector " \
                          + "are not of the same length")

     cdef np.ndarray raout=np.zeros(length,dtype=np.float64)
     cdef np.ndarray decout=np.zeros(length,dtype=np.float64)
     cdef double ra
     cdef double da
     cdef double camprms[21]
     for i in range(21):
          camprms[i] = amprms[i]
     for i in range(length):
         cpal.palMapqk( rm[i], dm[i], pr[i], pd[i], px[i], rv[i], camprms, &ra, &da )
         raout[i]=ra
         decout[i]=da

     return (raout, decout)

def mapqkz( double rm, double dm, np.ndarray[double, ndim=1] amprms not None):
     """
     (ra, da) = mapqkz( rm, dm, amprms )

     @palMapqkz@
     """
     cdef double ra
     cdef double da
     cdef double camprms[21]
     for i in range(21):
          camprms[i] = amprms[i]
     cpal.palMapqkz( rm, dm, camprms, &ra, &da )
     return (ra, da)


def mapqkzVector( np.ndarray[double, ndim=1] rm not None,
                  np.ndarray[double, ndim=1] dm not None,
                  np.ndarray[double, ndim=1] amprms not None):
     """
     (ra, da) = mapqkzVector( rm, dm, amprms )

     @palMapqkz@

     Raises
     ------
     ValueError if input arrays are not of the same length
     """
     cdef int i
     cdef int length=len(rm)

     if len(dm)!=length:
         raise ValueError("You did not pass as many RAs as Decs to " \
                           + "mapqkzVector")

     cdef np.ndarray raout=np.zeros(length,dtype=np.float64)
     cdef np.ndarray decout=np.zeros(length,dtype=np.float64)
     cdef double ra
     cdef double da
     cdef double camprms[21]
     for i in range(21):
          camprms[i] = amprms[i]
     for i in range(length):
          cpal.palMapqkz( rm[i], dm[i], camprms, &ra, &da )
          raout[i]=ra
          decout[i]=da
     return (raout, decout)

def nut( double date ):
     """
     rmatn = nut( date )

     @palNut@
     """
     cdef double crmatn[3][3]
     cpal.palNut( date, crmatn )
     cdef np.ndarray rmatn = np.zeros( [3,3], dtype=np.float64 )
     for i in range(3):
          for j in range(3):
               rmatn[i,j] = crmatn[i][j]
     return rmatn

def nutc( double date):
     """
     (dpsi, deps, eps0) = nutc( date )

     @palNutc@
     """
     cdef double dpsi
     cdef double deps
     cdef double eps0
     cpal.palNutc( date, &dpsi, &deps, &eps0 )
     return (dpsi, deps, eps0)

def oap( type, double ob1, double ob2, double date,
         double dut, double elongm, double phim, double hm,
         double xp, double yp, double tdk, double pmb,
         double rh, double wl, double tlr ):
     """
     (rap, dap) = oap( ob1, ob2, date, dut, elongm, phim,
                    hm, xp, yp, tdk, pmb, rh, wl, tlr )

     @palOap@
     """
     cdef double rap
     cdef double dap
     byte_string = type.encode('ascii')
     cdef char * ctype = byte_string
     cpal.palOap( ctype, ob1, ob2, date, dut, elongm, phim,
                  hm, xp, yp, tdk, pmb, rh, wl, tlr, &rap, &dap )
     return ( rap, dap )

def oapqk( type, double ob1, double ob2, np.ndarray[double, ndim=1] aoprms not None ):
     """
     (rap, dap) = oapqk( type, ob1, ob2, aoprms )

     @palOapqk@
     """
     cdef double rap
     cdef double dap
     byte_string = type.encode('ascii')
     cdef char * ctype = byte_string
     cdef double caoprms[14]
     for i in range(14):
          caoprms[i] = aoprms[i]
     cpal.palOapqk( ctype, ob1, ob2, caoprms, &rap, &dap )
     return ( rap, dap )

def oapqkVector(type,
                np.ndarray[double, ndim=1] ob1 not None,
                np.ndarray[double, ndim=1] ob2 not None,
                np.ndarray[double, ndim=1] aoprms not None):
    """
    (rap, dap) = oapqkVector( type, ob1, ob2, aoprms )

    where ob1, ob2, rap, dap are all numpy arrays

    @palOapqk@

    Raises
    ------
    ValueError if input coordinate arrays are of different length
    """

    cdef int length = len(ob1)
    if len(ob2)!=length:
        raise ValueError("The observed coordinate arrays you passed into " \
                         + "oapqkVector are of different lengths")

    cdef np.ndarray ra_out = np.empty(length, dtype=np.float64)
    cdef np.ndarray dec_out = np.empty(length, dtype=np.float64)

    byte_string = type.encode('ascii')

    cdef char * ctype = byte_string
    cdef double caoprms[14]

    for i in range(14):
        caoprms[i] = aoprms[i]

    cdef double rap
    cdef double dap

    for ii in range(length):
        cpal.palOapqk( ctype, ob1[ii], ob2[ii], caoprms, &rap, &dap)
        ra_out[ii] = rap
        dec_out[ii] = dap

    return (ra_out, dec_out)

# Numeric lookup is only useful when scanning through the
# list of telescopes. The python interface does not need this.
# Instead obs() returns a dict (which may be a bit less efficient
# than an iterable but it's easy) with a dict inside. No arguments.
def obs():
     """
     telescopes = pal.obs()

     Obtain telescope parameters.

     Returns a dict with keys corresponding to the short
     name of the telescope. The corresponding value is itself
     a dict with keys:

     name: long name of the telescope
     long: west longitude (radians)
     lat: geodetic latitude (radians)
     height: Height above sea level (metres)
     """
     cdef int n = 1
     cdef char c
     cdef char cident[11]
     cdef char cname[41]
     cdef double w
     cdef double p
     cdef double h
     cdef int retval = 0

     result = {}

     while True:
          retval = cpal.palObs( n, &c, cident, 11, cname, 41,
                                &w, &p, &h )
          n=n+1 # Next telescope

          if retval != 0:
               break
          newtel = { 'name': cname.decode('UTF-8'),
                     'long': w,
                     'lat': p,
                     'height': h }
          result[cident.decode('UTF-8')] = newtel

     return result

def pa( double ha, double dec, double phi):
     """
     a = pa( ha, dec, phi )

     @palPa@
     """
     return cpal.palPa( ha, dec, phi )

def paVector(np.ndarray[double, ndim=1] ha not None,
             np.ndarray[double, ndim=1] dec not None,
             double phi):
     """
     a = paVector( ha, dec, phi )

     where ha, dec, and a are all numpy arrays

     @palPa@

     Raises
     ------
     ValueError if input arrays are not of the same length
     """
     cdef int length = len(ha)

     if len(dec)!=length:
         raise ValueError("You did not pass as many Decs as " \
                          + "Hour Angles to paVector")

     cdef np.ndarray paOut = np.zeros(length, dtype=np.float64)
     cdef int i
     for i in range(length):
          paOut[i]= cpal.palPa(ha[i], dec[i], phi)
     return paOut

def pcd( double disco, double x, double y):
    """
    (cx, cy) = pcd( disco, x, y )

    Returns the distorted coordinates and does not
    modify the supplied arguments.

    @palPcd@
    """
    cpal.palPcd( disco, &x, &y )
    return (x, y)

def pcdVector( double disco, np.ndarray[ double, ndim=1] x_in not None,
               np.ndarray[ double, ndim=1] y_in not None):
    """
    (cx, cy) = pcd( disco, x, y )

    where x and y are numpy arrays of x, y tangent plane coordinates
    and cx, cy are numpy arrays of tangent plane coordinates with the
    pincushion distortion applied.

    @palPcd@

    Raises
    ------
    ValueError if input arrays have inconsistent lengths
    """

    cdef int length = len(x_in)
    if len(y_in) != length:
        raise ValueError("You did not pass the same number of x as y " \
                         + "coordinates to pcdVector")

    cdef np.ndarray x_out = np.zeros(length, dtype=np.float64)
    cdef np.ndarray y_out = np.zeros(length, dtype=np.float64)
    cdef double xx
    cdef double yy

    for i in range(length):
        xx = x_in[i]
        yy = y_in[i]
        cpal.palPcd( disco, &xx, &yy)
        x_out[i] = xx
        y_out[i] = yy

    return (x_out, y_out)

def polmo( double elongm, double phim, double xp, double yp):
    """
    (elong, phi, daz) = polmo( elongm, phim, xp, yp )

    @palPolmo@
    """
    cdef double elong
    cdef double phi
    cdef double daz

    cpal.palPolmo( elongm, phim, xp, yp, &elong, &phi, &daz )

    return (elong, phi, daz)

def pertel(int jform, double date0, double date1,
            double epoch0, double orbi0, double anode0,
            double perih0, double aorq0, double e0, double am0):
    """
    (epoch1, orbi1, anode1, perih1, aorq1, e1, am1) = pertel( jform, date0, date1,
                                                            epoch0, orbi0, anode0, perih0,
                                                            aorq0, e0, am0 )

    @palPertel@

    Raises
    ------
    ValueError: Illegal input arguments.
    ArithmeticError: Internal error.
    """
    cdef double epoch1
    cdef double orbi1
    cdef double anode1
    cdef double perih1
    cdef double aorq1
    cdef double e1
    cdef double am1
    cdef int jstat

    cpal.palPertel( jform, date0, date1,
                    epoch0, orbi0, anode0, perih0, aorq0, e0, am0,
                    &epoch1, &orbi1, &anode1, &perih1, &aorq1, &e1, &am1,
                    &jstat )
    if jstat == -1:
        raise ValueError( "Illegal jform" )
    elif jstat == -2:
        raise ValueError( "Illegal e0" )
    elif jstat == -3:
        raise ValueError( "Illegal aorq0" )
    elif jstat == -4:
        raise ArithmeticError( "Internal error" )
    elif jstat == -5:
        raise ArithmeticError( "Numerical error" )

    return ( epoch1, orbi1, anode1, perih1, aorq1, e1, am1 )

def pertue( double date, np.ndarray[double, ndim=1] u not None ):
    """
    u2 = pertue( date, u )

    @palPertue@

    Raises
    ------
    ArithmeticError: Numerical error
    """
    cdef double cu[13]
    cdef int jstat

    for i in range(13):
        cu[i] = u[i]

    cpal.palPertue( date, cu, &jstat )

    if jstat == -1:
        raise ArithmeticError( "Numerical error" )

    # We return the modified U and do not change in place
    cdef np.ndarray u2 = np.zeros( [13], dtype=np.float64)
    for i in range(13):
        u2[i] = cu[i]
    return u2

def planel( double date, int jform, double epoch, double orbinc,
            double anode, double perih, double aorq,  double e,
            double aorl, double dm ):
    """
    pv = planel( date, jform, epoch, orbinc, anode, perih, aorq, e, aorl, dm )

    @palPlanel@

    Raises
    ------
    ValueError: Illegal input arguments
    ArithmeticError: Numerical error
    """
    cdef int jstat
    cdef double cpv[6]
    cpal.palPlanel( date, jform, epoch, orbinc, anode, perih, aorq,
                    e, aorl, dm, cpv, &jstat )
    if jstat == -1:
        raise ValueError( "Illegal jform" )
    elif jstat == -2:
        raise ValueError( "Illegal e" )
    elif jstat == -3:
        raise ValueError( "Illegal aorq" )
    elif jstat == -4:
        raise ValueError( "Illegal dm" )
    elif jstat == -5:
        raise ArithmeticError( "Numerical error" )

    cdef np.ndarray pv = np.zeros( [6], dtype=np.float64 )
    for i in range(6):
        pv[i] = cpv[i]
    return pv

def planet( double date, int planetnum ):
    """
    pv = planet( date, planetnum )

    @palPlanet@

    Raises
    ------
    ValueError: Illegal planet number
    ArithmeticError: Solution did not converge
    """
    cdef int jstat
    cdef double cpv[6]
    cpal.palPlanet( date, planetnum, cpv, &jstat )
    if jstat == -2:
        raise ArithmeticError( "Solution didn't converge" )
    elif jstat == -1:
        raise ValueError( "Illegal planet number "+str(planetnum)+", must be in range (1-8)" )

    cdef np.ndarray pv = np.zeros( [6], dtype=np.float64 )
    for i in range(6):
        pv[i] = cpv[i]
    return pv

def plante( double date, double elong, double phi, int jform,
            double epoch, double orbinc, double anode, double perih,
            double aorq, double e, double aorl, double dm ):
    """
    (ra, dec, r) = plante( date, elong, phi, jform, epoch, orbinc, anode,
                           perih, aorq, e, aorl, dm )

    @palPlante@

    Raises
    ------
    ValueError: Illegal arguments
    ArithmeticError: Numerical error
    """
    cdef double ra
    cdef double dec
    cdef double r
    cdef int jstat

    cpal.palPlante(date, elong, phi, jform, epoch, orbinc, anode, perih, aorq,e, aorl, dm, &ra, &dec, &r, &jstat)
    if jstat == -1:
        raise ValueError( "Illegal jform" )
    elif jstat == -2:
        raise ValueError( "Illegal e" )
    elif jstat == -3:
        raise ValueError( "Illegal aorq" )
    elif jstat == -4:
        raise ValueError( "Illegal dm" )
    elif jstat == -5:
        raise ArithmeticError( "Numerical error" )
    return (ra, dec, r)

def plantu( double date, double elong, double phi, np.ndarray[double, ndim=1] u not None):
    """
    (ra, dec, r) = plantu( date, elong, phi, u )

    @palPlantu@

    Raises
    ------
    ValueError: Radius vector zero
    Arithmetic Error: Failed to converge
    """
    cdef double ra
    cdef double dec
    cdef double r
    cdef int jstat
    cdef double cu[13]

    for i in range(13):
        cu[i] = u[i]
    cpal.palPlantu( date, elong, phi, cu, &ra, &dec, &r, &jstat )
    if jstat == -1:
        raise ValueError( "Radius vector zero" )
    elif jstat == -2:
        raise ArithmeticError( "Failed to converge" )
    return (ra, dec, r )

def pm( double r0, double d0, double pr, double pd,
        double px, double rv, double ep0, double ep1 ):
     """
     (r1, d1) = pm( r0, d0, pr, pd, px, rv, ep0, ep1 )

     @palPm@
     """
     cdef double r1
     cdef double d1
     cpal.palPm( r0, d0, pr, pd, px, rv, ep0, ep1, &r1, &d1 )
     return (r1, d1)

def pmVector( np.ndarray[double, ndim=1] r0 not None,
              np.ndarray[double, ndim=1] d0 not None,
              np.ndarray[double, ndim=1] pr not None,
              np.ndarray[double, ndim=1] pd not None,
              np.ndarray[double, ndim=1] px not None,
              np.ndarray[double, ndim=1] rv not None,
              double ep0, double ep1):
     """
     (r1, d1) = pmVector( r0, d0, pr, pd, px, rv, ep0, ep1 )

     where r0, d0, pr, pd, px, rv, r1, and d1 are all
     numpy arrays

     @palPm@

     Raises
     ------
     ValueError if input arrays are not of the same length
     """

     cdef int i
     cdef int length=len(r0)

     if len(d0)!=length or len(pr)!=length or len(pd)!=length \
     or len(px)!=length or len(rv)!=length:

         raise ValueError("The arrays you passed into pmVector " \
                          + "are not all of the same length")

     cdef np.ndarray r1out=np.zeros(length,dtype=np.float64)
     cdef np.ndarray d1out=np.zeros(length,dtype=np.float64)
     cdef double r1
     cdef double d1

     for i in range(length):
          cpal.palPm(r0[i],d0[i],pr[i],pd[i],px[i],rv[i],ep0,ep1,&r1,&d1)
          r1out[i]=r1
          d1out[i]=d1

     return(r1out,d1out)

def prebn( double bep0, double bep1 ):
     """
     rmatp = prebn( bep0, bep1 )

     @palPrebn@
     """
     cdef double crmatp[3][3]
     cpal.palPrebn( bep0, bep1, crmatp )
     cdef np.ndarray rmatp = np.zeros( [3,3], dtype=np.float64 )
     for i in range(3):
          for j in range(3):
               rmatp[i,j] = crmatp[i][j]
     return rmatp

def prec( double ep0, double ep1 ):
     """
     rmatp = prec( ep0, ep1 )

     @palPrec@
     """
     cdef double crmatp[3][3]
     cpal.palPrec( ep0, ep1, crmatp )
     cdef np.ndarray rmatp = np.zeros( [3,3], dtype=np.float64 )
     for i in range(3):
          for j in range(3):
               rmatp[i,j] = crmatp[i][j]
     return rmatp

def preces( sys, double ep0, double ep1, double ra, double dc ):
     """
     (ra, dc) = preces( sys, ep0, ep1, ra, dec)

     @palPreces@
     """
     byte_string = sys.encode('ascii')
     cdef char * csys = byte_string
     cpal.palPreces( csys, ep0, ep1, &ra, &dc )
     return (ra, dc)

def prenut( double epoch, double date):
    """
    rmatpn = prenut( epoch, date )

    @palPrenut@
    """
    cdef double crmatpn[3][3]
    cpal.palPrenut( epoch, date, crmatpn )
    cdef np.ndarray rmatpn = np.zeros( [3,3], dtype=np.float64 )
    for i in range(3):
        for j in range(3):
            rmatpn[i,j]=crmatpn[i][j]
    return rmatpn

def pv2el(np.ndarray[double, ndim=1] pv not None, double date, double pmass, int jformr):
    """
    (jform, epoch, orbinc, anode, perih, aorq, e, aorl, dm) = pv2el( pv, date, pmass, jformr )

    @palPv2el@

    Raises
    ------
    ValueError: Illegal arguments
    """
    cdef int jform
    cdef double epoch
    cdef double orbinc
    cdef double anode
    cdef double perih
    cdef double aorq
    cdef double e
    cdef double aorl
    cdef double dm
    cdef int jstat
    cdef double cpv[6]

    for i in range(6):
        cpv[i] = pv[i]
    cpal.palPv2el( cpv, date, pmass, jformr, &jform, &epoch, &orbinc, &anode, &perih,
                   &aorq, &e, &aorl, &dm, &jstat )
    if jstat == -1:
        raise ValueError( "Illegal PMASS" )
    elif jstat == -2:
        raise ValueError( "Illegal JFORMR" )
    elif jstat == -3:
        raise ValueError( "Position/velocity out of range" )
    return (jform, epoch, orbinc, anode, perih, aorq, e, aorl, dm)

def pv2ue( np.ndarray[double, ndim=1] pv not None, double date, double pmass ):
    """
    u = pv2ue( pv, date, pmass )

    @palPv2ue@

    Raises
    ------
    ValueError: Illegal arguments
    """
    cdef int jstat
    cdef double cu[13]
    cdef double cpv[6]

    for i in range(6):
        cpv[i] = pv[i]
    cpal.palPv2ue( cpv, date, pmass, cu, &jstat )
    if jstat == -1:
        raise ValueError( "Illegal PMASS" )
    elif jstat == -2:
        raise ValueError( "Too close to Sun" )
    elif jstat == -3:
        raise ValueError( "Too slow" )

    cdef np.ndarray u = np.zeros( [13], dtype=np.float64 )
    for i in range(13):
        u[i] = cu[i]
    return u

def pvobs( double p, double h, double stl ):
     """
     pv = pvobs( p, h, stl )

     @palPvobs@
     """
     cdef double cpv[6]
     cpal.palPvobs( p, h, stl, cpv )
     cdef np.ndarray pv = np.zeros( [6], dtype=np.float64 )
     for i in range(6):
          pv[i] = cpv[i]
     return pv

def rdplan( double date, int np, double elong, double phi ):
    """
    (ra, dec, diam) = rdplan( date, np, elong, phi )

    @palRdplan@
    """
    cdef double ra
    cdef double dec
    cdef double diam
    cpal.palRdplan( date, np, elong, phi, &ra, &dec, &diam )
    return (ra, dec, diam)

def refco( double hm, double tdk, double pmb, double rh, double wl, double phi, double tlr, double eps):
    """
    (refa, refb) = refco( hm, tdk, pmb, rh, wl, phi, tlr, eps )

    @palRefco@
    """
    cdef double refa
    cdef double refb
    cpal.palRefco(hm,tdk,pmb,rh,wl,phi,tlr,eps,&refa,&refb)
    return (refa,refb)

def refcoq( double tdk, double pmb, double rh, double wl):
    """
    (refa, refb) = refcoq( tdk, pmb, rh, wl )

    @palRefcoq@
    """
    cdef double refa
    cdef double refb
    cpal.palRefcoq( tdk, pmb, rh, wl, &refa, &refb )
    return (refa, refb)

def refro( double zobs, double hm, double tdk, double pmb,
           double rh, double wl, double phi, double tlr, double eps):
    """
    ref = refro( zobs, hm, tdk, pmb, rh, wl, phi, tlr, eps )

    @palRefro@
    """
    cdef double ref
    cpal.palRefro(zobs, hm, tdk, pmb, rh, wl, phi, tlr, eps, &ref)
    return ref

def refroVector( np.ndarray[ double, ndim=1] zobs not None, double hm, double tdk,
                double pmb, double rh, double wl, double phi, double tlr,
                double eps):
    """
    ref = refroVector( zobs, hm, tdk, pmb, rh, wl, phi, tlr, eps )

    where zobs is a numpy array of zenith distances in radians.  ref will be
    a numpy array of refractions (zenith distance minus observed zenith distance
    in radians)

    @palRefro@
    """
    cdef int length = len(zobs)
    cdef double ref
    cdef np.ndarray refVector = np.zeros([length], dtype=np.float64 )
    for i in range(length):
        cpal.palRefro(zobs[i], hm, tdk, pmb, rh, wl, phi, tlr, eps, &ref)
        refVector[i] = ref
    return refVector

def refv( np.ndarray[ double, ndim=1] vu not None, double refa, double refb):
    """
    vr = refv( vu, refa, refb )

    @palRefv@
    """
    cdef double cvr[3]
    cdef double cvu[3]
    cdef np.ndarray vr = np.zeros( [3], dtype=np.float64 )
    for i in range(3):
        cvu[i] = vu[i]
    cpal.palRefv( cvu, refa, refb, cvr )
    for i in range(3):
        vr[i] = cvr[i]
    return vr

def refvVector(np.ndarray[double, ndim=2] vu not None, double refa, double refb):
    """
    vr = refvVector( vu, refa, refb )

    where vu is a 2-D numpy array in which each row is a different 3-D cartesian
    vector.  vr will be a 2-D numpy array in which each row is the
    refraction-corrected vector corresponding to the same row in vu

    @palRefv@
    """
    cdef double cvr[3]
    cdef double cvu[3]
    cdef np.ndarray vr = np.zeros( [vu.shape[0], vu.shape[1]], dtype=np.float64 )

    for row in range(vu.shape[0]):
        for i in range(3):
            cvu[i] = vu[row][i]
        cpal.palRefv( cvu, refa, refb, cvr )
        for i in range(3):
            vr[row][i] = cvr[i]
    return vr

def refz( double zu, double refa, double refb ):
    """
    zr = refz( zu, refa, refb )

    @palRefz@
    """
    cdef double zr
    cpal.palRefz(zu,refa,refb,&zr)
    return zr

def refzVector(np.ndarray[double, ndim=1] zu not None,
               double refa, double refb):
     """
     zr = refzVector( zu, refa, refb )

     where zu and zr are numpy arrays

     @palRefz@
     """
     cdef int length = len(zu)
     cdef np.ndarray zrout = np.zeros(length, dtype=np.float64)
     cdef double zr
     cdef int i
     for i in range(length):
          cpal.palRefz(zu[i], refa, refb, &zr)
          zrout[i] = zr
     return zrout

def rverot( double phi, double ra, double da, double st ):
     """
     r = rverot( phi, ra, da, st )

     @palRverot@
     """
     return cpal.palRverot( phi, ra, da, st )

def rvgalc( double r2000, double d2000 ):
     """
     rv = rvgalc( r2000, d2000 )

     @palRvgalc@
     """
     return cpal.palRvgalc( r2000, d2000 )

def rvlg( double r2000, double d2000 ):
     """
     rv = rvlg( r2000, d2000 )

     @palRvlg@
     """
     return cpal.palRvlg( r2000, d2000 )

def rvlsrd( double r2000, double d2000 ):
     """
     rv = rvlsrd( r2000, d2000 )

     @palRvlsrd@
     """
     return cpal.palRvlsrd( r2000, d2000 )

def rvlsrk( double r2000, double d2000 ):
     """
     rv = rvlsrk( r2000, d2000 )

     @palRvlsrk@
     """
     return cpal.palRvlsrk( r2000, d2000 )

def subet( double rc, double dc, double eq ):
     """
     (rm, dm) = subet( rc, dc, eq )

     @palSubet@
     """
     cdef double rm
     cdef double dm
     cpal.palSubet( rc, dc, eq, &rm, &dm )
     return ( rm, dm )

def supgal( double dsl, double dsb ):
     """
     (dl, db) = supgal( dsl, dsb )

     @palSupgal@
     """
     cdef double dl
     cdef double db
     cpal.palSupgal( dsl, dsb, &dl, &db )
     return (dl, db)

def supgalVector(np.ndarray[double, ndim=1] dsl not None,
                 np.ndarray[double, ndim=1] dsb not None):
    """
    (dl, db) = supgalVector( dsl, dsb )

    where dsl, dsb, dl, and db are all numpy arrays

    @palSupgal@

    Raises
    ------
    ValueError if input arrays are of different lengths
    """

    cdef int length=len(dsl)

    if len(dsb)!=length:
        raise ValueError("You did not pass the same number of " \
                         + "longitudes as latitudes into supgalVector")

    cdef np.ndarray dl_out = np.empty(length, dtype=np.float64)
    cdef np.ndarray db_out = np.empty(length, dtype=np.float64)

    cdef double dl
    cdef double db

    for ii in range(length):
        cpal.palSupgal(dsl[ii], dsb[ii], &dl, &db)
        dl_out[ii] = dl
        db_out[ii] = db

    return (dl_out, db_out)


def ue2el(np.ndarray[double, ndim=1] u not None, int jformr ):
    """
    (jform, epoch, orbinc, anode, perih, aorq, e, aorl, dm) = ue2el( u, jformr )

    @palUe2el@

    Raises
    ------
    ValueError: Illegal arguments
    """
    cdef int jform
    cdef double epoch
    cdef double orbinc
    cdef double anode
    cdef double perih
    cdef double aorq
    cdef double e
    cdef double aorl
    cdef double dm
    cdef int jstat
    cdef double cu[13]

    for i in range(13):
        cu[i] = u[i]
    cpal.palUe2el( cu, jformr, &jform, &epoch, &orbinc, &anode, &perih,
                   &aorq, &e, &aorl, &dm, &jstat )
    if jstat == -1:
        raise ValueError( "Illegal combined mass" )
    elif jstat == -2:
        raise ValueError( "Illegal jformr" )
    elif jstat == -3:
        raise ValueError( "Position/velocity out of range")
    return (jform, epoch, orbinc, anode, perih, aorq, e, aorl, dm)

#  Note that u is updated and returned
def ue2pv( double date, np.ndarray[double, ndim=1] u not None ):
    """
    (u2, pv) = ue2pv( date, u )

    @palUe2pv@

    Raises
    ------
    ValueError: Radius vector zero
    ArithmeticError: Failed to converge
    """
    cdef double cu[13]
    cdef double cpv[6]
    cdef int jstat

    for i in range(13):
        cu[i] = u[i]
    cpal.palUe2pv( date, cu, cpv, &jstat )
    if jstat == -1:
        raise ValueError( "Radius vector zero" )
    elif jstat == -2:
        raise ArithmeticError( "Failed to converge" )

    # We need to return a completely new updated U
    # rather than overwrite in place
    cdef np.ndarray u2 = np.zeros( [13], dtype=np.float64)
    cdef np.ndarray pv = np.zeros( [6], dtype=np.float64)
    for i in range(13):
        u2[i] = cu[i]
    for i in range(6):
        pv[i] = cpv[i]
    return (u2, pv)

def unpcd( double disco, double x, double y):
    """
    (tx, ty) = unpcd( disco, x, y )

    Returns the tangent plane coordinates and does not
    modify the supplied arguments.

    @palUnpcd@
    """
    cpal.palUnpcd( disco, &x, &y )
    return (x, y)

def unpcdVector( double disco, np.ndarray[ double, ndim=1 ] x_in not None,
                 np.ndarray[ double, ndim=1 ] y_in not None):
    """
    (tx, ty) = unpcdVector( disco, x, y)

    where x and y are numpy arrays of x, y tangent plane coordinates and
    tx, ty are numpy arrays of tangent plane coordinates with the
    pincushion distortion removed

    @palUnpcd@

    Raises
    ------
    ValueError if input arrays have inconsistent lengths
    """

    cdef int length=len(x_in)
    if len(y_in) != length:
        raise ValueError("You did not pass the same number of x as y " \
                         + "coordinates to unpcdVector")

    cdef np.ndarray x_out = np.zeros(length, dtype=np.float64)
    cdef np.ndarray y_out = np.zeros(length, dtype=np.float64)
    cdef double xx
    cdef double yy

    for i in range(length):
        xx = x_in[i]
        yy = y_in[i]
        cpal.palUnpcd(disco, &xx, &yy)
        x_out[i] = xx
        y_out[i] =yy

    return (x_out,y_out)