diff --git a/DESCRIPTION b/DESCRIPTION index a6e9931eb34f04583b0f7ecb9552608f5af073d9..19f96b4aefcba263c838500605ca02aa55c05030 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,9 +1,12 @@ Package: akima -Version: 0.4-4 +Version: 0.4-5 +Date: 2005-07-21 Title: Interpolation of irregularly spaced data Author: Fortran code by H. Akima R port by Albrecht Gebhardt aspline function by Thomas Petzoldt + enhancements and corrections by Martin Maechler Maintainer: Albrecht Gebhardt Description: Linear or cubic spline interpolation for irregular gridded data License: Fortran code: ACM, free for non-commercial use, R functions GPL + diff --git a/R/interp.R b/R/interp.R index bb2e4714ea5cfcf39e73a3c0cca0044e7d196812..ce1c98349221a6a4a43ac207d44720ad5eb369c9 100644 --- a/R/interp.R +++ b/R/interp.R @@ -1,10 +1,15 @@ -"interp"<-function(x, y, z, xo = seq(min(x), max(x), length = 40), - yo = seq(min(y), max(y), length = 40), - ncp = 0, extrap = FALSE, duplicate = "error", dupfun = NULL) +interp <- + function(x, y, z, + xo = seq(min(x), max(x), length = 40), + yo = seq(min(y), max(y), length = 40), ncp = 0, + extrap = FALSE, duplicate = "error", dupfun = NULL) { - if (ncp==0) - # use the old version for linear interpolation - interp.old(x, y, z, xo, yo, ncp, extrap, duplicate, dupfun) - else - interp.new(x, y, z, xo, yo, ncp, extrap, duplicate, dupfun) + if (ncp == 0) + ## use the old version for linear interpolation + interp.old(x, y, z, xo = xo, yo = yo, ncp = ncp, + extrap = extrap, duplicate = duplicate, dupfun = dupfun) + else ## use the new one + interp.new(x, y, z, xo = xo, yo = yo, linear = FALSE, + ncp = NULL,# not using 'ncp' argument + extrap = extrap, duplicate = duplicate, dupfun = dupfun) } diff --git a/R/interp.new.R b/R/interp.new.R index edd599abe20060f8db9e0e4d4d0c42345ce1a817..949bc9fa6f0d4e762d69855634c444014fdd07cf 100644 --- a/R/interp.new.R +++ b/R/interp.new.R @@ -1,12 +1,13 @@ -"interp.new"<-function(x, y, z, xo = seq(min(x), max(x), length = 40), - yo = seq(min(y), max(y), length = 40), linear=FALSE, - ncp = NULL, extrap = FALSE, duplicate = "error", - dupfun = NULL) +interp.new <- + function(x, y, z, + xo = seq(min(x), max(x), length = 40), + yo = seq(min(y), max(y), length = 40), linear = FALSE, + ncp = NULL, extrap = FALSE, duplicate = "error", dupfun = NULL) { if(!(all(is.finite(x)) && all(is.finite(y)) && all(is.finite(z)))) stop("missing values and Infs not allowed") - if(!is.null(ncp)){ - if(ncp!=0){ + if(!is.null(ncp)) { + if(ncp != 0) { cat("ncp not supported, it is automatically choosen by Fortran code\n") } else { @@ -14,15 +15,15 @@ stop("use interp.old().") } } - if(linear){ - cat("linear interpolation not yet implemented with interp.new().\n") - stop("use interp.old().") + if(linear) { + cat("linear interpolation not yet implemented with interp.new().\n") + stop("use interp.old().") } drx <- diff(range(x)) dry <- diff(range(y)) if(drx == 0 || dry == 0) - stop("all data collinear") # other cases caught in Fortran code + stop("all data collinear") # other cases caught in Fortran code if(drx/dry > 10000 || drx/dry < 0.0001) stop("scales of x and y are too dissimilar") n <- length(x) @@ -30,48 +31,42 @@ ny <- length(yo) if(length(y) != n || length(z) != n) stop("Lengths of x, y, and z do not match") - xy <- paste(x, y, sep =",") - i <- match(xy, xy) - if(duplicate=="user" && !is.function(dupfun)) - stop("duplicate=\"user\" requires dupfun to be set to a function") - if(duplicate!="error") - { - centre <- function(x) { + + xy <- paste(x, y, sep = ",")# trick for 'duplicated' (x,y)-pairs + if(duplicate == "error") { + if(any(duplicated(xy))) + stop("duplicate data points: need to set 'duplicate = ..' ") + } + else { ## duplicate != "error" + + i <- match(xy, xy) + if(duplicate == "user") + dupfun <- match.fun(dupfun)#> error if it fails + + ord <- !duplicated(xy) + if(duplicate != "strip") { + centre <- function(x) switch(duplicate, mean = mean(x), median = median(x), user = dupfun(x)) - } - if(duplicate!="strip"){ - z <- unlist(lapply(split(z,i), centre)) - ord <- !duplicated(xy) - x <- x[ord] - y <- y[ord] - n <- length(x) - } - else{ - ord <- (hist(i,plot=FALSE,freq=TRUE,breaks=seq(0.5,max(i)+0.5,1))$counts==1) - x <- x[ord] - y <- y[ord] - z <- z[ord] - n <- length(x) - } + z <- unlist(lapply(split(z,i), centre)) + } else { + z <- z[ord] } - else - if(any(duplicated(xy))) - stop("duplicate data points") + x <- x[ord] + y <- y[ord] + n <- length(x) + } + zo <- matrix(0, nx, ny) storage.mode(zo) <- "double" - miss <- !extrap #if not extrapolating use missing values - extrap <- matrix(TRUE, nx, ny) - if(!is.null(ncp)){ - if(extrap & ncp == 0) - warning("Cannot extrapolate with linear option") - } - else { - if(extrap & linear) + miss <- !extrap # if not extrapolating, set missing values + misso <- matrix(TRUE, nx, ny)# hmm, or rather 'miss' ?? + + if(extrap && if(is.null(ncp)) linear else (ncp == 0)) warning("Cannot extrapolate with linear option") - } + ans <- .Fortran("sdsf3p", as.integer(1), as.integer(n), @@ -86,13 +81,13 @@ ier = integer(1), double(36 * n), integer(25 * n), - extrap = as.logical(extrap), + extrap = as.logical(misso), near = integer(n), nxt = integer(n), dist = double(n), - PACKAGE = "akima") - temp <- ans[c("x", "y", "z", "extrap")] + PACKAGE = "akima")[c("x", "y", "z", "extrap")] if(miss) - temp$z[temp$extrap]<-NA - temp[c("x", "y", "z")] + ans$z[ans$extrap] <- NA + + ans[c("x", "y", "z")] } diff --git a/R/interp.old.R b/R/interp.old.R index 61e7953c245e4c4497ceb38bd1d62a1463d7ab15..433c304631bcecb8fc5b67652e365ba3b5b8b9ef 100644 --- a/R/interp.old.R +++ b/R/interp.old.R @@ -1,13 +1,14 @@ "interp.old"<-function(x, y, z, xo = seq(min(x), max(x), length = 40), - yo = seq(min(y), max(y), length = 40), + yo = seq(min(y), max(y), length = 40), ncp = 0, extrap = FALSE, duplicate = "error", dupfun = NULL) { if(!(all(is.finite(x)) && all(is.finite(y)) && all(is.finite(z)))) stop("missing values and Infs not allowed") - if(ncp>25){ + if(ncp > 25) { ncp <- 25 cat("ncp too large, using ncp=25\n") } + drx <- diff(range(x)) dry <- diff(range(y)) if(drx == 0 || dry == 0) @@ -19,36 +20,34 @@ ny <- length(yo) if(length(y) != n || length(z) != n) stop("Lengths of x, y, and z do not match") - xy <- paste(x, y, sep =",") - i <- match(xy, xy) - if(duplicate=="user" && !is.function(dupfun)) - stop("duplicate=\"user\" requires dupfun to be set to a function") - if(duplicate!="error") - { - centre <- function(x) { + + xy <- paste(x, y, sep = ",")# trick for 'duplicated' (x,y)-pairs + if(duplicate == "error") { + if(any(duplicated(xy))) + stop("duplicate data points: need to set 'duplicate = ..' ") + } + else { ## duplicate != "error" + + i <- match(xy, xy) + if(duplicate == "user") + dupfun <- match.fun(dupfun)#> error if it fails + + ord <- !duplicated(xy) + if(duplicate != "strip") { + centre <- function(x) switch(duplicate, mean = mean(x), median = median(x), user = dupfun(x)) - } - if(duplicate!="strip"){ - z <- unlist(lapply(split(z,i), centre)) - ord <- !duplicated(xy) - x <- x[ord] - y <- y[ord] - n <- length(x) - } - else{ - ord <- (hist(i,plot=FALSE,freq=TRUE,breaks=seq(0.5,max(i)+0.5,1))$counts==1) - x <- x[ord] - y <- y[ord] - z <- z[ord] - n <- length(x) - } + z <- unlist(lapply(split(z,i), centre)) + } else { + z <- z[ord] } - else - if(any(duplicated(xy))) - stop("duplicate data points") + x <- x[ord] + y <- y[ord] + n <- length(x) + } + zo <- matrix(0, nx, ny) storage.mode(zo) <- "double" miss <- !extrap #if not extrapolating use missing values @@ -70,8 +69,7 @@ integer((31 + ncp) * n + nx * ny), double(5 * n), misso = as.logical(misso), - PACKAGE = "akima") - temp <- ans[c("x", "y", "z", "misso")] - temp$z[temp$misso]<-NA - temp[c("x", "y", "z")] + PACKAGE = "akima")[c("x", "y", "z", "misso")] + ans$z[ans$misso] <- NA + ans[c("x", "y", "z")] } diff --git a/man/interp.Rd b/man/interp.Rd index abfca20489356f57d8e482f297e8e1b7ec2486e9..f6ad0d2765c5f12db6486c0e9057910e6dc167d5 100644 --- a/man/interp.Rd +++ b/man/interp.Rd @@ -1,14 +1,24 @@ \name{interp} -\title{ - Gridded Bivariate Interpolation for Irregular Data -} +\title{Gridded Bivariate Interpolation for Irregular Data} \alias{interp} \alias{interp.new} \alias{interp.old} +\description{ + These functions implement bivariate interpolation onto a grid + for irregularly spaced input data. Bilinear or bicubic spline + interpolation is applied using different versions of algorithms from + Akima. +} \usage{ -interp(x, y, z, xo=seq(min(x), max(x), length = 40), yo=seq(min(y), max(y), length = 40), ncp=0, extrap=FALSE, duplicate = "error", dupfun = NULL) -interp.old(x, y, z, xo= seq(min(x), max(x), length = 40), yo=seq(min(y), max(y), length = 40), ncp=0, extrap=FALSE, duplicate = "error", dupfun = NULL) -interp.new(x, y, z, xo= seq(min(x), max(x), length = 40), yo=seq(min(y), max(y), length = 40), linear=FALSE, ncp=NULL, extrap=FALSE, duplicate = "error", dupfun = NULL) +interp(x, y, z, xo=seq(min(x), max(x), length = 40), + yo=seq(min(y), max(y), length = 40), + ncp = 0, extrap=FALSE, duplicate = "error", dupfun = NULL) +interp.old(x, y, z, xo= seq(min(x), max(x), length = 40), + yo=seq(min(y), max(y), length = 40), ncp = 0, + extrap=FALSE, duplicate = "error", dupfun = NULL) +interp.new(x, y, z, xo = seq(min(x), max(x), length = 40), + yo = seq(min(y), max(y), length = 40), linear = FALSE, + ncp = NULL, extrap=FALSE, duplicate = "error", dupfun = NULL) } \arguments{ \item{x}{ @@ -22,72 +32,62 @@ interp.new(x, y, z, xo= seq(min(x), max(x), length = 40), yo=seq(min(y), max(y), \item{z}{ vector of z-coordinates of data points. Missing values are not accepted. - - - \code{x}, \code{y}, and \code{z} must be the same length and may contain no fewer - than four points. The points of \code{x} and \code{y} - cannot be collinear, i.e, they cannot fall on the same line (two vectors - \code{x} and \code{y} such that \code{y = ax + b} for some \code{a}, \code{b} will not be - accepted). \code{interp} is meant for cases in which you have \code{x}, \code{y} - values scattered over a plane and a \code{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 \code{outer()}. + + \code{x}, \code{y}, and \code{z} must be the same length and may + contain no fewer than four points. The points of \code{x} and + \code{y} cannot be collinear, i.e, they cannot fall on the same line + (two vectors \code{x} and \code{y} such that \code{y = ax + b} for + some \code{a}, \code{b} will not be accepted). \code{interp} is + meant for cases in which you have \code{x}, \code{y} values + scattered over a plane and a \code{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 \code{outer()}. } \item{xo}{ vector of x-coordinates of output grid. The default is 40 points evenly spaced over the range of \code{x}. If extrapolation is not being - used (\code{extrap=FALSE}, the default), \code{xo} should have a range that is - close to or inside of the range of \code{x} for the results to be meaningful. - } - \item{yo}{ - vector of y-coordinates of output grid. The default is 40 points - evenly spaced over the range of \code{y}. If extrapolation is not being - used (\code{extrap=FALSE}, the default), \code{yo} should have a range that is - close to or inside of the range of \code{y} for the results to be meaningful. + used (\code{extrap=FALSE}, the default), \code{xo} should have a + range that is close to or inside of the range of \code{x} for the + results to be meaningful. } - \item{linear}{logical, switch to linear interpolation in \code{interp.new}} + \item{yo}{vector of y-coordinates of output grid; analogous to + \code{xo}, see above.} + \item{linear}{logical, switch to linear interpolation in \code{interp.new}.} \item{ncp}{ number of additional points to be used in computing partial - derivatives at each data point. + derivatives at each data point. \code{ncp} must be either \code{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 \code{interp.old}. } \item{extrap}{ logical flag: should extrapolation be used outside of the - convex hull determined by the data points? -} -\item{duplicate}{ - indicates how to handle duplicate data points. Possible values are - \code{"error"} - produces an error message, \code{"strip"} - remove - duplicate z values, \code{"mean"},\code{"median"},\code{"user"} - - calculate mean , median or user defined function of duplicate z - values.} -\item{dupfun}{this function is applied to duplicate points if \code{duplicate="user"} -} - + convex hull determined by the data points?} + \item{duplicate}{character string indicating how to handle duplicate + data points. Possible values are + \describe{ + \item{\code{"error"}}{produces an error message,} + \item{\code{"strip"}}{remove duplicate z values,} + \item{ \code{"mean"},\code{"median"},\code{"user"}}{calculate + mean , median or user defined function (\code{dupfun}) of duplicate + z values.} + }} + \item{dupfun}{a function, applied to duplicate points if + \code{duplicate= "user"}.} } \value{ list with 3 components: - - \item{x}{ - vector of x-coordinates of output grid, the same as the input - argument \code{xo}, if present. Otherwise, a vector 40 points evenly spaced - over the range of the input \code{x}. - } - \item{y}{ - vector of y-coordinates of output grid, the same as the input - argument \code{yo}, if present. Otherwise, a vector 40 points evenly spaced - over the range of the input \code{x}. - } + \item{x,y}{ + vectors of x- and y- coordinates of output grid, the same as the input + argument \code{xo}, or \code{yo}, if present. Otherwise, their + default, a vector 40 points evenly spaced over the range of the + input \code{x}.} \item{z}{ matrix of fitted z-values. The value \code{z[i,j]} is computed - at the x,y point \code{x[i], y[j]}. \code{z} has - dimensions \code{length(x)} times \code{length(y)} (\code{length(xo)} times \code{length(yo)}). - - -}} + at the x,y point \code{xo[i], yo[j]}. \code{z} has + dimensions \code{length(xo)} times \code{length(yo)}.} +} \note{ \code{interp} is a wrapper for the two versions \code{interp.old} (it uses (almost) the same Fortran code from Akima 1978 as the S-Plus version) and @@ -95,7 +95,7 @@ interp.new(x, y, z, xo= seq(min(x), max(x), length = 40), yo=seq(min(y), max(y), (it is based on new Fortran code from Akima 1996). For linear interpolation the old version is choosen, but spline interpolation is done by the new version. - + At the moment \code{interp.new} ignores \code{ncp} and does only bicubic spline interpolation. @@ -103,18 +103,17 @@ interp.new(x, y, z, xo= seq(min(x), max(x), length = 40), yo=seq(min(y), max(y), functions \code{contour} and \code{image}. Check the requirements of these functions when choosing values for \code{xo} and \code{yo}. } -\description{ +\details{ If \code{ncp} is zero, linear - interpolation is used in the triangles bounded by data points. + interpolation is used in the triangles bounded by data points. Cubic interpolation is done if partial derivatives are used. If \code{extrap} is \code{FALSE}, z-values for points outside the - convex hull are returned as \code{NA}. + convex hull are returned as \code{NA}. No extrapolation can be performed if \code{ncp} is zero. - - - The \code{interp} function handles duplicate \code{(x,y)} points + + The \code{interp} function handles duplicate \code{(x,y)} points in different ways. As default it will stop with an error message. But - it can give duplicate points an unique \code{z} value according to the + it can give duplicate points an unique \code{z} value according to the parameter \code{duplicate} (\code{mean},\code{median} or any other user defined function). @@ -137,50 +136,75 @@ interp.new(x, y, z, xo= seq(min(x), max(x), length = 40), yo=seq(min(y), max(y), } \seealso{ - \code{\link[base]{contour}}, \code{\link[base]{image}}, \code{\link[base]{approx}}, \code{\link[base]{spline}}, \code{\link[base]{outer}}, \code{\link[base]{expand.grid}}. + \code{\link{contour}}, \code{\link{image}}, + \code{\link{approx}}, \code{\link{spline}}, + \code{\link{outer}}, \code{\link{expand.grid}}. } \examples{ data(akima) -# linear interpolation +plot(y ~ x, data = akima, main = "akima example data") +with(akima, text(x, y, formatC(z,dig=2), adj = -0.1)) + +## 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=TRUE) -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=TRUE) -points(akima$x,akima$y) -# use triangulation library to -# show underlying triangulation: +image (akima.li, add=TRUE) +contour(akima.li, add=TRUE) +points (akima, pch = 3) + +## increase smoothness (using finer grid): +akima.smooth <- + with(akima, interp(x, y, z, xo=seq(0,25, length=100), + yo=seq(0,20, length=100))) +image (akima.smooth) +contour(akima.smooth, add=TRUE) +points(akima, pch = 3, cex = 2, col = "blue") +# use triangulation package to show underlying triangulation: if(library(tripack, logical.return=TRUE)) - plot(tri.mesh(akima),add=TRUE,lty="dashed") + plot(tri.mesh(akima), add=TRUE, 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=TRUE) +akima.part <- with(akima, interp(x[1:15], y[1:15], z[1:15])) +image (akima.part) +contour(akima.part, add=TRUE) points(akima$x[1:15],akima$y[1:15]) -# spline interpolation, use 5 points to calculate derivatives -# interp gives `linear interpolation not yet implemented with interp.new()' -akima.spl <- interp.old(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=TRUE) -points(akima$x,akima$y) +## spline interpolation +## -------------------- +## "Old": use 5 points to calculate derivatives -> many NAs +akima.sO <- interp.old(akima$x, akima$y, akima$z, + xo=seq(0,25, length=100), yo=seq(0,20, length=100), ncp=5) +table(is.na(akima.sO$z)) ## 3990 NA's; = 40 \% +akima.sO <- with(akima, + interp.old(x,y,z, xo=seq(0,25, length=100), yo=seq(0,20, len=100), ncp = 4)) +sum(is.na(akima.sO$z)) ## still 3429 +image (akima.sO) # almost useless +contour(akima.sO, add = TRUE) + +## "New:" +akima.spl <- with(akima, interp.new(x,y,z, xo=seq(0,25, length=100), + yo=seq(0,20, length=100))) +contour(akima.spl) ; points(akima) + +full.pal <- function(n) hcl(h = seq(340, 20, length = n)) +cool.pal <- function(n) hcl(h = seq(120, 0, length = n) + 150) +warm.pal <- function(n) hcl(h = seq(120, 0, length = n) - 30) + +filled.contour(akima.spl, color.palette = full.pal, + plot.axes = { axis(1); axis(2); + points(akima, pch = 3, col= hcl(c=100, l = 20))}) +# no extrapolation! + +## example with duplicate points : -# example with duplicate points data(airquality) -air <- airquality[(!is.na(airquality$Temp) & - !is.na(airquality$Ozone) & - !is.na(airquality$Solar.R)),] -# gives an error: -\dontrun{air.ip <- interp(air$Temp,air$Solar.R,air$Ozone)} +air <- subset(airquality, + !is.na(Temp) & !is.na(Ozone) & !is.na(Solar.R)) +# gives an error {duplicate ..}: +try( air.ip <- interp.new(air$Temp,air$Solar.R,air$Ozone) ) # use mean of duplicate points: -air.ip <- interp(air$Temp,air$Solar.R,air$Ozone,duplicate="mean") +air.ip <- with(air, interp.new(Temp, Solar.R, log(Ozone), duplicate = "mean")) +image(air.ip, main = "Airquality: Ozone vs. Temp and Solar.R") +with(air, points(Temp, Solar.R)) } \keyword{dplot} -% Converted by Sd2Rd version 0.2-a3. + diff --git a/orig/akima.bib b/orig/akima.bib index 6358248813111e061aaefc10e9564e377c36bbd9..eb3bca69a85a16f6ef4156c1f0c39aa0d777a3f8 100644 --- a/orig/akima.bib +++ b/orig/akima.bib @@ -1,10 +1,3 @@ - -TOMS : Bibliographic record for "Akima:1996:ASS" - -

-Bibliographic record for "Akima:1996:ASS"

-

-
 @Article{Akima:1996:ASS,
   author =       "Hiroshi Akima",
   title =        "Algorithm 761: scattered-data surface fitting that has
@@ -30,8 +23,3 @@ Bibliographic record for "Akima:1996:ASS"
                  G.4}: Mathematics of Computing, MATHEMATICAL
                  SOFTWARE.",
 }
