C CHUONG TRINH NAY TINH PHO ENTROPI CUC DAI 
C CHO SO LIEU CACH NHAU MOT NAM (IYEAR=1) HOAC SL THANG (IYEAR=0)
C MODIFY ON 6/5/2000
      PROGRAM MEMSP
      REAL PM
        REAL DELTA_S
      INTEGER N,M,MM,MP,IYEAR, MONTH(12)   ! BIEN MONTH DE XAC DINH SE TINH CHO THANG NAO
      REAL TB(12)
C
      Character *50 FInp,FOut,IOFile
        REAL,ALLOCATABLE :: DAT1(:),COF(:),S(:),FREQ(:)
      REAL,ALLOCATABLE :: WK1(:),WK2(:),WKM(:)
      REAL,ALLOCATABLE :: DAT(:,:)       ! MANG CHUA SO LIEU BAN DAU
      INCLUDE 'MEMSP.PAR'
        ALLOCATE (DAT1(NMAX),COF(NMAX),S(NMAX),FREQ(NMAX))
      ALLOCATE (WK1(NMAX),WK2(NMAX),WKM(NMAX))
      
      MM = M  
C     CAN PHAI CO FILE SO THAM SO MEMSP.IO CHUA TEN FILE SO LIEU,...
C     IOFile='MEMSP.IOT':NHIET DO, ='MEMSP.IOM':MUA
      OPEN(2,FILE=IOFile,STATUS='OLD')  
C       
      IF (IYEAR==1) THEN   ! TINH VOI CHUOI SL NAM HOAC MUA     
C  -----------------------------------------------------------
 8888 READ(2,*,END=9999)
      READ(2,*)FInp
      WRITE(*,*)' FILE = ',FInp
      READ(2,*)FOut
C                  
      OPEN(1,FILE=FInp,STATUS='OLD')
      OPEN(3,FILE=FOut,STATUS='UNKNOWN')
      read(1,*)LPS
C      N=LPS*12 
        ALLOCATE (DAT(LPS,12))       ! MANG CHUA SO LIEU BAN DAU 
      
      DO I=1,LPS
         READ(1,*) NAM, (DAT(I,J),J=1,12)
      ENDDO
      N=LPS           
C                     GUI SO LIEU THANG MONTH SANG MANG DATA
      JMONTH=0
      DO J=1,12
         IF (MONTH(J).NE.0) JMONTH = JMONTH + 1
      ENDDO
      DO I=1,N
         DAT1(I) = 0.0
         DO J=1,12
            IF (MONTH(J).NE.0) DAT1(I) = DAT1(I) + DAT(I,J)
         ENDDO                
      IF (IYTO.EQ.2) DAT1(I) = DAT1(I)/FLOAT(JMONTH)
      ENDDO
            
C     NEU M=0 THI GNA M=N/4      
C      IF (MM.EQ.0) M = N/4 
C
C     GAN CO DINH M=N/5
      M = N/4      

C     QUI TAM SO LIEU
      TB1=0.0
      DO I=1,N 
          TB1=TB1+DAT1(I)
      ENDDO
      TB1=TB1/FLOAT(N)
      DO I=1,N
          DAT1(I)=DAT1(I)-TB1
      ENDDO          
C     CHUAN HOA SO LIEU
      PS=0.0
      DO I=1,N
          PS=PS+DAT1(I)**2
      ENDDO
      PS=SQRT(PS/FLOAT(N))
      DO I=1,N
          DAT1(I)=DAT1(I)/PS
      ENDDO
               
C       
      CALL MEMCOF(DAT1,N,M,PM,COF,WK1,WK2,WKM)
C                         
C     SO GIA TRI TAN SO VA MAT DO PHO DUOC TINH TU BUOC TAN SO = 0.0005
C      DELTA_S = 0.005
C       MP = 0.5/DELTA_S
C      MP = N
      
      Stb = 0.0
      DO I=1,MP+1
          FREQ(I)=(I-1)*DELTA_S
          S(I)=EVLMEM(FREQ(I),COF,M,PM) 
          Stb = Stb + S(I)
      ENDDO               

      Stb = Stb/FLOAT(MP+1)          
      DgrFree = (2*N-0.5*M)/FLOAT(M)
      IF (JMONTH.EQ.12) THEN
         WRITE(3,*)' TINH CHO CA NAM'
      ELSE
         WRITE(3,*)' TINH CHO THANG '
         WRITE(3,'(12I3)')(MONTH(J),J=1,12)
      ENDIF   
      WRITE(3,*) '    N     ','    M     ','    L     ','  Stb'
      WRITE(3,'(2I10,F10.2,F10.6)') N, M, DgrFree, Stb
      WRITE(3,*) '  Tan_So  ','Chuky_Nam ','  SPEC'
      DO I=1,MP+1
          IF (I.NE.1) THEN 
              PRI=1/FREQ(I)               
          ELSE
              PRI=999
          ENDIF    
          WRITE(3,'(F10.5,F10.4,F10.4)') FREQ(I),PRI,S(I)
      ENDDO      
      CLOSE(1)
      CLOSE(3)
C
      DEALLOCATE (DAT)
        GOTO 8888      
 9999 CONTINUE                 
C
C -----------------------------------------------------------
        ELSE  ! NEU KHONG THI TINH THEO CHUOI SL THANG
C                   
 7777 READ(2,*,END=6666)
      READ(2,*)FInp
      WRITE(*,*)' FILE = ',FInp
      READ(2,*)FOut
C                  
      OPEN(1,FILE=FInp,STATUS='OLD')
      OPEN(3,FILE=FOut,STATUS='UNKNOWN')
      read(1,*)LPS
