PROGRAM xevlmem
!     driver for routine evlmem
      INTEGER N,NFDT,M
!      PARAMETER(NN=2000,M=20,NFDT=512)
      PARAMETER(NN=2000,MM=300,NFDT=256)
!         PARAMETER(N=1000,M=10,NFDT=16)
      INTEGER i
      REAL evlmem,fdt,pm,dat(NN),cof(MM),X(NN),TEMP(120,12),TB(12)
          REAL OMEGA(5)
          DATA OMEGA /0.06,0.25,0.38,0.25,0.06/
      open(7,file='v.M',status='old')
!      read(7,*) (data(i),i=1,N)
      read(7,*)LPS
          N=1
  100 READ(7,*,END=101) NAM, (TEMP(N,J),J=1,12)
          N=N+1
          GOTO 100
  101 CONTINUE
          IF ((N-1).NE.LPS) THEN
             WRITE(*,*)' WRONG IN READING DATA FILE: N=',N,' LPS=',LPS
                 STOP
          ENDIF
      close(7)
! LOC BIEN TRINH NAM
          DO 102 J=1,12
                TB(J)=0.0
                DO I=1,LPS
                  TB(J)=TB(J)+TEMP(I,J)
            ENDDO
            TB(J)=TB(J)/FLOAT(LPS)
                DO I=1,LPS
                  TEMP(I,J)=(TEMP(I,J)-TB(J))/10.0
                ENDDO
  102 CONTINUE

          K=0
          DO I=1,LPS
                 DO J=1,12
            K=K+1
                    DAT(K)=TEMP(I,J)
         ENDDO
      ENDDO
          N=LPS*12 
!  LOC TRONG LUONG, HOAC TRUNG BINH TRUOT
!         J=0
!         DO I=3, N-2
!           J=J+1
!               X(J)=0.0
!               DO K=I-2,I+2
!                 X(J)=X(J)+DAT(K)
!               ENDDO
!               X(J)=X(J)/5.0
!         ENDDO
!         DAT=X
!         N=N-4

!     CHUAN HOA SO LIEU
      TB1=0.0
          PS=0.0
     DO I=1,N
          TB1=TB1+DAT(I)
                  PS=PS+DAT(I)**2
      ENDDO
          TB1=TB1/FLOAT(N)
      PS=SQRT(PS/FLOAT(N)-TB1*TB1)
      DO I=1,N
          DAT(I)=DAT(I)/PS
      ENDDO
      
!         M = N / 10
          M=30
      call memcof(dat,N,M,pm,cof)

          L = (2*N - 0.5*M)/M     ! SO BAC TU DO
!       TINH TRUNG BINH CUA M GIA TRI PHO 
          Stb = 0.0
      do 12 i=1,M
        fdt=0.5*i/M
!               FDT=1.0/FDT/12.0
                Stb = Stb + evlmem(fdt,cof,M,pm)
 12   continue
          STB=STB/FLOAT(M)

          OPEN(3,FILE='v.DAT',STATUS='UNKNOWN')
      write(3,*) 'Power spectrum estimate of DATA'
          WRITE(3,'(1X,3I12,F12.6)') N,M,L, STB
           
      write(3,'(1x,2(A12))') 'f*delta','power'
!      do 11 i=0,NFDT
      do 11 i=1,NFDT
        fdt=0.5*i/NFDT
                S = evlmem(fdt,cof,M,pm)
                FDT=1.0/FDT/12.0
        write(3,'(1x,2f12.6)') fdt, S
 11   continue
          CLOSE(3)
      END
! END MAIN PROGRAM
! ================

! THE FOLLOWING SUBROUTINES SHOULD NOT BE CHANGED

      SUBROUTINE memcof(data,n,m,xms,d)
      INTEGER m,n,MMAX,NMAX
      REAL xms,d(m),data(n)
      PARAMETER (MMAX=600,NMAX=2000)
      INTEGER i,j,k
      REAL denom,p,pneum,wk1(NMAX),wk2(NMAX),wkm(MMAX)
      if (m.gt.MMAX.or.n.gt.NMAX) pause 'workspace too small in memcof'
      p=0.
      do 11 j=1,n
        p=p+data(j)**2
11    continue
      xms=p/n
      wk1(1)=data(1)
      wk2(n-1)=data(n)
      do 12 j=2,n-1
        wk1(j)=data(j)
        wk2(j-1)=data(j)
12    continue
      do 17 k=1,m
        pneum=0.
        denom=0.
        do 13 j=1,n-k

          pneum=pneum+wk1(j)*wk2(j)
          denom=denom+wk1(j)**2+wk2(j)**2
13      continue
        d(k)=2.*pneum/denom
        xms=xms*(1.-d(k)**2)
        do 14 i=1,k-1
          d(i)=wkm(i)-d(k)*wkm(k-i)
14      continue
        if(k.eq.m)return
        do 15 i=1,k
          wkm(i)=d(i)
15      continue
        do 16 j=1,n-k-1
          wk1(j)=wk1(j)-wkm(k)*wk2(j)
          wk2(j)=wk2(j+1)-wkm(k)*wk1(j+1)
16      continue
17    continue
      pause 'never get here in memcof'
      END

      FUNCTION evlmem(fdt,d,m,xms)
      INTEGER m
      REAL evlmem,fdt,xms,d(m)
      INTEGER i
      REAL sumi,sumr
      DOUBLE PRECISION theta,wi,wpi,wpr,wr,wtemp
      theta=6.28318530717959d0*fdt
      wpr=cos(theta)
      wpi=sin(theta)
      wr=1.d0
      wi=0.d0
      sumr=1.
      sumi=0.
      do 11 i=1,m
        wtemp=wr
        wr=wr*wpr-wi*wpi
        wi=wi*wpr+wtemp*wpi
        sumr=sumr-d(i)*sngl(wr)
        sumi=sumi-d(i)*sngl(wi)
11    continue
      evlmem=xms/(sumr**2+sumi**2)
      return
      END