text:
"""
    eform(n::Ellipsoid)

Earth reference ellipsoids.

### Given ###

- `n`: Ellipsoid identifier (Note 1)

### Returned ###

- `a`: Equatorial radius (meters, Note 2)
- `f`: Flattening (Note 2)

### Notes ###

1. The identifier n is a number that specifies the choice of
   reference ellipsoid.  The following are supported:

        - `WGS84`
        - `GRS80`
        - `WGS72`

2. The ellipsoid parameters are returned in the form of equatorial
   radius in meters (a) and flattening (f).  The latter is a number
   around 0.00335, i.e. around 1/298.

3. For the case where an unsupported n value is supplied, zero a and
   f are returned, as well as error status.

### References ###

- Department of Defense World Geodetic System 1984, National
    Imagery and Mapping Agency Technical Report 8350.2, Third
    Edition, p3-2.

- Moritz, H., Bull. Geodesique 66-2, 187 (1992).

- The Department of Defense World Geodetic System 1972, World
    Geodetic System Committee, May 1974.

- Explanatory Supplement to the Astronomical Almanac,
    P. Kenneth Seidelmann (ed), University Science Books (1992),
    p220.

"""
function eform(n::Ellipsoid)
    a = Ref(0.0)
    f = Ref(0.0)
    i = ccall((:eraEform, liberfa), Cint,
              (Cint, Ref{Cdouble}, Ref{Cdouble}),
              n, a, f)
    if i == -1
        throw(ERFAException("illegal identifier"))
    end
    a[], f[]
end

"""
    eors(rnpb, s)

Equation of the origins, given the classical NPB matrix and the
quantity s.

### Given ###

- `rnpb`: Classical nutation x precession x bias matrix
- `s`: The quantity s (the CIO locator)

### Returned ###

- The equation of the origins in radians.

### Notes ###

1.  The equation of the origins is the distance between the true
    equinox and the celestial intermediate origin and, equivalently,
    the difference between Earth rotation angle and Greenwich
    apparent sidereal time (ERA-GST).  It comprises the precession
    (since J2000.0) in right ascension plus the equation of the
    equinoxes (including the small correction terms).

2.  The algorithm is from Wallace & Capitaine (2006).

### References ###

- Capitaine, N. & Wallace, P.T., 2006, Astron.Astrophys. 450, 855

- Wallace, P. & Capitaine, N., 2006, Astron.Astrophys. 459, 981

"""
function eors(rnpb, s)
    ccall((:eraEors, liberfa), Cdouble,
          (Ptr{Cdouble}, Cdouble),
          rnpb, s)
end

"""
    epv00(date1, date2)

Earth position and velocity, heliocentric and barycentric, with
respect to the Barycentric Celestial Reference System.

### Given ###

- `date1`, `date2`: TDB date (Note 1)

### Returned ###

- `pvh`: Heliocentric Earth position/velocity
- `pvb`: Barycentric Earth position/velocity

### Notes ###

1. The TDB date date1+date2 is a Julian Date, apportioned in any
   convenient way between the two arguments.  For example,
   JD(TDB)=2450123.7 could be expressed in any of these ways, among
   others:

   | `date1`   |     `date2` |                    |
   |:----------|:------------|:-------------------|
   | 2450123.7 |         0.0 | JD method          |
   | 2451545.0 |     -1421.3 | J2000 method       |
   | 2400000.5 |     50123.2 | MJD method         |
   | 2450123.5 |         0.2 | date & time method |

   The JD method is the most natural and convenient to use in cases
   where the loss of several decimal digits of resolution is
   acceptable.  The J2000 method is best matched to the way the
   argument is handled internally and will deliver the optimum
   resolution.  The MJD method and the date & time methods are both
   good compromises between resolution and convenience.  However,
   the accuracy of the result is more likely to be limited by the
   algorithm itself than the way the date has been expressed.

   n.b. TT can be used instead of TDB in most applications.

2. On return, the arrays pvh and pvb contain the following:

      pvh[0][0]  x       }
      pvh[0][1]  y       } heliocentric position, au
      pvh[0][2]  z       }

      pvh[1][0]  xdot    }
      pvh[1][1]  ydot    } heliocentric velocity, au/d
      pvh[1][2]  zdot    }

      pvb[0][0]  x       }
      pvb[0][1]  y       } barycentric position, au
      pvb[0][2]  z       }

      pvb[1][0]  xdot    }
      pvb[1][1]  ydot    } barycentric velocity, au/d
      pvb[1][2]  zdot    }

   The vectors are with respect to the Barycentric Celestial
   Reference System.  The time unit is one day in TDB.

3. The function is a SIMPLIFIED SOLUTION from the planetary theory
   VSOP2000 (X. Moisson, P. Bretagnon, 2001, Celes. Mechanics &
   Dyn. Astron., 80, 3/4, 205-213) and is an adaptation of original
   Fortran code supplied by P. Bretagnon (private comm., 2000).

4. Comparisons over the time span 1900-2100 with this simplified
   solution and the JPL DE405 ephemeris give the following results:

                              RMS    max
         Heliocentric:
            position error    3.7   11.2   km
            velocity error    1.4    5.0   mm/s

         Barycentric:
            position error    4.6   13.4   km
            velocity error    1.4    4.9   mm/s

   Comparisons with the JPL DE406 ephemeris show that by 1800 and
   2200 the position errors are approximately double their 1900-2100
   size.  By 1500 and 2500 the deterioration is a factor of 10 and
   by 1000 and 3000 a factor of 60.  The velocity accuracy falls off
   at about half that rate.

5. It is permissible to use the same array for pvh and pvb, which
   will receive the barycentric values.

"""
function epv00(date1, date2)
    pvh = zeros((2, 3))
    pvb = zeros((2, 3))
    i = ccall((:eraEpv00, liberfa),
              Cint,
              (Cdouble, Cdouble, Ptr{Cdouble}, Ptr{Cdouble}),
              date1, date2, pvh, pvb)
    if i == 1
        @warn "date outside the range 1900-2100 AD"
    end
    pvh, pvb
end

"""
    eceq06(date1, date2, dl, db)

Transformation from ecliptic coordinates (mean equinox and ecliptic
of date) to ICRS RA,Dec, using the IAU 2006 precession model.

### Given ###

- `date1`, `date2`: TT as a 2-part Julian date (Note 1)
- `dl`, `db`: Ecliptic longitude and latitude (radians)

### Returned ###

- `dr`, `dd`: ICRS right ascension and declination (radians)

1. The TT date date1+date2 is a Julian Date, apportioned in any
   convenient way between the two arguments.  For example,
   JD(TT)=2450123.7 could be expressed in any of these ways,
   among others:

   | `date1`   |     `date2` |                    |
   |:----------|:------------|:-------------------|
   | 2450123.7 |         0.0 | JD method          |
   | 2451545.0 |     -1421.3 | J2000 method       |
   | 2400000.5 |     50123.2 | MJD method         |
   | 2450123.5 |         0.2 | date & time method |

   The JD method is the most natural and convenient to use in
   cases where the loss of several decimal digits of resolution
   is acceptable.  The J2000 method is best matched to the way
   the argument is handled internally and will deliver the
   optimum resolution.  The MJD method and the date & time methods
   are both good compromises between resolution and convenience.

2. No assumptions are made about whether the coordinates represent
   starlight and embody astrometric effects such as parallax or
   aberration.

3. The transformation is approximately that from ecliptic longitude
   and latitude (mean equinox and ecliptic of date) to mean J2000.0
   right ascension and declination, with only frame bias (always
   less than 25 mas) to disturb this classical picture.

### Called ###

- `eraS2c`: spherical coordinates to unit vector
- `eraEcm06`: J2000.0 to ecliptic rotation matrix, IAU 2006
- `eraTrxp`: product of transpose of r-matrix and p-vector
- `eraC2s`: unit vector to spherical coordinates
- `eraAnp`: normalize angle into range 0 to 2pi
- `eraAnpm`: normalize angle into range +/- pi

"""
eceq06

"""
    eqec06(date1, date2, dr, dd)

Transformation from ICRS equatorial coordinates to ecliptic
coordinates (mean equinox and ecliptic of date) using IAU 2006
precession model.

### Given ###

- `date1`, `date2`: TT as a 2-part Julian date (Note 1)
- `dr`, `dd`: ICRS right ascension and declination (radians)

### Returned ###

- `dl`, `db`: Ecliptic longitude and latitude (radians)

1. The TT date date1+date2 is a Julian Date, apportioned in any
   convenient way between the two arguments.  For example,
   JD(TT)=2450123.7 could be expressed in any of these ways,
   among others:

   | `date1`   |     `date2` |                    |
   |:----------|:------------|:-------------------|
   | 2450123.7 |         0.0 | JD method          |
   | 2451545.0 |     -1421.3 | J2000 method       |
   | 2400000.5 |     50123.2 | MJD method         |
   | 2450123.5 |         0.2 | date & time method |

   The JD method is the most natural and convenient to use in
   cases where the loss of several decimal digits of resolution
   is acceptable.  The J2000 method is best matched to the way
   the argument is handled internally and will deliver the
   optimum resolution.  The MJD method and the date & time methods
   are both good compromises between resolution and convenience.

2. No assumptions are made about whether the coordinates represent
   starlight and embody astrometric effects such as parallax or
   aberration.

3. The transformation is approximately that from mean J2000.0 right
   ascension and declination to ecliptic longitude and latitude
   (mean equinox and ecliptic of date), with only frame bias (always
   less than 25 mas) to disturb this classical picture.

### Called ###

- `eraS2c`: spherical coordinates to unit vector
- `eraEcm06`: J2000.0 to ecliptic rotation matrix, IAU 2006
- `eraRxp`: product of r-matrix and p-vector
- `eraC2s`: unit vector to spherical coordinates
- `eraAnp`: normalize angle into range 0 to 2pi
- `eraAnpm`: normalize angle into range +/- pi

"""
eqec06

for name in ("eceq06",
             "eqec06")
    f = Symbol(name)
    fc = "era" * uppercasefirst(name)
    @eval begin
        function ($f)(date1, date2, d1, d2)
            r1 = [0.0]
            r2 = [0.0]
            ccall(($fc, liberfa), Cvoid,
                  (Cdouble, Cdouble, Cdouble, Cdouble, Ptr{Cdouble}, Ptr{Cdouble}),
                  date1, date2, d1, d2, r1, r2)
            r1[1], r2[1]
        end
    end
end

"""
    epb2jd(epj)

Besselian Epoch to Julian Date.

### Given ###

- `epb`: Besselian Epoch (e.g. 1957.3)

### Returned ###

- `djm0`: MJD zero-point: always 2400000.5
- `djm`: Modified Julian Date

### Note ###

   The Julian Date is returned in two pieces, in the usual ERFA
   manner, which is designed to preserve time resolution.  The
   Julian Date is available as a single number by adding djm0 and
   djm.

### Reference ###

- Lieske, J.H., 1979, Astron.Astrophys. 73, 282.

"""
epb2jd

"""
    epj2jd(epj)

Julian Epoch to Julian Date.

### Given ###

- `epj`: Julian Epoch (e.g. 1996.8)

### Returned ###

- `djm0`: MJD zero-point: always 2400000.5
- `djm`: Modified Julian Date

### Note ###

   The Julian Date is returned in two pieces, in the usual ERFA
   manner, which is designed to preserve time resolution.  The
   Julian Date is available as a single number by adding djm0 and
   djm.

### Reference ###

- Lieske, J.H., 1979, Astron.Astrophys. 73, 282.

"""
epj2jd

for name in ("epb2jd",
             "epj2jd")
    f = Symbol(name)
    fc = "era" * uppercasefirst(name)
    @eval begin
        function ($f)(d)
            r1 = Ref(0.0)
            r2 = Ref(0.0)
            ccall(($fc, liberfa), Cvoid,
                  (Cdouble, Ref{Cdouble}, Ref{Cdouble}),
                  d, r1, r2)
            r1[], r2[]
        end
    end
end

