onedint.f File Reference

Go to the source code of this file.

Functions/Subroutines
subroutine	onedint (XE, YE, NE, XA, YA, NA, IART, IEXP, IER)

Function/Subroutine Documentation

onedint()

subroutine onedint	(	real*8, dimension(ne)	XE,
		real*8, dimension(ne)	YE,
		integer	NE,
		real*8, dimension(na)	XA,
		real*8, dimension(na)	YA,
		integer	NA,
		integer	IART,
		integer	IEXP,
		integer	IER
	)

      implicit none
      INTEGER ne,na,na1,ne1,ig,ier,ia,iart,ie2,i,iexp,ie1,l
      REAL*8 xe(ne),ye(ne),xa(na),ya(na),zw1,zw2,xo,yo,rab,xd,yd,
     &  xz,yz,xu,yu,eq,eqd,x
C
C INTERPOLATION FUNCTION
C ------------------------
      eq(x) = yu + yu * (x-xu) / xu +
     1        ((yz-yu)/(xz-xu) - yu/xu) * (x-xu) * x / xz
      eqd(x) = yz * x / xz +
     1        (yd / xd - yz / xz) * x * (x - xz) / (xd - xz)
C
C INPUT/DATA TEST,INTERPOLATION DIVERGENCE,EXTRAPOLATION LIMIT
C----------------------------------------------------------------
      na1 = na - 1
      IF (na .LE. 0) GO TO 900
      ne1 = ne - 1
      IF (ne1.lt.0) then
         go to 900
      elseif(ne1.eq.0) then
         go to 22
      else
         go to 18
      endif
   18 DO 20 l = 1,ne1
   20 IF ((xe(l+1)-xe(l)) .LE. 0) GO TO 900
   22 ie1 = iexp / 10
      ie2 = iexp - 10*ie1
      ia = iart
      IF (ne1 .LT. ia)   ia = ne1
      IF (ia .LT. ie1)   ie1 = ia
      IF (ia .LT. ie2)   ie2 = ia
C
C SUCCESSIVE PROCESSING THE INTERPOLATION EXIGENCES
C-------------------------------------------------------
C
C     ZUR ERHOEHUNG DER NUMERISCHEN GENAUIGKEIT WIRD EINE
C     TRANSLATION VON (XO,YO) IN (0,0) DURCHGEFUEHRT. DIES
C     BEWIRKT AUSSERDEM  EINE BESCHLEUNIGUNG DES VERFAHRENS.
C
      DO 100 i = 1,na
      DO 24 l = 1,ne
      IF (xa(i) .LT. xe(l)) GO TO 30
   24 CONTINUE
      l = ne
      IF ((ie2 - 1).lt.0) then
         go to 50
      elseif((ie2-1).eq.0) then
         go to 35
      else
         go to 70
      endif
   30 IF (l .GT. 1) GO TO 40
      IF ((ie1 - 1).lt.0) then
         go to 50
      elseif((ie1-1).eq.0) then
         go to 25
      else
         go to 70
      endif
   40 IF ((ia-1).lt.0) then
         go to 45
      elseif((ia-1).eq.0) then
         go to 60
      else
         go to 70
      endif
C
C CONSTANT INTERPOLATION
C -----------------------
   45 l = l - 1
   50 ya(i) = ye(l)
      GO TO 100
C
C LINEAR EXTRAPOLATION
C ------------------------------
   25 IF (ia .EQ. 1) GO TO 60
      xo = xe(2)
      xu = xe(1) - xo
      yo = ye(2)
      yu = ye(1) - yo
      xz = xe(3) - xo
      yz = ye(3) - yo
      GO TO 38
   35 IF (ia .EQ. 1) GO TO 60
      xo = xe(ne1)
      xz = xe(ne1-1) - xo
      xu = xe(ne) - xo
      yo = ye(ne1)
      yz = ye(ne1-1) - yo
      yu = ye(ne) - yo
C
C LINEAR EXTRAPOLATION WITH QUADRATIC INTERPOLATION
C -----------------------------------------------------
   38 rab = yu / xu + xu * ((yz-yu) / (xz-xu) - yu/xu) / xz
      ya(i) = yu + yo + (xa(i) -xu-xo)*rab
      GO TO 100
C
C LINEAR INTERPOLATION
C ---------------------
  60  ig = l - 1
      IF (ig .LT. 1) ig = 1
      ya(i) = ye(ig) + (xa(i)-xe(ig))*(ye(ig+1)-ye(ig))
     1        / (xe(ig+1)-xe(ig))
      GO TO 100
   70 IF (l .GT. 2) GO TO 80
      xo = xe(2)
      xu = xe(1) - xo
      yo = ye(2)
      yu = ye(1) - yo
      xz = xe(3) - xo
      yz = ye(3) - yo
      GO TO 85
   80 IF (l .LT. ne) GO TO 90
      xo = xe(ne1)
      xu = xe(ne1-1) - xo
      xz = xe(ne) - xo
      yo = ye(ne1)
      yu = ye(ne1-1) - yo
      yz = ye(ne) - yo
   85 ya(i) = eq(xa(i)-xo) + yo
      GO TO 100
C
C DOUBLE QUADRATIC INTERPOLATION
C ----------------------------------
   90 xo = xe(l-1)
      xu = xe(l-2) - xo
      xz = xe(l) - xo
      xd = xe(l+1) - xo
      yo = ye(l-1)
      yu = ye(l-2) - yo
      yz = ye(l) - yo
      yd = ye(l+1) - yo
      zw1 = eq(xa(i)-xo)
      zw2 = eqd(xa(i)-xo)
      ya(i) = zw1 + (zw2 - zw1) * (xa(i) - xo)/xz + yo
  100 CONTINUE
C
C RETURN BY NORMAL PROCEEDING
C -------------------------------
      ier = 0
      RETURN
C
C ERROR RETURN
C ------------
  900 ier = -1
      RETURN