C
C *** gramp.f is a FORTRAN code to obtain Fourier amplitudes,
C     which are computed from grpoint.f and grfault.f
C
      PROGRAM gramp
C
      PARAMETER(MWV=4100)
      CHARACTER INFILE1*20
      REAL TIME(MWV)
      REAL XR1(MWV),XI1(MWV),YR1(MWV),YI1(MWV),ZR1(MWV),ZI1(MWV)
      COMPLEX*8 FX1(MWV),FY1(MWV),FZ1(MWV)
C
      WRITE(*,*)
      WRITE(*,*)
      WRITE(*,*) 
     *'"grfft" - To obtain seismograms for grflt3.f and grpoint.f'
      WRITE(*,*)'           at a single observation point.'
      WRITE(*,*)
C
C * Read data *
C
      WRITE(*,1010) 
 1010 FORMAT('Choose the input datafile for FFT (Enter Number 1 to 6):'/ 
     *' 1) grfault.dxyz (x,y,z-comp. of disp. from seismic sources)'/
     *' 2) grfault.dqvw (u,v,w-comp. of disp. from seismic sources)'/
     *' 3) grpoint.dxyz (x,y,z-comp. of disp. from point sources)'/
     *' 4) grpoint.sxyz (x,y,z-comp. of stress from point sources)'/
     *' 5) grpoint.dqvw (u,v,w-comp. of disp.  from point sources)'/
     *' 6) grpoint.sqvw (u,v,w-comp. of stress from point sources)')
      READ(*,*) IDF
C    
      IF(IDF.EQ.1)THEN
        INFILE1='grfault.dxyz'
      ELSE 
        IF(IDF.EQ.2)THEN
          INFILE1='grfault.dqvw'
        ELSE 
          IF(IDF.EQ.3)THEN
            INFILE1='grpoint.dxyz'
          ELSE
            IF(IDF.EQ.4)THEN
              INFILE1='grpoint.sxyz'
            ELSE 
              IF(IDF.EQ.5)THEN
                INFILE1='grpoint.dqvw'
              ELSE 
                IF(IDF.EQ.6)THEN
                  INFILE1='grpoint.sqvw'
                END IF
              END IF
            END IF
          END IF
        END IF
      END IF
c
      OPEN(11,FILE=INFILE1,STATUS='OLD')
c
      READ(11,1100)DOM,IOM,NDATA
 1100 FORMAT(1X,F10.0,2I10)
      WRITE(*,*)
      WRITE(*,*)' DOM,IOM,NDATA =',DOM,IOM,NDATA
      NDATA=NDATA+1
      WRITE(*,*)' New NDATA including OM=0 =',NDATA
      READ(11,1110)NOB
 1110 FORMAT(1X,I10///)
      WRITE(*,*)
      WRITE(*,*)' Total number of receivers =',NOB
      WRITE(*,1120) 
 1120 FORMAT(1X,' Enter the receiver number to output : ')
      READ(*,*)NOB1
      WRITE(*,*)
      WRITE(*,*)' XR,XI,YR,YI,ZR,ZI='
      DO 150 J=1,NOB1
      DO 150 I = IOM+1,NDATA
	 READ(11,1130) XR1(I),XI1(I),YR1(I),YI1(I),ZR1(I),ZI1(I)
         IF(J.NE.NOB1)GOTO 150
c	 WRITE(*,1130) XR1(I),XI1(I),YR1(I),YI1(I),ZR1(I),ZI1(I)
  150 CONTINUE
 1130 FORMAT(1X,6E12.5)
C
      CLOSE(11)
C     
C * DIFFERENTIATION *
C
      WRITE(*,*)
      WRITE(*,*)
     * ' You can transfer the displacements to vel. or acc.;'
      WRITE(*,*)
     * '   Do not transfer                              (Enter 0),'
      WRITE(*,*)
     * '   Transfer from displacements to velocities    (Enter 1),'
      WRITE(*,*)
     * '   Transfer from displacements to accelerations (Enter 2):'
      READ(*,*)NYES
C     
C * BAND PASS *
C
      WRITE(*,*)
      WRITE(*,*)' Band Pass Filter? (Yes=1 or No=0) :'
      READ(*,*)NBD
C
      OML1=0.0
      OMH2=99999.9
      IF(NBD.EQ.1)THEN
        WRITE(*,*)'   Data for TRAPEZOIDAL BAND PASS FILTER:'
        WRITE(*,*)'   INITIAL AND FINAL OMs =',DOM*IOM,DOM*(NDATA-1)
        WRITE(*,*)
        WRITE(*,*)
     *  '   Enter OML1,OML2,OMH1 and OMH2 (OML1<OML2<OMH1<OMH2): '
        READ(*,*)OML1,OML2,OMH1,OMH2
      END IF
C
C * Write data to a complex array *
C
      FX1(1) = CMPLX(0.0,0.0)
      FY1(1) = CMPLX(0.0,0.0)
      FZ1(1) = CMPLX(0.0,0.0)
      DO 200 I = 2,NDATA
           OMEGA=FLOAT(I-1)*DOM
           IF(OMEGA.LE.OML1)GOTO 200
           IF(OMEGA.GE.OMH2)GOTO 200
           AA=1.0
           IF(NBD.EQ.1)THEN
             IF(OMEGA.LT.OML2)AA=(OMEGA-OML1)/(OML2-OML1)
             IF(OMEGA.GT.OMH1)AA=(OMH2-OMEGA)/(OMH2-OMH1)
           END IF
	   FX1(I) = AA*CMPLX(XR1(I),XI1(I))
	   FY1(I) = AA*CMPLX(YR1(I),YI1(I))
	   FZ1(I) = AA*CMPLX(ZR1(I),ZI1(I))
           IF(NYES.EQ.1)THEN
             FX1(I)=FX1(I)*OMEGA*(0.0,1.0)
             FY1(I)=FY1(I)*OMEGA*(0.0,1.0)
             FZ1(I)=FZ1(I)*OMEGA*(0.0,1.0)
           END IF
           IF(NYES.EQ.2)THEN
             FX1(I)=-FX1(I)*OMEGA*OMEGA
             FY1(I)=-FY1(I)*OMEGA*OMEGA
             FZ1(I)=-FZ1(I)*OMEGA*OMEGA
           END IF
  200 CONTINUE
C      
C * Output data file  *
C
      WRITE(*,*)
      WRITE(*,*)'Output data are stored in "amp_dat"'
      WRITE(*,*)'You can plot the waveform using "amp.m" for Matlab'
      WRITE(*,*)
      OPEN(11,FILE='amp_dat',STATUS='unknown')
      do i=1,ndata
      write(11,2000)dom*(i-1),fx1(I),fy1(I),fz1(I)
      end do
      CLOSE(11)
 2000 FORMAT(F7.4,6E12.4)
c
      STOP
      END
