C  * * * * * * * * * * * * * * * * * * * * * * * *                      
C     SUBROUTINE FOR SIMULATED WAVE GENERATION                          
C  * * * * * * * * * * * * * * * * * * * * * * * *                      
C                                                                       
C                            CODED BY Y.OHSAKI                          
C                                                                       
C     PURPOSE                                                           
C       TO GENERATE, FOR GIVEN MAGNITUDE AND EPICENTRAL DISTANCE OF     
C       AN EARTHQUAKE, AN ACCELERATION TIME-HISTORY CONSISTENT WITH     
C       OHSAKI'S VELOCITY RESPONSE SPECTRUM OF 5-PERCENT DAMPING        
C                                                                       
C     USAGE                                                             
C       CALL WGEN(EM,R,NN,IR,ACC,ND,DT,AMAX,VMAX,MXCYCL,ERR,UW1,UW2)    
C                                                                       
C     DESCRIPTION OF ARGUMENTS                                          
C       EM      - MAGNITUDE OF EARTHQUAKE                               
C       R       - EPICENTRAL DISTANCE IN KILOMETERS                     
C       NN      - TOTAL NUMBER OF DATA IN TIME-HISTORY                  
C       IR      - ANY INTEGER TO INITIATE THE SEQUENCE FOR RANDOM       
C                 NUMBER GENERATION                                     
C       ACC(ND) - ACCELERATION TIME-HISTORY IN GALS                     
C       ND      - DIMENSION OF ACC,UW1,UW2 IN CALLING PROGRAM           
C       DT      - TIME INTERVAL IN TIME-HISTORY IN SEC                  
C       AMAX    - MAX. ACCELERATION OF TIME-HISTORY IN GALS             
C       VMAX    - MAX. VELOCITY OF TIME-HISTORY IN KINES                
C       MXCYCL  - MAX. NUMBER OF ITERATION                              
C       ERR     - ROOT-MEAN-SQUARE ERROR FOR CONVERGENCE IN DECIMAL     
C                 FRACTION                                              
C       UW1(ND) - WORKING AREA                                          
C       UW2(ND) - WORKING AREA                                          
C                                                                       
C     REMARKS                                                           
C       NN MUST BE EQUAL TO POWER OF 2 NOT LARGER THAN 4096             
C                                                                       
C     SUBROUTINES AND FUNCTION SUBPROGRAM REQUIRED                      
C       ENVL OHSP VELK FAST CRAC IACC ERES RAND                         
C                                                                       
      SUBROUTINE WGEN(EM,R,NN,IR,ACC,ND,DT,AMAX,VMAX,MXCYCL,ERR,UW1,    
     *               UW2)                                               
C                                                                       
      COMPLEX   C(4096)                                                 
      DIMENSION ACC(ND),UW1(ND),UW2(ND)                                 
      DIMENSION E(33),X(33),EE(33)                                      
      DIMENSION PDIF(2046),PHI(2049),F(2049),T(2049),SV(2049),H(1),     
     *          RES(2049,1),RR(2049)                                    
      PARAMETER (PI2=6.283185)                                          
      DATA      DX/0.03125/,H0/0./,H/0.05/                              
C                                                                       
C     PHASE DIFFERENCES                                                 
C                                                                       
      CALL ENVL(EM,TB,TC,TD,33,E,33)                                    
      DT=TD/REAL(NN)                                                    
      NN2=NN/2                                                          
      NFOLD=NN2+1                                                       
      X(1)=0.                                                           
      EE(1)=0.                                                          
      DO 110 M=2,33                                                     
      X(M)=REAL(M-1)*DX                                                 
      EE(M)=EE(M-1)+E(M)                                                
  110 CONTINUE                                                          
      DO 120 M=2,33                                                     
      EE(M)=EE(M)/EE(33)                                                
  120 CONTINUE                                                          
      DO 150 K=1,NN2-2                                                  
      P=RAN2(IR)                                                        
      DO 130 J=2,33                                                     
      IF(P.LE.EE(J)) GO TO 140                                          
  130 CONTINUE                                                          
  140 PDIF(K)=-(X(J-1)+(P-EE(J-1))/(EE(J)-EE(J-1))*DX)*PI2              
  150 CONTINUE                                                          
