

   GGrriiddddeedd BBiivvaarriiaattee IInntteerrppoollaattiioonn ffoorr IIrrrreegguullaarr DDaattaa

        interp(x, y, z, xo=<<see below>>, yo=<<see below>>, ncp=0, extrap=F)
        interp.old(x, y, z, xo=<<see below>>, yo=<<see below>>, ncp=0, extrap=F)
        interp.new(x, y, z, xo=<<see below>>, yo=<<see below>>, ncp=0, extrap=F)

   AArrgguummeennttss::

          x: vector of x-coordinates of data points.  Missing
             values are not accepted.

          y: vector of y-coordinates of data points.  Missing
             values are not accepted.

          z: vector of z-coordinates of data points.  Missing
             values are not accepted.

             `x', `y', and `z' must be the same length  and may
             contain no fewer than four points. The points of
             `x' and `y' cannot be collinear, i.e, they cannot
             fall on the same line (two vectors `x' and `y'
             such that `y = ax + b' for some `a', `b' will not
             be accepted). `interp' is meant for cases in which
             you have `x', `y' values scattered over a plane
             and a `z' value for each.  If, instead, you are
             trying to evaluate a mathematical function, or get
             a graphical interpretation of relationships that
             can be described by a polynomial, try
             `outer()'.

         xo: vector of x-coordinates of output grid.  The
             default is 40 points evenly spaced over the range
             of `x'.  If extrapolation is not being used
             (`extrap=F', the default), `xo' should have a
             range that is close to or inside of the range of
             `x' for the results to be meaningful.

         yo: vector of y-coordinates of output grid.  The
             default is 40 points evenly spaced over the range
             of `y'.  If extrapolation is not being used
             (`extrap=F', the default), `yo' should have a
             range that is close to or inside of the range of
             `y' for the results to be meaningful.

        ncp: number of additional points to be used in comput-
             ing partial derivatives at each data point.  `ncp'
             must be either `0' (partial derivatives are not
             used), or at least 2 but smaller than the number
             of data points (and smaller than 25). This option
             is only supported by `interp.old'.

     extrap: logical flag: should extrapolation be used outside
             of the convex hull determined by the data points?

   duplicate: indicates how to handle duplicate data points.
             Possible values are `"error"' - produces an error
             message, `"strip"' - remove duplicate z values,
             `"mean"',`"median"',`"user"'  - calculate  mean ,
             median or user defined function of duplicate z
             values.

     dupfun: this function is applied to duplicate points if
             `duplicate="user"'

   DDeessccrriippttiioonn::

        If `ncp' is zero, linear interpolation is used in the
        triangles bounded by data points.  Cubic interpolation
        is done if partial derivatives are used.  If `extrap'
        is `FALSE', z-values for points outside the convex hull
        are returned as `NA'.  No extrapolation can be per-
        formed if `ncp' is zero.

        The `interp' function handles duplicate `(x,y)' points
        in different ways. As default it will stop with an
        error message. But it can give duplicate points an
        unique `z' value according to the parameter `duplicate'
        (`mean',`median' or any other user defined function).

        The triangulation scheme used by `interp' works well if
        `x' and `y' have similar scales but will appear
        stretched if they have very different scales.  The
        spreads of `x' and `y' must be within four orders of
        magnitude of each other for `interp' to work.

   VVaalluuee::

        list with 3 components:

          x: vector of x-coordinates of output grid, the same
             as the input argument `xo', if present.  Other-
             wise, a vector 40 points evenly spaced over the
             range of the input `x'.

          y: vector of y-coordinates of output grid, the same
             as the input argument `yo', if present.  Other-
             wise, a vector 40 points evenly spaced over the
             range of the input `x'.

          z: matrix of fitted z-values.  The value `z[i,j]' is
             computed at the x,y point `x[i], y[j]'. `z' has
             dimensions `length(x)' times `length(y)'
             (`length(xo)' times `length(yo)').

   NNoottee::

        `interp' is a wrapper for the two versions `interp.old'
        (it uses (almost) the same Fortran code from Akima 1978
        as the S-Plus version) and `interp.new' (it is based on
        new Fortran code from Akima 1996). For linear interpo-
        lation the old version is choosen, but spline interpo-
        lation is done by the new version.

        At the moment `interp.new' ignores `ncp' and does only
        bicubic spline interpolation.

        The resulting structure is suitable for input to the
        functions `contour' and `image'.  Check the require-
        ments of these functions when choosing values  for `xo'
        and `yo'.

   RReeffeerreenncceess::

        Akima, H. (1978). A Method of Bivariate Interpolation
        and Smooth Surface Fitting for Irregularly Distributed
        Data Points.  ACM Transactions on Mathematical Soft-
        ware, 4, 148-164.

        Akima, H. (1996). Algorithm 761: scattered-data surface
        fitting that has the accuracy of a cubic polynomial.
        ACM Transactions on Mathematical Software, 22, 362-371

   SSeeee AAllssoo::

        `contour', `image', `approx', `spline', `outer',
        `expand.grid'.

   EExxaammpplleess::

        data(akima)
        # linear interpolation
        akima.li <- interp(akima$x, akima$y, akima$z)
        image(akima.li$x,akima.li$y,akima.li$z)
        contour(akima.li$x,akima.li$y,akima.li$z,add=T)
        points(akima$x,akima$y)

        # increase smoothness
        akima.smooth <- interp(akima$x, akima$y, akima$z,
              xo=seq(0,25, length=100),  yo=seq(0,20, length=100))
        image(akima.smooth$x,akima.smooth$y,akima.smooth$z)
        contour(akima.smooth$x,akima.smooth$y,akima.smooth$z,add=T)
        points(akima$x,akima$y)
        # use triangulation library to
        # show underlying triangulation:
        library(tripack)
        plot(tri.mesh(akima),add=T,lty="dashed")

        # use only 15 points (interpolation only within convex hull!)
        akima.part <- interp(akima$x[1:15],akima$y[1:15],akima$z[1:15])
        image(akima.part$x,akima.part$y,akima.part$z)
        contour(akima.part$x,akima.part$y,akima.part$z,add=T)
        points(akima$x[1:15],akima$y[1:15])

        # spline interpolation, use 5 points to calculate derivatives
        akima.spl <- interp(akima$x, akima$y, akima$z,
              xo=seq(0,25, length=100),  yo=seq(0,20, length=100),ncp=5)
        image(akima.spl$x,akima.spl$y,akima.spl$z)
        contour(akima.spl$x,akima.spl$y,akima.spl$z,add=T)
        points(akima$x,akima$y)

        # example with duplicate points
        data(airquality)
        air <- airquality[(!is.na(airquality$Temp) &
                           !is.na(airquality$Ozone) &
                           !is.na(airquality$Solar.R)),]
        # gives an error:
        air.ip <- interp(air$Temp,air$Solar.R,air$Ozone)
        # use mean of duplicate points:
        air.ip <- interp(air$Temp,air$Solar.R,air$Ozone,duplicate="mean")