"""
    ee00a(dj1, dj2)

Equation of the equinoxes, compatible with IAU 2000 resolutions.

### Given ###

- `date1`, `date2`: TT as a 2-part Julian Date (Note 1)

### Returned ###

- Equation of the equinoxes (Note 2)

### Notes ###

1. The TT date date1+date2 is a Julian Date, apportioned in any
   convenient way between the two arguments.  For example,
   JD(TT)=2450123.7 could be expressed in any of these ways,
   among others:

   | `date1`   |     `date2` |                    |
   |:----------|:------------|:-------------------|
   | 2450123.7 |         0.0 | JD method          |
   | 2451545.0 |     -1421.3 | J2000 method       |
   | 2400000.5 |     50123.2 | MJD method         |
   | 2450123.5 |         0.2 | date & time method |

   The JD method is the most natural and convenient to use in
   cases where the loss of several decimal digits of resolution
   is acceptable.  The J2000 method is best matched to the way
   the argument is handled internally and will deliver the
   optimum resolution.  The MJD method and the date & time methods
   are both good compromises between resolution and convenience.

2. The result, which is in radians, operates in the following sense:

      Greenwich apparent ST = GMST + equation of the equinoxes

3. The result is compatible with the IAU 2000 resolutions.  For
   further details, see IERS Conventions 2003 and Capitaine et al.
   (2002).

### Called ###

- `eraPr00`: IAU 2000 precession adjustments
- `eraObl80`: mean obliquity, IAU 1980
- `eraNut00a`: nutation, IAU 2000A
- `eraEe00`: equation of the equinoxes, IAU 2000

### References ###

- Capitaine, N., Wallace, P.T. and McCarthy, D.D., "Expressions to
    implement the IAU 2000 definition of UT1", Astronomy &
    Astrophysics, 406, 1135-1149 (2003).

- McCarthy, D. D., Petit, G. (eds.), IERS Conventions (2003),
    IERS Technical Note No. 32, BKG (2004).

"""
ee00a

"""
    ee00b(dj1, dj2)

Equation of the equinoxes, compatible with IAU 2000 resolutions but
using the truncated nutation model IAU 2000B.

### Given ###

- `date1`, `date2`: TT as a 2-part Julian Date (Note 1)

### Returned ###

- Equation of the equinoxes (Note 2)

### Notes ###

1. The TT date date1+date2 is a Julian Date, apportioned in any
   convenient way between the two arguments.  For example,
   JD(TT)=2450123.7 could be expressed in any of these ways,
   among others:

   | `date1`   |     `date2` |                    |
   |:----------|:------------|:-------------------|
   | 2450123.7 |         0.0 | JD method          |
   | 2451545.0 |     -1421.3 | J2000 method       |
   | 2400000.5 |     50123.2 | MJD method         |
   | 2450123.5 |         0.2 | date & time method |

   The JD method is the most natural and convenient to use in
   cases where the loss of several decimal digits of resolution
   is acceptable.  The J2000 method is best matched to the way
   the argument is handled internally and will deliver the
   optimum resolution.  The MJD method and the date & time methods
   are both good compromises between resolution and convenience.

2. The result, which is in radians, operates in the following sense:

      Greenwich apparent ST = GMST + equation of the equinoxes

3. The result is compatible with the IAU 2000 resolutions except
   that accuracy has been compromised for the sake of speed.  For
   further details, see McCarthy & Luzum (2001), IERS Conventions
   2003 and Capitaine et al. (2003).

### Called ###

- `eraPr00`: IAU 2000 precession adjustments
- `eraObl80`: mean obliquity, IAU 1980
- `eraNut00b`: nutation, IAU 2000B
- `eraEe00`: equation of the equinoxes, IAU 2000

### References ###

- Capitaine, N., Wallace, P.T. and McCarthy, D.D., "Expressions to
    implement the IAU 2000 definition of UT1", Astronomy &
    Astrophysics, 406, 1135-1149 (2003)

- McCarthy, D.D. & Luzum, B.J., "An abridged model of the
    precession-nutation of the celestial pole", Celestial Mechanics &
    Dynamical Astronomy, 85, 37-49 (2003)

- McCarthy, D. D., Petit, G. (eds.), IERS Conventions (2003),
    IERS Technical Note No. 32, BKG (2004)

"""
ee00b

"""
    ee06a(dj1, dj2)

Equation of the equinoxes, compatible with IAU 2000 resolutions and
IAU 2006/2000A precession-nutation.

### Given ###

- `date1`, `date2`: TT as a 2-part Julian Date (Note 1)

### Returned ###

- Equation of the equinoxes (Note 2)

### Notes ###

1. The TT date date1+date2 is a Julian Date, apportioned in any
   convenient way between the two arguments.  For example,
   JD(TT)=2450123.7 could be expressed in any of these ways,
   among others:

   | `date1`   |     `date2` |                    |
   |:----------|:------------|:-------------------|
   | 2450123.7 |         0.0 | JD method          |
   | 2451545.0 |     -1421.3 | J2000 method       |
   | 2400000.5 |     50123.2 | MJD method         |
   | 2450123.5 |         0.2 | date & time method |

   The JD method is the most natural and convenient to use in
   cases where the loss of several decimal digits of resolution
   is acceptable.  The J2000 method is best matched to the way
   the argument is handled internally and will deliver the
   optimum resolution.  The MJD method and the date & time methods
   are both good compromises between resolution and convenience.

2. The result, which is in radians, operates in the following sense:

      Greenwich apparent ST = GMST + equation of the equinoxes

### Called ###

- `eraAnpm`: normalize angle into range +/- pi
- `eraGst06a`: Greenwich apparent sidereal time, IAU 2006/2000A
- `eraGmst06`: Greenwich mean sidereal time, IAU 2006

### Reference ###

- McCarthy, D. D., Petit, G. (eds.), 2004, IERS Conventions (2003),
    IERS Technical Note No. 32, BKG

"""
ee06a

"""
    eect00(date1, date2)

Equation of the equinoxes complementary terms, consistent with
IAU 2000 resolutions.

### Given ###

- `date1`, `date2`: TT as a 2-part Julian Date (Note 1)

### Returned ###

- Complementary terms (Note 2)

### Notes ###

1. The TT date date1+date2 is a Julian Date, apportioned in any
   convenient way between the two arguments.  For example,
   JD(TT)=2450123.7 could be expressed in any of these ways,
   among others:

   | `date1`   |     `date2` |                    |
   |:----------|:------------|:-------------------|
   | 2450123.7 |         0.0 | JD method          |
   | 2451545.0 |     -1421.3 | J2000 method       |
   | 2400000.5 |     50123.2 | MJD method         |
   | 2450123.5 |         0.2 | date & time method |

   The JD method is the most natural and convenient to use in
   cases where the loss of several decimal digits of resolution
   is acceptable.  The J2000 method is best matched to the way
   the argument is handled internally and will deliver the
   optimum resolution.  The MJD method and the date & time methods
   are both good compromises between resolution and convenience.

2. The "complementary terms" are part of the equation of the
   equinoxes (EE), classically the difference between apparent and
   mean Sidereal Time:

      GAST = GMST + EE

   with:

      EE = dpsi * cos(eps)

   where dpsi is the nutation in longitude and eps is the obliquity
   of date.  However, if the rotation of the Earth were constant in
   an inertial frame the classical formulation would lead to
   apparent irregularities in the UT1 timescale traceable to side-
   effects of precession-nutation.  In order to eliminate these
   effects from UT1, "complementary terms" were introduced in 1994
   (IAU, 1994) and took effect from 1997 (Capitaine and Gontier,
   1993):

      GAST = GMST + CT + EE

   By convention, the complementary terms are included as part of
   the equation of the equinoxes rather than as part of the mean
   Sidereal Time.  This slightly compromises the "geometrical"
   interpretation of mean sidereal time but is otherwise
   inconsequential.

   The present function computes CT in the above expression,
   compatible with IAU 2000 resolutions (Capitaine et al., 2002, and
   IERS Conventions 2003).

### Called ###

- `eraFal03`: mean anomaly of the Moon
- `eraFalp03`: mean anomaly of the Sun
- `eraFaf03`: mean argument of the latitude of the Moon
- `eraFad03`: mean elongation of the Moon from the Sun
- `eraFaom03`: mean longitude of the Moon's ascending node
- `eraFave03`: mean longitude of Venus
- `eraFae03`: mean longitude of Earth
- `eraFapa03`: general accumulated precession in longitude

### References ###

- Capitaine, N. & Gontier, A.-M., Astron. Astrophys., 275,
    645-650 (1993)

- Capitaine, N., Wallace, P.T. and McCarthy, D.D., "Expressions to
    implement the IAU 2000 definition of UT1", Astronomy &
    Astrophysics, 406, 1135-1149 (2003)

- IAU Resolution C7, Recommendation 3 (1994)

- McCarthy, D. D., Petit, G. (eds.), IERS Conventions (2003),
    IERS Technical Note No. 32, BKG (2004)

"""
eect00

"""
    eo06a(date1, date2)

Equation of the origins, IAU 2006 precession and IAU 2000A nutation.

### Given ###

- `date1`, `date2`: TT as a 2-part Julian Date (Note 1)

### Returned ###

- Equation of the origins in radians

### Notes ###

1. The TT date date1+date2 is a Julian Date, apportioned in any
   convenient way between the two arguments.  For example,
   JD(TT)=2450123.7 could be expressed in any of these ways,
   among others:

   | `date1`   |     `date2` |                    |
   |:----------|:------------|:-------------------|
   | 2450123.7 |         0.0 | JD method          |
   | 2451545.0 |     -1421.3 | J2000 method       |
   | 2400000.5 |     50123.2 | MJD method         |
   | 2450123.5 |         0.2 | date & time method |

   The JD method is the most natural and convenient to use in
   cases where the loss of several decimal digits of resolution
   is acceptable.  The J2000 method is best matched to the way
   the argument is handled internally and will deliver the
   optimum resolution.  The MJD method and the date & time methods
   are both good compromises between resolution and convenience.

2. The equation of the origins is the distance between the true
   equinox and the celestial intermediate origin and, equivalently,
   the difference between Earth rotation angle and Greenwich
   apparent sidereal time (ERA-GST).  It comprises the precession
   (since J2000.0) in right ascension plus the equation of the
   equinoxes (including the small correction terms).

### Called ###

- `eraPnm06a`: classical NPB matrix, IAU 2006/2000A
- `eraBpn2xy`: extract CIP X,Y coordinates from NPB matrix
- `eraS06`: the CIO locator s, given X,Y, IAU 2006
- `eraEors`: equation of the origins, given NPB matrix and s

### References ###

- Capitaine, N. & Wallace, P.T., 2006, Astron.Astrophys. 450, 855

- Wallace, P.T. & Capitaine, N., 2006, Astron.Astrophys. 459, 981

"""
eo06a

"""
    epb(dj1, dj2)

Julian Date to Besselian Epoch.

### Given ###

- `dj1`, `dj2`: Julian Date (see note)

### Returned ###

- Besselian Epoch.

### Note ###

   The Julian Date is supplied in two pieces, in the usual ERFA
   manner, which is designed to preserve time resolution.  The
   Julian Date is available as a single number by adding dj1 and
   dj2.  The maximum resolution is achieved if dj1 is 2451545.0
   (J2000.0).

### Reference ###

- Lieske, J.H., 1979. Astron.Astrophys., 73, 282.

"""
epb

"""
    epj(dj1, dj2)

Julian Date to Julian Epoch.

### Given ###

- `dj1`, `dj2`: Julian Date (see note)

### Returned ###

- Julian Epoch

### Note ###

   The Julian Date is supplied in two pieces, in the usual ERFA
   manner, which is designed to preserve time resolution.  The
   Julian Date is available as a single number by adding dj1 and
   dj2.  The maximum resolution is achieved if dj1 is 2451545.0
   (J2000.0).

### Reference ###

- Lieske, J.H., 1979, Astron.Astrophys. 73, 282.

"""
epj