-
-
-

- - diff --git a/src/idptip.f b/src/idptip.f index 5ffc09c2075fe1cf72e73df0d12c480e2b229d78..3013b84f306290e879180f4eb8a220d582345042 100644 --- a/src/idptip.f +++ b/src/idptip.f @@ -3,6 +3,7 @@ 1 ZII,MISSII) C THIS SUBROUTINE PERFORMS PUNCTUAL INTERPOLATION OR EXTRAPOLA- C TION, I.E., DETERMINES THE Z VALUE AT A POINT. + C THE INPUT PARAMETERS ARE C XD,YD,ZD = ARRAYS OF DIMENSION NDP CONTAINING THE X, C Y, AND Z COORDINATES OF THE DATA POINTS, WHERE @@ -22,9 +23,11 @@ C THE POINT FOR WHICH INTERPOLATION IS TO BE C PERFORMED, C XII,YII = X AND Y COORDINATES OF THE POINT FOR WHICH C INTERPOLATION IS TO BE PERFORMED. -C THE OUTPUT PARAMETER IS -C ZII = INTERPOLATED Z VALUE. + +C THE OUTPUT PARAMETERs are +C ZII = INTERPOLATED Z VALUE. C MISSII = LOCIGAL INDICATING MISSING VALUE + C DECLARATION STATEMENTS IMPLICIT DOUBLE PRECISION (A-D,P-Z) LOGICAL MISSII @@ -35,6 +38,7 @@ C DECLARATION STATEMENTS 1 ZU(3),ZV(3),ZUU(3),ZUV(3),ZVV(3) DOUBLE PRECISION LU,LV EQUIVALENCE (P5,P50) + C PRELIMINARY PROCESSING 10 IT0=ITI NTL=NT+NL @@ -170,6 +174,7 @@ C EVALUATES THE POLYNOMIAL. ZII=P0+U*(P1+U*(P2+U*(P3+U*(P4+U*P5)))) MISSII=.FALSE. RETURN + C CALCULATION OF ZII BY EXTRAPOLATION IN THE RECTANGLE. C CHECKS IF THE NECESSARY COEFFICIENTS HAVE BEEN CALCULATED. 40 IF(IT0.EQ.ITPV) GO TO 50