C 
C Program SS (Synthetic Seismograms) to read or generate and filter the
C source time function.  It may store the filtered source time function
C and its Hilbert transform in the GSE data exchange format, or read the
C frequency-domain response function and generate synthetic seismograms
C in the GSE data exchange format.
C
C Version: 5.10
C Date: 1997, September 26
C
C Coded by: Ludek Klimes
C     Department of Geophysics, Charles University Prague,
C     Ke Karlovu 3, 121 16 Praha 2, Czech Republic,
C     E-mail: klimes@seis.karlov.mff.cuni.cz
C
C.......................................................................
C
C This Fortran77 file consists of the following external procedures:
C     SS...   Main program.
C             SS
C     SIGNAL..Subroutine to generate Gabor's or Mueller's signal of
C             given parameters.
C             SIGNAL
C     PLSIG...Subroutine to determine the maximum amplitude of the given
C             function, detect zeros beyond and behind the signal, write
C             the parameters of the signal to the output log file, and
C             eventually plot the signal.
C             PLSIG
C     PLOPN...Simple subroutine to initiate plotting.
C             PLOPN
C     PLTIM...Simple subroutine to supplement the signal plots with the
C             time labels.
C             PLTIM
C     FCOOLR..Subroutine for the fast Fourier transform of N=2**K
C             complex data points.
C             FCOOLR
C
C Other external procedures required:
C     WGSE1,WGSE2,WGSE3... Subroutines of the Fortran 77 file 'gse.for'
C             (package MODEL), designed to write seismograms in the GSE
C             data exchange format.
C     PLOTS,PLOT,SYMBOL,NUMBER... CALCOMP plotting subroutines.  For
C             example, Fortran 77 routines of file 'calcops.for'
C             (package MODEL) may be used to generate seismogram plots
C             in the PostScript files.
C
C.......................................................................
C
C                                                    
C Description of the data files:
C
C The data are read in by the list directed input (free format).
C In the description of data files, each numbered paragraph indicates
C the beginning of a new input operation (new READ statement).
C If the symbolic name of the input variable is enclosed in apostrophes,
C the corresponding value in input data is of the type CHARACTER, i.e.
C it should be a character string enclosed in apostrophes.  If the first
C letter of the symbolic name is I-N, the corresponding value is of the
C type INTEGER.  Otherwise, the input parameter is of the type REAL and
C may or may not contain a decimal point.
C
C Input data read from the * external unit:
C     The interactive * external unit may also be redirected to the file
C     containing the relevant data.
C (1) 'SSDAT','RF','SSGSE','SSLOG',KCOMP,/
C     'SSDAT'... String with the name of the input data file containing
C             the data describing the parameters of the source time
C             function and frequency-domain filter.
C             Description of file SSDAT
C     'RF'... String with the name of the input data file with the
C             frequency-domain response functions at individual
C             receivers.  Must be specified for positive KCOMP, but is
C             not taken into account if KCOMP=0.
C             Description of file RF
C     'SSGSE'... String with the name of the output data file in the GSE
C             data exchange format, containing the seismograms at the
C             receivers (for KCOMP positive), or the filtered source
C             time function and its Hilbert transform (if KCOMP=0).
C             Description of file SSGSE
C     'SSLOG'... String with the name of the output log file.
C             Description of file SSLOG
C     KCOMP...KCOMP=0: Output file 'SSGSE' will contain the filtered
C               source time function and its Hilbert transform, no
C               synthetic seismograms are generated.
C             KCOMP=1: Output file 'SSGSE' will contain the X1 component
C               of the synthetic seismograms.
C             KCOMP=2: Output file 'SSGSE' will contain the X2 component
C               of the synthetic seismograms.
C             KCOMP=3: Output file 'SSGSE' will contain the X3 component
C               of the synthetic seismograms.
C     Default: 'SSDAT'='ss.dat', 'RF'='rf.out', 'SSGSE'='ss.gse',
C             'SSLOG'='sslog.out', KCOMP=0.
C
C                                                   
C Input file SSDAT with the data describing the source time function:
C (1) KSGNL,KEXP,MPTS,NPTS,DER,HILB,/
C     KSGNL...Type of the source time function:
C             KSGNL=0: Digitized time function specified point by point.
C             KSGNL=1: Gabor signal.
C             KSGNL=2: Mueller signal.
C     KEXP... Exponent controlling the frequency-domain filter, see (2).
C     MPTS... Number of points of the time functions at the output check
C             plot.  The corresponding MPTS-1 time intervals are all
C             together scaled to 10.23cm.  The signal is approximately
C             centred.  MPTS does not apply to the spectra at the output
C             check plot, NPTS/2 points of each spectrum is plotted.
C             MPTS=0: No output check plot is generated.
C     NPTS... Number of the time samples for the fast Fourier transform.
C             NPTS is rounded up to the nearest power of 2.
C     DER,HILB... The source time function is DER-th derivative and
C             HILB-th Hilbert transform of the given signal.
C     Default: KEXP=1, MPTS=0, NPTS=MSS, DER=0, HILB=0,
C             where MSS is the array dimension declared in the code.
C (2) FMIN,FL,FR,FMAX,TL,TD,TRED,VRED,/
C     FMIN,FL,FR,FMAX... Parameters of the frequency filter to be
C             applied to the source time function.  The filter is zero
C             for frequencies F smaller than FMIN or greater than FMAX.
C             The filter is 1 between FL and FR.
C             Between FMIN and FL, cosine tappering
C               ( 0.5-0.5*cos(pi*(F-FMIN)/(FL-FMIN) )**KEXP
C             is used.
C             Between FR and FMAX, cosine tappering
C               ( 0.5-0.5*cos(pi*(F-FMAX)/(FR-FMAX) )**KEXP
C             is used.
C     TL...   Reference time of the given signal.
C     TD...   Time interval to digitize the source time function and
C             seismograms.
C     TRED,VRED... Specification of the time window for synthetic
C             seismograms.
C             VRED.EQ.0: Seismogram is centred in the time interval of
C               length (NPTS-1)*TD according to the travel times of the
C               first and last arrivals at the receiver.
C             VRED.NE.0: Time of the first sample of the time window is
C               T=TRED+R/VRED, where R is the hypocentral distance.
C     Default: FMIN=0, FL=0, FR=FMAX, FMAX=0.5/TD, TL=0, TRED=0, VRED=0.
C (3) PAR1,PAR2,PAR3,PAR4,...,/
C     Parameters of the given signal.
C     KSGNL=0: Digitized time function specified point by point:
C             PAR1,PAR2,PAR3,PAR4,PAR5,PAR6,... ...time function
C               digitized with step TD.  PAR1 corresponds to time TL.
C     KSGNL=1: Gabor signal:
C               S(t)=PAR4*exp(-(2*pi*PAR1*(t-TL)/PAR2)**2)
C                        *cos(2*pi*PAR1*(t-TL)+PAR3)
C             PAR1... Prevailing frequency.
C             PAR2... Relative width of the signal.
C             PAR3... Phase.
C             PAR4... Amplitude of the envelope.
C     KSGNL=2: Mueller signal:
C               For 0.LE.(t-TL).LE.PAR2/(2*PAR1):
C               S(t)=A*(sin(2*pi*PAR1*(t-TL))
C                      -sin(2*pi*PAR1*(t-TL)*(PAR2+2)/PAR2)
C                        * PAR2 / (PAR2+2)
C                      -(1-cos(2*pi*PAR1*(t-TL)/PAR2))
C                        * sin(pi*PAR2) / (PAR2+2) )
C               Otherwise:
C               S(t)=0
C             PAR1... Reference frequency.
C             PAR2... Relative width of the signal.  If PAR2 is integer,
C                     the number of local extrema of the signal.
C             PAR3... No meaning.
C             PAR4... Maximum amplitude of the signal.  Determines A.
C Example of data SSDAT
C
C                                                      
C Input file 'RF' with the response functions:
C (1) 'TEXT1',/
C     'TEXT1'... Text string in apostrophes.  The first 34 characters
C             will be passed to the header of the output GSE file.
C (2) XS1,XS2,XS3,/
C     XS1,XS2,XS3... Coordinates of the hypocentre.
C (3) FMINIM,FD,NF,/
C     FMINIM..The lowest frequency at the digitized response function.
C     FD...   The frequency step.
C     NF...   The number of discrete frequencies.
C (4) For each receiver (4.1) and (4.2):
C (4.1) X1,X2,X3,TMIN,TMAX,AMAX,/
C     X1,X2,X3... Coordinates of the receiver.
C     TMIN,TMAX...Travel times of the first and last arrivals at the
C             receiver.
C     AMAX... Maximum absolute value over the real and imaginary part of
C             all 3 components of the response function.
C (4.2) Written only if TMIN.LE.TMAX (Attention: Formatted input!):
C     6*NF integer numbers, FORMAT(12I6):
C     (IPR(I,IF),I=1,6,IF=1,NF)
C     IPR...  IPR(I,IF)=IFIX(99999.1*SPR(I,IF)/AMAX), where
C             SPR(*,IF) is the response function at the frequency
C               F=FMIN+(IF-1)*FD.
C             SPR(1,IF) is the real part of the X1 component.
C             SPR(2,IF) is the imaginary part of the X1 component.
C             SPR(3,IF) is the real part of the X2 component.
C             SPR(4,IF) is the imaginary part of the X2 component.
C             SPR(5,IF) is the real part of the X3 component.
C             SPR(6,IF) is the imaginary part of the X3 component.
C (5) / or end of file.
C
C                                                   
C Output file 'SSGSE' with the seismograms or source time function:
C     File in the GSE data exchange format, see the description in file
C     'gse.for'.
C     Description of format GSE
C
C                                                   
C Output log file 'SSLOG':
C     Contains information on the input data, source time function,
C     synthetic seismograms.  Often may be discarded.
C
C                                                    
C Output check plot:
C     If MPTS.GT.0, contains plots of:
C             1 Input signal,
C             2 Spectrum of the input signal,
C             3 Spectrum of the filter,
C             4 Spectrum of the filtered signal,
C             5 Filtered  signal,
C             6 Hilbert transform of the filtered signal.
C     If KCOMP.GT.0, contains plots of the synthetic seismograms at the
C             receivers.
C     Horizontal size of each function is 10.23cm, vertical scaling is
C             1cm per the maximum absolute amplitude of the function.
C     MPTS points are plotted for Input signal, Filtered  signal and its
C             Hilbert transform.
C     NPTS/2 points are plotted for the spectra.
C     NPTS points are plotted for the synthetic seismograms.
C
C=======================================================================
C
C     
C
C Common block /RAMC/:
      INCLUDE 'ram.inc'
C     ram.inc
C
C     Allocation of working arrays:
      INTEGER MSS
      PARAMETER (MSS=MRAM/5)
      REAL S1(2,MSS),S2(2,MSS),SS(MSS)
      EQUIVALENCE (S1,RAM         )
      EQUIVALENCE (S2,RAM(2*MSS+1))
      EQUIVALENCE (SS,RAM(4*MSS+1))
C
C-----------------------------------------------------------------------
C
C     Filenames:
      CHARACTER*80 FSSDAT,FILERF,FSSGSE,FSSLOG
C
      PARAMETER (UNDEF=-999999.)
      CHARACTER*80 TEXT1
      REAL PSGNL(10)
C
C     Source coordinates transferred through the GSE file:
      INTEGER NCOM
      PARAMETER (NCOM=3)
      CHARACTER*80 LINE
      CHARACTER*3  TCOM(NCOM)
      REAL         VCOM(NCOM)
      DATA TCOM/'XS1','XS2','XS3'/
C
C.......................................................................
C
C     Opening data files:
      FSSDAT='ss.dat'
      FILERF='rf.out'
      FSSGSE='ss.gse'
      FSSLOG='sslog.out'
      KCOMP=0
      WRITE(*,'(A)')
     *        ' Enter 4 filenames (ss.dat, rf.out, ss.gse, sslog.out): '
      READ(*,*) FSSDAT,FILERF,FSSGSE,FSSLOG,KCOMP
      OPEN(5,FILE=FSSDAT,STATUS='OLD')
      OPEN(6,FILE=FSSLOG)
C     Input data for the source signal
      WRITE(6,'(A)') ' Input data for source signal'
      KEXP=1
      MPTS=0
      NPTS=MSS
      DER=0.
      HILB=0.
      READ (5,*) KSGNL,KEXP,MPTS,NPTS,DER,HILB
      WRITE(6,'(A,5I4,2F8.3)')
     *    ' ***',KSGNL,KEXP,KCOMP,MPTS,NPTS,DER,HILB
      FMIN=0.
      FL=0.
      FR=UNDEF
      FMAX=UNDEF
      TL=0.
      TRED=0.
      VRED=0.
      READ (5,*) FMIN,FL,FR,FMAX,TL,TD,TRED,VRED
      IF(FMAX.EQ.UNDEF) THEN
        FMAX=0.5/TD
      END IF
      IF(FR.EQ.UNDEF) THEN
        FR=FMAX
      END IF
      WRITE(6,'(A,10F8.3)') ' ***',FMIN,FL,FR,FMAX,TL,TD,TRED,VRED
      N   = NPTS-1
      NPTS= 1
      DO 1 KPTS=1,15
        NPTS= NPTS+NPTS
        IF (N-1.LE.0) GO TO 10
        N= N/2
    1 CONTINUE
   10 CONTINUE
      IF (NPTS.GT.MSS) THEN
C       SS-01
        PAUSE 'Error SS-01: Array dimension MSS less than NPTS.'
        STOP
      END IF
      IF (KSGNL.LE.0) THEN
        DO 12 I=1,NPTS
          S1(1,I)=0.
 12     CONTINUE
        READ (5,*) (S1(1,I),I=1,NPTS)
        DO 13 I=NPTS,1,-1
          IF(S1(1,I).NE.0.) GO TO 14
          S1(1,I)=0.
 13     CONTINUE
 14     CONTINUE
        N = I
        N1= (NPTS-N)/2
        N2= NPTS-N-N1
        J = NPTS
        DO 15 I=1,N2
          S1(1,J)= 0.
          J = J-1
   15   CONTINUE
        DO 16 I=1,N
          K = J-N1
          S1(1,J)= S1(1,K)
          J = J-1
   16   CONTINUE
        DO 17 I=1,N1
          S1(1,I)= 0.
   17   CONTINUE
        TL= TL-TD*FLOAT(N1)
      ELSE
        READ (5,*) (PSGNL(I),I=1,10)
        WRITE(6,'(A,10F8.3)') ' ***',(PSGNL(I),I=1,10)
        CALL SIGNAL(KSGNL,NPTS,TL,TD,S1,PSGNL)
      END IF
      DO 20 I=1,NPTS
        S1(2,I)= 0.
   20 CONTINUE
C
      IF(KCOMP.GT.0) THEN
        OPEN(4,FILE=FILERF,STATUS='OLD')
      END IF
      OPEN(7,FILE=FSSGSE)
C
C     Plotting the input signal:
      WRITE(6,'(/A,I2)') ' Source signal No.',KSGNL
      WRITE(6,'(2A/2A)') '   *   Left-hand  Left-hand Right-hand',
     *                   ' Right-hand  Non-zero   Maximum ',
     *                   '          tip     hill-side  hill-side',
     *                   '     tip       range   amplitude'
      CALL PLSIG(MPTS,1,1,MPTS,NPTS,TL,TD,S1,A,I,J)
C
C     Spectrum of the input signal:
      CALL FCOOLR(KPTS,S1,1.)
      FD= 1./TD/FLOAT(NPTS)
C
C     Amplitude spectrum, frequency window:
      A = 2./FLOAT(NPTS)
      DO 38 I=1,NPTS/2
      S2(1,I)= SQRT(S1(1,I)*S1(1,I)+S1(2,I)*S1(2,I))
      F = FD*FLOAT(I-1)
      IF (F-FMIN) 31,31,32
   31 S2(2,I)= 0.
      GO TO 38
   32 IF (F-FL)   33,34,34
   33 S2(2,I)= A*(.5+.5*COS(3.14159*(F-FL)/(FMIN-FL)))**KEXP
      GO TO 38
   34 IF (F-FR)   35,35,36
   35 S2(2,I)= A
      GO TO 38
   36 IF (F-FMAX) 37,31,31
   37 S2(2,I)= A*(.5+.5*COS(3.14159*(F-FR)/(FMAX-FR)))**KEXP
   38 CONTINUE
      IF (DER) 39,41,39
   39 FDA= 6.283185*FD
      F  = 0.
      DO 40 I=2,NPTS/2
      F  = F+FDA
   40 S2(2,I)= S2(2,I)*F**DER
   41 CONTINUE
      DO 42 I=NPTS/2+1,NPTS
      S2(1,I)= 0.
   42 S2(2,I)= 0.
      CALL PLSIG(MPTS,2,2,NPTS/2,NPTS/2,0.,FD,S2,A,I,J)
      CALL PLSIG(MPTS,3,3,NPTS/2,NPTS/2,0.,FD,S2(2,1),A,I,J)
C
C     Filtration of the input signal:
      A = 1.570796*(DER+HILB)
      C =  COS(A)
      S = -SIN(A)
      DO 47 I=1,NPTS
      A      = S1(1,I)*C-S1(2,I)*S
      S1(2,I)= S1(1,I)*S+S1(2,I)*C
   47 S1(1,I)= A
   48 DO 49 I=1,NPTS
      S1(1,I)= S1(1,I)*S2(2,I)
      S1(2,I)= S1(2,I)*S2(2,I)
   49 S2(1,I)= S2(1,I)*S2(2,I)
      CALL PLSIG(MPTS,4,4,NPTS/2,NPTS/2,0.,FD,S2,A,NF1,NF2)
      DO 50 I=1,NPTS
      S2(1,I)= S1(1,I)
   50 S2(2,I)= S1(2,I)
      CALL FCOOLR(KPTS,S2,-1.)
      CALL PLSIG(MPTS,5,5,MPTS,NPTS,TL,TD,S2,A,N1,N2)
      CALL PLSIG(MPTS,6,6,MPTS,NPTS,TL,TD,S2(2,1),A,N,J)
C     Legend:
      CALL SYMBOL(-1.4,-3.,0.8,'1',0.,1)
      CALL SYMBOL( 0.0,-3.,0.3,'INPUT SIGNAL',0.,12)
      CALL SYMBOL( 5.6,-3.,0.3,'RIGHT:',0.,6)
      CALL SYMBOL( 7.6,-3.,0.4,'MAXIMUM AMPLITUDE',0.,17)
      CALL SYMBOL(-1.4,-4.,0.8,'2',0.,1)
      CALL SYMBOL( 0.0,-4.,0.3,'SPECTRUM OF THE INPUT SIGNAL',0.,28)
      CALL SYMBOL(-1.4,-5.,0.8,'3',0.,1)
      CALL SYMBOL( 0.0,-5.,0.3,'SPECTRUM OF THE FILTER',0.,22)
      CALL SYMBOL(-1.4,-6.,0.8,'4',0.,1)
      CALL SYMBOL( 0.0,-6.,0.3,'SPECTRUM OF THE FILTERED SIGNAL',0.,31)
      CALL SYMBOL(-1.4,-7.,0.8,'5',0.,1)
      CALL SYMBOL( 0.0,-7.,0.3,'FILTERED  SIGNAL',0.,16)
      CALL SYMBOL(-1.4,-8.,0.8,'6',0.,1)
      CALL SYMBOL( 0.0,-8.,0.3,
     *                 'HILBERT TRANSFORM OF THE FILTERED SIGNAL',0.,40)
      IF(KCOMP.LE.0) THEN
C       Writing the source time function and its Hilbert transform:
        CALL WGSE1(7,' ')
        DO 51 I=1,NPTS
          SS(I)=S2(1,I)
   51   CONTINUE
        CALL WGSE2(7,' ',' ',0,0.,0.,0.,TL,TD,NPTS,SS)
        DO 52 I=1,NPTS
          SS(I)=S2(2,I)
   52   CONTINUE
        CALL WGSE2(7,' ',' ',0,0.,0.,0.,TL,TD,NPTS,SS)
        CALL WGSE3(7)
        WRITE(*,'(A)') '+SS: Source time function generated.'
        STOP
      END IF
      IF(MPTS.GT.0) CALL PLOT(0.,0.,999)
      IF (N1.GE.NPTS.OR.N.GE.NPTS) THEN
C       SS-02
        PAUSE
     *     'Error SS-02: Too small number NPTS of time samples for FFT.'
        STOP
      END IF
      N1= MIN0(N1,I)
      N2= MIN0(N2,J)
      IF (KCOMP.GT.3) THEN
C       SS-03
        PAUSE 'Error SS-03: KOMP greater than 3.'
        STOP
      END IF
C
C.......................................................................
C
C     Headlines of files:
      WRITE(6,'(/A)')
     *            ' Beginning of the input file with frequency response'
      READ (4,*) TEXT1
      WRITE(6,'(2A)') ' ***',TEXT1
      READ (4,*) (VCOM(I),I=1,3)
      WRITE(6,'(A,10F8.3)') ' ***',(VCOM(I),I=1,3)
      READ (4,*) FMINIM,A,NF
      WRITE(6,'(A,2E12.5,I4)') ' ***',FMINIM,A,NF
      IF (NPTS.NE.INT(1./A/TD+.5)) THEN
C       SS-04
        PAUSE 'Error SS-04: Inconsistent time and frequency steps.'
        STOP
      END IF
      MINIM= INT(FMINIM/FD+1.5)
      MAXIM= MINIM+NF-1
      IF (NF1+1.LT.MINIM) THEN
C       SS-05
        PAUSE
     *     'Error SS-05: Missing low frequencies in response function.'
        STOP
      END IF
      IF (MAXIM+NF2.LT.NPTS/2) THEN
C       SS-06
        PAUSE
     *     'Error SS-06: Missing high frequencies in response function.'
        STOP
      END IF
      CALL WGSE1(7,TEXT1)
      WRITE(6,'(/A)') ' Synthetic sections at receivers'
      WRITE(6,'(2A/2A)') '   *     Coordinates of the receiver ',
     *                   '     First       Last     Upper  ',
     *                   '           X          Y          Z   ',
     *                   '    arrival    arrival  amplitude'
      WRITE(6,'(2A/2A)') '   *   Left-hand  Left-hand Right-han',
     *                   'd Right-hand  Non-zero   Maximum ',
     *                   '          tip     hill-side  hill-sid',
     *                   'e     tip       range   amplitude'
C
C     Preparing source coordinates for output in the GSE file:
      DO 55 I=1,NCOM
        CALL WSEPR(LINE,TCOM(I),VCOM(I))
        CALL WGSE2C(LINE)
   55 CONTINUE
C
C     Loop for the receivers:
      NUMS= 1
      DO 79 IREC=1,999999
        X=UNDEF
        TMIN= 999999.
        TMAX=-999999.
        AMAX=0.
        READ (4,*,END=90) X,Y,Z,TMIN,TMAX,AMAX
        IF(X.EQ.UNDEF) THEN
          GO TO 90
        END IF
        WRITE(6,'(I4,5F11.3,E11.3)') IREC,X,Y,Z,TMIN,TMAX,AMAX
        IF(TMIN.LE.TMAX) THEN
C         Zero range in frequency domain:
          N = MINIM-1
          DO 58 I=1,N
            S2(1,I)= 0.
            S2(2,I)= 0.
   58     CONTINUE
          N = MIN0(NPTS,MAXIM+1)
          DO 59 I=N,NPTS
            S2(1,I)= 0.
            S2(2,I)= 0.
   59     CONTINUE
C
          IF(KCOMP.EQ.1) THEN
            READ (4,'(12F6.0)')
     *                   (S2(1,I),S2(2,I),AUX,AUX,AUX,AUX,I=MINIM,MAXIM)
          ELSE IF(KCOMP.EQ.2) THEN
            READ (4,'(12F6.0)')
     *                   (AUX,AUX,S2(1,I),S2(2,I),AUX,AUX,I=MINIM,MAXIM)
          ELSE IF(KCOMP.EQ.3) THEN
            READ (4,'(12F6.0)')
     *                   (AUX,AUX,AUX,AUX,S2(1,I),S2(2,I),I=MINIM,MAXIM)
          END IF
          A = AMAX/99999.
          DO 65 I=MINIM,MAXIM
            B      = A*(S1(1,I)*S2(1,I)-S1(2,I)*S2(2,I))
            S2(2,I)= A*(S1(1,I)*S2(2,I)+S1(2,I)*S2(1,I))
            S2(1,I)= B
   65     CONTINUE
          CALL FCOOLR(KPTS,S2,-1.)
          N = INT((TMAX+TMIN)/(TD+TD))
          IF(VRED.GT.0) N=NPTS/2+
     *      INT((TRED+SQRT((VCOM(1)-X)**2+(VCOM(2)-Y)**2
     *                                   +(VCOM(3)-Z)**2)/VRED)/TD)
          TMIN= TL+TD*FLOAT(N)
          N = MOD(N,NPTS)
          IF(N.NE.0) THEN
            DO 66 I=1,N
              J = NPTS-N+I
              S2(2,J)= S2(1,I)
   66       CONTINUE
            K = NPTS-N
          END IF
          DO 68 I=1,K
            J = N+I
            S2(2,I)= S2(1,J)
   68     CONTINUE
          CALL PLSIG
     *           (MPTS+1,NUMS,IREC,NPTS,NPTS,TMIN,TD,S2(2,1),AMAX,N1,N2)
          IF(N1.NE.NPTS) THEN
            NUMS= NUMS+1
            TMIN= TMIN+TD*FLOAT(N1)
            N2= NPTS-N2
            N = N2-N1
            DO 69 I=1,N
              N1= N1+1
              SS(I)=S2(2,N1)
   69       CONTINUE
            CALL WGSE2(7,' ',' ',KCOMP,X,Y,Z,TMIN,TD,N,SS)
          ELSE
            CALL WGSE2(7,' ',' ',KCOMP,X,Y,Z,0.,TD,0,SS)
          END IF
        ELSE
          CALL WGSE2(7,' ',' ',KCOMP,X,Y,Z,0.,TD,0,SS)
        END IF
   79 CONTINUE
C
C     End of computation:
   90 IF(MPTS.GT.-1) CALL PLOT(0.,0.,999)
      CALL WGSE3(7)
      WRITE(*,'(A)') '+SS: done.                                       '
      STOP
      END
C
C=======================================================================
C
C     
C
      SUBROUTINE SIGNAL(KSGNL,NPTS,TL,TD,S,PAR)
C
      REAL S(2,NPTS),PAR(*)
C
      GO TO (10,20),KSGNL
C
C     Gabor's signal
   10 T = -TD*FLOAT(NPTS/2)
      TL= TL+T
      A = 6.283185*PAR(1)
      B = A*A/PAR(2)/PAR(2)
      DO 11 I=1,NPTS
      S(1,I)=0.
      IF(B*T*T.LT.70.) S(1,I)= COS(A*T+PAR(3))*EXP(-B*T*T)
      IF(PAR(4).NE.0.) S(1,I)=S(1,I)*PAR(4)
   11 T = T+TD
      RETURN
C
C     Mueller's signal:
   20 N2= IFIX(PAR(2)/PAR(1)/TD/2.)+1
      N1= (NPTS-N2)/2
      TL= TL-TD*FLOAT(N1)
      A = 6.283185*PAR(1)
      B = PAR(2)/(2.+PAR(2))
      C = A/B
      D = SIN(3.141593*PAR(2))/(2.+PAR(2))
      E = A/PAR(2)
      DO 21 I=1,N1
   21 S(1,I)= 0.
      T = 0.
      F = 0.
      N2= N1+N2
      N1= N1+1
      DO 22 I=N1,N2
      S(1,I)= SIN(A*T)-B*SIN(C*T)+D*COS(E*T)-D
      F = AMAX1(F,ABS(S(1,I)))
   22 T = T+TD
      IF(PAR(4).NE.0.) F=F/PAR(4)
      DO 23 I=N1,N2
   23 S(1,I)= S(1,I)/F
      N2= N2+1
      DO 24 I=N2,NPTS
   24 S(1,I)= 0.
      RETURN
      END
C
C=======================================================================
C
C     
C
      SUBROUTINE PLSIG(KPLOT,NUMS,NUM,MPTS,NPTS,TL,TD,S,AMP,N1,N2)
C
      REAL S(2,NPTS)
C
      IF(MPTS.GT.NPTS) THEN
C       SS-07
        WRITE(*,'(2(A,I6))') ' MPTS=',MPTS,' NPTS=',NPTS
        PAUSE 'Error SS-07: MPTS greater than NPTS'
        STOP
      END IF
C
C     Maximum amplitude:
      AMP= 0.000001
      DO 1 I=1,NPTS
    1 AMP= AMAX1(AMP,ABS(S(1,I)))
      EPS= 0.002*AMP
C
C     Zeros beyond and behind the signal:
      DO 2 N1=1,NPTS
      IF(ABS(S(1,N1)).GT.EPS) GO TO 3
    2 CONTINUE
      N1= NPTS
      RETURN
    3 N1= N1-1
      DO 4 N2=1,NPTS
      I = NPTS-N2+1
      IF(ABS(S(1,I)).GT.EPS) GO TO 5
    4 CONTINUE
    5 N2= N2-1
C
C     Writing the parameters of the signal:
      N3= (NPTS-MPTS)/2+1
      N4= N3+MPTS-1
      T1= TL+TD*FLOAT(N1)
      T2= TL+TD*FLOAT(NPTS-N2-1)
      T3= TL+TD*FLOAT(N3-1)
      T4= TL+TD*FLOAT(N4-1)
      A = T2-T1
      WRITE(6,'(I4,5F11.3,E11.3)') NUM,T3,T1,T2,T4,A,AMP
      IF(KPLOT.LE.0) RETURN
C
C     Plotting the signal:
      NUM1= MOD(NUMS-1,14)+1
      IF(NUM1.EQ.1.AND.NUMS.NE.1) CALL PLOT(0.,0.,999)
      IF(NUM1.EQ.1) CALL PLOPN
      CALL PLOT(0.,-2.,-3)
ccc   A = -0.8*FLOAT(NUM/10)-1.428
      A = -0.8*AINT(ALOG10(FLOAT(NUM)+.5)+1.)-1.428
      CALL NUMBER(A,-0.4,0.8,FLOAT(NUM),0.,-1)
      CALL PLTIM(T3,T4,T3,-.3)
      CALL PLTIM(T3,T4,T1,-.5)
      CALL PLTIM(T3,T4,T2,-.5)
      CALL PLTIM(T3,T4,T4,-.3)
      CALL NUMBER(11.016,-0.200,0.4,AMP,0.,6)
      CALL PLOT(10.23,0.00,3)
      CALL PLOT( 0.00,0.00,2)
      A = 10.23/FLOAT(MPTS-1)
      X = 0.
      DO 11 I=N3,N4
      Y = S(1,I)/AMP
      CALL PLOT(X,Y,2)
   11 X = X+A
      RETURN
      END
C
C=======================================================================
C
C     
C
      SUBROUTINE PLOPN
C
      CALL PLOTS(0,0,0)
*     CALL PLOT( 0. , 0. ,3)
*     CALL PLOT(21.0, 0. ,2)
*     CALL PLOT(21.0,29.7,2)
*     CALL PLOT( 0. ,29.7,2)
*     CALL PLOT( 0. , 0. ,2)
      CALL PLOT(5.38,29.7,-3)
      RETURN
      END
C
C=======================================================================
C
C     
C
      SUBROUTINE PLTIM(T3,T4,T,B)
C
      A = (T-T3)/(T4-T3)
      IF(A.LT.-0.01) RETURN
      IF(A.GT. 1.01) RETURN
      A = 10.23*A
      CALL PLOT(A, 0.2,3)
      CALL PLOT(A,-0.2,2)
      A = A-0.457
      CALL NUMBER(A,B-0.1,0.2,T,0.,2)
      RETURN
      END
C
C=======================================================================
C
C       
C
C       Fast Fourier transform of N = 2**K complex data points
C
        SUBROUTINE FCOOLR(K,D,SN)
C
        DIMENSION INU(15),D(*)
C
        LX=2**K
        Q1=LX
        IL=LX
        SH=SN*6.28318530718/Q1
        DO 10 I=1,K
        IL=IL/2
10      INU(I)=IL
        NKK=1
        DO 40 LA=1,K
        NCK=NKK
        NKK=NKK+NKK
        LCK=LX/NCK
        L2K=LCK+LCK
        NW=0
        DO 40 ICK=1,NCK
        FNW=NW
        AA=SH*FNW
        W1=COS(AA)
        W2=SIN(AA)
        LS=L2K*(ICK-1)
        DO 20 I=2,LCK,2
        J1=I+LS
        J=J1-1
        JH=J+LCK
        JH1=JH+1
        Q1=D(JH)*W1-D(JH1)*W2
        Q2=D(JH)*W2+D(JH1)*W1
        D(JH)=D(J)-Q1
        D(JH1)=D(J1)-Q2
        D(J)=D(J)+Q1
20      D(J1)=D(J1)+Q2
        DO 29 I=2,K
        ID=INU(I)
        IL=ID+ID
        IF(NW-ID-IL*(NW/IL)) 40,30,30
30      NW=NW-ID
29      CONTINUE
40      NW=NW+ID
        NW=0
        DO 6 J=1,LX
        IF(NW-J) 8,7,7
7       JJ=NW+NW+1
        J1=JJ+1
        JH1=J+J
        JH=JH1-1
        Q1=D(JJ)
        D(JJ)=D(JH)
        D(JH)=Q1
        Q1=D(J1)
        D(J1)=D(JH1)
        D(JH1)=Q1
8       DO 9 I=1,K
        ID=INU(I)
        IL=ID+ID
        IF(NW-ID-IL*(NW/IL)) 6,5,5
5       NW=NW-ID
9       CONTINUE
6       NW=NW+ID
        RETURN
        END
C
C=======================================================================
C
      INCLUDE 'sep.for'
C     sep.for
      INCLUDE 'gse.for'
C     gse.for
      INCLUDE 'length.for'
C     length.for
      INCLUDE 'calcops.for'
C     calcops.for
C
C=======================================================================
C