"""
    eqeq94(date1, date2)

Equation of the equinoxes, IAU 1994 model.

### Given ###

- `date1`, `date2`: TDB date (Note 1)

### Returned ###

- Equation of the equinoxes (Note 2)

### Notes ###

1. The date date1+date2 is a Julian Date, apportioned in any
   convenient way between the two arguments.  For example,
   JD(TT)=2450123.7 could be expressed in any of these ways,
   among others:

   | `date1`   |     `date2` |                    |
   |:----------|:------------|:-------------------|
   | 2450123.7 |         0.0 | JD method          |
   | 2451545.0 |     -1421.3 | J2000 method       |
   | 2400000.5 |     50123.2 | MJD method         |
   | 2450123.5 |         0.2 | date & time method |

   The JD method is the most natural and convenient to use in
   cases where the loss of several decimal digits of resolution
   is acceptable.  The J2000 method is best matched to the way
   the argument is handled internally and will deliver the
   optimum resolution.  The MJD method and the date & time methods
   are both good compromises between resolution and convenience.

2. The result, which is in radians, operates in the following sense:

      Greenwich apparent ST = GMST + equation of the equinoxes

### Called ###

- `eraAnpm`: normalize angle into range +/- pi
- `eraNut80`: nutation, IAU 1980
- `eraObl80`: mean obliquity, IAU 1980

### References ###

- IAU Resolution C7, Recommendation 3 (1994).

- Capitaine, N. & Gontier, A.-M., 1993, Astron. Astrophys., 275,
    645-650.

"""
eqeq94

"""
    era00(dj1, dj2)

Earth rotation angle (IAU 2000 model).

### Given ###

- `dj1`, `dj2`: UT1 as a 2-part Julian Date (see note)

### Returned ###

- Earth rotation angle (radians), range 0-2pi

### Notes ###

1. The UT1 date dj1+dj2 is a Julian Date, apportioned in any
   convenient way between the arguments dj1 and dj2.  For example,
   JD(UT1)=2450123.7 could be expressed in any of these ways,
   among others:

           dj1            dj2

       2450123.7           0.0       (JD method)
       2451545.0       -1421.3       (J2000 method)
       2400000.5       50123.2       (MJD method)
       2450123.5           0.2       (date & time method)

   The JD method is the most natural and convenient to use in
   cases where the loss of several decimal digits of resolution
   is acceptable.  The J2000 and MJD methods are good compromises
   between resolution and convenience.  The date & time method is
   best matched to the algorithm used:  maximum precision is
   delivered when the dj1 argument is for 0hrs UT1 on the day in
   question and the dj2 argument lies in the range 0 to 1, or vice
   versa.

2. The algorithm is adapted from Expression 22 of Capitaine et al.
   2000.  The time argument has been expressed in days directly,
   and, to retain precision, integer contributions have been
   eliminated.  The same formulation is given in IERS Conventions
   (2003), Chap. 5, Eq. 14.

### Called ###

- `eraAnp`: normalize angle into range 0 to 2pi

### References ###

- Capitaine N., Guinot B. and McCarthy D.D, 2000, Astron.
    Astrophys., 355, 398-405.

- McCarthy, D. D., Petit, G. (eds.), IERS Conventions (2003),
    IERS Technical Note No. 32, BKG (2004)

"""
era00

for name in ("ee00a",
             "ee00b",
             "ee06a",
             "eect00",
             "eo06a",
             "epb",
             "epj",
             "eqeq94",
             "era00")
    f = Symbol(name)
    fc = "era" * uppercasefirst(name)
    @eval ($f)(d1, d2) = ccall(($fc, liberfa), Cdouble, (Cdouble, Cdouble), d1, d2)
end

"""
    ee00(date1, date2, epsa, dpsi)

The equation of the equinoxes, compatible with IAU 2000 resolutions,
given the nutation in longitude and the mean obliquity.

### Given ###

- `date1`, `date2`: TT as a 2-part Julian Date (Note 1)
- `epsa`: Mean obliquity (Note 2)
- `dpsi`: Nutation in longitude (Note 3)

### Returned ###

- Equation of the equinoxes (Note 4)

### Notes ###

1. The TT date date1+date2 is a Julian Date, apportioned in any
   convenient way between the two arguments.  For example,
   JD(TT)=2450123.7 could be expressed in any of these ways,
   among others:

   | `date1`   |     `date2` |                    |
   |:----------|:------------|:-------------------|
   | 2450123.7 |         0.0 | JD method          |
   | 2451545.0 |     -1421.3 | J2000 method       |
   | 2400000.5 |     50123.2 | MJD method         |
   | 2450123.5 |         0.2 | date & time method |

   The JD method is the most natural and convenient to use in
   cases where the loss of several decimal digits of resolution
   is acceptable.  The J2000 method is best matched to the way
   the argument is handled internally and will deliver the
   optimum resolution.  The MJD method and the date & time methods
   are both good compromises between resolution and convenience.

2. The obliquity, in radians, is mean of date.

3. The result, which is in radians, operates in the following sense:

      Greenwich apparent ST = GMST + equation of the equinoxes

4. The result is compatible with the IAU 2000 resolutions.  For
   further details, see IERS Conventions 2003 and Capitaine et al.
   (2002).

### Called ###

- `eraEect00`: equation of the equinoxes complementary terms

### References ###

- Capitaine, N., Wallace, P.T. and McCarthy, D.D., "Expressions to
    implement the IAU 2000 definition of UT1", Astronomy &
    Astrophysics, 406, 1135-1149 (2003)

- McCarthy, D. D., Petit, G. (eds.), IERS Conventions (2003),
    IERS Technical Note No. 32, BKG (2004)

"""
ee00(date1, date2, epsa, dpsi) = ccall((:eraEe00, liberfa), Cdouble, (Cdouble, Cdouble, Cdouble, Cdouble), date1, date2, epsa, dpsi)

"""
    ecm06(date1, date2)

ICRS equatorial to ecliptic rotation matrix, IAU 2006.

### Given ###

- `date1`, `date2`: TT as a 2-part Julian date (Note 1)

### Returned ###

- `rm`: ICRS to ecliptic rotation matrix

### Notes ###

1. The TT date date1+date2 is a Julian Date, apportioned in any
   convenient way between the two arguments.  For example,
   JD(TT)=2450123.7 could be expressed in any of these ways,
   among others:

   | `date1`   |     `date2` |                    |
   |:----------|:------------|:-------------------|
   | 2450123.7 |         0.0 | JD method          |
   | 2451545.0 |     -1421.3 | J2000 method       |
   | 2400000.5 |     50123.2 | MJD method         |
   | 2450123.5 |         0.2 | date & time method |

   The JD method is the most natural and convenient to use in
   cases where the loss of several decimal digits of resolution
   is acceptable.  The J2000 method is best matched to the way
   the argument is handled internally and will deliver the
   optimum resolution.  The MJD method and the date & time methods
   are both good compromises between resolution and convenience.

1. The matrix is in the sense

      E_ep = rm x P_ICRS,

   where P_ICRS is a vector with respect to ICRS right ascension
   and declination axes and E_ep is the same vector with respect to
   the (inertial) ecliptic and equinox of date.

2. P_ICRS is a free vector, merely a direction, typically of unit
   magnitude, and not bound to any particular spatial origin, such
   as the Earth, Sun or SSB.  No assumptions are made about whether
   it represents starlight and embodies astrometric effects such as
   parallax or aberration.  The transformation is approximately that
   between mean J2000.0 right ascension and declination and ecliptic
   longitude and latitude, with only frame bias (always less than
   25 mas) to disturb this classical picture.

### Called ###

- `eraObl06`: mean obliquity, IAU 2006
- `eraPmat06`: PB matrix, IAU 2006
- `eraIr`: initialize r-matrix to identity
- `eraRx`: rotate around X-axis
- `eraRxr`: product of two r-matrices

"""
function ecm06(date1, date2)
    r = zeros((3, 3))
    ccall((:eraEcm06, liberfa), Cvoid,
            (Cdouble, Cdouble, Ptr{Cdouble}),
            date1, date2, r)
    r
end


text:
using Trilinos
using Test
using MPI

if Test.get_testset_depth() == 0
  MPI.Init()
end

ccomm = MPI.CComm(MPI.COMM_WORLD)
comm = Teuchos.MpiComm(ccomm)

@test Teuchos.getRank(comm) == MPI.Comm_rank(MPI.COMM_WORLD)
@test Teuchos.getSize(comm) == MPI.Comm_size(MPI.COMM_WORLD)

if Test.get_testset_depth() == 0
  MPI.Finalize()
end

ext: jl


text:
using DataFrames

#VelCopy(du,uInit,tmp,3600.0)
#VelComp(du,uInit,uvetc,3600.0)

"""
    VelComp(du,u,p::Dict,tim)

Interpolate velocity from gridded fields and return position increment `du`
"""
function VelComp(du::Array{Float64,1},u::Array{Float64,1},p::Dict,tim)
    #compute positions in index units
    dt=(tim-p["t0"])/(p["t1"]-p["t0"])
    #
    x,y = u[1:2]
    nx,ny=p["u0"].grid.ioSize
    x,y=[mod(x,nx),mod(y,ny)]
    #
    dx,dy=[x - floor(x),y - floor(y)]
    i_c,j_c = Int32.(floor.([x y])) .+ 1
    i_c==0 ? i_c=nx : nothing
    j_c==0 ? j_c=ny : nothing
    #
    i_w,i_e=[i_c i_c+1]
    x>=nx-1 ? (i_w,i_e)=(nx,1) : nothing
    j_s,j_n=[j_c j_c+1]
    y>=ny-1 ? (j_s,j_n)=(ny,1) : nothing
    #debugging stuff
    if false
        println((x,y,i_c,j_c))
        println((i_w,i_e,j_s,j_n))
        println((dx,dy,dt))
        #println(du)
    end
    #interpolate u to position and time
    du[1]=(1.0-dx)*(1.0-dt)*p["u0"].f[1][i_w,j_c]+
    dx*(1.0-dt)*p["u0"].f[1][i_e,j_c]+
    (1.0-dx)*dt*p["u1"].f[1][i_w,j_c]+
    dx*dt*p["u1"].f[1][i_e,j_c]
    #interpolate v to position and time
    du[2]=(1.0-dy)*(1.0-dt)*p["v0"].f[1][i_c,j_s]+
    dy*(1.0-dt)*p["v0"].f[1][i_c,j_n]+
    (1.0-dy)*dt*p["v1"].f[1][i_c,j_s]+
    dy*dt*p["v1"].f[1][i_c,j_n]
    #
    return du
end

function VelComp(du::Array{Float64,2},u::Array{Float64,2},p::Dict,tim)
    for i=1:size(u,2)
        tmpdu=du[1:2,i]
        VelComp(tmpdu,u[1:2,i],p,tim)
        du[1:2,i]=tmpdu
    end
    return du
end

# VelCopy.jl

"""
    VelCopy(du,u,p::DataFrame,t)

Interpolate velocity from MITgcm float_trajectories output and return
position increment `du`.
"""
function VelCopy(du,u,p::DataFrame,t)
    tt=t/3600.0
    tt0=Int32(floor(tt))
    w=tt-tt0
    du[1]=(1.0-w)*p[tt0+1,:uVel]+w*p[tt0+2,:uVel]
    du[2]=(1.0-w)*p[tt0+1,:vVel]+w*p[tt0+2,:vVel]
end


text:
# Copyright 2016 Martin Holters
# See accompanying license file.

export SimpleSolver, HomotopySolver, CachingSolver

type ParametricNonLinEq
    func::Function
    res::Vector{Float64}
    Jp::Matrix{Float64}
    J::Matrix{Float64}
    function ParametricNonLinEq(func::Function, nn::Integer, np::Integer)
        res = zeros(nn)
        Jp = zeros(nn, np)
        J = zeros(nn, nn)
        return new(func, res, Jp, J)
    end
end

nn(nleq::ParametricNonLinEq) = length(nleq.res)
np(nleq::ParametricNonLinEq) = size(nleq.Jp, 2)

evaluate!(nleq::ParametricNonLinEq, p, z) =
    nleq.func(nleq.res, nleq.J, nleq.Jp, p, z)

