C$TEST PLYA
C TO RUN AS A MAIN PROGRAM REMOVE NEXT LINE
      SUBROUTINE PLYA
C***********************************************************************
C
C  TEST OF THE PORT PROGRAM TRIGONOMETRIC POLYNOMIAL ROUTINES
C
C***********************************************************************
      INTEGER IWRITE
      REAL COEF1(1001),COEF2(1001),X,STEP
      DOUBLE PRECISION X0,X1,Y0,Y1,Z0,Z1,T0(2),T1(2),T2(2)
      IWRITE = I1MACH(2)
      WRITE(IWRITE,10)
 10   FORMAT(3H  X,14X,3HSUM,25X,10HTRIGP A.E.,5X,10HTRIGP R.E.,5X,
     111HHORNER A.E.,4X,11HHORNER R.E.)
      I = 1000
      STEP = 1.0/1024.0
      PI2 = 3.141592653/2.0
      X=STEP
23000 IF(.NOT.(X.LE.PI2))GOTO 23002
      X0 = X
      T0(1) = 1.0D0 - DCOS(FLOAT(I)*X0)
      T0(2) = -DSIN(FLOAT(I)*X0)
      T1(1) = 1.0D0 - DCOS(X0)
      T1(2) = -DSIN(X0)
      IF(.NOT.(T1(1) .EQ. 0.0D0 .AND. T1(2) .EQ. 0.0D0))GOTO 23003
      T2(1) =FLOAT(I)
      T2(2) = 0.0D0
      GOTO 23004
23003 CONTINUE
      CALL CDDIV(T0,T1,T2)
23004 CONTINUE
      DO 23005 J=1,I
      COEF1(J) = 1.0
      COEF2(J) = 0.0
23005 CONTINUE
      X0 = TRIGP(I-1,COEF1,COEF2,X)
      X1 = TRIG01(COEF1,COEF2,I,X)
      Y0 = DABS(T2(1) - X0)
      Y1 = DABS(T2(1) - X1)
      Z0 = Y0/DABS(T2(1))
      Z1 = Y1/DABS(T2(1))
      WRITE(IWRITE,20) X,T2(1),Y0,Z0,Y1,Z1
20    FORMAT(1PE15.6,D25.17,4D15.4)
      X0 = TRIGP(I-1,COEF2,COEF1,X)
      X1 = TRIG01(COEF2,COEF1,I,X)
      Y0 = DABS(T2(2) - X0)
      Y1 = DABS(T2(2) - X1)
      Z0 = Y0/DABS(T2(2))
      Z1 = Y1/DABS(T2(2))
      WRITE(IWRITE,30) T2(2),Y0,Z0,Y1,Z1
30    FORMAT(15X,1PD25.17,4D15.4)
       X=X+STEP
      GOTO 23000
23002 CONTINUE
      STOP
      END
      FUNCTION TRIG01(ALFA,BETA,NP1,THETA)
      INTEGERNP1,K
      REAL ALFA(NP1),BETA(NP1),THETA,C,S,MU,GK,GK1,DK,DK1
      IF(.NOT.(NP1 .LE. 0))GOTO 23007
C/6S
      CALL SETERR(26HTRIG01 - INVALID DIMENSION,26,1,2)
C/7S
C     CALL SETERR('TRIG01 - INVALID DIMENSION',26,1,2)
C/
      GOTO 23008
23007 CONTINUE
      IF(.NOT.(NP1 .EQ. 1))GOTO 23009
      TRIG01 = ALFA(1)
      GOTO 23010
23009 CONTINUE
      C = COS(THETA)
      S = SIN(THETA)
      GK1 = ALFA(NP1)
      DK1 = BETA(NP1)
      K=NP1
23011 IF(.NOT.(K.GE.2))GOTO 23013
      GK = C*GK1 + S*DK1
      DK = -S*GK1 + C*DK1
      GK1 = GK + ALFA(K-1)
      DK1 = DK + BETA(K-1)
       K=K-1
      GOTO 23011
23013 CONTINUE
      TRIG01 = ALFA(1) + GK
23010 CONTINUE
23008 CONTINUE
      RETURN
      END
