Skip to content

GitLab

Commit 43507cd4 authored by agebhard's avatar agebhard
Browse files

upload

parent ca58ed3a
Show whitespace changes
Inline Side-by-side
DESCRIPTION
View file @ 43507cd4
Package: akima
Version: 0.5-11
Date: 2013-01-19
Title: Interpolation of irregularly spaced data
Version: 0.5-12
Date: 2015-08-26
Title: Interpolation of Irregularly and Regularly Spaced Data
Author: Fortran code by H. Akima
R port by Albrecht Gebhardt <albrecht.gebhardt@uni-klu.ac.at>
R port by Albrecht Gebhardt <albrecht.gebhardt@aau.at>
aspline function by Thomas Petzoldt <thomas.petzoldt@tu-dresden.de>
interp2xyz, enhancements and corrections by Martin Maechler <maechler@stat.math.ethz.ch>
Maintainer: Albrecht Gebhardt <albrecht.gebhardt@uni-klu.ac.at>
Description: Linear or cubic spline interpolation for irregular gridded data
Description: Several cubic spline interpolation methods of H. Akima for irregular and
regular gridded data are available through this package, both for the bivariate case
(irregular data: ACM 526 and ACM 761, regular data: ACM 433) and univariate case (ACM 697).
Linear interpolation of irregular gridded data is also covered by reusing D. J. Renkas
triangulation code which is part of Akimas Fortran code.
License: ACM | file LICENSE
Depends: R (>= 2.0.0)
Imports: sp
NAMESPACE
View file @ 43507cd4
useDynLib(akima)
export(aspline,interp,interpp,interp.old,interp.new,interp2xyz,bicubic,bicubic.grid)
importFrom("sp","coordinates")
importFrom("sp","coordinates<-")
importFrom("sp","gridded<-")
importFrom("grDevices", "xy.coords")
importFrom("graphics", "hist")
importFrom("stats", "median")
PORTING deleted 100644 → 0
View file @ ca58ed3a
This library contains an R implementation of the S-Plus function interp().
See ChangeLog for details about further versions
Version 0.1-1
I used the S-Plus version of interp() (S-Plus 3.3 f. Win 3.11) as starting
point. Then I searched for Akimas idsfft function, which is referenced in
interp().
I found the appropriate Fortran code (1978) in the ACM Collected
Algorithms archive under
http://www.netlib.org/toms/526
(also included here) and splitted it into the Fortran files under src/ .
The test driver ttidbs is also included and can be compiled with "make test".
However, it seems that this code differs a little bit from that used in S-Plus:
It implements "only" bicubic spline interpolation, no linear interpolation as
in S-Plus.
So I modified IDSFFT and added a subroutine IDPTLI which does linear
interpolation within the triangles, generated by IDTANG.
Further changes are
+ REAL -> DOUBLE PRECISION
+ static DIMENSIONs replaced with dynamic
+ option to toggle extrapolation outside of convex hull added in IDSFFT and
IDBVIP.
Because I don't know how to generate NAs in Fortran (I use g77 on Linux),
I added a logical array MISSI, that indicates if a returned value should
be NA. These values will be set to NA after the Fortran call.
+ option to handle duplicate data points added (according to an example in the
S-Plus help page)
+ man pages converted and rewritten
+ data set akima (from S-Plus) added
+ function interpp() added, it evaluates the interpolated function at
arbitraryly choosen points and generates no regular grid as interp() does.
There where some problems with interpp() when using the Fortran version:
- it crashes when compiled with "g77 -O2 -fpic" and called with more than
one output point.
- compilation with "g77 -g -fpic" fails (see src/Makefile for details)
- compilation with "g77 -g" works (and no crashes occur)
These problems do not occur in the C Version (generated by f2c), so it would
be better to use only the src-c tree.
After I finished the above steps I found a more recent version (ACM 761, 1996)
of Akimas interpolation code which uses the tripack package (also available at
ACM as algorithm no. 751) for triangulation and now I'm trying to use it for
the next version of interp().
------------------------------------------------------------------
Albrecht Gebhardt email: albrecht.gebhardt@uni-klu.ac.at
Institut fuer Mathematik Tel. : (++43 463) 2700/837
Universitaet Klagenfurt Fax : (++43 463) 2700/834
Villacher Str. 161
A-9020 Klagenfurt, Austria
------------------------------------------------------------------
R/interp.R
View file @ 43507cd4
interp <-
function(x, y, z,
xo = seq(min(x), max(x), length = 40),
yo = seq(min(y), max(y), length = 40), linear=TRUE,
extrap = FALSE, duplicate = "error", dupfun = NULL, ncp=NULL)
function(x, y=NULL, z,
xo = seq(min(x), max(x), length = nx),
yo = seq(min(y), max(y), length = ny), linear=TRUE,
extrap = FALSE, duplicate = "error", dupfun = NULL, ncp=NULL,
nx=40, ny=40)
{
# for backward compatibility
## for backward compatibility
if(!is.null(ncp)){
warning('use of \'ncp\' parameter is deprecated!')
if(ncp==0)
......@@ -15,12 +16,49 @@ interp <-
else
stop('ncp < 0 ?')
}
## handle sp data, save coordinate and value names
is.sp <- FALSE
sp.coord <- NULL
sp.z <- NULL
sp.proj4string <- NULL
if(is.null(y)&&is.character(z)){
if(class(x)=="SpatialPointsDataFrame"){
sp.coord <- dimnames(coordinates(x))[[2]]
sp.z <- z
sp.proj4string <- x@proj4string
z <- x@data[,z]
y <- coordinates(x)[,2]
x <- coordinates(x)[,1]
is.sp <- TRUE
xo = seq(min(x), max(x), length = nx)
yo = seq(min(y), max(y), length = ny)
} else
stop("either x,y,z are numerical or x is SpatialPointsDataFrame and z a name of a data column in x")
}
if(linear)
## use the old version for linear interpolation
interp.old(x, y, z, xo = xo, yo = yo, ncp = 0,
extrap = extrap, duplicate = duplicate, dupfun = dupfun)
ret <- interp.old(x, y, z, xo = xo, yo = yo, ncp = 0,
extrap = extrap, duplicate = duplicate,
dupfun = dupfun)
else ## use the new one
interp.new(x, y, z, xo = xo, yo = yo, linear = FALSE,
ret <- interp.new(x, y, z, xo = xo, yo = yo, linear = FALSE,
ncp = NULL,# not using 'ncp' argument
extrap = extrap, duplicate = duplicate, dupfun = dupfun)
extrap = extrap, duplicate = duplicate,
dupfun = dupfun)
if(is.sp){
zm <- dim(ret$z)[1]
zn <- dim(ret$z)[2]
zvec <- c(ret$z)
xvec <- c(matrix(rep(ret$x,zn),nrow=zm,ncol=zn,byrow=FALSE))
yvec <- c(matrix(rep(ret$y,zm),nrow=zm,ncol=zn,byrow=TRUE))
nona <- !is.na(zvec)
ret <- data.frame(xvec[nona],yvec[nona],zvec[nona])
names(ret) <- c(sp.coord[1],sp.coord[2],sp.z)
coordinates(ret) <- sp.coord
ret@proj4string <- sp.proj4string
gridded(ret) <- TRUE
}
ret
}
R/interpp.R
View file @ 43507cd4
"interpp"<-function(x, y, z, xo, yo, linear = TRUE, extrap = FALSE,
"interpp"<-function(x, y=NULL, z, xo, yo=NULL, linear = TRUE, extrap = FALSE,
duplicate = "error", dupfun = NULL, ncp=NULL)
{
......@@ -12,11 +12,43 @@
else
stop('ncp < 0 ?')
}
## handle sp data, save coordinate and value names
is.sp <- FALSE
sp.coord <- NULL
sp.z <- NULL
sp.proj4string <- NULL
if(is.null(y)&&is.character(z)){
if(class(xo)=="SpatialPointsDataFrame"){
yo <- coordinates(xo)[,2]
xo <- coordinates(xo)[,1]
} else
stop("either x,y,z,xo,yo have to be numeric vectors or
both x and xo have to be SpatialPointsDataFrames and z a name of a data column in x")
if(class(x)=="SpatialPointsDataFrame"){
sp.coord <- dimnames(coordinates(x))[[2]]
sp.z <- z
sp.proj4string <- x@proj4string
z <- x@data[,z]
y <- coordinates(x)[,2]
x <- coordinates(x)[,1]
is.sp <- TRUE
} else
stop("either x,y,z,xo,yo have to be numeric vectors or
both x and xo have to be SpatialPointsDataFrames and z a name of a data column in x")
}
if(linear)
# the old Akima code:
interpp.old(x, y, z, xo, yo, ncp=0, extrap, duplicate, dupfun)
## the old Akima code:
ret <- interpp.old(x, y, z, xo, yo, ncp=0, extrap, duplicate, dupfun)
else
# new code for splines
interpp.new(x, y, z, xo, yo, extrap, duplicate, dupfun)
## new code for splines
ret <- interpp.new(x, y, z, xo, yo, extrap, duplicate, dupfun)
if(is.sp){
nona <- !is.na(ret$z)
ret <- data.frame(ret$x[nona],ret$y[nona],ret$z[nona])
names(ret) <- c(sp.coord[1],sp.coord[2],sp.z)
coordinates(ret) <- sp.coord
ret@proj4string <- sp.proj4string
}
ret
}
README
View file @ 43507cd4
......@@ -21,11 +21,83 @@ Thomas Petzoldt <petzoldt@rcs.urz.tu-dresden.de> provided an R interface to
Akimas univariate spline interpolation functions (aspline) which has been
included into this library.
The original papers from H. Akima are available (access is not free, depends on your
ACM subscription):
TOMS 433: http://dl.acm.org/citation.cfm?id=355604.355605
TOMS 526: http://dl.acm.org/citation.cfm?id=355780.355787
TOMS 697: http://dl.acm.org/citation.cfm?id=114697.116808
TOMS 761: http://dl.acm.org/citation.cfm?id=232826.232856
Original Fortran code is freely available at Netlib/TOMS -- Collected Algorithms
of the ACM at:
http://www.netlib.org/toms/433
http://www.netlib.org/toms/526
http://www.netlib.org/toms/697
http://www.netlib.org/toms/761
------------------------------------------------------------------------------
Albrecht Gebhardt email: albrecht.gebhardt@uni-klu.ac.at
Albrecht Gebhardt email: albrecht.gebhardt@aau.at
Institut fuer Mathematik Tel. : (++43 463) 2700/3118
Universitaet Klagenfurt Fax : (++43 463) 2700/3198
Universitaetsstr. 65-67 WWW : http://www.math.uni-klu.ac.at/~agebhard
Universitaetsstr. 65-67 WWW : http://www.stat.aau.at/~agebhard
A-9020 Klagenfurt, Austria
------------------------------------------------------------------------------
Details:
See ChangeLog for details about further versions
Version 0.1-1
I used the S-Plus version of interp() (S-Plus 3.3 f. Win 3.11) as starting
point. Then I searched for Akimas idsfft function, which is referenced in
interp().
I found the appropriate Fortran code (1978) in the ACM Collected
Algorithms archive under
http://www.netlib.org/toms/526
(also included here) and splitted it into the Fortran files under src/ .
The test driver ttidbs was also included in older versions and could be
compiled with "make test", now it is removed here.
However, it seems that this code differs a little bit from that used in S-Plus:
It implements "only" bicubic spline interpolation, no linear interpolation as
in S-Plus.
So I modified IDSFFT and added a subroutine IDPTLI which does linear
interpolation within the triangles, generated by IDTANG.
Further changes are
+ REAL -> DOUBLE PRECISION
+ static DIMENSIONs replaced with dynamic
+ option to toggle extrapolation outside of convex hull added in IDSFFT and
IDBVIP.
Because I don't know how to generate NAs in Fortran (I use g77 on Linux),
I added a logical array MISSI, that indicates if a returned value should
be NA. These values will be set to NA after the Fortran call.
+ option to handle duplicate data points added (according to an example in the
S-Plus help page)
+ man pages converted and rewritten
+ data set akima (from S-Plus) added
+ function interpp() added, it evaluates the interpolated function at
arbitraryly choosen points and generates no regular grid as interp() does.
There where some problems with interpp() when using the Fortran version:
- it crashes when compiled with "g77 -O2 -fpic" and called with more than
one output point.
- compilation with "g77 -g -fpic" fails (see src/Makefile for details)
- compilation with "g77 -g" works (and no crashes occur)
These problems do not occur in the C Version (generated by f2c), so it would
be better to use only the src-c tree.
After I finished the above steps I found a more recent version (ACM 761, 1996)
of Akimas interpolation code which uses the tripack package (also available at
ACM as algorithm no. 751) for triangulation and now I'm trying to use it for
the next version of interp().
TITLE deleted 100644 → 0
View file @ ca58ed3a
akima Interpolation of irregularly spaced data
cleanup
View file @ 43507cd4
......@@ -4,4 +4,5 @@ rm -f src/Makevars
rm -f src/*.o
rm -f src/*.so
rm -f src/*.dll
rm -f config.*
man/interp.Rd
View file @ 43507cd4
......@@ -10,9 +10,10 @@
Akima.
}
\usage{
interp(x, y, z, xo=seq(min(x), max(x), length = 40),
yo=seq(min(y), max(y), length = 40),
linear = TRUE, extrap=FALSE, duplicate = "error", dupfun = NULL, ncp = NULL)
interp(x, y=NULL, z, xo=seq(min(x), max(x), length = nx),
yo=seq(min(y), max(y), length = ny),
linear = TRUE, extrap=FALSE, duplicate = "error", dupfun = NULL,
ncp = NULL, nx = 40, ny = 40)
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)
......@@ -22,18 +23,27 @@ interp.new(x, y, z, xo = seq(min(x), max(x), length = 40),
}
\arguments{
\item{x}{
vector of x-coordinates of data points.
vector of x-coordinates of data points or a
\code{SpatialPointsDataFrame} object.
Missing values are not accepted.
}
\item{y}{
vector of y-coordinates of data points.
Missing values are not accepted.
If left as NULL indicates that \code{x} should be a
\code{SpatialPointsDataFrame} and \code{z} names the variable of
interest in this dataframe.
}
\item{z}{
vector of z-coordinates of data points.
vector of z-coordinates of data points or a character variable
naming the variable of interest in the
\code{SpatialPointsDataFrame} \code{x}.
Missing values are not accepted.
\code{x}, \code{y}, and \code{z} must be the same length and may
\code{x}, \code{y}, and \code{z} must be the same length
(execpt if \code{x} is a \code{SpatialPointsDataFrame}) 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
......@@ -59,7 +69,7 @@ interp.new(x, y, z, xo = seq(min(x), max(x), length = 40),
deprecated, use parameter \code{linear}. Now only used by
\code{interp.old()}.
meaning was:
old meaning was:
number of additional points to be used in computing partial
derivatives at each data point.
\code{ncp} must be either \code{0} (partial derivatives are not used), or at
......@@ -80,6 +90,8 @@ interp.new(x, y, z, xo = seq(min(x), max(x), length = 40),
}}
\item{dupfun}{a function, applied to duplicate points if
\code{duplicate= "user"}.}
\item{nx}{dimension of output grid in x direction}
\item{ny}{dimension of output grid in y direction}
}
\value{
list with 3 components:
......@@ -92,6 +104,9 @@ interp.new(x, y, z, xo = seq(min(x), max(x), length = 40),
matrix of fitted z-values. The value \code{z[i,j]} is computed
at the x,y point \code{xo[i], yo[j]}. \code{z} has
dimensions \code{length(xo)} times \code{length(yo)}.}
If input is a \code{SpatialPointsDataFrame} a
\code{SpatialPixelssDataFrame} is returned.
}
\note{
\code{interp} is a wrapper for the two versions \code{interp.old} (it
......@@ -163,12 +178,12 @@ akima.smooth <-
image (akima.smooth, main = "interp(<akima data>, *) on finer grid")
contour(akima.smooth, add = TRUE, col = "thistle")
points(akima, pch = 3, cex = 2, col = "blue")
# use triangulation package to show underlying triangulation:
## use triangulation package to show underlying triangulation:
\dontrun{
if(library(tripack, logical.return=TRUE))
plot(tri.mesh(akima), add=TRUE, lty="dashed")
}
# use only 15 points (interpolation only within convex hull!)
## use only 15 points (interpolation only within convex hull!)
akima.part <- with(akima, interp(x[1:15], y[1:15], z[1:15]))
image(akima.part)
title("interp() on subset of only 15 points")
......@@ -206,20 +221,32 @@ filled.contour(akima.spl, color.palette = full.pal,
plot.axes = { axis(1); axis(2);
title("smooth interp(*, linear = FALSE)");
points(akima, pch = 3, col= hcl(c=100, l = 20))})
# no extrapolation!
## no extrapolation!
## example with duplicate points :
data(airquality)
air <- subset(airquality,
!is.na(Temp) & !is.na(Ozone) & !is.na(Solar.R))
# gives an error {duplicate ..}:
## gives an error {duplicate ..}:
try( air.ip <- interp(air$Temp,air$Solar.R,air$Ozone, linear=FALSE) )
# use mean of duplicate points:
## use mean of duplicate points:
air.ip <- with(air, interp(Temp, Solar.R, log(Ozone), duplicate = "mean",
linear = FALSE))
image(air.ip, main = "Airquality: Ozone vs. Temp and Solar.R")
with(air, points(Temp, Solar.R))
\dontrun{
## interp can handle spatial point dataframes created by the sp package:
library(sp)
data(meuse)
coordinates(meuse) <- ~x+y
## argument z has to be named, y has to be omitted!
z <- interp(meuse,z="zinc",nx=100,ny=150)
spplot(z,"zinc")
z <- interp(meuse,z="zinc",nx=100,ny=150,linear=FALSE)
spplot(z,"zinc")
}
}
\keyword{dplot}
man/interpp.Rd
View file @ 43507cd4
......@@ -6,23 +6,32 @@
\alias{interpp.old}
\alias{interpp.new}
\usage{
interpp(x, y, z, xo, yo, linear=TRUE, extrap=FALSE, duplicate = "error",
dupfun = NULL, ncp)
interpp(x, y=NULL, z, xo, yo=NULL, linear=TRUE, extrap=FALSE,
duplicate = "error", dupfun = NULL, ncp)
}
\arguments{
\item{x}{
vector of x-coordinates of data points.
vector of x-coordinates of data points or a
\code{SpatialPointsDataFrame} object.
Missing values are not accepted.
}
\item{y}{
vector of y-coordinates of data points.
Missing values are not accepted.
If left as NULL indicates that \code{x} should be a
\code{SpatialPointsDataFrame} and \code{z} names the variable of
interest in this dataframe.
}
\item{z}{
vector of z-coordinates of data points.
vector of z-coordinates of data points or a character variable
naming the variable of interest in the
\code{SpatialPointsDataFrame} \code{x}.
Missing values are not accepted.
\code{x}, \code{y}, and \code{z} must be the same length and may contain no fewer
\code{x}, \code{y}, and \code{z} must be the same length
(execpt if \code{x} is a \code{SpatialPointsDataFrame}) 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
......@@ -30,10 +39,15 @@ dupfun = NULL, ncp)
}
\item{xo}{
vector of x-coordinates of points at which to evaluate the interpolating
function.}
function. If \code{x} is a \code{SpatialPointsDataFrame} this has
also to be a \code{SpatialPointsDataFrame}.
}
\item{yo}{
vector of y-coordinates of points at which to evaluate the interpolating
function.}
function.
If operating on \code{SpatialPointsDataFrame}s this is left as \code{NULL}
}
\item{linear}{logical -- indicating wether linear or spline
interpolation should be used. supersedes old \code{ncp} parameter}
\item{ncp}{
......@@ -74,7 +88,9 @@ dupfun = NULL, ncp)
\item{z}{
fitted z-values. The value \code{z[i]} is computed
at the x,y point \code{x[i], y[i]}.
}
}
If input is \code{SpatialPointsDataFrame} than an according
\code{SpatialPointsDataFrame} is returned.
}
\section{NOTE}{
Use \code{interp} if interpolation on a regular grid is wanted.
......@@ -134,5 +150,19 @@ akima.lip$z
akima.sip<-interpp(akima$x, akima$y, akima$z,c(1,5,10),c(2,6,12),
linear=FALSE)
akima.sip$z
\dontrun{
## interaction with sp objects:
library(sp)
## take 30 sample points out of meuse grid:
data(meuse.grid)
m0 <- meuse.grid[sample(1:3103,30),]
coordinates(m0) <- ~x+y
## interpolate on this 30 points:
## note: both "meuse" and "m0" are sp objects
## (SpatialPointsDataFrame) !!
## arguments z and xo have to named, y has to be omitted!
ipp <- interpp(meuse,z="zinc",xo=m0)
spplot(ipp)
}
}
\keyword{dplot}
src/ttidbs.f deleted 100644 → 0
View file @ ca58ed3a
PROGRAM TTIDBS
C PROGRAM TTIDBS(OUTPUT,TAPE6=OUTPUT) ID000070
C THIS PROGRAM IS A TEST PROGRAM FOR THE IDBVIP/IDSFFT SUBPRO- ID000080
C GRAM PACKAGE. ALL ELEMENTS OF RESULTING DZI1 AND DZI2 ARRAYS ID000090
C ARE EXPECTED TO BE ZERO. ID000100
C THE LUN CONSTANT IN THE LAST DATA INITIALIZATION STATEMENT IS ID000110
C THE LOGICAL UNIT NUMBER OF THE STANDARD OUTPUT UNIT AND IS, ID000120
C THEREFORE, SYSTEM DEPENDENT. ID000130
C DECLARATION STATEMENTS ID000140
IMPLICIT DOUBLE PRECISION (A-D,P-Z)
DIMENSION XD(30),YD(30),ZD(30), ID000150
1 XI(6),YI(5),ZI(6,5),MISSI(6,5), ID000160
2 ZI1(6,5),ZI2(6,5),DZI1(6,5),DZI2(6,5), ID000170
3 IWK(1030),WK(240) ID000180
LOGICAL MISSI
DATA NCP/4/ ID000190
DATA NDP/30/ ID000200
DATA XD(1), XD(2), XD(3), XD(4), XD(5), XD(6), ID000210
1 XD(7), XD(8), XD(9), XD(10),XD(11),XD(12), ID000220
2 XD(13),XD(14),XD(15),XD(16),XD(17),XD(18), ID000230
3 XD(19),XD(20),XD(21),XD(22),XD(23),XD(24), ID000240
4 XD(25),XD(26),XD(27),XD(28),XD(29),XD(30)/ ID000250
5 11.16, 24.20, 19.85, 10.35, 19.72, 0.00, ID000260
6 20.87, 19.99, 10.28, 4.51, 0.00, 16.70, ID000270
7 6.08, 25.00, 14.90, 0.00, 9.66, 5.22, ID000280
8 11.77, 15.10, 25.00, 25.00, 14.59, 15.20, ID000290
9 5.23, 2.14, 0.51, 25.00, 21.67, 3.31/ ID000300
DATA YD(1), YD(2), YD(3), YD(4), YD(5), YD(6), ID000310
1 YD(7), YD(8), YD(9), YD(10),YD(11),YD(12), ID000320
2 YD(13),YD(14),YD(15),YD(16),YD(17),YD(18), ID000330
3 YD(19),YD(20),YD(21),YD(22),YD(23),YD(24), ID000340
4 YD(25),YD(26),YD(27),YD(28),YD(29),YD(30)/ ID000350
5 1.24, 16.23, 10.72, 4.11, 1.39, 20.00, ID000360
6 20.00, 4.62, 15.16, 20.00, 4.48, 19.65, ID000370
7 4.58, 11.87, 3.12, 0.00, 20.00, 14.66, ID000380
8 10.47, 17.19, 3.87, 0.00, 8.71, 0.00, ID000390
9 10.72, 15.03, 8.37, 20.00, 14.36, 0.13/ ID000400
DATA ZD(1), ZD(2), ZD(3), ZD(4), ZD(5), ZD(6), ID000410
1 ZD(7), ZD(8), ZD(9), ZD(10),ZD(11),ZD(12), ID000420
2 ZD(13),ZD(14),ZD(15),ZD(16),ZD(17),ZD(18), ID000430
3 ZD(19),ZD(20),ZD(21),ZD(22),ZD(23),ZD(24), ID000440
4 ZD(25),ZD(26),ZD(27),ZD(28),ZD(29),ZD(30)/ ID000450
5 22.15, 2.83, 7.97, 22.33, 16.83, 34.60, ID000460
6 5.74, 14.72, 21.59, 15.61, 61.77, 6.31, ID000470
7 35.74, 4.40, 21.70, 58.20, 4.73, 40.36, ID000480
8 13.62, 12.57, 8.74, 12.00, 14.81, 21.60, ID000490
9 26.50, 53.10, 49.43, 0.60, 5.52, 44.08/ ID000500
DATA NXI/6/, NYI/5/ ID000510
DATA XI(1), XI(2), XI(3), XI(4), XI(5), XI(6)/ ID000520
1 0.00, 5.00, 10.00, 15.00, 20.00, 25.00/ ID000530
DATA YI(1), YI(2), YI(3), YI(4), YI(5)/ ID000540
1 0.00, 5.00, 10.00, 15.00, 20.00/ ID000550
DATA ZI(1,1),ZI(2,1),ZI(3,1),ZI(4,1),ZI(5,1),ZI(6,1), ID000560
1 ZI(1,2),ZI(2,2),ZI(3,2),ZI(4,2),ZI(5,2),ZI(6,2), ID000570
2 ZI(1,3),ZI(2,3),ZI(3,3),ZI(4,3),ZI(5,3),ZI(6,3), ID000580
3 ZI(1,4),ZI(2,4),ZI(3,4),ZI(4,4),ZI(5,4),ZI(6,4), ID000590
4 ZI(1,5),ZI(2,5),ZI(3,5),ZI(4,5),ZI(5,5),ZI(6,5)/ ID000600
5 58.20, 39.55, 26.90, 21.71, 17.68, 12.00, ID000610
6 61.58, 39.39, 22.04, 21.29, 14.36, 8.04, ID000620
7 59.18, 27.39, 16.78, 13.25, 8.59, 5.36, ID000630
8 52.82, 40.27, 22.76, 16.61, 7.40, 2.88, ID000640
9 34.60, 14.05, 4.12, 3.17, 6.31, 0.60/ ID000650
DATA LUN/6/ ID000660
C CALCULATION ID000670
10 MD=1 ID000680
DO 12 IYI=1,NYI ID000690
DO 11 IXI=1,NXI ID000700
IF(IXI.NE.1.OR.IYI.NE.1) MD=2 ID000710
MISSI(IXI,IYI)=.FALSE.
CALL IDBVIP(MD,NCP,NDP,XD,YD,ZD,1,XI(IXI),YI(IYI), ID000720
1 ZI1(IXI,IYI),IWK,WK,MISSI) ID000730
11 CONTINUE ID000740
12 CONTINUE ID000750
15 CALL IDSFFT(1,NCP,NDP,XD,YD,ZD,NXI,NYI,XI,YI,ZI2,IWK,WK,MISSI) ID000760
DO 17 IYI=1,NYI ID000770
DO 16 IXI=1,NXI ID000780
DZI1(IXI,IYI)=ABS(ZI1(IXI,IYI)-ZI(IXI,IYI)) ID000790
DZI2(IXI,IYI)=ABS(ZI2(IXI,IYI)-ZI(IXI,IYI)) ID000800
16 CONTINUE ID000810
17 CONTINUE ID000820
C PRINTING OF INPUT DATA ID000830
C WRITE (LUN,2020) NDP ID000840
DO 23 IDP=1,NDP ID000850
C WRITE (LUN,2021) ID000860
C WRITE (LUN,2022) IDP,XD(IDP),YD(IDP),ZD(IDP) ID000870
23 CONTINUE ID000880