C
C *** WAVE READ, FOURIER TRANS., MULTIPLE REFLECTION METHOD, 
C     SPECTRUM PLOT, INVERSE FOURIER TRANS., AND WAVE PLOT
C
      PROGRAM MWVSHML
C
C PARAMETER
C *** MWV    : WAVE AMPLITUDE DATA IN MAIN PROGRAM DIMENSION
C            : IF MWV > 2**13(=8192), CAN NOT USE IN SUBROUTINE !!!
C READ DATA FROM FILE
C *** TITL   : TITLE OF EACH CHANNEL
C *** WV     : WAVE AMPLITUDE DATA
C *** NWV    : WAVE DATA NUMBER
C *** DT     : DELTA TIME ( SEC )
C DATA FOT PLOTTER
C *** FOUR(I): FOURIER COEFFICIENTS OR SPECTRUM
C            : FOURIER SPECTRUM = FOURIER COEFFICIENTS * DURATION/2.0
C
      PARAMETER(ML=30,MWV=9000)
      CHARACTER TITL*50,FNAME*30
      DIMENSION DNS(ML),DEP(ML),VS(ML),QF(ML),WV(MWV)
      COMPLEX*8 FOUR(MWV)
C
C *** WAVE DATA READ AND FOURIER COEFF. TRANS. ***
C
      CALL WVREAD(TITL,WV,NWV,MWV,DT)
      CALL SOILIN(FOM,DOM,NOM,DNS,VS,DEP,QF,NQF,NL,ML,THETA,NSF)
      CALL WVFOUR(WV,NWV,MWV,FOUR,NFAST)
C
      DUR=FLOAT(NFAST)*DT
      FOM=0.0
      DOM=6.2831853070/DUR
      NOM=NFAST/2+1
      CALL WVSHML(FOM,DOM,NOM,MWV,DNS,VS,DEP,QF,NQF
     *           ,NL,ML,THETA,NSF,FOUR)
C      CALL SPCOUT(FOM,DOM,NOM,MWV,FOUR)
C
C *** OUTPUT OF SPEC DATA TO FILE FOR LAZER SHOT ***
C
      WRITE(*,*)
   10 WRITE(*,1000)
      READ(*,'(A)')FNAME
      IF(FNAME(1:1).EQ.' ')GOTO 10
 1000 FORMAT( ' From MWVSHML : '
     *       /' Enter spec file name for Lazer Shot : ',\)
C      
      OPEN(8,FILE=FNAME)
        WRITE(8,1100)DOM,NOM-1
        WRITE(8,1110)
        DO 100 K=2,NOM,5
          WRITE(8,1120)(FOUR(K+J)*DUR,J=0,4)
  100   CONTINUE
      CLOSE(8)
C
 1100 FORMAT(F15.8,I15)
 1110 FORMAT(///)
 1120 FORMAT(10F8.2)
C
C *** TRANS. FROM FOURIER COEFF. TO AMPLITUDE SPECTRUM ***
C
      DO 200 I=1,NOM
      WV(I)=CABS(FOUR(I))*DUR
  200 CONTINUE
C
C *** SPECTRUM PLOTTER DATA OUTPUT ***
C
      CALL SPAPLT(FOM,DOM,NOM,MWV,WV)
C
C *** INVERSE FOURIER TRANSFORM ***
C
      DO 300 I=2,NFAST/2
      FOUR(NFAST+2-I)=CMPLX(REAL(FOUR(I)),-IMAG(FOUR(I)))
  300 CONTINUE
C
      CALL WVFAST(NFAST,FOUR,MWV,1)
C
      DO 400 I=1,NFAST
      WV(I)=REAL(FOUR(I))
  400 CONTINUE
C
      CALL WVOTS(WV,NWV,MWV,DT)
      CALL WVPLOT(TITL,WV,NWV,MWV,DT)
C
      STOP
      END
C
C
C
C
C *** OUTPUT OF WAVE DATA TO FILE ***
C
      SUBROUTINE WVOTS(WV,NWV,MWV,DT)
C
C WV(I) : WAVE DATA
C NWV   : NUMBER OF WAVE DATA
C MWV   : NWV IN MAIN
C DT    : DELTA TIME
C
      DIMENSION WV(MWV)
      CHARACTER FNAME*30
C
C *** FILE OUTPUT TO FILE DATA ***
C
      WRITE(*,1000)
      WRITE(*,1010)
      READ(*,*)IYES
      IF(IYES.NE.0)GOTO 20
 1000 FORMAT(/' From SUBROUTINE WVOT : wave data to file')
 1010 FORMAT(/' File out (yes=0)? :',\)
C
   10 WRITE(*,1020)
      READ(*,'(A)',ERR=10)FNAME
      IF(FNAME(1:1).EQ.' ')GOTO 10
 1020 FORMAT(/' Enter file name :',\)
C
      OPEN(7,FILE=FNAME,ERR=10)
      WRITE(7,2000)NWV,DT
      DO 100 IWV=1,NWV,8
      WRITE(7,2010)(WV((IWV-1)+I),I=1,8)
  100 CONTINUE
      CLOSE(7)
 2000 FORMAT(I5,3F10.4)
 2010 FORMAT(8F10.4)
C
   20 WRITE(*,*)
C
      RETURN
      END