C                                                                       
C     PHASE ANGLES                                                      
C                                                                       
      PHI(2)=0.                                                         
      DO 160 K=1,NN2-2                                                  
      PHI(K+2)=AMOD(PHI(K+1)+PDIF(K),PI2)                               
  160 CONTINUE                                                          
C                                                                       
C     TARGET SPECTRUM AND FIRST FOURIER APPROXIMATION                   
C                                                                       
      T(1)=TD*2.                                                        
      DO 170 K=2,NFOLD                                                  
      T(K)=TD/REAL(K-1)                                                 
  170 CONTINUE                                                          
      DO 180 K=2,NFOLD                                                  
      IF(T(K).LE.2.) GO TO 190                                          
  180 CONTINUE                                                          
  190 KMIN=K                                                            
      DO 200 K=KMIN,NFOLD                                               
      IF(T(K).LT.0.02) GO TO 210                                        
  200 CONTINUE                                                          
      KMAX=NFOLD                                                        
      GO TO 220                                                         
  210 KMAX=K-1                                                          
  220 NE=KMAX-KMIN+1                                                    
      DO 230 K=2,NFOLD                                                  
      CALL OHSP(EM,R,H0,T(K),SV0,K-2)                                   
      IF(K.LT.KMIN.OR.K.GT.KMAX) SV0=0.                                 
      F(K)=SV0/TD                                                       
  230 CONTINUE                                                          
      DO 240 K=2,NFOLD                                                  
      CALL OHSP(EM,R,H(1),T(K),SV(K),K-2)                               
      IF(K.LT.KMIN.OR.K.GT.KMAX) SV(K)=0.                               
  240 CONTINUE                                                          
C                                                                       
C     ITERATIVE COMPUTATION                                             
C                                                                       
      ENN=1./REAL(NN)                                                   
      NCYCL=0                                                           
  250 NCYCL=NCYCL+1                                                     
      C(1)=(0.,0.)                                                      
      DO 260 K=2,NN2                                                    
      C(K)=F(K)*CMPLX(COS(PHI(K)),SIN(PHI(K)))                          
      C(NN+2-K)=CONJG(C(K))                                             
  260 CONTINUE                                                          
      C(NFOLD)=F(NFOLD)*(1.,0.)                                         
      CALL FAST(NN,C,4096,+1)                                           
      AMAX=0.                                                           
      DO 270 M=1,NN                                                     
      ACC(M)=REAL(C(M))                                                 
      AMAX=AMAX1(AMAX,ABS(ACC(M)))                                      
  270 CONTINUE                                                          
C                                                                       
C     BASE LINE CORRECTION                                              
C                                                                       
      CALL CRAC(DT,NN,AMAX,ACC,ND,UW1,UW2)                              
      DO 280 M=1,NN                                                     
      C(M)=CMPLX(ACC(M),0.)                                             
  280 CONTINUE                                                          
      CALL FAST(NN,C,4096,-1)                                           
      DO 290 K=2,NFOLD                                                  
      F(K)=CABS(C(K))*ENN                                               
  290 CONTINUE                                                          
C                                                                       
C     MODIFICATION OF FOURIER AMPLITUDES                                
C                                                                       
      CALL ERES(1,H,1,NFOLD,T,2049,DT,NN,ACC,ND,2,VMAX,RES)             
      DO 300 K=2,NFOLD                                                  
      RR(K)=SV(K)/RES(K,1)                                              
      F(K)=F(K)*RR(K)                                                   
  300 CONTINUE                                                          
C                                                                       
C     ERROR FOR CONVERGENCE                                             
C                                                                       
      EPS=0.                                                            
      DO 310 K=KMIN,KMAX                                                
      EPS=EPS+(1.-RR(K))**2                                             
  310 CONTINUE                                                          
      EPS=SQRT(EPS/REAL(NE))                                            
C                                                                       
      IF(EPS.LE.ERR) GO TO 320                                          
      IF(NCYCL.EQ.MXCYCL) GO TO 330                                     
      GO TO 250                                                         
  320 MXCYCL=NCYCL                                                      
  330 ERR=EPS                                                           
      RETURN                                                            
      END                                                               