type SimpleSolver
    nleq::ParametricNonLinEq
    z::Vector{Float64}
    last_z::Vector{Float64}
    last_p::Vector{Float64}
    last_Jp::Matrix{Float64}
    JLU::Base.LU{Float64,Matrix{Float64}}
    iters::Int
    ressumabs2::Float64
    tol::Float64
    tmp_nn::Vector{Float64}
    tmp_np::Vector{Float64}
    function SimpleSolver(nleq::ParametricNonLinEq, initial_p::Vector{Float64},
                          initial_z::Vector{Float64})
        z = zeros(nn(nleq))
        last_z = zeros(nn(nleq))
        last_p = zeros(np(nleq))
        last_Jp = zeros(nn(nleq), np(nleq))
        JLU = lufact(eye(nn(nleq)))
        tmp_nn = zeros(nn(nleq))
        tmp_np = zeros(np(nleq))
        solver = new(nleq, z, last_z, last_p, last_Jp, JLU, 0, 0.0, 1e-20,
                     tmp_nn, tmp_np)
        set_extrapolation_origin(solver, initial_p, initial_z)
        return solver
    end
end

set_resabs2tol!(solver::SimpleSolver, tol) = solver.tol = tol

function set_extrapolation_origin(solver::SimpleSolver, p, z)
    evaluate!(solver.nleq, p, z)
    JLU = lufact!(solver.nleq.J)
    set_extrapolation_origin(solver, p, z, solver.nleq.Jp, JLU)
end

function set_extrapolation_origin(solver::SimpleSolver, p, z, Jp, JLU)
    copy!(solver.JLU.factors, JLU.factors)
    solver.JLU = Base.LU{Float64, Matrix{Float64}}(solver.JLU.factors, JLU.ipiv,
                                                   JLU.info)
    copy!(solver.last_Jp, Jp)
    copy!(solver.last_p, p)
    copy!(solver.last_z, z)
end

get_extrapolation_origin(solver::SimpleSolver) = solver.last_p, solver.last_z

hasconverged(solver::SimpleSolver) = solver.ressumabs2 < solver.tol

needediterations(solver::SimpleSolver) = solver.iters

function solve(solver::SimpleSolver, p::AbstractVector{Float64}, maxiter=500)
    #solver.z = solver.last_z - solver.JLU\(solver.last_Jp * (p-solver.last_p))
    copy!(solver.tmp_np, p)
    BLAS.axpy!(-1.0, solver.last_p, solver.tmp_np)
    BLAS.gemv!('N', 1.,solver.last_Jp, solver.tmp_np, 0., solver.tmp_nn)
    A_ldiv_B!(solver.JLU, solver.tmp_nn)
    copy!(solver.z, solver.last_z)
    BLAS.axpy!(-1.0, solver.tmp_nn, solver.z)
    local JLU
    for solver.iters=1:maxiter
        evaluate!(solver.nleq, p, solver.z)
        solver.ressumabs2 = sumabs2(solver.nleq.res)
        if ~isfinite(solver.ressumabs2) || ~all(isfinite(solver.nleq.J))
            return solver.z
        end
        JLU = lufact!(solver.nleq.J)
        if JLU.info > 0 # J was singular
            return solver.z
        end
        hasconverged(solver) && break
        #solver.z -= JLU\solver.nleq.res
        copy!(solver.tmp_nn, solver.nleq.res)
        A_ldiv_B!(JLU, solver.tmp_nn)
        BLAS.axpy!(-1.0, solver.tmp_nn, solver.z)
    end
    if hasconverged(solver)
        set_extrapolation_origin(solver, p, solver.z, solver.nleq.Jp, JLU)
    end
    return solver.z
end

type HomotopySolver{BaseSolver}
    basesolver::BaseSolver
    start_p::Vector{Float64}
    iters::Int
    function HomotopySolver(nleq::ParametricNonLinEq,
                            initial_p::Vector{Float64},
                            initial_z::Vector{Float64})
        basesolver = BaseSolver(nleq, initial_p, initial_z)
        return new(basesolver, zeros(np(nleq)), 0)
    end
end

set_resabs2tol!(solver::HomotopySolver, tol) =
    set_resabs2tol!(solver.basesolver, tol)

set_extrapolation_origin(solver::HomotopySolver, p, z) =
    set_extrapolation_origin(solver.basesolver, p, z)

function solve(solver::HomotopySolver, p)
    z = solve(solver.basesolver, p)
    solver.iters = needediterations(solver.basesolver)
    if ~hasconverged(solver)
        a = 0.5
        best_a = 0.0
        copy!(solver.start_p, get_extrapolation_origin(solver.basesolver)[1])
        while best_a < 1
            pa = (1-a) * solver.start_p + a * p
            z = solve(solver.basesolver, pa)
            if hasconverged(solver)
                best_a = a
                a = 1.0
            else
                new_a = (a + best_a) / 2
                if !(best_a < new_a < a)
                    # no floating point value inbetween best_a and a
                    break
                end
                a = new_a
            end
        end
    end
    return z
end

hasconverged(solver::HomotopySolver) = hasconverged(solver.basesolver)
needediterations(solver::HomotopySolver) = solver.iters

type CachingSolver{BaseSolver}
    basesolver::BaseSolver
    ps_tree::KDTree{Vector{Float64}, Matrix{Float64}}
    zs::Matrix{Float64}
    num_ps::Int
    new_count::Int
    new_count_limit::Int
    function CachingSolver(nleq::ParametricNonLinEq, initial_p::Vector{Float64},
                          initial_z::Vector{Float64})
        basesolver = BaseSolver(nleq, initial_p, initial_z)
        ps_tree = KDTree(hcat(initial_p))
        zs = reshape(copy(initial_z), nn(nleq), 1)
        return new(basesolver, ps_tree, zs, 1, 0, 2)
    end
end

set_resabs2tol!(solver::CachingSolver, tol) =
    set_resabs2tol!(solver.basesolver, tol)

hasconverged(solver::CachingSolver) = hasconverged(solver.basesolver)
needediterations(solver::CachingSolver) = needediterations(solver.basesolver)

function solve(solver::CachingSolver, p)
    best_diff = sumabs2(p - get_extrapolation_origin(solver.basesolver)[1])
    idx = 0
    for i in (solver.num_ps-solver.new_count+1):solver.num_ps
        diff = 0.
        for j in 1:size(solver.ps_tree.ps, 1)
            diff += abs2(solver.ps_tree.ps[j,i] - p[j])
        end
        if diff < best_diff
            best_diff = diff
            idx = i
        end
    end

    idx = indnearest(solver.ps_tree, p,
                     Alts([AltEntry(1, zeros(p), 0.0)], best_diff, idx))[1]

    if idx ≠ 0
        set_extrapolation_origin(solver.basesolver,
                                 solver.ps_tree.ps[:,idx], solver.zs[:,idx])
    end

    z = solve(solver.basesolver, p)
    if needediterations(solver.basesolver) > 5 && hasconverged(solver.basesolver)
        solver.num_ps += 1
        if solver.num_ps > size(solver.ps_tree.ps, 2)
            solver.ps_tree.ps =
                copy!(zeros(size(solver.ps_tree.ps, 1), 2solver.num_ps),
                      solver.ps_tree.ps)
            solver.zs =
                copy!(zeros(size(solver.zs, 1), 2solver.num_ps), solver.zs)
        end
        solver.ps_tree.ps[:,solver.num_ps] = p
        solver.zs[:,solver.num_ps] = z
        solver.new_count += 1
    end
    if solver.new_count > 0
        solver.new_count_limit -= 1
    end
    if solver.new_count > solver.new_count_limit
        solver.ps_tree = KDTree(solver.ps_tree.ps, solver.num_ps)
        solver.new_count = 0
        solver.new_count_limit = 2size(solver.ps_tree.ps, 2)
    end
    return z
end

get_extrapolation_origin(solver::CachingSolver) =
    get_extrapolation_origin(solver.basesolver)


text:
using Measurements
using Base.Test

test_approx_eq(a::Measurement, b::Measurement) =
    (@test_approx_eq(a.val, b.val) ; @test_approx_eq(a.err, b.err))
test_approx_eq{T1<:AbstractFloat,T2<:AbstractFloat}(a::Complex{Measurement{T1}}, b::Complex{Measurement{T2}}) =
    (@test_approx_eq(real(a), real(b)) ; @test_approx_eq(imag(a), imag(b)))
test_approx_eq{T1<:AbstractFloat,T2<:AbstractFloat}(a::Complex{Measurement{T1}}, b::Measurement{T2}) =
    test_approx_eq(a, complex(b))
test_approx_eq{T1<:AbstractFloat,T2<:AbstractFloat}(a::Measurement{T1}, b::Complex{Measurement{T2}}) =
    test_approx_eq(complex(a), b)
test_approx_eq_eps(a::Measurement, b::Measurement, tol::Real) =
    (@test_approx_eq_eps(a.val, b.val, tol) ; @test_approx_eq_eps(a.err, b.err, tol))
test_approx_eq_eps{T1<:AbstractFloat,T2<:AbstractFloat}(a::Complex{Measurement{T1}}, b::Complex{Measurement{T2}}, tol::Real) =
    (@test_approx_eq_eps(real(a), real(b), tol) ; @test_approx_eq_eps(imag(a), imag(b), tol))

w = -0.5 ± 0.03
x = 3 ± 0.1
y = 4 ± 0.2
z = complex(x)

# Standard Score
test_approx_eq(stdscore(x, y), -10 ± 2)
@test_approx_eq stdscore(w, x.val) -350/3

# Vectorial version of Measurement
@test measurement(1:2:5, 2:2:6) == [1±2, 3±4, 5±6]

# Weighted Average with Inverse-Variance Weighting
test_approx_eq(weightedmean((w, x, y)),
               measurement(-0.12584269662921355, 0.028442727788398632))

# Derivative
@test_approx_eq(Measurements.derivative(3*x^2, (x.val, x.err, x.tag)), 18)
@test_approx_eq(Measurements.derivative(3*x^2, x), 18)
@test_approx_eq(Measurements.gradient(2x + y - w, [x, y, w]), [2, 1, -1])

