SUBROUTINE INTERP(X,Y,XX,YY,NX,NXX,NPOL,ILOGX,ILOGY)
C     ====================================================
C
C     General interpolation procedure of the (NPOL-1)-th order
C
C     for  ILOGX = 1  logarithmic interpolation in X
C     for  ILOGY = 1  logarithmic interpolation in Y
C
C     Input:
C      X    - array of original x-coordinates
C      Y    - array of corresponding functional values Y=y(X)
C      NX   - number of elements in arrays X or Y
C      XX   - array of new x-coordinates (to which is to be 
C             interpolated
C      NXX  - number of elements in array XX
C     Output:
C      YY   - interpolated functional values YY=y(XX)
C
      INCLUDE 'PARAMS.FOR'
      DIMENSION X(1),Y(1),XX(1),YY(1)
      EXP10(X0)=EXP(X0*2.30258509299405D0)
      IF(NPOL.LE.0.OR.NX.LE.0) GO TO 200
      IF(ILOGX.NE.0) THEN
         DO I=1,NX
            X(I)=LOG10(X(I))
         END DO
         DO I=1,NXX
            XX(I)=LOG10(XX(I))
         END DO
      END IF
      IF(ILOGY.NE.0) THEN
         DO I=1,NX
            Y(I)=LOG10(Y(I))
         END DO
      END IF
      NM=(NPOL+1)/2      
      NM1=NM+1
      NUP=NX+NM1-NPOL
      DO ID=1,NXX
         XXX=XX(ID)
         DO I=NM1,NUP
            IF(XXX.LE.X(I)) GO TO 70
         END DO
         I=NUP
   70    J=I-NM
         JJ=J+NPOL-1
         YYY=0.
         DO K=J,JJ
            T=1.
            DO 80 M=J,JJ
               IF(K.EQ.M) GO TO 80
               T=T*(XXX-X(M))/(X(K)-X(M))
   80       CONTINUE
            YYY=Y(K)*T+YYY
         END DO
         YY(ID)=YYY
      END DO
      IF(ILOGX.NE.0) THEN
         DO I=1,NX
            X(I)=EXP10(X(I))
         END DO
         DO I=1,NXX
            XX(I)=EXP10(XX(I))
         END DO
      END IF
      IF(ILOGY.NE.0) THEN
         DO I=1,NX
            Y(I)=EXP10(Y(I))
         END DO
         DO I=1,NXX
           YY(I)=EXP10(YY(I))
         END DO
      END IF
      RETURN
  200 N=NX
      IF(NXX.GE.NX) N=NXX
      DO I=1,N
         XX(I)=X(I)
         YY(I)=Y(I)
      END DO
      RETURN
      END