Introduction
The group of
points shown to the left becomes a function through the process of
interpolation. In general the region in
which the interpolating function is most valid is that between points 12 and
13. The fact that the polynomial goes through points 12 and 13 give rise to the
hope that it is reasonably accurate between them. The fact that it also goes through points 11
and 14 implies that the first derivatives should be accurate at each side of
the (12,13) region. The set of points
will in general be spaced such that there are more points in the regions where
the function is rapidly changing.
The Fortran test code is in Flocate.zip TLOCATE.WPJ.
The very first problem is counting the lines in the data file. The main code opens a file and calls
FUNCTION NLINES(NUNIT)
N=0 ß
This could be NLINES, but other languages object
REWIND(NUNIT)
5 READ(NUNIT,*,ERR=10,END=10) ß
the end= and err= could be
N=N+1
ß separated and the whle loop coulc be do
while
GOTO 5
ß Watcom’s version of do while differs from
F90
10 REWIND(NUNIT)
NLINES=N
RETURN
END
There are no data conversions; FORTRAN uses line ends to know when to stop. This routine is a minimal overhead paid for not having the number of items in the file. There is a problem. If for some reason the end and err are not triggered, the code will continue reading through the entire hard drive. Fortunately the code does not write anything, so it will not destroy the computer software. It will eventually stop or overflow the integers. A do while does nothing to help this.
Tlocate.for
IMPLICIT REAL*8 (A-H,O-Z)
DIMENSION XI(:),FI(:) ß
dynamic arrays
CI ALLOCATABLE XI,FI ß
necessary in F90, error in Watcom
OPEN(1,FILE='LORENT.OUT') ß
a specific file
N=NLINES(1) ß
the function above
ALLOCATE(XI(N),FI(N))
DO I=1,N
READ(1,*)XI(I),FI(I) ß
the * tells the code to figure it out
ENDDO
PRINT*,'LORENT.OUT HAS',N,' DATA POINTS'
15 PRINT*,'XI(1),XI(N)',XI(1),XI(N),' ENTER X
OR -1000 TO STOP'
READ(*,*)X
IF(X.EQ.-1000)STOP
CALL LOCATE(X,XI,N,J)
IF(J.EQ.0)THEN
PRINT*,' X IS LESS THAN ',XI(1)
ELSEIF(J.EQ.N)THEN
PRINT*,' x IS GREATER THAN last point
',XI(N)
ELSE
PRINT*,' J= ',J
PRINT*,XI(J),X,XI(J+1)
ENDIF
GOTO 15
END
Locate.for
SUBROUTINE LOCATE(X,R,NMAX,J)
C
RETURNS A VALUE J SUCH THAT R(J)<=X<=R(J+1)
IMPLICIT REAL*8 (A-H,O-Z)
DIMENSION R(NMAX) ß
compilers gave up and allow us to use * or 1
JL=0 ß
but NMAX is more likely to find errors
JU=NMAX+1
10 IF(JU-JL.LE.1)GOTO 20
JM=(JU+JL)/2
IF(X.GT.R(JM))THEN
JL=JM
ELSE
JU=JM
ENDIF
GOTO 10
20 J=JL
RETURN
END
The tlocate source in C
#include <stdio.h>
#include <io.h>
void *calloc( size_t n, size_t size ); ß these help in debugging
int fgetc( FILE *fp );
void rewind( FILE *fp );
int linecount( FILE *fp ) ß there is a small advantage to having linecount in this file
{int c,n1;
n1=0;
while( (c = fgetc( fp )) != EOF ){ ß these two lines do the work. It is slightly more if(c == '\n')++n1;} ß than the FORTRAN file since it reads every character.
rewind(fp);
return n1;}
int locate(double x,double r[],int n); ß a needed header
void main(void)
{
double x,*xi,*fi;
int n,i;
FILE *fp;
fp = fopen("loren.out","r");
if(fp == NULL){printf(" cannot open file loren.out \n");
return;}
n=linecount(fp);
xi = (double *) calloc(n, sizeof(double));
fi = (double *) calloc(n, sizeof(double));
for (i=0;i<n;i++)fscanf(fp, "%lg %lg", &xi[i],&fi[i]); ß more compact than FORTRAN
fclose(fp);
while (x != -10000)
{printf("xi[1] xi[n-1] %lg %lg enter x or -10000 to quit \n",xi[0],xi[n-1]);
scanf("%lf",&x);
i=locate(x,xi,n);
printf(" after locate i = %d \n",i);
if(i >= 0 && i < n-1)printf(" %f %f %f \n",xi[i],x,xi[i+1]);
else if (i == -1 )printf(" %f %f \n",x,xi[i+1]);
else if (i == n-1 )printf(" %f %f \n",xi[i],x);}
return; }
The locate source in C
int locate(double x,double r[],int n)
{ int jl = -1;
int ju = n;
int jm;
while (ju-jl>1)
{jm=(ju+jl) >> 1; ß integer division by 2 is a right shift
if(x > r[jm])jl=jm;
else ju=jm; }
return jl;}