C....*...1.........2.........3.........4.........5.........6.........7.*.......8
C     SPEC     2/18/75
C
C     PURPOSE
C     COMPUTE THE ESTIMATED SPECTRAL DENSITY FROM ESTIMATED
C     AUTOCOVARIANCES AT A GIVEN FREQUENCY.
C
C     USAGE
C     F=SPEC(FREQ,R,LAG,ISW,WORK)
C
C     ARGUMENTS
C     FREQ - FREQUENCY IN RADIANS.
C            REAL*8
C     R    - INPUT VECTOR OF LENGTH LAG+1 CONTAINING ESTIMATED
C            AUTOCOVARIANCES FOR LAGS O,1,...,LAG.
C            REAL*8
C     LAG  - ORDER OF LAG TO BE USED TO COMPUTE THE SPECTRAL DENSITY
C            REAL*8
C     ISW  - INPUT SWITCH USED TO SELECT SMOOTHING WEIGHTS:
C            ISW=1  RECTANGULAR LAG WINDOW
C            ISW=2  BARTLETT LAG WINDOW
C            ISW=3  TUKEY LAG WINDOW
C            ISW=4  PARZEN LAG WINDOW
C            INTEGER*4
C     WORK - WORK VECTOR OF LENGTH LAG+1.  CONTAINS WEIGHTS ON RETURN.
C            REAL*8
C
C     REMARK
C     THE PARAMETRIC QUANTITY ESTIMATED IS
C     F=(1/(2*PI))*(G(0)+2*SUM(G(I)*COS(I*FREQ):I=1,2,...))
C     WHERE THE G(I) ARE THE AUTOCOVARIANCES OF THE TIME SERIES.
C
C     REFERENCE
C     JENKINS, G. M. & WATTS, D. G.  SPECTRAL ANALYSIS AND ITS
C     APPLICATIONS.  SAN FRANCISCO: HOLDEN-DAY,INC. 1968.
C
C
      DOUBLE PRECISION FUNCTION SPEC(FREQ,R,LAG,ISW,WORK)
      IMPLICIT REAL*8 (A-H,O-Z)
      save
      REAL*8 R(1),WORK(1)
      PI=4.D0*DATAN(1.D0)
      L=LAG+1
      IF(ISW.EQ.1) GO TO 10
      IF(ISW.EQ.2) GO TO 20
      IF(ISW.EQ.4) GO TO 40
      IF(ISW.EQ.3) GO TO 30
      RETURN
C
C     RECTANGULAR WEIGHTS
C
10    DO 15 I=1,L
15    WORK(I)=1.D0
      GO TO 100
C
C     BARTLETT WEIGHTS
C
20    DO 25 I=1,L
25    WORK(I)=1.D0-DABS(DFLOAT(I-1))/DFLOAT(LAG)
      GO TO 100
C
C     TUKEY WEIGHTS
C
30    DO 35 I=1,L
35    WORK(I)=.5D0+.5D0*DCOS(PI*DFLOAT(I-1)/DFLOAT(LAG))
      GO TO 100
C
C     PARZEN WEIGHTS
C
40    M=LAG/2
      DO 45 I=1,L
      U=DFLOAT(I-1)/DFLOAT(LAG)
      IF(I-1.LE.M) WORK(I)=1.D0+6.D0*U*U*(U-1.D0)
45    IF(I-1.GT.M) WORK(I)=2.D0*((1.D0-U)**3)
      GO TO 100
C
C     COMPUTE ESTIMATED SPECTRUM
C
100   CONTINUE
      SUM=WORK(1)*R(1)
      IF(LAG.EQ.0) GO TO 120
      DO 110 I=1,LAG
110   SUM=SUM+2.D0*DCOS(FREQ*DFLOAT(I))*R(I+1)*WORK(I+1)
120   SPEC=SUM/(2.D0*PI)
      RETURN
      END