C      N=LPS*12 
        ALLOCATE (DAT(LPS,12))       ! MANG CHUA SO LIEU BAN DAU 
      
      DO I=1,LPS
         READ(1,*) NAM, (DAT(I,J),J=1,12)
      ENDDO
      N=LPS 
C  TINH BIEN TRINH TRUNG BINH NAM       
      DO J=1,12
                TB(J)=0.0
            DO I=1,N
                        TB(J)=TB(J)+DAT(I,J)
            ENDDO
            TB(J)=TB(J)/FLOAT(N)
        ENDDO
C  QUI TAM SO LIEU DE LOC BIEN TRINH NAM
        DO I=1,N
                DO J=1,12
                    DAT(I,J)=DAT(I,J)-TB(J)
            END DO
        ENDDO
C  CHUYEN THANH MANG MOT CHIEU
        K=0
        DO I=1,N
                DO J=1,12
                        K=K+1
                        DAT1(K)=DAT(I,J)
                ENDDO
        ENDDO 
      N=N*12      

C     QUI TAM SO LIEU LAN 2 
      TB1=0.0
      DO I=1,N 
          TB1=TB1+DAT1(I)
      ENDDO
      TB1=TB1/FLOAT(N)
      DO I=1,N
          DAT1(I)=DAT1(I)-TB1
      ENDDO          
C     CHUAN HOA SO LIEU
      PS=0.0
      DO I=1,N
          PS=PS+DAT1(I)**2
      ENDDO
      PS=SQRT(PS/FLOAT(N))
      DO I=1,N
          DAT1(I)=DAT1(I)/PS
      ENDDO
               
C       
C     NEU M=0 THI GNA M=N/4      
C      IF (MM.EQ.0) M = N/4 
C
C     GAN CO DINH M=N/5
      M = N/4      
      CALL MEMCOF(DAT1,N,M,PM,COF,WK1,WK2,WKM)
C                         
C     SO GIA TRI TAN SO VA MAT DO PHO DUOC TINH TU BUOC TAN SO = 0.0005
C      DELTA_S = 0.005/12.0
C       MP = 0.5/DELTA_S
C      MP = N
      
      Stb = 0.0
      DO I=1,MP+1
          FREQ(I)=(I-1)*DELTA_S
          S(I)=EVLMEM(FREQ(I),COF,M,PM) 
          Stb = Stb + S(I)
      ENDDO               

      Stb = Stb/FLOAT(MP+1)          
      DgrFree = (2*N-0.5*M)/FLOAT(M)
      
        WRITE(3,*)' TINH THEO SO LIEU THANG'
      
      WRITE(3,*) '    N     ','    M     ','    L     ','  Stb'
      WRITE(3,'(2I10,F10.2,F10.6)') N, M, DgrFree, Stb
      WRITE(3,*) '  Tan_So  ','Chuky_Nam ','  SPEC'
      DO I=1,MP+1
          IF (I.NE.1) THEN 
              PRI=1./FREQ(I)/12.
          ELSE
              PRI=999
          ENDIF    
          WRITE(3,'(F10.5,F10.4,F10.4)') FREQ(I),PRI,S(I)
      ENDDO      
      CLOSE(1)
      CLOSE(3)
C
      DEALLOCATE (DAT)
        GOTO 7777      
 6666 CONTINUE      
        END IF  ! IYEAR=1           
C  ----------------------------------------------------
      STOP
      END   
C      
C   
               

      SUBROUTINE MEMCOF(DATA1,N,M,PM,COF,WK1,WK2,WKM)
      DIMENSION DATA1(N),COF(M),WK1(N),WK2(N),WKM(M)
C
      P=0.0
      DO 11 J=1,N
          P=P+DATA1(J)**2
   11 CONTINUE
      PM=P/N
      WK1(1)=DATA1(1)
      WK2(N-1)=DATA1(N)
C      
      DO 12 J=2,N-1
          WK1(J)=DATA1(J)
          WK2(J-1)=DATA1(J)
   12 CONTINUE
C
      DO 17 K=1,M  
          PNEUM=0.0
          DENOM=0.0
C
          DO 13 J=1,N-K
              PNEUM=PNEUM+WK1(J)*WK2(J)
              DENOM=DENOM+WK1(J)**2+WK2(J)**2
   13     CONTINUE
          COF(K)=2.0*PNEUM/DENOM
          PM=PM*(1.0-COF(K)**2)
C
          DO 14 I=1,K-1
              COF(I)=WKM(I)-COF(K)*WKM(K-I)
   14     CONTINUE
          IF (K.EQ.M) RETURN
C
          DO 15 II=1,K
              WKM(II)=COF(II)
   15     CONTINUE
C
          DO 16 JJ=1,N-K-1
              WK1(JJ)=WK1(JJ)-WKM(K)*WK2(JJ)
              WK2(JJ)=WK2(JJ+1)-WKM(K)*WK1(JJ+1)
   16     CONTINUE
   17 CONTINUE
      END
C      
C
C
      FUNCTION EVLMEM(FDT,COF,M,PM)
      DIMENSION COF(M)
      REAL*8 WR,WI,WPR,WPI,WTEMP,THETA
      THETA=6.28318530717959D0*FDT
      WPR=DCOS(THETA)
      WPI=DSIN(THETA)
      WR=1.D0
      WI=0.D0
      SUMR=1.0
      SUMI=0.0
      DO 11 I=1,M
          WTEMP=WR
          WR=WR*WPR-WI*WPI
          WI=WI*WPR+WTEMP*WPI
          SUMR=SUMR-COF(I)*SNGL(WR)
          SUMI=SUMI-COF(I)*SNGL(WI)
   11 CONTINUE
      EVLMEM=PM/(SUMR**2+SUMI**2)
      RETURN
      END