C  * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *        
C     SUBROUTINE FOR SMOOTHED SPECTRA BY PARZEN'S LAG WINDOW            
C  * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *        
C                                                                       
C                                          CODED BY Y.OHSAKI            
C                                                                       
C     PURPOSE                                                           
C       TO COMPUTE FOURIER AND POWER SPECTRA OF A SERIES OF EQUI-       
C       SPACED DATA SMOOTHED BY APPLICATION OF PARZEN'S LAG WINDOW      
C       WITH SPECIFIED BANDWIDTH IN FREQUENCY DOMAIN                    
C                                                                       
C     USAGE                                                             
C       CALL LWIN(N,X,ND1,DT,F,G,ND2,NFOLD,DF,BAND)                     
C                                                                       
C     DESCRIPTION OF ARGUMENTS                                          
C       N      - TOTAL NUMBER OF DATA                                   
C       X(ND1) - EQUI-SPACED DATA                                       
C       ND1    - DIMENSION OF X IN CALLING PROGRAM  ND1.LE.8192         
C       DT     - TIME INCREMENT IN DATA IN SEC                          
C       F(ND2) - SMOOTHED FOURIER SPECTRUM                              
C       G(ND2) - SMOOTHED POWER SPECTRUM                                
C       ND2    - DIMENSION OF F,G IN CALLING PROGRAM                    
C       NFOLD  - NUMBER OF DATA IN THE SPECTRA                          
C       DF     - FREQUENCY INCREMENT IN THE SPECTRA IN HZ               
C       BAND   - SPECIFIED BANDWIDTH IN HZ                              
C                                                                       
C     REMARKS                                                           
C      (1) ND2.GE.NT/2+1, WHERE NT IS THE POWER OF 2 EQUAL TO N OR      
C          MINIMUM LARGER THAN N                                        
C      (2) BAND/DF MUST BE LARGER THAN 560/151 (APPROXIMATELY 3.71)     
C      (3) IF THE ARGUMENT BAND=0., NO SMOOTHING IS MADE                
C                                                                       
C     SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED                     
C       FAST                                                            
C                                                                       
      SUBROUTINE LWIN(N,X,ND1,DT,F,G,ND2,NFOLD,DF,BAND)                 
C                                                                       
      COMPLEX   A(8192)                                                 
      DIMENSION X(ND1),F(ND2),G(ND2)                                    
C                                                                       
C     INITIALIZATION                                                    
C                                                                       
      DO 110 M=1,N                                                      
      A(M)=CMPLX(X(M),0.)                                               
  110 CONTINUE                                                          
      NT=2                                                              
  120 IF(NT.GE.N) GO TO 130                                             
      NT=NT*2                                                           
      GO TO 120                                                         
  130 IF(NT.EQ.N) GO TO 150                                             
      DO 140 M=N+1,NT                                                   
      A(M)=(0.,0.)                                                      
  140 CONTINUE                                                          
  150 NFOLD=NT/2+1                                                      
      T=REAL(NT)*DT                                                     
      DF=1./T                                                           
C                                                                       
C     AUTOCORRELATION                                                   
C                                                                       
      CALL FAST(NT,A,8192,-1)                                           
      IF(BAND.NE.0.) GO TO 170                                          
      DO 160 K=1,NFOLD                                                  
      F(K)=CABS(A(K))*DT                                                
  160 CONTINUE                                                          
      GO TO 240                                                         
  170 DO 180 K=1,NT                                                     
      A(K)=A(K)*CONJG(A(K))/REAL(NT)**2                                 
  180 CONTINUE                                                          
      CALL FAST(NT,A,8192,+1)                                           
C                                                                       
C     LAG WINDOW                                                        
C                                                                       
      U=3.708609/BAND*DF                                                
      IF(U.GT.1.) GO TO 260                                             
      U=REAL(NFOLD-1)*U                                                 
      DO 210 J=2,NFOLD                                                  
      TAU=REAL(J-1)/U                                                   
      IF(TAU.GT.0.5) GO TO 190                                          
      A(J)=A(J)*(1.-6.*TAU**2*(1.-TAU))                                 
      GO TO 210                                                         
  190 IF(TAU.GT.1.) GO TO 200                                           
      A(J)=A(J)*2.*(1.-TAU)**3                                          
      GO TO 210                                                         
  200 A(J)=(0.,0.)                                                      
  210 CONTINUE                                                          
      DO 220 J=2,NFOLD-1                                                
      A(NT+2-J)=A(J)                                                    
  220 CONTINUE                                                          
C                                                                       
C     SMOOTHED SPECTRA                                                  
C                                                                       
      CALL FAST(NT,A,8192,-1)                                           
      DO 230 K=1,NFOLD                                                  
      F(K)=SQRT(ABS(REAL(A(K)))*REAL(NT))*DT                            
  230 CONTINUE                                                          
  240 G(1)=F(1)**2/T                                                    
      DO 250 K=2,NFOLD-1                                                
      G(K)=2.*F(K)**2/T                                                 
  250 CONTINUE                                                          
      G(NFOLD)=F(NFOLD)**2/T                                            
      RETURN                                                            
C                                                                       
  260 WRITE(6,601)                                                      
      STOP                                                              
C                                                                       
C     FORMAT STATEMENT                                                  
C                                                                       
  601 FORMAT('BANDWIDTH IS TOO NARROW')                                 
      END                                                               
