Commit 8c525f1e by agebhard

### add aspline interface to Akimas univariate interpolation functions,

provided by Thomas Petzoldt <petzoldt@rcs.urz.tu-dresden.de>.
parent 0ed973ce
 Package: akima Version: 0.3-4 Version: 0.4-1 Title: Interpolation of irregularly spaced data Author: Fortran code by H. Akima R port by Albrecht Gebhardt aspline function by Thomas Petzoldt 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
 akima Waveform Distortion Data for Bivariate Interpolation aspline Univariate Akima interpolation interp Gridded Bivariate Interpolation for Irregular Data interpp Pointwise Bivariate Interpolation for ... ...
R/aspline.R 0 → 100644
 aspline <- function(x, y, xout, n = 50, ties=mean, algorithm="original", np=3) { if (! algorithm %in% c("original", "improved")) stop(paste("unknown algorithm:", algorithm)) x <- xy.coords(x, y) # -> (x,y) numeric of same length y <- x$y x <- x$x nx <- length(x) if(any(na <- is.na(x) | is.na(y))) { ok <- !na x <- x[ok] y <- y[ok] nx <- length(x) } if (!identical(ties, "ordered")) { if (length(ux <- unique(x)) < nx) { if (missing(ties)) warning("Collapsing to unique x values") y <- as.vector(tapply(y,x,ties))# as.v: drop dim & dimn. x <- sort(ux) nx <- length(x) } else { o <- order(x) x <- x[o] y <- y[o] } } if (nx <= 1) stop("need at least two non-NA values to interpolate") if (missing(xout)) { if (n <= 0) stop("approx requires n >= 1") xout <- seq(x[1], x[nx], length = n) } nout <- length(xout) yout <- numeric(nout) if (algorithm == "improved") { y <- .Fortran("uvip3p", as.integer(np), as.integer(nx), as.double(x), as.double(y), as.integer(nout), as.double(xout), yout=as.double(yout), PACKAGE="akima")$yout } else { y <- .Fortran("intrpl", as.integer(nx), as.double(x), as.double(y), as.integer(nout), as.double(xout), yout=as.double(yout), PACKAGE="akima")$yout } list(x = xout, y = y) }
 ... ... @@ -15,11 +15,17 @@ interp() actually calls interp.new() or interp.old() depending on its arguments, linear interpolation is done by interp.old(), spline interpolation by interp.new . UPDATE (13.12.2004): Thomas Petzoldt provided an R interface to Akimas univariate spline interpolation functions (aspline) which has been included into this library. ------------------------------------------------------------------------------ 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 WWW : http://pc02-stat.sci.uni-klu.ac.at/~agebhard 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 A-9020 Klagenfurt, Austria ------------------------------------------------------------------------------
man/aspline.Rd 0 → 100644
 \name{aspline} \alias{aspline} %- Also NEED an '\alias' for EACH other topic documented here. \title{Univariate Akima interpolation} \description{ The function returns a list of points which smoothly interpolate given data points, similar to a curve drawn by hand.} \usage{ aspline(x, y, xout, n = 50, ties = mean, algorithm="original", np=3) } %- maybe also 'usage' for other objects documented here. \arguments{ \item{x, y}{vectors giving the coordinates of the points to be interpolated. Alternatively a single plotting structure can be specified: see \code{\link{xy.coords}}.} \item{xout}{an optional set of values specifying where interpolation is to take place.} \item{n}{If \code{xout} is not specified, interpolation takes place at \code{n} equally spaced points spanning the interval [\code{min(x)}, \code{max(x)}].} \item{ties}{Handling of tied \code{x} values. Either a function with a single vector argument returning a single number result or the string \code{"ordered"}.} \item{algorithm}{either \code{"original"} method after Akima (1970) or \code{"improved"} method after Akima (1991)} \item{np}{im improved algorithm is selected: degree of the polynomials for the interpolating function} } \details{ The original algorithm is based on a piecewise function composed of a set of polynomials, each of degree three, at most, and applicable to successive interval of the given points. In this method, the slope of the curve is determined at each given point locally, and each polynomial representing a portion of the curve between a pair of given points is determined by the coordinates of and the slopes at the points. } \value{ A list with components \code{x} and \code{y}, containing \code{n} coordinates which interpolate the given data points. } \references{ Akima, H. (1970) A new method of interpolation and smooth curve fitting based on local procedures, J. ACM 17(4), 589-602 Akima, H. (1991) A Method of Univariate Interpolation that Has the Accuracy of a Third-degree Polynomial. ACM Transactions on Mathematical Software, 17(3), 341-366. } \author{Thomas Petzoldt} %%\note{ ~~further notes~~ } %% ~Make other sections like Warning with \section{Warning }{....} ~ \seealso{\code{\link{approx}}, \code{\link{spline}}} \examples{ ## regular spaced data x <- 1:10 y <- c(rnorm(5), c(1,1,1,1,3)) xnew <- seq(-1, 11, 0.1) plot(x, y, ylim=c(-3, 3), xlim=range(xnew)) lines(spline(x, y, xmin=min(xnew), xmax=max(xnew), n=200), col="blue") lines(aspline(x, y, xnew), col="red") lines(aspline(x, y, xnew, algorithm="improved"), col="black", lty="dotted") lines(aspline(x, y, xnew, algorithm="improved", np=10), col="green", lty="dashed") ## irregular spaced data x <- sort(runif(10, max=10)) y <- c(rnorm(5), c(1,1,1,1,3)) xnew <- seq(-1, 11, 0.1) plot(x, y, ylim=c(-3, 3), xlim=range(xnew)) lines(spline(x, y, xmin=min(xnew), xmax=max(xnew), n=200), col="blue") lines(aspline(x, y, xnew), col="red") lines(aspline(x, y, xnew, algorithm="improved"), col="black", lty="dotted") lines(aspline(x, y, xnew, algorithm="improved", np=10), col="green", lty="dashed") ## an example of Akima, 1991 x <- c(-3, -2, -1, 0, 1, 2, 2.5, 3) y <- c( 0, 0, 0, 0, -1, -1, 0, 2) plot(x, y, ylim=c(-3, 3)) lines(spline(x, y, n=200), col="blue") lines(aspline(x, y, n=200), col="red") lines(aspline(x, y, n=200, algorithm="improved"), col="black", lty="dotted") lines(aspline(x, y, n=200, algorithm="improved", np=10), col="green", lty="dashed") } \keyword{arith} \keyword{dplot} \ No newline at end of file
orig/433 0 → 100644
 C ALGORITHM 433 COLLECTED ALGORITHMS FROM ACM. C ALGORITHM APPEARED IN COMM. ACM, VOL. 15, NO. 10, C P. 914. SUBROUTINE INTRPL(IU,L,X,Y,N,U,V) INTR 10 C INTERPOLATION OF A SINGLE-VALUED FUNCTION C THIS SUBROUTINE INTERPOLATES, FROM VALUES OF THE FUNCTION C GIVEN AS ORDINATES OF INPUT DATA POINTS IN AN X-Y PLANE C AND FOR A GIVEN SET OF X VALUES (ABSCISSAS), THE VALUES OF C A SINGLE-VALUED FUNCTION Y = Y(X). C THE INPUT PARAMETERS ARE C IU = LOGICAL UNIT NUMBER OF STANDARD OUTPUT UNIT C L = NUMBER OF INPUT DATA POINTS C (MUST BE 2 OR GREATER) C X = ARRAY OF DIMENSION L STORING THE X VALUES C (ABSCISSAS) OF INPUT DATA POINTS C (IN ASCENDING ORDER) C Y = ARRAY OF DIMENSION L STORING THE Y VALUES C (ORDINATES) OF INPUT DATA POINTS C N = NUMBER OF POINTS AT WHICH INTERPOLATION OF THE C Y VALUE (ORDINATE) IS DESIRED C (MUST BE 1 OR GREATER) C U = ARRAY OF DIMENSION N STORING THE X VALUES C (ABSCISSAS) OF DESIRED POINTS C THE OUTPUT PARAMETER IS C V = ARRAY OF DIMENSION N WHERE THE INTERPOLATED Y C VALUES (ORDINATES) ARE TO BE DISPLAYED C DECLARATION STATEMENTS DIMENSION X(L),Y(L),U(N),V(N) EQUIVALENCE (P0,X3),(Q0,Y3),(Q1,T3) REAL M1,M2,M3,M4,M5 EQUIVALENCE (UK,DX),(IMN,X2,A1,M1),(IMX,X5,A5,M5), 1 (J,SW,SA),(Y2,W2,W4,Q2),(Y5,W3,Q3) C PRELIMINARY PROCESSING 10 L0=L LM1=L0-1 LM2=LM1-1 LP1=L0+1 N0=N IF(LM2.LT.0) GO TO 90 IF(N0.LE.0) GO TO 91 DO 11 I=2,L0 IF(X(I-1)-X(I)) 11,95,96 11 CONTINUE IPV=0 C MAIN DO-LOOP DO 80 K=1,N0 UK=U(K) C ROUTINE TO LOCATE THE DESIRED POINT 20 IF(LM2.EQ.0) GO TO 27 IF(UK.GE.X(L0)) GO TO 26 IF(UK.LT.X(1)) GO TO 25 IMN=2 IMX=L0 21 I=(IMN+IMX)/2 IF(UK.GE.X(I)) GO TO 23 22 IMX=I GO TO 24 23 IMN=I+1 24 IF(IMX.GT.IMN) GO TO 21 I=IMX GO TO 30 25 I=1 GO TO 30 26 I=LP1 GO TO 30 27 I=2 C CHECK IF I = IPV 30 IF(I.EQ.IPV) GO TO 70 IPV=I C ROUTINES TO PICK UP NECESSARY X AND Y VALUES AND C TO ESTIMATE THEM IF NECESSARY 40 J=I IF(J.EQ.1) J=2 IF(J.EQ.LP1) J=L0 X3=X(J-1) Y3=Y(J-1) X4=X(J) Y4=Y(J) A3=X4-X3 M3=(Y4-Y3)/A3 IF(LM2.EQ.0) GO TO 43 IF(J.EQ.2) GO TO 41 X2=X(J-2) Y2=Y(J-2) A2=X3-X2 M2=(Y3-Y2)/A2 IF(J.EQ.L0) GO TO 42 41 X5=X(J+1) Y5=Y(J+1) A4=X5-X4 M4=(Y5-Y4)/A4 IF(J.EQ.2) M2=M3+M3-M4 GO TO 45 42 M4=M3+M3-M2 GO TO 45 43 M2=M3 M4=M3 45 IF(J.LE.3) GO TO 46 A1=X2-X(J-3) M1=(Y2-Y(J-3))/A1 GO TO 47 46 M1=M2+M2-M3 47 IF(J.GE.LM1) GO TO 48 A5=X(J+2)-X5 M5=(Y(J+2)-Y5)/A5 GO TO 50 48 M5=M4+M4-M3 C NUMERICAL DIFFERENTIATION 50 IF(I.EQ.LP1) GO TO 52 W2=ABS(M4-M3) W3=ABS(M2-M1) SW=W2+W3 IF(SW.NE.0.0) GO TO 51 W2=0.5 W3=0.5 SW=1.0 51 T3=(W2*M2+W3*M3)/SW IF(I.EQ.1) GO TO 54 52 W3=ABS(M5-M4) W4=ABS(M3-M2) SW=W3+W4 IF(SW.NE.0.0) GO TO 53 W3=0.5 W4=0.5 SW=1.0 53 T4=(W3*M3+W4*M4)/SW IF(I.NE.LP1) GO TO 60 T3=T4 SA=A2+A3 T4=0.5*(M4+M5-A2*(A2-A3)*(M2-M3)/(SA*SA)) X3=X4 Y3=Y4 A3=A2 M3=M4 GO TO 60 54 T4=T3 SA=A3+A4 T3=0.5*(M1+M2-A4*(A3-A4)*(M3-M4)/(SA*SA)) X3=X3-A4 Y3=Y3-M2*A4 A3=A4 M3=M2 C DETERMINATION OF THE COEFFICIENTS 60 Q2=(2.0*(M3-T3)+M3-T4)/A3 Q3=(-M3-M3+T3+T4)/(A3*A3) C COMPUTATION OF THE POLYNOMIAL 70 DX=UK-P0 80 V(K)=Q0+DX*(Q1+DX*(Q2+DX*Q3)) RETURN C ERROR EXIT 90 WRITE (IU,2090) GO TO 99 91 WRITE (IU,2091) GO TO 99 95 WRITE (IU,2095) GO TO 97 96 WRITE (IU,2096) 97 WRITE (IU,2097) I,X(I) 99 WRITE (IU,2099) L0,N0 RETURN C FORMAT STATEMENTS 2090 FORMAT(1X/22H *** L = 1 OR LESS./) 2091 FORMAT(1X/22H *** N = 0 OR LESS./) 2095 FORMAT(1X/27H *** IDENTICAL X VALUES./) 2096 FORMAT(1X/33H *** X VALUES OUT OF SEQUENCE./) 2097 FORMAT(6H I =,I7,10X,6HX(I) =,E12.3) 2099 FORMAT(6H L =,I7,10X,3HN =,I7/ 1 36H ERROR DETECTED IN ROUTINE INTRPL) END SUBROUTINE CRVFIT(IU,MD,L,X,Y,M,N,U,V) CRVF1670 C SMOOTH CURVE FITTING C THIS SUBROUTINE FITS A SMOOTH CURVE TO A GIVEN SET OF IN- C PUT DATA POINTS IN AN X-Y PLANE. IT INTERPOLATES POINTS C IN EACH INTERVAL BETWEEN A PAIR OF DATA POINTS AND GENER- C ATES A SET OF OUTPUT POINTS CONSISTING OF THE INPUT DATA C POINTS AND THE INTERPOLATED POINTS. IT CAN PROCESS EITHER C A SINGLE-VALUED FUNCTION OR A MULTIPLE-VALUED FUNCTION. C THE INPUT PARAMETERS ARE C IU = LOGICAL UNIT NUMBER OF STANDARD OUTPUT UNIT C MD = MODE OF THE CURVE (MUST BE 1 OR 2) C = 1 FOR A SINGLE-VALUED FUNCTION C = 2 FOR A MULTIPLE-VALUED FUNCTION C L = NUMBER OF INPUT DATA POINTS C (MUST BE 2 OR GREATER) C X = ARRAY OF DIMENSION L STORING THE ABSCISSAS OF C INPUT DATA POINTS (IN ASCENDING OR DESCENDING C ORDER FOR MD = 1) C Y = ARRAY OF DIMENSION L STORING THE ORDINATES OF C INPUT DATA POINTS C M = NUMBER OF SUBINTERVALS BETWEEN EACH PAIR OF C INPUT DATA POINTS (MUST BE 2 OR GREATER) C N = NUMBER OF OUTPUT POINTS C = (L-1)*M+1 C THE OUTPUT PARAMETERS ARE C U = ARRAY OF DIMENSION N WHERE THE ABSCISSAS OF C OUTPUT POINTS ARE TO BE DISPLAYED C V = ARRAY OF DIMENSION N WHERE THE ORDINATES OF C OUTPUT POINTS ARE TO BE DISPLAYED C DECLARATION STATEMENTS DIMENSION X(L),Y(L),U(N),V(N) EQUIVALENCE (M1,B1),(M2,B2),(M3,B3),(M4,B4), 1 (X2,P0),(Y2,Q0),(T2,Q1) REAL M1,M2,M3,M4 EQUIVALENCE (W2,Q2),(W3,Q3),(A1,P2),(B1,P3), 1 (A2,DZ),(SW,R,Z) C PRELIMINARY PROCESSING 10 MD0=MD MDM1=MD0-1 L0=L LM1=L0-1 M0=M MM1=M0-1 N0=N IF(MD0.LE.0) GO TO 90 IF(MD0.GE.3) GO TO 90 IF(LM1.LE.0) GO TO 91 IF(MM1.LE.0) GO TO 92 IF(N0.NE.LM1*M0+1) GO TO 93 GO TO (11,16), MD0 11 I=2 IF(X(1)-X(2)) 12,95,14 12 DO 13 I=3,L0 IF(X(I-1)-X(I)) 13,95,96 13 CONTINUE GO TO 18 14 DO 15 I=3,L0 IF(X(I-1)-X(I)) 96,95,15 15 CONTINUE GO TO 18 16 DO 17 I=2,L0 IF(X(I-1).NE.X(I)) GO TO 17 IF(Y(I-1).EQ.Y(I)) GO TO 97 17 CONTINUE 18 K=N0+M0 I=L0+1 DO 19 J=1,L0 K=K-M0 I=I-1 U(K)=X(I) 19 V(K)=Y(I) RM=M0 RM=1.0/RM C MAIN DO-LOOP 20 K5=M0+1 DO 80 I=1,L0 C ROUTINES TO PICK UP NECESSARY X AND Y VALUES AND C TO ESTIMATE THEM IF NECESSARY IF(I.GT.1) GO TO 40 30 X3=U(1) Y3=V(1) X4=U(M0+1) Y4=V(M0+1) A3=X4-X3 B3=Y4-Y3 IF(MDM1.EQ.0) M3=B3/A3 IF(L0.NE.2) GO TO 41 A4=A3 B4=B3 31 GO TO (33,32), MD0 32 A2=A3+A3-A4 A1=A2+A2-A3 33 B2=B3+B3-B4 B1=B2+B2-B3 GO TO (51,56), MD0 40 X2=X3 Y2=Y3 X3=X4 Y3=Y4 X4=X5 Y4=Y5 A1=A2 B1=B2 A2=A3 B2=B3 A3=A4 B3=B4 IF(I.GE.LM1) GO TO 42 41 K5=K5+M0 X5=U(K5) Y5=V(K5) A4=X5-X4 B4=Y5-Y4 IF(MDM1.EQ.0) M4=B4/A4 GO TO 43 42 IF(MDM1.NE.0) A4=A3+A3-A2 B4=B3+B3-B2 43 IF(I.EQ.1) GO TO 31 GO TO (50,55), MD0 C NUMERICAL DIFFERENTIATION 50 T2=T3 51 W2=ABS(M4-M3) W3=ABS(M2-M1) SW=W2+W3 IF(SW.NE.0.0) GO TO 52 W2=0.5 W3=0.5 SW=1.0 52 T3=(W2*M2+W3*M3)/SW IF(I-1) 80,80,60 55 COS2=COS3 SIN2=SIN3 56 W2=ABS(A3*B4-A4*B3) W3=ABS(A1*B2-A2*B1) IF(W2+W3.NE.0.0) GO TO 57 W2=SQRT(A3*A3+B3*B3) W3=SQRT(A2*A2+B2*B2) 57 COS3=W2*A2+W3*A3 SIN3=W2*B2+W3*B3 R=COS3*COS3+SIN3*SIN3 IF(R.EQ.0.0) GO TO 58 R=SQRT(R) COS3=COS3/R SIN3=SIN3/R 58 IF(I-1) 80,80,65 C DETERMINATION OF THE COEFFICIENTS 60 Q2=(2.0*(M2-T2)+M2-T3)/A2 Q3=(-M2-M2+T2+T3)/(A2*A2) GO TO 70 65 R=SQRT(A2*A2+B2*B2) P1=R*COS2 P2=3.0*A2-R*(COS2+COS2+COS3) P3=A2-P1-P2 Q1=R*SIN2 Q2=3.0*B2-R*(SIN2+SIN2+SIN3) Q3=B2-Q1-Q2 GO TO 75 C COMPUTATION OF THE POLYNOMIALS 70 DZ=A2*RM Z=0.0 DO 71 J=1,MM1 K=K+1 Z=Z+DZ U(K)=P0+Z 71 V(K)=Q0+Z*(Q1+Z*(Q2+Z*Q3)) GO TO 79 75 Z=0.0 DO 76 J=1,MM1 K=K+1 Z=Z+RM U(K)=P0+Z*(P1+Z*(P2+Z*P3)) 76 V(K)=Q0+Z*(Q1+Z*(Q2+Z*Q3)) 79 K=K+1 80 CONTINUE RETURN C ERROR EXIT 90 WRITE (IU,2090) GO TO 99 91 WRITE (IU,2091) GO TO 99 92 WRITE (IU,2092) GO TO 99 93 WRITE (IU,2093) GO TO 99 95 WRITE (IU,2095) GO TO 98 96 WRITE (IU,2096) GO TO 98 97 WRITE (IU,2097) 98 WRITE (IU,2098) I,X(I),Y(I) 99 WRITE (IU,2099) MD0,L0,M0,N0 RETURN C FORMAT STATEMENTS 2090 FORMAT(1X/31H *** MD OUT OF PROPER RANGE./) 2091 FORMAT(1X/22H *** L = 1 OR LESS./) 2092 FORMAT(1X/22H *** M = 1 OR LESS./) 2093 FORMAT(1X/25H *** IMPROPER N VALUE./) 2095 FORMAT(1X/27H *** IDENTICAL X VALUES./) 2096 FORMAT(1X/33H *** X VALUES OUT OF SEQUENCE./) 2097 FORMAT(1X/33H *** IDENTICAL X AND Y VALUES./) 2098 FORMAT(7H I =,I4,10X,6HX(I) =,E12.3, 1 10X,6HY(I) =,E12.3) 2099 FORMAT(7H MD =,I4,8X,3HL =,I5,8X, 1 3HM =,I5,8X,3HN =,I5/ 2 36H ERROR DETECTED IN ROUTINE CRVFIT) END
orig/697 0 → 100644
This diff is collapsed.
src/akima433.f 0 → 100644
 C Modified after ACM Algorithm 433 by ThPe C C ALGORITHM 433 COLLECTED ALGORITHMS FROM ACM. C ALGORITHM APPEARED IN COMM. ACM, VOL. 15, NO. 10, C P. 914. SUBROUTINE INTRPL(L,X,Y,N,U,V) C INTERPOLATION OF A SINGLE-VALUED FUNCTION C THIS SUBROUTINE INTERPOLATES, FROM VALUES OF THE FUNCTION C GIVEN AS ORDINATES OF INPUT DATA POINTS IN AN X-Y PLANE C AND FOR A GIVEN SET OF X VALUES (ABSCISSAS), THE VALUES OF C A SINGLE-VALUED FUNCTION Y = Y(X). C THE INPUT PARAMETERS ARE C L = NUMBER OF INPUT DATA POINTS C (MUST BE 2 OR GREATER) C X = ARRAY OF DIMENSION L STORING THE X VALUES C (ABSCISSAS) OF INPUT DATA POINTS C (IN ASCENDING ORDER) C Y = ARRAY OF DIMENSION L STORING THE Y VALUES C (ORDINATES) OF INPUT DATA POINTS C N = NUMBER OF POINTS AT WHICH INTERPOLATION OF THE C Y VALUE (ORDINATE) IS DESIRED C (MUST BE 1 OR GREATER) C U = ARRAY OF DIMENSION N STORING THE X VALUES C (ABSCISSAS) OF DESIRED POINTS C THE OUTPUT PARAMETER IS C V = ARRAY OF DIMENSION N WHERE THE INTERPOLATED Y C VALUES (ORDINATES) ARE TO BE DISPLAYED C DECLARATION STATEMENTS INTEGER L, N DOUBLE PRECISION X(L),Y(L),U(N),V(N) EQUIVALENCE (P0,X3),(Q0,Y3),(Q1,T3) DOUBLE PRECISION M1,M2,M3,M4,M5 DOUBLE PRECISION A2,A4 EQUIVALENCE (UK,DX),(IMN,X2,A1,M1),(IMX,X5,A5,M5), 1 (J,SW,SA),(Y2,W2,W4,Q2),(Y5,W3,Q3) C PRELIMINARY PROCESSING 10 L0=L LM1=L0-1 LM2=LM1-1 LP1=L0+1 N0=N M2=0 A2=0 A4=0 IF(LM2.LT.0) GO TO 90 IF(N0.LE.0) GO TO 91 DO 11 I=2,L0 IF(X(I-1)-X(I)) 11,95,96 11 CONTINUE IPV=0 C MAIN DO-LOOP DO 80 K=1,N0 UK=U(K) C ROUTINE TO LOCATE THE DESIRED POINT 20 IF(LM2.EQ.0) GO TO 27 IF(UK.GE.X(L0)) GO TO 26 IF(UK.LT.X(1)) GO TO 25 IMN=2 IMX=L0 21 I=(IMN+IMX)/2 IF(UK.GE.X(I)) GO TO 23 22 IMX=I GO TO 24 23 IMN=I+1 24 IF(IMX.GT.IMN) GO TO 21 I=IMX GO TO 30 25 I=1 GO TO 30 26 I=LP1 GO TO 30 27 I=2