# Conversion and Promotion
@test convert(Measurement{Float64}, pi) == pi ± 0
@test convert(Measurement{Float64}, 1//2) == 0.5 ± 0
@test convert(Measurement{Float64}, 3) == 3.0 ± 0.0
@test convert(Measurement{Float64}, 3 ± 1) == 3.0 ± 1.0
@test convert(Measurement, x) === x
@test convert(Measurement, pi) == pi ± 0
@test convert(Measurement, 1//2) == 0.5 ± 0
@test convert(Measurement, 3) == 3 ± 0
@test convert(Signed, x) == 3
@test float(3 ± 1) == 3.0 ± 1.0
@test float(x) === x
@test promote(Measurement{Float32}(3 ± 0.5), Measurement{Float64}(4 ± 0)) ==
    (Measurement{Float64}(3 ± 0.5), Measurement{Float64}(4 ± 0))
@test promote(x, complex(7)) == (complex(3.0 ± 0.1),
                                  complex(measurement(7.0)))
@test promote(complex(0, 1 ± 0), 2.1 ± 0.2) ==
    (complex(0, 1.0 ± 0), complex(2.1 ± 0.2))

##### Comparisons and Tests
@test z == x != y
@test x == x
@test -2 < w <= x < y < 5
@test 3 == x
@test y == 4 != w
@test measurement(0.8, 0.01) == 4//5
@test measurement(π) == π
@test e == measurement(e)
@test 4//5 == measurement(0.8, 0.03)
@test isnan(x) == false
@test isfinite(y) == true && isfinite(measurement(Inf)) == false
@test isinf(measurement(Inf)) == true && isinf(x) == false
@test (isinteger(x) == true && isinteger(w) == false)

##### Mathematical Operations
# Addition
test_approx_eq(+x, x)
test_approx_eq(x + y, measurement(7, 0.22360679774997896))
test_approx_eq(x + y, y + x) # Commutativity
test_approx_eq(2 + x, measurement(5, 0.1))
test_approx_eq(x + 5//2, measurement(5.5, 0.1))
test_approx_eq(x + true, measurement(4, 0.1))
test_approx_eq(x + 2, measurement(5, 0.1))
# Test correlation
for a in (w, x, y); test_approx_eq(a + a + a, 3a); end

# Subtraction
test_approx_eq(-x, measurement(-3, 0.1))
test_approx_eq(x - y, measurement(-1, 0.22360679774997896))
test_approx_eq(2 - x, measurement(-1, 0.1))
test_approx_eq(x - 2, measurement(1, 0.1))
# Test correlation
for a in (w, x, y); test_approx_eq(a - a, measurement(0)); end

# Multiplication
test_approx_eq(x*y, measurement(12, 0.7211102550927979))
test_approx_eq(x*y, y*x) # Commutativity
test_approx_eq(2x, measurement(6, 0.2))
test_approx_eq(x*3, measurement(9, 0.3))
test_approx_eq(w*0, measurement(0))
test_approx_eq(true*x, x)
test_approx_eq(y*false, measurement(0))
for a in (w, x, y); test_approx_eq(a*(0 ± 1), measurement(0, abs(a.val))); end
# Test correlation
for a in (x, y); test_approx_eq(a*a*a, a^3); end
for f in (muladd, fma)
    test_approx_eq(f(w, x, y), w*x + y)
    test_approx_eq(f(w, x, 1), w*x + 1)
    test_approx_eq(f(w, 2, y), w*2 + y)
    test_approx_eq(f(w, 2, 1), w*2 + 1)
    test_approx_eq(f(3, x, y), 3*x + y)
    test_approx_eq(f(3, x, 1), 3*x + 1)
    test_approx_eq(f(3, 2, y), 3*2 + y)
end

# Division
test_approx_eq(x/y, measurement(0.75, 0.04506939094329987))
test_approx_eq(x/10, measurement(0.3, 0.01))
test_approx_eq(1/y, measurement(0.25, 0.0125))
test_approx_eq(0/x, measurement(0))
test_approx_eq(div(1.2*x, w), measurement(-7))
test_approx_eq(div(x, 1.2), measurement(2))
test_approx_eq(div(9.4, y), measurement(2))
test_approx_eq(fld(1.2*x, w), measurement(-8))
test_approx_eq(fld(x, 1.2), measurement(2))
test_approx_eq(fld(9.4, y), measurement(2))
test_approx_eq(cld(1.2*x, w), measurement(-7))
test_approx_eq(cld(x, 1.2), measurement(3))
test_approx_eq(cld(9.4, y), measurement(3))
for a in (w, x, y); test_approx_eq((0 ± 1)/a, measurement(0, 1/abs(a.val))); end
# Test correlation
for a in (w, x, y); test_approx_eq(a/a, measurement(1)); end
# Test derivatives of "div", "fld", and "cld".  They're defined to be exactly 0.
# Should you discover this is not correct, update the test accordingly.
for f in (div, fld, cld), a in (w, x, y), b in (w, x, y)
    Measurements.gradient(f(a, b), [a, b]) == [0.0, 0.0]
end

# Inverse
for a in (w, x, y); test_approx_eq(inv(a), 1/a); end

# signbit
@test signbit(x) == false
@test signbit(w) == true

# Power
test_approx_eq(x^y, measurement(81, 20.818061515800505))
for a in (w, x, y); test_approx_eq(a^(-1), inv(a)); end
test_approx_eq(x^2, measurement(9, 0.6))
for a in (w, x, y); test_approx_eq(a^2, a^2.0); end
test_approx_eq(y^(1//2), measurement(2, 0.05))
for a in (w, x, y); test_approx_eq(a^(4//2), a*a); end
test_approx_eq(x^(pi), measurement(31.54428070019754, 3.3033093503504967))
test_approx_eq(2^x, measurement(8, 0.5545177444479562))
for a in (w, x, y); test_approx_eq(2^a, 2.0^a); end
test_approx_eq(pi^x, measurement(31.006276680299816, 3.5493811564854525))
for val in (w, x, y); test_approx_eq(e^val, exp(val)); end
for val in (w, x, y); test_approx_eq(exp2(val), 2^val); end
test_approx_eq(z^2.5, x^2.5)
test_approx_eq(z^3, x^3)
# Make sure "p ± 0" behaves like "p", in particular with regard to the
# uncertainty.
for p in (-3, 0, 3); test_approx_eq((0 ± 0.1)^(p ± 0), (0 ± 0.1)^p); end

# rad2deg
test_approx_eq(rad2deg(x), measurement(171.88733853924697, 5.729577951308232))

# Cosine
test_approx_eq(cos(x), measurement(-0.9899924966004454, 0.014112000805986721))
test_approx_eq(cosd(x), measurement(0.9986295347545738, 9.134347536190512e-5))
test_approx_eq(cosh(x), measurement(10.067661995777765, 1.0017874927409902))

# Sine
test_approx_eq(sin(y), measurement(-0.7568024953079282, 0.13072872417272238))
test_approx_eq(sind(y), measurement(0.0697564737441253, 0.0034821554353128255))
test_approx_eq(sinh(y), measurement(27.28991719712775, 5.461646567203298))

#Tangent
for val in (w, x, y); test_approx_eq(tan(val), sin(val)/cos(val)); end
for val in (w, x, y); test_approx_eq(tand(val), sind(val)/cosd(val)); end
for val in (w, x, y); test_approx_eq(tanh(val), sinh(val)/cosh(val)); end

# Inverse trig functions
test_approx_eq(acos(w), measurement(2.0943951023931957, 0.034641016151377546))
test_approx_eq(acosd(w), measurement(120, 1.9847840235184515))
test_approx_eq(acosh(x), measurement(1.7627471740390859, 0.035355339059327376))
test_approx_eq(cos(acos(w)), w)
test_approx_eq(cosd(acosd(w)), w)
test_approx_eq(cosh(acosh(x)), x)
test_approx_eq(asin(w), measurement(-0.5235987755982989, 0.034641016151377546))
test_approx_eq(asind(w), measurement(-30, 1.9847840235184515))
test_approx_eq(asinh(x), measurement(1.8184464592320668, 0.0316227766016838))
test_approx_eq(sin(asin(w)), w)
test_approx_eq(sind(asind(w)), w)
test_approx_eq(sinh(asinh(x)), x)
test_approx_eq(atan(w), measurement(-0.4636476090008061, 0.024))
test_approx_eq(atand(w), measurement(-26.56505117707799, 1.3750987083139758))
test_approx_eq(atanh(w), measurement(-0.5493061443340548, 0.04))
for a in (w, x, y); test_approx_eq(tan(atan(a)), a); end
for a in (w, x, y); test_approx_eq(tand(atand(a)), a); end
test_approx_eq(tanh(atanh(w)), w)
test_approx_eq(atan2(x, y), measurement(0.6435011087932844, 0.028844410203711916))
test_approx_eq(atan2(x, 5), measurement(0.5404195002705842, 0.014705882352941178))
test_approx_eq(atan2(-3, y), measurement(-0.6435011087932844, 0.024))

# Reciprocal trig functions
test_approx_eq(csc(y), measurement(-1.3213487088109024, 0.228247438348944))
test_approx_eq(cscd(y), measurement(14.335587026203676, 0.7156144742178738))
test_approx_eq(csch(y), measurement(0.03664357032586561, 0.007333632734561779))
for a in (w, x, y); test_approx_eq(csc(a),  1/sin(a)); end
for a in (w, x, y); test_approx_eq(cscd(a), 1/sin(deg2rad(a))); end
for a in (w, x, y); test_approx_eq(csch(a), 1/sinh(a)); end
test_approx_eq(sec(w), measurement(1.139493927324549, 0.018675251089778414))
test_approx_eq(secd(w), measurement(1.000038078385737, 4.5695512845605615e-6))
test_approx_eq(sech(w), measurement(0.886818883970074, 0.012294426649942352))
for a in (w, x, y); test_approx_eq(sec(a),  1/cos(a)); end
for a in (w, x, y); test_approx_eq(secd(a), 1/cos(deg2rad(a))); end
for a in (w, x, y); test_approx_eq(sech(a), 1/cosh(a)); end
test_approx_eq(cot(x), measurement(-7.015252551434534, 5.021376836040872))
test_approx_eq(cotd(x), measurement(19.08113668772821, 0.6372018679183983))
test_approx_eq(coth(x), measurement(1.0049698233136892, 0.000996434577114765))
for a in (w, x, y); test_approx_eq(cot(a),  1/tan(a)); end
for a in (w, x, y); test_approx_eq(cotd(a), 1/tan(deg2rad(a))); end
for a in (w, x, y); test_approx_eq(coth(a), 1/tanh(a)); end

# Exponentials
test_approx_eq(exp(x), measurement(20.085536923187668, 2.008553692318767))
for a in (w, x, y); test_approx_eq(expm1(a), exp(a) - 1); end
for a in (w, x, y); test_approx_eq(exp10(a), 10^a); end
for a in (w, 3//5*w, x/10, x, y/50, y); test_approx_eq(ldexp(frexp(a)...), a); end

# Logarithm
test_approx_eq(log(x, y), measurement(1.261859507142915, 0.059474298734200806))
test_approx_eq(log(y), measurement(1.3862943611198906, 0.05))
test_approx_eq(log(y, 4), measurement(1, 0.03606737602222409))
for a in (x, y); test_approx_eq(log(e, a), log(a)); end
for a in (x, y); test_approx_eq(log(2, a), log2(a)); end
for a in (x, y); test_approx_eq(log(10, a), log10(a)); end
for a in (x, y); test_approx_eq(log1p(a), log(1 + a)); end
test_approx_eq(log(pi, x), measurement(0.9597131185693899, 0.029118950894341064))
for a in (w, x, y); test_approx_eq(log(exp(a)), a); end
for a in (x, y); test_approx_eq(exp(log(a)), a); end

# Hypotenus
test_approx_eq(hypot(x, y), measurement(5, 0.17088007490635065))
test_approx_eq(hypot(x, y), hypot(y, x)) # Commutativity
test_approx_eq(hypot(x, 4), measurement(5, 0.06))
test_approx_eq(hypot(3, y), measurement(5, 0.16))

# Square root
test_approx_eq(sqrt(y), measurement(2, 0.05))
for a in (x, y); test_approx_eq(sqrt(a), a^0.5); end
for a in (x, y); test_approx_eq(sqrt(a)*sqrt(a), a); end
# Derivative of sqrt diverges in 0, but if the measurement is exact (like "a-a"
# is) also the resulting quantity must have 0 uncertainty.
for a in (w, x, y); test_approx_eq(sqrt(a - a), measurement(0)); end

# Cube root
test_approx_eq(cbrt(x), measurement(1.4422495703074083, 0.01602499522563787))
for a in (x, y); test_approx_eq(cbrt(a), a^(1/3)); end
for a in (x, y); test_approx_eq(cbrt(a)*cbrt(a)^2, a); end

# Absolute value
test_approx_eq(abs(-x), x)
test_approx_eq(abs(measurement(0, 2)), measurement(0, 2))
for a in (w, x, y); test_approx_eq(abs2(a), a*a); end

# Zero
test_approx_eq(zero(x), measurement(0))

# Sign
test_approx_eq(sign(x), measurement(1))
test_approx_eq(sign(-y), measurement(-1))
test_approx_eq(sign(measurement(0, 5)), measurement(0))
@test copysign(x, -5) == -x
@test copysign(x, w) == -x
@test copysign(5, w) == -5
@test copysign(1//2, w) == -1//2
@test copysign(w, w) == w
@test copysign(Float32(5), w) == -5
@test copysign(Float64(5), w) == -5
@test copysign(pi, w) == -pi
@test flipsign(x, -5) == -x
@test flipsign(w, w) == -w
@test flipsign(-5, w) == 5
@test flipsign(Float32(-5), w) == 5
@test flipsign(Float64(-5), w) == 5
@test flipsign(pi, w) == -pi

# One
test_approx_eq(one(y), measurement(1))

# Error function
test_approx_eq(erf(x), measurement(0.9999779095030014, 1.3925305194674787e-5))
test_approx_eq(erfinv(erf(w)), w)
test_approx_eq(erfc(w), 1 - erf(w))
for a in (w, x, y); test_approx_eq(erfcinv(erfc(a)), a); end
for a in (w, x, y); test_approx_eq(erfcx(a), exp(a^2)*erfc(a)); end
test_approx_eq(erfi(x), 1629.9946226015657 ± 914.3351093102547)
for a in (w, x, y); test_approx_eq(dawson(a), 0.5*exp(-a^2)*erfi(a)*sqrt(pi)); end

# Factorial and gamma
test_approx_eq(factorial(x), measurement(6, 0.7536706010590813))
for a in (w, x, y); test_approx_eq(gamma(a), factorial(a - 1)); end
for a in (w, x, y); test_approx_eq(gamma(a + 1), factorial(a)); end
for a in (x, y); test_approx_eq(lgamma(a), log(gamma(a))); end
test_approx_eq(digamma(y), 1.256117668431802 ± 0.056764591147422994)
for a in (w, x, y); test_approx_eq(digamma(a), polygamma(0, a)); end
for a in (w, x, y); test_approx_eq(digamma(invdigamma(a)), a); end
test_approx_eq(polygamma(3, w), 193.40909103400242 ± 0.10422749480000776)
for a in (w, x, y); test_approx_eq(trigamma(a), polygamma(1, a)); end

# Beta function
for a in (w, x, y); test_approx_eq(beta(a, x), gamma(a)*gamma(x)/gamma(a + x)); end
for a in (w, x, y); test_approx_eq(beta(a, pi), gamma(a)*gamma(pi)/gamma(a + pi)); end
for a in (w, x, y); test_approx_eq(beta(e, a), gamma(e)*gamma(a)/gamma(e + a)); end
for a in (x, y); test_approx_eq(lbeta(a, x), log(beta(a, x))); end
for a in (x, y); test_approx_eq(lbeta(a, pi), log(beta(a, pi))); end
for a in (x, y); test_approx_eq(lbeta(e, a), log(beta(e, a))); end

# Airy functions
test_approx_eq(airy(0, x), 0.006591139357460721 ± 0.0011912976705951322)
test_approx_eq(airy(1, x), -0.01191297670595132 ± 0.0019773418072382165)
test_approx_eq(airy(2, x), 14.037328963730136 ± 2.2922214966382017)
test_approx_eq(airy(3, x), 22.922214966382015 ± 4.211198689119041)
for a in (w, x, y); test_approx_eq(airy(a), airy(0, a)); end
for a in (w, x, y); test_approx_eq(airyai(a), airy(0, a)); end
for a in (w, x, y); test_approx_eq(airyprime(a), airy(1, a)); end
for a in (w, x, y); test_approx_eq(airyaiprime(a), airy(1, a)); end
for a in (w, x, y); test_approx_eq(airybi(a), airy(2, a)); end
for a in (w, x, y); test_approx_eq(airybiprime(a), airy(3, a)); end

# Bessel functions
for a in (x, y); test_approx_eq(2*besselj1(a)/a, besselj0(a) + besselj(2, a)); end
for a in (x, y); test_approx_eq(besselj(-1/2, a), sqrt(2/(pi*a))*cos(a)); end
for a in (x, y); test_approx_eq(besselj(1/2, a), sqrt(2/(pi*a))*sin(a)); end
test_approx_eq(bessely0(x), 0.3768500100127904 ± 0.03246744247917999)
for a in (x, y); test_approx_eq(bessely1(a), -bessely(-1, a)); end
for a in (x, y); test_approx_eq(bessely(5/2, a),
                                (besselj(5/2, a)*cos(2.5pi) -
                                 besselj(-5/2, a))/sin(2.5pi)); end
for a in (x, y), k in (1, 2), nu in -1:1
    sgn = k == 1 ? +1 : -1
    test_approx_eq(besselh(nu, k, a),
                   besselj(nu, a) + sgn*im*bessely(nu, a))
end
test_approx_eq(besseli(5//2, y), 4.757626874823528 ± 1.0398232869843944)
for a in (x, y); test_approx_eq(besselix(3, a), besseli(3, a)*exp(-abs(a))); end
test_approx_eq(besselk(7//3, x), 0.07521953258226349 ± 0.010340691203742959)
for a in (x, y); test_approx_eq(besselkx(3//4, a), besselk(3//4, a)*exp(a)); end

# Modulo
let
    frac, int = modf(x + w)
    test_approx_eq(frac, x + w - int)
    @test_approx_eq int floor((x + w).val)
end

for a in (x, y) # Test property of "mod" function
    b = 10.0*w
    test_approx_eq(mod(a, b), a - fld(a, b)*b)
end
test_approx_eq(mod(y, -3), y + fld(y, -3)*3)
test_approx_eq(mod(-5.8, x), 0.2 ± 0.2)

for a in (x, y) # Test property of "rem" function
    b = 10.0*w
    test_approx_eq(rem(a, b), a - div(a, b)*b)
end
test_approx_eq(rem(y, -3), y + div(y, -3)*3)
test_approx_eq(rem(-5.8, x), -2.8 ± 0.1)
test_approx_eq(mod2pi(pi*x), measurement(pi, 0.1*pi))

# Machine precisionx
@test_approx_eq eps(Measurement{Float64}) eps(Float64)
@test_approx_eq eps(x) eps(x.val)
@test_approx_eq nextfloat(x) nextfloat(x.val)
@test_approx_eq maxintfloat(Measurement{Float64}) maxintfloat(Float64)

# Rounding
@test_approx_eq round(w) round(w.val)
@test_approx_eq round(Int, w) round(Int, w.val)
@test_approx_eq floor(w) floor(w.val)
@test_approx_eq floor(Int, w) floor(Int, w.val)
@test_approx_eq ceil(w) ceil(w.val)
@test_approx_eq ceil(Int, w) ceil(Int, w.val)
@test_approx_eq trunc(w) trunc(w.val)
@test_approx_eq trunc(Int, w) trunc(Int, w.val)

# Widening
@test widen(Measurement{Float32}) == Measurement{Float64}

# To big float and widemul
test_approx_eq(widemul(x, y), big(x)*big(y))

# Dummy call to show
show(DevNull, x)

# Sum
test_approx_eq(sum((x, w, x, y, -w)), 2x + y)
test_approx_eq(sum([x, w, x, y, -w]), 2x + y)

# Prod
test_approx_eq(prod((w, x, y, -w)), -w^2 * x * y)
test_approx_eq(prod([w, x, y, -w]), -w^2 * x * y)
test_approx_eq(prod((w, x, y, -w, 0)), 0 ± 0)
test_approx_eq(prod([w, x, y, -w, 0]), 0 ± 0)

##### Some functions not defined in Measurements.jl that work out-of-the-box.
# Mean
test_approx_eq(mean((x, w, x, y, -w)), (2x + y)/5)
test_approx_eq(mean([x, w, x, y, -w]), (2x + y)/5)

# min, max, extrema
@test min(w, x, y) === w
@test max(w, x, y) === y
@test extrema([w, x, y]) == (w, y)

# sort
@test sort([y, w, x]) == [w, x, y]

##### Calculations with NaNs
# NaN as nominal value
@test isequal(2*(NaN ± 3), NaN ± 6)
# NaN as uncertainty
@test isequal(value(2*(3 ± NaN)), 6) && isequal(uncertainty(2*(3 ± NaN)), NaN)
# Both
@test isequal(value(2*(NaN ± NaN)), NaN) && isequal(uncertainty(2*(NaN ± NaN)), NaN)

##### Test @uncertain macro
test_approx_eq_eps(@uncertain(tan(x)), tan(x), 2e-11)
test_approx_eq_eps(@uncertain((a -> a + a + a)(x)), 3x, 3e-12)
test_approx_eq(@uncertain(zeta(x)),
               measurement(1.2020569031595951, 0.019812624290876782))
for f in (log, hypot, atan2); test_approx_eq_eps(@uncertain(f(x, y)), f(x, y), 2e-12); end
test_approx_eq_eps(@uncertain(((a,b,c,d,e,f) -> a+b+c+d+e+f)(x, 2x, y, log(y), -w, w^2)),
                   3x + y + log(y) - w + w^2, 7e-12)

let
    # Test with a "ccall"
    f(x) = x*x
    ptr = cfunction(f, Cdouble, (Cdouble,))
    g(x) = ccall(ptr, Cdouble, (Cdouble,), x)*x
    test_approx_eq_eps(@uncertain(g(x)), x^3, 4e-11)
end

# Test getindex with Derivatives type
@test_throws KeyError getindex(x.der, 0)

##### value, uncertainty
@test value([w, x, y]) == [-0.5, 3.0, 4.0]
@test value([complex(w, x)]) == [complex(-0.5, 3.0)]
@test uncertainty([w, x, y]) == [0.03, 0.1, 0.2]
@test uncertainty([complex(w, x)]) == [complex(0.03, 0.1)]

##### Test `length' method
@test length((w + w + 2x + y).der) == 3

##### Big Floats
a = big"3.00000001" ± big"1e-17"
b = big"4.00000001" ± big"1e-17"
c = big"5.00000001" ± big"1e-15"
for X in (a, b, c); test_approx_eq(sin(X*X + X*X)/cos(X*X + X*X), tan(2X^2)); end
test_approx_eq(a*b,
               big"12.0000000700000001" ±
               big"5.000000014000000000399999998880000003119999991353600023834879934652928178154734e-17")
test_approx_eq(a/x + w/b,
               big"8.750000036458332770409319189914605364095140370607375793457032012939451217651416e-01" ±
               big"3.416666677095189699499391052146002403356891669706809079678073000173723744722607e-02")

include("complex.jl")


text:
###############################
# Trait-functions with dispatch
###############################
#
# Provides @traitfn macro which can be used like so:
#
# @traitfn f1{X,Y<:Integer; D1{Y}, D4{X,Y}}(x::X,y::Y) = x + sin(y)
# @traitfn f1{S,T<:Integer; D1{S}, D1{T}  }(s::S,t::T) = sin(s) - sin(t)
# @traitfn f1{X,Y<:AbstractFloat; D1{X}, D1{Y}  }(x::X,y::Y) = cos(x) - cos(y)

typealias FName Union{Symbol,Expr}

# generates: X1, X2,... or x1, x2.... (just symbols not actual TypeVar)
type GenerateTypeVars{CASE}
end
Base.start(::GenerateTypeVars) = 1
Base.next(::GenerateTypeVars{:upcase}, state) = (symbol("X$state"), state+1) # X1,..
Base.next(::GenerateTypeVars{:lcase}, state) = (symbol("x$state"), state+1)  # x1,...
Base.done(::GenerateTypeVars, state) = false

# Type to hold parsed function defs:
type ParsedFn  # (probably should adapt MetaTools.jl...)
    name::FName  # f1
    fun::Expr # :(f1{X<:Int,Y})
    typs::Vector{Any} # [:(X<:Int),:Y]
    sig::Vector{Any} # [:(x::X), :(y::Y)]
    traits::Vector{Any} # [D1{X}, D2{X,Y}]
    body::Expr # quote ... end
end

# Parsing:
function parsetraitfn_head(head::Expr)
    # Transforms
    # f1{X<:Int,Y; D1{X}, D2{X,Y}}(x::X,y::Y)
    #
    # into a ParsedFn

    nametyp = head.args[1]
    sig = head.args[2:end]
    name = nametyp.args[1]
    fun = Expr(:curly, deepcopy(nametyp.args[[1;3:end]])...)
    typs = fun.args[2:end]
    traits = nametyp.args[2].args
    return ParsedFn(name, fun, typs, sig, traits, :())
end

gettypesymbol(x::Expr) = x.args[1] # :(X1<:Int)
gettypesymbol(x::Symbol) = x
function translate_head(fn::ParsedFn)
    # Takes output from parsetraitfn_head and
    # renames sig and TypeVar:
    # f1{X,Y; D1{X}, D2{X,Y}}(x::X,y::Y)
    # ->
    # f1{X1,X2; D1{X1}, D2{X1,X2}}(x::X1,y::X2)
    #
    # Returns translated ParsedFn

    function make_trans(sig, typs)
        # makes two dictionaries with keys the old typevars, and
        # values the lowercase and uppercase variables

        # make variable-symbol translation map:
        trans_var = Dict{Symbol,Symbol}()
        for (i,si) in enumerate(GenerateTypeVars{:lcase}())
            if length(sig)<i;  break end
            trans_var[sig[i].args[1]] = si
        end

        trans_Tvar = Dict{Symbol,Symbol}()
        # make type-symbol translation map:
        for (i,tv) in enumerate(GenerateTypeVars{:upcase}())
            if length(typs)<i;  break end
            trans_Tvar[gettypesymbol(typs[i])] = tv
        end
        return trans_var, trans_Tvar
    end
    trans_var, trans_Tvar = make_trans(fn.sig, fn.typs)

    # do the translations:
    fnt = deepcopy(fn)
    for i in 2:length(fnt.fun.args)
        tv = fnt.fun.args[i]
        if isa(tv, Expr)
            # :(X<:Int)
            tv.args[1] = trans_Tvar[tv.args[1]]
        else
            # :X
            fnt.fun.args[i] = trans_Tvar[tv]
        end
    end
    fnt.typs = fnt.fun.args[2:end]
    for s in fnt.sig
        s.args[1] = trans_var[s.args[1]]
        if isa(s, Expr)
            t = s.args[2]
            if isa(t, Symbol)
                if haskey(trans_Tvar, s.args[2])
                    s.args[2] = trans_Tvar[s.args[2]]
                end
            elseif t.head==:curly
                for i=2:length(t.args)
                    if haskey(trans_Tvar, t.args[i])
                        t.args[i] = trans_Tvar[t.args[i]]
                    end
                end
            else
                error("Cannot parse $t in $(fnt.sig)")
            end
        end
    end
    for t in fnt.traits
        t.args[2:end] = map(x->trans_Tvar[x], t.args[2:end])
    end

    return fnt
end

function parsetraitfn(fndef::Expr)
    # Transforms
    # f1{X,Y; D1{X}, D2{X,Y}}(x::X,y::Y) = sin(x) + y
    #
    # into a ParsedFn

    # checks
    length(fndef.args)==2 || throw(TraitException("Something is wrong with $fndef"))
    # parse
    head = fndef.args[1]
    body = fndef.args[2]
    fn = parsetraitfn_head(head)
    fn.body = body
    fnt = translate_head(fn)

    return fn, fnt
end

# Macro building blocks:
function makefnhead(fname, typs, sig)
    # Makes a Expr for a function head:
    # f1, {X,Y}, (x::X, r::Y) -> f1{X,Y}(x::X, r::Y)
    fntyp = Expr(:curly)
    push!(fntyp.args, fname)
    append!(fntyp.args, typs)
    outfn = Expr(:call)
    push!(outfn.args, fntyp)
    append!(outfn.args, sig)
    return outfn
end
function strip_typeasserts(sig)
    # [:(x::X), :(r::R)] -> [:x,:r]
    sig = deepcopy(sig)
    for i = 1:length(sig)
        if isa(sig[i], Expr) && sig[i].head==:(::)
            # :(x::X)
            sig[i] = sig[i].args[1]
        end
    end
    return sig
end
function makefncall(fname, sig)
    # Makes an Expr to call function fname:
    # f1, (x::X, r::R) -> f1(x, r)
    # (i.e. strips the extraneous type-assertions)
    return :( $fname( $(strip_typeasserts(sig)...) ) )
end
function get_concrete_type_symb(typs)
    # Returns the most general type satisfying typs as a list of symbols:
    # [:(X<:Int), :Y] -> [:Int, :Any]
    out = Any[]
    for t in typs
        if isa(t, Symbol)
            push!(out, :Any)
        elseif  t.head==:.
            push!(out, t)
        else
            push!(out, t.args[2])
        end
    end
    return out
end
function get_concrete_type_Typetuple(typs)
    # Returns the most general type satisfying typs as a tuple of
    # actual types:
    # [:(X<:Int), :Y] -> Tuple{Int, Any}
    out = get_concrete_type_symb(typs)
    for i in 1:length(out)
        out[i] = Type{eval_curmod(out[i])}
    end
    return Tuple{out...}
end

function make_Type_sig(typs)
    # [:(X<:Int), :Y] -> [:(::Type{X}), :(::Type{Y})]
    out = Any[]
    for t in typs
        push!(out, :(::Type{$t}))
        # if isa(t, Symbol) || t.head==:.

        # else
        #     push!(out, :(::Type{$(t.args[1])}))
        # end
    end
    return out
end
function has_only_one_method(fname::Symbol, typs)
    # Checks whether fname has one and only one method for types in
    # typs.
    if isdefined(current_module(), fname)
        1 == length( methods(eval_curmod(fname),
                             get_concrete_type_Typetuple(typs)) )
    else
        false
    end
end
function has_only_one_method(fname::Expr, typs)
    # Checks whether fname has one and only one method for types in
    # typs.
    if isdefined(eval_curmod(fname.args[1]), fname.args[2].args[1])
        1 == length( methods(eval_curmod(fname),
                             get_concrete_type_Typetuple(typs)) )
    else
        false
    end
end

@doc """The heart, the trait-dispatch function.

     Trait-function (TF) dispatch works like:

     - first dispatch on the normal types

     Then dispatch on traits using the following rules, terminating
     when only one or zero possibilities are left

     - find all matching traits
     - discriminate using subtraits, i.e. a subtrait will win over its supertrait
     - score all traits according to:
       1 point for all single parameter traits,
       2 points for all two parameter traits,
       etc.
       Now pick the highest scoring method.
     - if still ambiguous throw an error
     """->
function traitdispatch(traittypes, fname)
    poss = Any[]
    for Tr in traittypes
        if istrait(Tr)
            push!(poss, Tr)
        end
    end
    # discriminate using subtraits
    if length(poss)>1
        topurge = []
        for (i,p1) in enumerate(poss)
            for p2 in poss
                if length(p1)==length(p2)
                    checks = true
                    for j=1:length(p1)
                        checks = checks && issubtrait(p2[j], p1[j]) && p2[j]!=p1[j]
                    end
                    if checks
                        push!(topurge, i)
                    end
                end
            end
        end
        topurge = sort(unique(topurge))
        deleteat!(poss, topurge)
    end
    # check whether we can discriminate on more specific definitions:
    if length(poss)>1
        # strategy:
        # - each single parameter trait gives 1 point, 2-parameter 2 points, etc
        # - pick method with most points
        #
        # This is not the end of the story but better...
        score = zeros(Int, length(poss))
        for (i,p1) in enumerate(poss)
            for t in p1
                score[i] += length(t.parameters)
            end
        end
        poss = poss[find(maximum(score).==score)]
    end
    # check we got a unique trait for dispatch
    if length(poss)==0
        throw(TraitException("No matching trait found for function $(string(fname))"))
    elseif length(poss)>1
        throw(TraitException("For function $(string(fname)) there are several matching traits:\n $poss"))
    end
    return poss[1]
end

@doc """Functions which also dispatch using traits can be defined with
     `@traitfn`.  Once defined they can be used like ordinary
     functions.  For instance a library user would not need to treat
     them specifically.

     Example continuing from the documentation of `@traitdef` and `@traitimpl`:

     ```
     @traitfn tf{X,Y; MyArith{X,Y}}(x::X, y::Y) = x/y+y
     tf(5, Int8(4)) # -> 5.25
     tf(A(2), AB(4)) # -> A(4.5)
     tf(5, Uint8(4)) # -> ERROR: TraitException("No matching trait found for function tf")
     ```

     - all the arguments of `tf` which participate in trait-dispatch
       need to be parameterized.
     - trait-methods and non-trait methods can be mixed.
     - for details on trait dispatch see doc of Traits.traitdispatch
     """ ->
macro traitfn(fndef)
    fn, fnt = parsetraitfn(fndef)
    if length(unique(fnt.traits))!=length(fnt.traits)
        throw(TraitException(
          "There are repeated traits in the trait signature of $(fndef.args[1])"))
    end

    ## make primary function: f
    #### tf(x, y)
    # (Just overwrite definitions of f if they exists already,
    # generates warnings though...)

    # definition head: fn{X,Y}(x::X,y::Y)
    f = makefnhead(fn.name, fn.typs, fn.sig)
    # definition body: _trait_fn(_trait_type_f1(x,y) ), x, y)
    args1 = Any[:(Traits._TraitDispatch), fn.sig...]
    args2 = Any[makefncall(fn.name, args1), fn.sig...]
    body = makefncall(fn.name, args2)
    f = :($f = $body)

    ## make function containing the logic: trait_f
    #### tf(::Type{Tuple{Traits...}}, x, y)
    # 1) make the traits-type
    trait_typ = :(Tuple{})
    append!(trait_typ.args, fn.traits)
    trait_typ = :(::Type{$trait_typ})
    args = Any[trait_typ, fn.sig...]
    trait_f = makefnhead(fn.name, fn.typs, args)
    trait_f = :($trait_f = $(fn.body))
    # add @inline
    trait_fn = Expr(:macrocall, symbol("@inline"), trait_f)

    ## Make function storing the trait-types
    #### tf(Traits._TraitStorage, sig...)
    # This function will return all defined Trait-tuples for a certain
    # signature. Example:
    #
    # julia> yt1(Traits._TraitStorage, Int, Int)
    # (Any[ Tuple{Traits.Arith{Int64,Int64}} ], Any[:( Tuple{Arith{X1,X2}} )], Any[:X1,:X2])

    ## 1) Get the existing traits out of the trait_type_f:
    #    These can be retrieved with the call:
    #    trait_type_f(Traits._TraitStorage, ::Type{X}, ::Type{Y}...) for suitable X, Y...
    args1 = Any[:(Traits._TraitStorage), get_concrete_type_symb(fn.typs)...]
    trait_type_f_store_call = makefncall(fn.name, args1)

    args2 = Any[:(Traits._TraitStorage), fn.typs...]
    if has_only_one_method(fn.name, args2)
        traittypes = eval_curmod(trait_type_f_store_call)[2]
    else
        traittypes = Any[]
    end

    ## 2) update old_traittypes with the new ones
    newtrait = :( Tuple{$(fnt.traits...)} )
    if !(newtrait in traittypes)
        push!(traittypes, newtrait)
    end

    ## 3) make new trait-type storage function
    #     tf(::Type{Traits._TraitStorage}, ::Type{X}, ::Type{Y}...)
    sig2typs(sig) = [s.args[2] for s in fnt.sig]
    sig = make_Type_sig([:(Traits._TraitStorage), sig2typs(fnt.sig)...])
    trait_type_f_store_head = makefnhead(fn.name,
                                          fnt.typs, sig)
    trait_type_f_store = :(
                            $trait_type_f_store_head = (Any[$(traittypes...)], )
                            )
    push!(trait_type_f_store.args[2].args, traittypes)
    push!(trait_type_f_store.args[2].args, fnt.typs)

    ## make trait-dispatch stagedfunction: trait_type_f
    #### tf(Traits._TraitDispatch, sig...)

    sig = make_Type_sig([:(Traits._TraitDispatch)])
    append!(sig,fnt.sig)
    trait_type_f = Expr(:stagedfunction, makefnhead(fn.name, fnt.typs, sig))
    args = Any[:(Traits._TraitStorage), fnt.sig...]
    body = quote
        # get the stored trait-types:
        traittypes = $(makefncall(fn.name, args))[1]
        traittyp = Traits.traitdispatch(traittypes, $(fn.name))
        # construct function from traittyp
        out = :(Tuple{})
        for s in traittyp
            push!(out.args, s)
        end
        return out
    end
    push!(trait_type_f.args, body)

    ## now put all together
    ####
    out = quote
        $trait_fn
        $trait_type_f_store
        $trait_type_f
        $f
    end
    return esc(out)
end

##########
# Helper functions
##########

function traitmethods(f::Function, nsig=Tuple{Vararg{Any}}; print=false)
    out = Any[]
    for m in methods(f, Tuple{Any, nsig...})
        if isa(m.sig[1], Type) && m.sig[1].parameters[1]<:Tuple{Vararg{Traits.Trait}}
            push!(out, m)
        end
    end
    if print # pretty print it
        if length(out)>0
            println("Function $f has the following trait methods:")
        else
            println("Function $f has no trait methods.")
            return nothing
        end
        for m in out
            showtraitmethod(m)
        end
        return nothing
    else
        return out
    end
end

# adapted from methodshow.jl:
showtraitmethod(m::Method) = showtraitmethod(STDOUT, m)
function showtraitmethod(io::IO, m::Method)
    print(io, m.func.code.name)
    tv, decls, file, line = Base.arg_decl_parts(m)
    if isempty(tv)
        error("Not a trait-method")
    end
    tv = string(tv)[2:end-1]
    traits = strip(decls[1][2][7:end-2], ',')
    sig = [isempty(d[2]) ? d[1] : d[1]*"::"*d[2] for d in decls[2:end]]
    out = "{$tv; $traits}($(sig...))\n"
    print(io, out)
end


text:
### FieldTuple types

# FieldTuple is a thin wrapper around a Tuple or NamedTuple holding some Fields
# and behaving like a Field itself
struct FieldTuple{B<:Basis,FS<:Union{Tuple,NamedTuple},T} <: Field{B,Spin,Pix,T}
    fs::FS
    # the constructor for FieldTuples is a bit complex because there's alot of
    # different ways we want to allow constructing them. see the
    # FieldTuples/Constructors section of runtests.jl for a list. the use of
    # @generated is an easier way to deal with all the different cases and avoid
    # ambiguities vs. trying to do this via dispatch.
    @generated function (::Type{FT})(fs::Union{Tuple,NamedTuple}) where {FT<:FieldTuple}
        T =  (eltype(FT)==Any) ? :(promote_type(map(eltype,values(fs))...))     : eltype(FT)
        B = (basis(FT)==Basis) ? :(BasisTuple{Tuple{map(basis,values(fs))...}}) : basis(FT)
        getnames(::Type{FT}) where {Names, FT<:FieldTuple{<:Any, <:NamedTuple{Names}}} = Names
        getnames(::Any)  = nothing
        Names = getnames(FT)
        fs′ = (Names==nothing) ? :fs : :(NamedTuple{$Names}(fs))
        quote
            fs′ = $fs′
            $(Expr(:new, :(FieldTuple{$B,typeof(fs′),$T}), fs′))
        end
    end
end
# FieldTuple(args...) or FieldTuple(;kwargs...) calls the inner constructor
# which takes a single Tuple/NamedTuple:
(::Type{FT})(;kwargs...) where {FT<:FieldTuple} = FT((;kwargs...))
(::Type{FT})(f1::Field,f2::Field,fs::Field...) where {FT<:FieldTuple} = FT((f1,f2,fs...))
# converting a FieldTuple to a different kind of FieldTuple just changes the
# basis then asserts we end up with the right type:
(::Type{FT})(ft::FieldTuple) where {B<:Basis,FT<:FieldTuple{B}} = FieldTuple{B}(ft.fs) :: FT

### printing
getindex(f::FieldTuple,::Colon) = vcat(getindex.(values(f.fs),:)...)[:]
getindex(D::DiagOp{<:FieldTuple}, i::Int, j::Int) = (i==j) ? D.diag[:][i] : diagzero(D, i, j)
@show_datatype show_datatype(io::IO, t::Type{FT}) where {B,Names,T,FS,FT<:FieldTuple{B,NamedTuple{Names,FS},T}} =
    print(io, "Field$(tuple_type_len(FS))Tuple{$(Names), $(B), $(T)}")
@show_datatype show_datatype(io::IO, t::Type{FT}) where {B,T,FS<:Tuple,FT<:FieldTuple{B,FS,T}} =
    print(io, "Field$(tuple_type_len(FS))Tuple{$(B), $(T)}")

### array interface
size(f::FieldTuple) = (sum(map(length, f.fs)),)
copyto!(dest::FT, src::FT) where {FT<:FieldTuple} = (map(copyto!,dest.fs,src.fs); dest)
iterate(ft::FieldTuple, args...) = iterate(ft.fs, args...)
getindex(f::FieldTuple, i::Union{Int,UnitRange}) = getindex(f.fs, i)
fill!(ft::FieldTuple, x) = (map(f->fill!(f,x), ft.fs); ft)
adapt_structure(to, f::FieldTuple{B}) where {B} = FieldTuple{B}(map(f->adapt(to,f),f.fs))
similar(ft::FT) where {FT<:FieldTuple} = FT(map(f->similar(f),ft.fs))
function similar(ft::FT, ::Type{T}, dims::Dims) where {T<:Number, B, FT<:FieldTuple{B}}
    @assert size(ft)==dims "Tried to make a field similar to $F but dims should have been $(size(f)), not $dims."
    FieldTuple{B}(map(f->similar(f,T),ft.fs))
end
mapreduce(func, op, ft::FieldTuple; kw...) = mapreduce(f->mapreduce(func, op, f; kw...), op, ft.fs; kw...)

### broadcasting
struct FieldTupleStyle{B,Names,FS<:Tuple} <: AbstractArrayStyle{1} end
(::Type{FTS})(::Val{1}) where {FTS<:FieldTupleStyle} = FTS()
BroadcastStyle(::FieldTupleStyle{B1}, ::FieldTupleStyle{B2}) where {B1,B2} = invalid_broadcast_error(B1,B2)
BroadcastStyle(::Type{FT}) where {B,FS<:Tuple,FT<:FieldTuple{B,FS}} = FieldTupleStyle{B,Nothing,Tuple{map_tupleargs(typeof∘BroadcastStyle,FS)...}}()
BroadcastStyle(::Type{FT}) where {B,Names,FS,NT<:NamedTuple{Names,FS},FT<:FieldTuple{B,NT}} = FieldTupleStyle{B,Names,Tuple{map_tupleargs(typeof∘BroadcastStyle,FS)...}}()
similar(::Broadcasted{FTS}, ::Type{T}) where {T, FTS<:FieldTupleStyle} = similar(FTS,T)
similar(::Type{FieldTupleStyle{B,Nothing,FS}}, ::Type{T}) where {B,FS,T} = FieldTuple{B}(map_tupleargs(F->similar(F,T), FS))
similar(::Type{FieldTupleStyle{B,Names,FS}}, ::Type{T}) where {B,Names,FS,T} = FieldTuple{B}(NamedTuple{Names}(map_tupleargs(F->similar(F,T), FS)))
instantiate(bc::Broadcasted{<:FieldTupleStyle}) = bc
fieldtuple_data(f::FieldTuple, i) = f.fs[i]
fieldtuple_data(f::Field, i) = f
fieldtuple_data(x, i) = x
function copyto!(dest::FieldTuple, bc::Broadcasted{Nothing})
    for (i,d) in enumerate(dest.fs)
        copyto!(d, map_bc_args(arg->fieldtuple_data(arg,i), bc))
    end
    dest
end

### promotion
function promote(ft1::FieldTuple, ft2::FieldTuple)
    fts = map(promote,ft1.fs,ft2.fs)
    FieldTuple(map(first,fts)), FieldTuple(map(last,fts))
end

### conversion
# The basis we're converting to is always B′. The FieldTuple's basis is B (if
# its different). Each of these might be a concrete basis or a BasisTuple, and
# the FieldTuple might be named or not. And we have Basis(f) which should be a
# no-op. This giant matrix of possibilities presents some ambiguity problems,
# hence why the rules below are so lengthy. Perhaps there's a more succinct
# way to do it, but for now this works.
#
#
# cases where no conversion is needed
Basis(f::FieldTuple{<:Basis}) = f
Basis(f::FieldTuple{<:BasisTuple}) = f
(::Type{B′})(f::F)  where {B′<:Basis,      F<:FieldTuple{B′}} = f
(::Type{B′})(f::F)  where {B′<:BasisTuple, F<:FieldTuple{B′,<:Tuple}} = f
(::Type{B′})(f::F)  where {B′<:BasisTuple, Names,F<:FieldTuple{B′,<:NamedTuple{Names}}} = f
# cases where FieldTuple is in BasisTuple
(::Type{B′})(f::F) where {B′<:BasisTuple,B<:BasisTuple,F<:FieldTuple{B,<:Tuple}} =
    FieldTuple(map((B,f)->B(f), tuple(B′.parameters[1].parameters...), f.fs))
(::Type{B′})(f::F) where {B′<:BasisTuple,B<:BasisTuple,Names,F<:FieldTuple{B,<:NamedTuple{Names}}} =
    FieldTuple(NamedTuple{Names}(map((B,f)->B(f), tuple(B′.parameters[1].parameters...), values(f.fs))))
(::Type{B′})(f::F) where {B′<:Basis,     B<:BasisTuple,F<:FieldTuple{B}} = FieldTuple(map(B′,f.fs))
(::Type{B′})(f::F) where {B′<:Basislike, B<:BasisTuple,F<:FieldTuple{B}} = FieldTuple(map(B′,f.fs))
# cases FieldTuple is in a concrete basis
(::Type{B′})(f::F) where {B′<:Basis,     B<:Basis,     F<:FieldTuple{B}} = FieldTuple(map(B′,f.fs))
(::Type{B′})(f::F) where {B′<:Basislike, B<:Basis,     F<:FieldTuple{B}} = B′(F)(f)

# in-place conversions
(::Type{B′})(f′::FieldTuple, f::FieldTuple) where {B′<:BasisTuple} =
    (tmap((B,f′,f)->B(f′,f), tuple(B′.parameters[1].parameters...), f′.fs, f.fs); f′)

### properties
getproperty(f::FieldTuple, s::Symbol) = getproperty(f, Val(s))
getproperty(f::FieldTuple, ::Val{:fs}) = getfield(f,:fs)
getproperty(f::FieldTuple, ::Val{s}) where {s} = getproperty(getfield(f,:fs),s)
propertynames(f::FieldTuple) = (:fs, propertynames(f.fs)...)

### simulation
white_noise(::Type{<:FieldTuple{B,FS}}) where {B,FS<:Tuple} =
    FieldTuple(map(white_noise, tuple(FS.parameters...)))
white_noise(::Type{<:FieldTuple{B,NamedTuple{Names,FS}}}) where {B,Names,FS<:Tuple} =
    FieldTuple(NamedTuple{Names}(map(white_noise, tuple(FS.parameters...))))

# generic AbstractVector inv/pinv don't work with FieldTuples because those
# implementations depends on get/setindex which we don't implement for FieldTuples
for func in [:inv, :pinv]
    @eval $(func)(D::DiagOp{FT}) where {FT<:FieldTuple} =
        Diagonal(FT(map(firstfield, map($(func), map(Diagonal,D.diag.fs)))))
end

# promote before recursing for these
≈(a::FieldTuple, b::FieldTuple) = all(map(≈, getfield.(promote(a,b),:fs)...))
dot(a::FieldTuple, b::FieldTuple) = sum(map(dot, getfield.(promote(a,b),:fs)...))
hash(ft::FieldTuple, h::UInt) = hash(ft.fs, h)

function ud_grade(f::FieldTuple, args...; kwargs...)
    FieldTuple(map(f->ud_grade(f,args...; kwargs...), f.fs))
end

### adjoint tuples

# represents a field which is adjoint over just the "tuple" indices. multiplying
# such a field by a non-adjointed one should be the inner product over just the
# tuple indices, and hence return a tuple-less, i.e a spin-0, field.
# note: these are really only lightly used in one place in LenseFlow, so they
# have almost no real functionality, the code here is in fact all there is.
struct TupleAdjoint{T<:Field}
    f :: T
end
tuple_adjoint(f::Field) = TupleAdjoint(f)

mul!(dst::Field{<:Any,S0}, a::TupleAdjoint{FT}, b::FT) where {FT<:Field{<:Any,S0}} = dst .= a.f .* b
mul!(dst::Field{<:Any,S0}, a::TupleAdjoint{FT}, b::FT) where {FT<:FieldTuple{<:Any,<:NamedTuple{<:Any,NTuple{2}}}} =
    (@. dst = a.f[1]*b[1] + a.f[2]*b[2])
# todo: make this generic case efficient:
mul!(dst::Field{<:Any,S0}, a::TupleAdjoint{FT}, b::FT) where {FT<:FieldTuple} =
    dst .= sum(map((a,b)->mul!(copy(dst),tuple_adjoint(a),b), a.f.fs, b.fs))


text:
using Dates: DateTime
using Test: @test, @testset

using BrokenRecord: HTTP, configure!, playback
using JSON3: JSON3

using Discord: Discord
const D = Discord

configure!(
    path=joinpath(@__DIR__, "fixtures"),
    ignore_headers=["Authorization", "User-Agent"],
)

@testset "Discord.jl" begin
    sets = ("snowflake", "objects", "routes")
    @testset "$set" for set in sets
        include("$set.jl")
    end
end