C  * * * * * * * * * * * * * * * * * * * * *                            
C     SUBROUTINE FOR PERIOD DISTRIBUTION                                
C  * * * * * * * * * * * * * * * * * * * * *                            
C                                                                       
C                      CODED BY Y.OHSAKI                                
C                                                                       
C     PURPOSE                                                           
C       TO COMPUTE THE PERIOD-FREQUENCY DISTRIBUTION OF A TIME SERIES   
C       OF EQUI-SPACED DATA BY MEANS OF ZERO-CROSSING OR PEAK METHOD    
C                                                                       
C     USAGE                                                             
C       CALL PERD(N,X,ND,DT,IND,T,NFREQ,RFREQ,EPS)                      
C                                                                       
C     DESCRIPTION OF ARGUMENTS                                          
C       N         - TOTAL NUMBER OF DATA                                
C       X(ND)     - EQUI-SPACED DATA                                    
C       ND        - DIMENSION OF X IN CALLING PROGRAM                   
C       DT        - TIME INCREMENT OF DATA IN SEC                       
C       IND       - IF IND.EQ.0, ZERO-CROSSING METHOD                   
C                   IF IND.EQ.1, PEAK METHOD                            
C       T(20)     - PERIODS IN SEC                                      
C       NFREQ(20) - FREQUENCY IN EACH CLASS OF PERIOD                   
C       RFREQ(20) - RELATIVE FREQUENCY IN EACH CLASS OF PERIOD IN       
C                   PERCENT                                             
C       EPS       - IRREGULARITY INDEX                                  
C                                                                       
C     SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED                     
C       NONE                                                            
C                                                                       
      SUBROUTINE PERD(N,X,ND,DT,IND,T,NFREQ,RFREQ,EPS)                  
C                                                                       
      DIMENSION X(ND),T(20),NFREQ(20),RFREQ(20)                         
      DIMENSION BOUND(21)                                               
      DATA      BOUND/0.05,0.06,0.07,0.08,0.10,0.12,0.15,0.18,0.22,     
     *                0.27,0.32,0.40,0.50,0.60,0.75,0.90,1.10,1.30,     
     *                1.60,2.00,2.50/                                   
C                                                                       
C     INITIALIZATION                                                    
C                                                                       
      DO 110 I=1,20                                                     
      T(I)=(BOUND(I)+BOUND(I+1))/2.                                     
      NFREQ(I)=0                                                        
  110 CONTINUE                                                          
      XMIN=99999.                                                       
      DO 120 M=1,N                                                      
      IF(X(M).EQ.0.) GO TO 120                                          
      XMIN=AMIN1(XMIN,ABS(X(M)))                                        
  120 CONTINUE                                                          
      ZERO=XMIN/1000.                                                   
      IF(X(1).EQ.0.0) X(2)=X(2)+ZERO                                    
      DO 130 M=2,N-1                                                    
      IF(ABS(X(M)-X(M+1)).GT.ZERO) GO TO 130                            
      X(M+1)=X(M)+SIGN(ZERO,X(M)-X(M-1))                                
  130 CONTINUE                                                          
      N0=0                                                              
      DO 140 M=1,N-1                                                    
      IF(X(M).EQ.0..OR.X(M)*X(M+1).LT.0.) GO TO 150                     
  140 CONTINUE                                                          
  150 TZ1=(REAL(M-1)+ABS(X(M)/(X(M)-X(M+1))))*DT                        
      NZ=1                                                              
      IF(X(M+1).GT.0.) N0=N0+1                                          
      MZ1=M+1                                                           
      DO 160 M=2,N-1                                                    
      IF(X(M)-X(M-1).LT.0..OR.X(M+1)-X(M).GT.0.) GO TO 160              
      GO TO 170                                                         
  160 CONTINUE                                                          
  170 TPP1=REAL(M-1)*DT                                                 
      NP=1                                                              
      MPP1=M+1                                                          
      IF(IND.EQ.0) GO TO 200                                            
      DO 180 M=2,N-1                                                    
      IF(X(M)-X(M-1).GT.0..OR.X(M+1)-X(M).LT.0.) GO TO 180              
      GO TO 190                                                         
  180 CONTINUE                                                          
  190 TPM1=REAL(M-1)*DT                                                 
      MPM1=M+1                                                          
C                                                                       
C     ZERO-CROSSING METHOD                                              
C                                                                       
  200 DO 260 M=MZ1,N                                                    
      IF(M.EQ.N) GO TO 210                                              
      IF(X(M)*X(M+1).GT.0.) GO TO 260                                   
      IF(X(M).EQ.0..AND.X(M-1)*X(M+1).GT.0..OR.X(M+1).EQ.0.) GO TO 260  
      TZ2=(REAL(M-1)+ABS(X(M)/(X(M)-X(M+1))))*DT                        
      GO TO 220                                                         
  210 IF(X(N).NE.0.) GO TO 260                                          
      TZ2=REAL(N-1)*DT                                                  
  220 TT=(TZ2-TZ1)*2.                                                   
      TZ1=TZ2                                                           
C                                                                       
C     FREQUENCY COUNT                                                   
C                                                                       
      IF(TT.LE.BOUND(1).OR.TT.GT.BOUND(21)) GO TO 260                   
      IF(IND.EQ.1) GO TO 250                                            
      DO 230 I=1,20                                                     
      IF(TT.GT.BOUND(I+1)) GO TO 230                                    
      NFREQ(I)=NFREQ(I)+1                                               
      GO TO 240                                                         
  230 CONTINUE                                                          
  240 NZ=NZ+1                                                           
  250 IF(X(M-1).LT.0.) N0=N0+1                                          
  260 CONTINUE                                                          
      TOTAL=REAL(NZ-1)                                                  
C                                                                       
C     PEAK METHOD                                                       
C                                                                       
      GRAD=1.                                                           
      MP1=MPP1                                                          
      TP1=TPP1                                                          
  270 DO 300 M=MP1,N-1                                                  
      IF((X(M)-X(M-1))*GRAD.LT.0..OR.(X(M+1)-X(M))*GRAD.GT.0.)          
     * GO TO 300                                                        
      TP2=REAL(M-1)*DT                                                  
      TT=TP2-TP1                                                        
      TP1=TP2                                                           
      IF(TT.LE.BOUND(1).OR.TT.GT.BOUND(21)) GO TO 300                   
      IF(IND.EQ.0) GO TO 290                                            
      DO 280 I=1,20                                                     
      IF(TT.GT.BOUND(I+1)) GO TO 280                                    
      NFREQ(I)=NFREQ(I)+1                                               
      GO TO 290                                                         
  280 CONTINUE                                                          
  290 NP=NP+1                                                           
  300 CONTINUE                                                          
      IF(GRAD.GT.0.) NM=NP                                              
      IF(IND.EQ.0) GO TO 310                                            
      GRAD=GRAD-2.                                                      
      MP1=MPM1                                                          
      TP1=TPM1                                                          
      IF(GRAD.GT.-2.) GO TO 270                                         
      TOTAL=REAL(NP-1)                                                  
C                                                                       
C     RELATIVE FREQ. AND IRREGULARITY INDEX                             
C                                                                       
  310 DO 320 I=1,20                                                     
      RFREQ(I)=REAL(NFREQ(I))/TOTAL*100.                                
  320 CONTINUE                                                          
      IF(NM.LE.N0) GO TO 330                                            
      EPS=SQRT(1.-(REAL(N0)/REAL(NM))**2)                               
      RETURN                                                            
  330 EPS=0.                                                            
      RETURN                                                            
      END                                                               
