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: 6.00
C Date: 2006, June 9
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 Description of data files:
C
C Input data read from the standard input device (*):
C     The data are read by the list directed input (free format) and
C     consist of a single string 'SEP':
C     'SEP'...String in apostrophes containing the name of the input
C             SEP parameter or history file with the input data.
C     No default, 'SEP' must be specified and cannot be blank.
C
C                                                     
C Input data file 'SEP':
C     File 'SEP' has the form of the SEP
C     parameter file.  The parameters, which do not differ from their
C     defaults, need not be specified in file 'SEP'.
C Name of the input file:
C     RF='string'... String with the name of the input data file with
C             the frequency-domain response functions at individual
C             receivers.  Is not taken into account if SS=' '.
C             The file is generated by program 'greenss.for'.
C             Description of file RF
C             Default: RF='rf.out' (mostly sufficient)
C Names of output files:
C     SIGGSE='string'... String with the name of the output data file in
C             the GSE data exchange format, containing the filtered
C             source time function and its Hilbert transform.
C             Generated just if SIGGSE.NE.' '.
C             For the GSE data exchange format, refer to the description
C             in file 'gse.for'.
C             Default: SIGGSE=' '
C     SS='string'... String with the name of the output data file in the
C             GSE data exchange format, containing the seismograms at
C             the receivers.
C             Generated just if SS.NE.' '.
C             For the GSE data exchange format, refer to the description
C             in file 'gse.for'.
C             Default: SS='ss.gse'
C     SSLOG='string'... String with the name of the output log file of
C             program SS.
C             Do not specify blank name.
C             Description of file SSLOG
C             Default: SSLOG='sslog.out' (mostly sufficient)
C     SIGPLOT='string'... String with the name of the output PostScript
C             file with the sketches 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             Generated just if SIGPLOT is not blank and MPTS is
C             positive.
C             Description of check plots
C             Default: SIGPLOT=' '
C     SSPLOT='string'... String with the name of the output PostScript
C             file with the sketches of the synthetic seismograms at the
C             receivers.
C             Generated just if both SSPLOT and SS are not blank.
C             Description of check plots
C             Default: SSPLOT=' '
C Time step and time interval for the Fast Fourier Transform:
C     DT=real... Time interval to digitize the source time function and
C             seismograms.
C             Default: DT=1.
C     NFFT=integer... Number of the time samples for the fast Fourier
C             transform.
C             NFFT must be an integer power of 2 (NFFT is rounded up to
C             the nearest power of 2 in this program but not in other
C             programs which may work with it).
C             Default: NFFT=MSS, where MSS is the array dimension
C               declared in the code.
C     TRED=real, SSVRED=real... Specification of the time window for
C             the calculation of synthetic seismograms.  Usually need
C             not be specified.
C             SSVRED.EQ.0: Seismogram is centred in the time interval of
C               length (NFFT-1)*DT according to the travel times of the
C               first and last arrivals at the receiver (default).
C             SSVRED.NE.0: Time of the first sample of the time window
C               is T=TRED+R/SSVRED, where R is the hypocentral distance.
C             Default: TRED=0., SSVRED=0. (mostly sufficient)
C Data describing the filtration of the source time function:
C     DER=real, HILB=real... The source time function is DER-th
C             derivative and HILB-th Hilbert transform of the given
C             signal.
C             Default: DER=0.
C             Default: HILB=0.
C     FMIN=real, FLOW=real, FHIGH=real, FMAX=real... Parameters of the
C             frequency filter to be applied to the source time
C             function.  The filter is zero
C             for frequencies F smaller than FMIN or greater than FMAX.
C             The filter is 1 between FLOW and FHIGH.
C             Between FMIN and FLOW, cosine tapering
C               ( 0.5-0.5*cos(pi*(F-FMIN)/(FLOW-FMIN)) )**KEXP
C             is used.
C             Between FMIN and FLOW, cosine tapering
C               ( 0.5-0.5*cos(pi*(F-FMAX)/(FHIGH-FMAX)) )**KEXP
C             is used.
C             Default: FMIN=0.
C             Default: FLOW=0.
C             Default: FMAX=0.5/DT
C             Default: FHIGH=FMAX
C     KEXP=real... Exponent controlling the cosine frequency-domain
C             filter.  Usually need not be specified because the default
C             is the most common option.
C             Default: KEXP=1 (mostly sufficient)
C Input data to control optional plotting:
C     MPTS=integer... Number of points of the time functions at the
C             output check plot.  The corresponding MPTS-1 time
C             intervals are all together scaled to 10.23cm.  The signal
C             is approximately centred.  MPTS does not apply to the
C             spectra at the output check plot, NFFT/2 points of each
C             spectrum is plotted.
C             MPTS=0: No output check plot of the source time function
C                     is generated.
C             MPTS.LT.0: No output plot is generated, including the
C                     seismograms.
C             Default: MPTS=0
C     SMALL=positive real... Amplitudes in absolute value not exceeding
C             SMALL times the maximum absolute amplitude of the
C             seismogram are assumed zero.
C             Default: SMALL=0.002
C     CALCOPS='string'... String with the PostScript instructions, see
C             file calcops.for.
C Data describing the source time function:
C     KSIG=integer... Type of the source time function:
C             KSIG=0: Digitized time function specified point by point.
C             KSIG=1: Gabor signal.
C             * KSIG=2: Hermite-Gaussian (Ricker) signal.
C             KSIG=3: Kuepper (Mueller) signal.
C             * KSIG=4: Rayleigh signal.
C             * KSIG=5: Generalized Berlage signal.
C             * KSIG=6: ?
C             * Only KSIG=0,1,3 is implemented in this version.
C             Default: KSIG=0
C     SIGDIG='string'... Name of the file containing the digitized time
C             function specified point by point and terminated by a
C             slash.  The file is read by the list-directed input.
C             Required just for KSIG=0
C             Default: SIGDIG=' '.
C     SIGT=real... Reference time of the given signal.
C             Used for all signals.
C             Default: SIGT=0. (mostly sufficient)
C     SIGF=positive real... Reference frequency.
C             Required for all analytical signals (KSIG=1,2,3,4,5,6).
C             Default: SIGF=1.
C     SIGW=positive real... Relative width of the signal.
C             Often the width of the signal envelope expressed in the
C             reference half-periods 1/(2*SIGF).
C             Very roughly defined for non-causal functions.
C             Required for all analytical signals (KSIG=1,2,3,4,5,6).
C             Default: SIGW=4.
C     SIGA=real... Amplitude of the maximum of the signal or its
C             envelope.
C             Used for all analytical signals (KSIG=1,2,3,4,5,6).
C             Default: SIGA=1. (mostly sufficient)
C     SIGPH=real... Phase in radians.
C             Used for KSIG=1,3,4,5,6.
C             Default: SIGPH=0.
C     PAR5=real...
C             Used for KSIG=5,6.
C     PAR6=real...
C             Used for KSIG=5.
C Detailed description of the source time functions:
C     KSIG=0: Time function digitized with step DT, specified point by
C             point and terminated by a slash is read from the file
C             given by input parameter SIGDIG.  The first sample
C             corresponds to time SIGT.
C     KSIG=1: Gabor signal:
C               S(t)=SIGA*exp(-(2*pi*SIGF*(t-SIGT)/SIGW)**2)
C                        *cos(2*pi*SIGF*(t-SIGT)+SIGPH)
C             Gabor signal is not causal but has all derivatives
C               continuous.
C             SIGF... Prevailing frequency.
C             SIGW... Relative width of the signal.  In the interval of
C                     SIGW half-periods 1/(2*SIGF) the signal envelope
C                     exceeds 8.48 per cent of its maximum.
C             SIGA... Amplitude of the envelope.
C             SIGPH...Phase.
C     KSIG=2: Ricker signal (according to Sheriff: Encyclopedic
C             Dictionary of Applied Geophysics)
C              f(t)=(1-2*A*t'*t')*exp(-A*t'*t')
C               where A=pi*pi*SIGF*SIGF and t'=t-SIGT-0.5/SIGF.
C     KSIG=7: Hermite-Gaussian (Ricker) signal:
C             * for future extension, not implemented in this version *
C              S(t)=A*(-1)**n*Hn(SIGF*(t-SIGT))*exp(-(SIGF*(t-SIGT))**2)
C               where Hn(x) is the Hermite polynomial of order n=SIGW,
C                 (-1)**n*Hn(x)*exp(-x**2)=(d/dx)**n[exp(-x**2)]
C               and the scaling factor is
C                 A=SIGA*(n/2)!/n!
C             Then
C               S(t)=A*(1/SIGF*d/dt)**n[exp(-(SIGF*(t-SIGT))**2)]
C             Note that Ricker suggested the Hermite-Gaussian signal
C               with n=1, 2 and 3 for the approximation of particle
C               diplacement, velocity and acceleration, respectively,
C               at the receivers very distant from a point source.
C             Most widely used is the special case for SIGW=2,
C               originally suggested by Ricker for particle velocity,
C               S(t)=SIGA*(2*(SIGF*(t-SIGT))**2-1)
C                    *exp(-(SIGF*(t-SIGT))**2)
C             Hermite-Gaussian signal is not causal but has all
C               derivatives continuous.
C             SIGF... Reference frequency.  The mean frequency is
C                     F0=SIGF*(SIGW/2)!/(((SIGW-1)/2))!/pi and the
C                     spectrum has maximum at FM=SIGF*SQRT(SIGW/2)/pi.
C                     Note that for even SIGW
C                       (SIGW/2)!/(((SIGW-1)/2))!
C                            =2*4*...*SIGW/(1*3*...*(SIGW-1))/SQRT(pi)
C                     and for odd SIGW
C                       (SIGW/2)!/(((SIGW-1)/2))!
C                            =3*5*...*SIGW/(2*4*...*(SIGW-1))*SQRT(pi)/2
C             SIGW... Order n of the Hermite-Gaussian signal.  The best
C                     known option is SIGW=2 (Ricker signal for
C                     far-field particle velocity).
C                     SIGW is rounded to the nearest integer.
C     KSIG=3: Kuepper (Mueller) signal:
C             For 0.LE.(t-SIGT).LE.SIGW/(2*SIGF):
C               S(t)=SIGA*A*(sin(2*pi*SIGF*(t-SIGT))
C                             -sin(2*pi*SIGF*(t-SIGT)*(SIGW+2)/SIGW)
C                                * SIGW / (SIGW+2)
C                             -(1-cos(2*pi*SIGF*(t-SIGT)/SIGW))
C                               * sin(pi*SIGW) / (SIGW+2) )
C               with A=1/SMAX, where SMAX is the maximum absolute value
C               of the signal without multiplication factors SIGA*A.
C             Outside the interval:
C               S(t)=0
C             Note that the standard Kuepper signal is defined just for
C               integer SIGW and consists only of the two sine terms.
C               The last term (with cosine) has been added to make the
C               signal continuous also for non-integer SIGW.
C             Kuepper signal is causal, has continuous the first
C               derivative, and for integer SIGW also the second
C               derivative.
C             SIGF... Reference frequency.
C             SIGW... Relative width of the signal, i.e. the duration of
C                     the signal in half-periods 1/(2*SIGF).
C                     SIGW rounded to the next greater integer is the
C                     number of local extrema of the signal.
C             SIGA... Maximum amplitude of the signal.  Determines A.
C             SIGPH...Not applicable.
C     KSIG=4: Rayleigh signal:
C             * for future extension, not implemented in this version *
C               S(t)=SIGA*(COS(SIGPH)-SIN(SIGPH)*SIGF*(t-SIGT))
C                        /(1+(SIGF*(t-SIGT))**2)
C             Note that the Rayleigh signal is the Dirac distribution
C               D(t-SIGT) shifted by i/SIGF in the complex plane,
C                 S(t)=Re(exp(i*SIGPH)*D(t-SIGT+i/SIGF))
C             Rayleigh signal is not causal but has all derivatives
C               continuous.
C             SIGF... Reference frequency.
C             SIGW... Not applicable.
C             SIGA... Amplitude of the envelope.
C             SIGPH...Phase.
C     KSIG=5: Generalized Berlage signal:
C             * for future extension, not implemented in this version *
C             For t.LE.SIGT:
C               S(t)=0
C             For t.GT.SIGT:
C               S(t)=SIGA*A(t)*(t-SIGT)**SIGW*exp(-PAR5*(t-SIGT))
C                                         *sin(2*pi*SIGF*(t-SIGT)+SIGPH)
C               where A(t) normalizes the maximum of the envelope to 1
C               and also enables to generalize the Berlage signal by
C               means of choosing PAR6 positive,
C                 A(t)=(SIGW/(PAR5*(1+SIGW*PAR6*(t-SIGT))*TMAX**2)
C                                          *exp(PAR5*TMAX)
C               where time t=SIGT+TMAX corresponds to the maximum of the
C               envelope,
C                 TMAX=(SQRT(1+4*SIGW*SIGW*PAR6/PAR5)-1)/(2*SIGW*PAR6)
C             Note that the limiting case for PAR6=0 is the standard
C               Berlage signal, with A(t) constant,
C                 A(t)=PAR5/SIGW*exp(SIGW)
C                 TMAX=SIGW/PAR5
C             Note that the limiting case for SIGW=+infinity is given by
C               S(t)=SIGA
C                 *exp(-1/PAR6/(t-SIGT)+2*SQRT(PAR5/PAR6)-PAR5*(t-SIGT))
C                 *sin(2*pi*SIGF*(t-SIGT)+SIGPH)
C               with maximum of the envelope in time t=SIGT+TMAX,
C                 TMAX=1/SQRT(PAR5*PAR6)
C             The generalized Berlage signal is causal, has continuous
C               derivatives of orders less than INT(SIGW),
C               and for SIGPH=0 or SIGPH=pi even of orders less than
C               INT(SIGW)+1.
C             SIGF... Prevailing frequency.
C             SIGW... Controls the onset of the signal.  Derivatives of
C                     orders less than INT(SIGW), and for SIGPH=0 or
C                     SIGPH=pi even of orders less than INT(SIGW)+1, are
C                     continuous.  For PAR6=0 (standard Berlage),
C                     distance TMAX beetween the begining of the signal
C                     and the maximum of the envelope is SIGW/PAR5.
C                     SIGW=999 or greater is understood as infinity.
C             SIGA... Amplitude of the envelope.
C             SIGPH...Phase.
C             PAR5... Controls the effective duration of the signal.
C             PAR6... Zero for standard Berlage signal.  Positive values
C                     enable to decrease distance TMAX beetween the
C                     begining of the signal and the maximum of the
C                     envelope.  To decrease TMAX from SIGW/PAR5 to
C                     K/PAR5, select PAR6=PAR5*(1/K-1/SIGW)/K.
C     KSIG=6: ?
C             * for future extension, not implemented in this version *
C               For 0.LE.(t-SIGT):
C               S(t)=SIGA*exp(-(1/TRED-2+TRED)*pi*PAR5/SIGW)
C                        *sin(2*pi*SIGF*(t-SIGT)+SIGPH)
C               Where:
C                        TRED=4*SIGF*(t-SIGT)/PAR5
C               Otherwise:
C               S(t)=0
C             SIGF... Prevailing frequency.
C             SIGW... Controls the relative width of the signal.
C             SIGA... Amplitude of the envelope.
C             SIGPH...Phase.
C             PAR5... Distance beetween the begining of the signal and
C                     the maximum of the envelope is SIGW/4 prevailing
C                     periods 1/SIGF.
C                     Reasonable value for SIGPH=0 is PAR5=3.
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) X1SRC,X2SRC,X3SRC,/
C     X1SRC,X2SRC,X3SRC... 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) 'REC',X1,X2,X3,TMIN,TMAX,AMAX,/
C     'REC'...Name of the receiver (6 characters).
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 plots:
C     File SIGPLOT (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     File SSPLOT (if SS.NE.' ') contains plots of the synthetic
C             seismograms at the receivers.
C     Horizontal size of each function is 10.23cm, vertical scaling is
C             1cm per the maximum absolute amplitude of the function.
C     Endpoints of the time (or frequency) interval plotted and the
C             leftmost and rightmost samples of amplitudes in absolute
C             value grater than SMALL times the maximum absolute
C             amplitude of the seismogram are labelled with the
C             corresponding time (or frequency) values.
C     MPTS points are plotted for input signal, filtered  signal and its
C             Hilbert transform.
C     NFFT/2 points are plotted for the spectra.
C     NFFT points are plotted for the synthetic seismograms.
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 or Mueller 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 FORMS), 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 FORMS) may be used to generate seismogram plots
C             in the PostScript files.
C
C=======================================================================
C
C     
C
      REAL UARRAY
      EXTERNAL UARRAY
      EXTERNAL RSEP1,RSEP3T,RSEP3I,RSEP3R,WSEPR,WGSE1,WGSE2,WGSE2C,WGSE3
      EXTERNAL ERROR,SIGNAL,PLSIG,SYMBOL,FCOOLR
C
C Common block /RAMC/:
      INCLUDE 'ram.inc'
C     ram.inc
C
C     Allocation of working arrays:
      INTEGER MSS
      PARAMETER (MSS=MRAM/9)
      REAL S1(2,MSS),S2(2,MSS),S3(2,MSS),S4(2,MSS),SS(MSS)
      EQUIVALENCE (S1,RAM         )
      EQUIVALENCE (S2,RAM(2*MSS+1))
      EQUIVALENCE (S3,RAM(4*MSS+1))
      EQUIVALENCE (S4,RAM(6*MSS+1))
      EQUIVALENCE (SS,RAM(8*MSS+1))
C
C-----------------------------------------------------------------------
C
C     Filenames:
      CHARACTER*80 FILDAT,FILLOG,FILRF,FILGSE,FILSIG,FILPS
      INTEGER LU4,LU5,LU6,LU7
      PARAMETER (LU4=1,LU5=2,LU6=3,LU7=4)
C
      INTEGER KSGNL,MPTS,NPTS,KPTS,N,I,N1,N2,J,K,NF,NF1,NF2,IREC,NUMS,
     *        MINIM,MAXIM
      REAL UNDEF,PSGNL(10),DER,HILB,DT,FMIN,FR,FMAX,FL,SIGT,TRED,VRED,
     *     A,F,FD,FDA,C,S,FMINIM,X,Y,Z,TMIN,TMAX,AMAX,B,T0
      CHARACTER*80 TEXT1
C
C     Receiver name:
      CHARACTER*6 REC
C
C     Source coordinates transferred through the GSE file:
      INTEGER NCOM
      PARAMETER (NCOM=3)
      CHARACTER*80 LINE
      CHARACTER*5  TCOM(NCOM)
      REAL         VCOM(NCOM),KEXP
      DATA TCOM/'X1SRC','X2SRC','X3SRC'/
C
      UNDEF=UARRAY()
C
C-----------------------------------------------------------------------
C
C     Reading name of SEP file with input data:
      WRITE(*,'(A)') '+SS: Enter input filename: '
      FILDAT=' '
      READ(*,*) FILDAT
C
C     Reading all data from the SEP file into the memory:
      IF (FILDAT.NE.' ') THEN
        CALL RSEP1(LU5,FILDAT)
      ELSE
C       SS-12
        CALL ERROR('SS-12: SEP file not given')
C       Input file in the form of the SEP (Stanford Exploration Project)
C       parameter or history file must be specified.
C       There is no default filename.
      ENDIF
C
      WRITE(*,'(A)') '+SS: Working...            '
C
      CALL RSEP3T('SSLOG',FILLOG,'sslog.out')
      OPEN(LU6,FILE=FILLOG)
C
C     Input data for the source signal
      CALL RSEP3I('KSIG  ',KSGNL ,0)
      CALL RSEP3I('MPTS  ',MPTS  ,0)
      CALL RSEP3I('NFFT  ',NPTS  ,MSS)
      CALL RSEP3R('DER   ',DER   ,0.)
      CALL RSEP3R('HILB  ',HILB  ,0.)
      CALL RSEP3R('DT   ' ,DT    ,1.)
      CALL RSEP3R('FMIN ' ,FMIN  ,0.)
      CALL RSEP3R('FLOW ' ,FL    ,0.)
      CALL RSEP3R('FMAX ' ,FMAX  ,0.5/DT)
      CALL RSEP3R('FHIGH' ,FR    ,FMAX)
      CALL RSEP3R('KEXP ' ,KEXP  ,1.0)
      CALL RSEP3R('SIGT ' ,SIGT  ,0.)
      CALL RSEP3R('TRED ' ,TRED  ,0.)
      CALL RSEP3R('SSVRED',VRED  ,0.)
      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
C       SS-05
        CALL ERROR
     *          ('SS-05: Too large number NFFT of time samples for FFT')
C       NFFT should not exceed 2**15.  Check the input data.
C       To use greater NFFT, edit subroutine FCOOLR at the end of this
C       file.
   10 CONTINUE
      IF (NPTS.GT.MSS) THEN
C       SS-01
        CALL ERROR('SS-01: Array dimension MSS less than NPTS.')
      END IF
      IF (KSGNL.LE.0) THEN
        CALL RSEP3T('SIGDIG',FILSIG,' ')
        IF (FILSIG.EQ.' ') THEN
C         SS-02
          CALL ERROR('SS-02: No filename with source time function')
C         For input parameter KSIG=0 (default value), input parameter
C         SIGDIG must contain the name of the file with the digitized
C         source time function.  Refer to the description of the
C         input data.
        END IF
        OPEN(LU5,FILE=FILSIG,STATUS='OLD')
        DO 12 I=1,NPTS
          S1(1,I)=0.
 12     CONTINUE
        READ (LU5,*) (S1(1,I),I=1,NPTS)
        CLOSE(LU5)
        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
        SIGT= SIGT-DT*FLOAT(N1)
      ELSE
        CALL RSEP3R('SIGF' ,PSGNL(1),1.)
        CALL RSEP3R('SIGW' ,PSGNL(2),4.)
        CALL RSEP3R('SIGPH',PSGNL(3),0.)
        CALL RSEP3R('SIGA' ,PSGNL(4),1.)
        CALL RSEP3R('SIG5' ,PSGNL(5),0.)
        CALL RSEP3R('SIG6' ,PSGNL(6),0.)
        CALL SIGNAL(KSGNL,NPTS,SIGT,DT,S1,PSGNL)
      END IF
      DO 20 I=1,NPTS
        S1(2,I)= 0.
   20 CONTINUE
C
C     Plotting the input signal:
      CALL RSEP3T('SIGPLOT',FILPS,' ')
      WRITE(LU6,'(/A,I2)') ' Source signal No.',KSGNL
      WRITE(LU6,'(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,SIGT,DT,S1,A,I,J,FILPS)
C
C     Spectrum of the input signal:
      CALL FCOOLR(KPTS,S1,1.)
      FD= 1./DT/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,FILPS)
      CALL PLSIG(MPTS,3,3,NPTS/2,NPTS/2,0.,FD,S2(2,1),A,I,J,FILPS)
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,FILPS)
      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,SIGT,DT,S2,A,N1,N2,FILPS)
      CALL PLSIG(MPTS,6,6,MPTS,NPTS,SIGT,DT,S2(2,1),A,N,J,FILPS)
C     Legend:
      IF(FILPS.NE.' '.AND.MPTS.GT.0) THEN
        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)
      END IF
      CALL RSEP3T('SIGGSE',FILSIG,' ')
      IF(FILSIG.NE.' ') THEN
C       Writing the source time function and its Hilbert transform:
        OPEN(LU7,FILE=FILSIG)
        CALL WGSE1(LU7,' ')
        DO 51 I=1,NPTS
          SS(I)=S2(1,I)
   51   CONTINUE
        CALL WGSE2(LU7,' ',' ',0,0.,0.,0.,SIGT,DT,NPTS,SS)
        DO 52 I=1,NPTS
          SS(I)=S2(2,I)
   52   CONTINUE
        CALL WGSE2(LU7,' ',' ',0,0.,0.,0.,SIGT,DT,NPTS,SS)
        CALL WGSE3(LU7)
        CLOSE(LU7)
        WRITE(*,'(A)') '+SS: Source time function generated.'
        STOP
      END IF
      IF(FILPS.NE.' '.AND.MPTS.GT.0) CALL PLOT(0.,0.,999)
      IF (N1.GE.NPTS.OR.N.GE.NPTS) THEN
C       SS-04
        CALL ERROR
     *          ('SS-04: Too small number NFFT of time samples for FFT')
      END IF
      N1= MIN0(N1,I)
      N2= MIN0(N2,J)
C
C.......................................................................
C
C     Opening output file with seismograms:
      CALL RSEP3T('SS',FILGSE,'ss.gse')
      IF (FILGSE.EQ.' ') THEN
        IF(FILSIG.EQ.' ') THEN
          IF(FILPS.EQ.' ') THEN
            WRITE(*,'(A)') '+SS: Nothing to do.             '
          ELSE
            WRITE(*,'(A)') '+SS: Signal plotted.            '
          END IF
        ELSE
          WRITE(*,'(A)') '+SS: Signal generated.            '
        END IF
        STOP
      END IF
      OPEN(LU7,FILE=FILGSE)
      CALL RSEP3T('SSPLOT',FILPS,' ')
C
C     Opening input files with the response function:
      CALL RSEP3T('RF',FILRF,'rf.out')
      IF (FILRF.EQ.' ') THEN
C       SS-03
        CALL ERROR('SS-03: No file with response function specified')
      END IF
      OPEN(LU4,FILE=FILRF,STATUS='OLD')
C
C     Headlines of files:
      WRITE(LU6,'(/A)')
     *            ' Beginning of the input file with frequency response'
      TEXT1=' '
      READ (LU4,*) TEXT1
      WRITE(LU6,'(2A)') ' ***',TEXT1
      READ (LU4,*) (VCOM(I),I=1,3)
      WRITE(LU6,'(A,10F8.3)') ' ***',(VCOM(I),I=1,3)
      READ (LU4,*) FMINIM,A,NF
      WRITE(LU6,'(A,2E12.5,I4)') ' ***',FMINIM,A,NF
      IF (NPTS.NE.INT(1./A/DT+.5)) THEN
C       SS-06
        CALL ERROR('SS-06: Inconsistent time and frequency steps.')
      END IF
      MINIM= INT(FMINIM/FD+1.5)
      MAXIM= MINIM+NF-1
      IF (NF1+1.LT.MINIM) THEN
C       SS-07
        CALL ERROR
     *          ('SS-07: Missing low frequencies in response function.')
      END IF
      IF (MAXIM+NF2.LT.NPTS/2) THEN
C       SS-08
        CALL ERROR
     *         ('SS-08: Missing high frequencies in response function.')
      END IF
      CALL WGSE1(LU7,TEXT1)
      WRITE(LU6,'(/A)') ' Synthetic sections at receivers'
      WRITE(LU6,'(2A/2A)') '   *     Coordinates of the receiver ',
     *                   '     First       Last     Upper  ',
     *                   '           X          Y          Z   ',
     *                   '    arrival    arrival  amplitude'
      WRITE(LU6,'(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 (LU4,*,END=90) REC,X,Y,Z,TMIN,TMAX,AMAX
        IF(X.EQ.UNDEF) THEN
          GO TO 90
        END IF
        WRITE(LU6,'(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.
            S3(1,I)= 0.
            S3(2,I)= 0.
            S4(1,I)= 0.
            S4(2,I)= 0.
   58     CONTINUE
          N = MIN0(NPTS,MAXIM+1)
          DO 59 I=N,NPTS
            S2(1,I)= 0.
            S2(2,I)= 0.
            S3(1,I)= 0.
            S3(2,I)= 0.
            S4(1,I)= 0.
            S4(2,I)= 0.
   59     CONTINUE
C
          READ(LU4,'(12F6.0)') (S2(1,I),S2(2,I),S3(1,I),S3(2,I),
     *                          S4(1,I),S4(2,I),I=MINIM,MAXIM)
          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
            B      = A*(S1(1,I)*S3(1,I)-S1(2,I)*S3(2,I))
            S3(2,I)= A*(S1(1,I)*S3(2,I)+S1(2,I)*S3(1,I))
            S3(1,I)= B
            B      = A*(S1(1,I)*S4(1,I)-S1(2,I)*S4(2,I))
            S4(2,I)= A*(S1(1,I)*S4(2,I)+S1(2,I)*S4(1,I))
            S4(1,I)= B
   65     CONTINUE
          CALL FCOOLR(KPTS,S2,-1.)
          CALL FCOOLR(KPTS,S3,-1.)
          CALL FCOOLR(KPTS,S4,-1.)
          N = INT((TMAX+TMIN)/(DT+DT))
          IF(VRED.GT.0) N=NPTS/2+
     *      INT((TRED+SQRT((VCOM(1)-X)**2+(VCOM(2)-Y)**2
     *                                   +(VCOM(3)-Z)**2)/VRED)/DT)
          TMIN= SIGT+DT*FLOAT(N)
          N = MOD(N,NPTS)
          IF(N.LT.0) THEN
            N=N+NPTS
          END IF
          DO 66 I=1,N
            J = NPTS-N+I
            S2(2,J)= S2(1,I)
            S3(2,J)= S3(1,I)
            S4(2,J)= S4(1,I)
   66     CONTINUE
          K = NPTS-N
          DO 68 I=1,K
            J = N+I
            S2(2,I)= S2(1,J)
            S3(2,I)= S3(1,J)
            S4(2,I)= S4(1,J)
   68     CONTINUE
          CALL PLSIG
     *     (MPTS+1,NUMS,IREC,NPTS,NPTS,TMIN,DT,S2(2,1),AMAX,N1,N2,FILPS)
          IF(N1.LT.NPTS) THEN
C           Non-zero signal
            NUMS= NUMS+1
            T0= TMIN+DT*FLOAT(N1)
            N2= NPTS-N2
            N = N2-N1
            DO 71 I=1,N
              N1= N1+1
              SS(I)=S2(2,N1)
   71       CONTINUE
            CALL WGSE2(LU7,REC,' ',1,X,Y,Z,T0,DT,N,SS)
          ELSE
C           Zero signal
            CALL WGSE2(LU7,REC,' ',1,X,Y,Z,0.,DT,0,SS)
          END IF
          CALL PLSIG
     *     (MPTS+1,NUMS,IREC,NPTS,NPTS,TMIN,DT,S3(2,1),AMAX,N1,N2,FILPS)
          IF(N1.LT.NPTS) THEN
C           Non-zero signal
            NUMS= NUMS+1
            T0= TMIN+DT*FLOAT(N1)
            N2= NPTS-N2
            N = N2-N1
            DO 72 I=1,N
              N1= N1+1
              SS(I)=S3(2,N1)
   72       CONTINUE
            CALL WGSE2(LU7,REC,' ',2,X,Y,Z,T0,DT,N,SS)
          ELSE
C           Zero signal
            CALL WGSE2(LU7,REC,' ',2,X,Y,Z,0.,DT,0,SS)
          END IF
          CALL PLSIG
     *     (MPTS+1,NUMS,IREC,NPTS,NPTS,TMIN,DT,S4(2,1),AMAX,N1,N2,FILPS)
          IF(N1.LT.NPTS) THEN
C           Non-zero signal
            NUMS= NUMS+1
            T0= TMIN+DT*FLOAT(N1)
            N2= NPTS-N2
            N = N2-N1
            DO 73 I=1,N
              N1= N1+1
              SS(I)=S4(2,N1)
   73       CONTINUE
            CALL WGSE2(LU7,REC,' ',3,X,Y,Z,T0,DT,N,SS)
          ELSE
C           Zero signal
            CALL WGSE2(LU7,REC,' ',3,X,Y,Z,0.,DT,0,SS)
          END IF
        ELSE
          CALL WGSE2(LU7,REC,' ',1,X,Y,Z,0.,DT,0,SS)
          CALL WGSE2(LU7,REC,' ',2,X,Y,Z,0.,DT,0,SS)
          CALL WGSE2(LU7,REC,' ',3,X,Y,Z,0.,DT,0,SS)
        END IF
   79 CONTINUE
C
C     End of computation:
   90 IF(FILPS.NE.' '.AND.MPTS.GT.-1.AND.NUMS.GT.1) CALL PLOT(0.,0.,999)
      CALL WGSE3(LU7)
      CLOSE(LU7)
      CLOSE(LU6)
      CLOSE(LU4)
      WRITE(*,'(A)') '+SS: Done.                            '
      STOP
      END
C
C=======================================================================
C
C     
C
      SUBROUTINE SIGNAL(KSGNL,NPTS,SIGT,DT,S,PAR)
      INTEGER NPTS,KSGNL
      REAL S(2,NPTS),PAR(*),SIGT,DT
C
      EXTERNAL ERROR
      REAL T,A,B,C,D,E,F,TMAX,TRED
      INTEGER I,N1,N2
C
      GO TO (10,20,30,1,1,1) KSGNL
    1 CONTINUE
C     SS-09
      CALL ERROR('SS-09: Only signals KSIG=1,3 allowed')
C
C     Gabor signal
   10 CONTINUE
      T = -DT*FLOAT(NPTS/2)
      SIGT= SIGT+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)
        T = T+DT
   11 CONTINUE
      RETURN
C
C     Ricker signal (according to Sheriff: Encyclopedic
C     Dictionary of Applied Geophysics)
   20 CONTINUE
      T= -DT*FLOAT(NPTS/2)
      SIGT= SIGT+T+1./PAR(1)/2.
C     SIGT= SIGT+T
      A= 3.141593*3.141593*PAR(1)*PAR(1)
      DO 21 I=1,NPTS
        S(1,I)=0.
        S(1,I)= (1.-2*A*T*T)*EXP(-A*T*T)
        T = T+DT
   21 CONTINUE
      RETURN
C
C     Ricker signal (according to specification in this code)
C  20 CONTINUE
C     T = -DT*FLOAT(NPTS/2)
C     SIGT= SIGT+T
C     A = PAR(1)*PAR(1)
C     DO 21 I=1,NPTS
C       S(1,I)=0.
C       S(1,I)= -1.*PAR4*(2*A*T*T-1)*EXP(-A*T*T)
C       T = T+DT
C  21 CONTINUE
C     RETURN
C
C     Kuepper (Mueller) signal:
   30 CONTINUE
      N2= IFIX(PAR(2)/PAR(1)/DT/2.)+1
      N1= (NPTS-N2)/2
      SIGT= SIGT-DT*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 31 I=1,N1
        S(1,I)= 0.
   31 CONTINUE
      T = 0.
      F = 0.
      N2= N1+N2
      N1= N1+1
      DO 32 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)))
        T = T+DT
   32 CONTINUE
      IF(PAR(4).NE.0.) F=F/PAR(4)
      DO 33 I=N1,N2
        S(1,I)= S(1,I)/F
   33 CONTINUE
      N2= N2+1
      DO 34 I=N2,NPTS
        S(1,I)= 0.
   34 CONTINUE
      RETURN
C
C     Generalized Berlage signal:
   50 CONTINUE
      N2= IFIX(ALOG(1000.)/PAR(5)/DT)+1
      N1= (NPTS-N2)/2
      SIGT= SIGT-DT*FLOAT(N1)
      A = 6.283185*PAR(1)
      B=2.*SQRT(PAR(5)/PAR(6))
      TMAX = 2.*PAR(2)*PAR(2)*PAR(6)/PAR(5)
      IF(TMAX.LT.0.000001) THEN
        TMAX=1.
      ELSE
        TMAX=(SQRT(1.+2.*B)-1.)/B
      ENDIF
      TMAX=TMAX*PAR(2)/PAR(5)
      DO 51 I=1,N1+1
        S(1,I)= 0.
   51 CONTINUE
      T = 0.
      F = 0.
      N1= N1+1
      DO 52 I=N1+1,NPTS
        T=T+DT
        S(1,I)=0.
        IF(PAR(2).LE.998.) THEN
          TRED=1.+PAR(2)*PAR(6)*T*TMAX*TMAX
          TRED=PAR(2)*T/(PAR(5)*TRED)
          C=PAR(5)*T
          IF(C.LT.70.) S(1,I)= SIN(A*T+PAR(3))*EXP(-C)*TRED**PAR(2)
        ELSE
          IF(PAR(6).LE.0.) THEN
C           SS-10
            CALL ERROR('SS-10: Signal parameter PAR6 not positive')
          ENDIF
          B=2.*SQRT(PAR(5)/PAR(6))
          C=1./(PAR(6)*T)-B+PAR(5)*T
          IF(C.LT.70.) S(1,I)= SIN(A*T+PAR(3))*EXP(-C)
        ENDIF
        IF(PAR(4).NE.0.) S(1,I)=S(1,I)*PAR(4)
   52 CONTINUE
      RETURN
C
C     Signal No.6:
   60 CONTINUE
      N2= IFIX(PAR(2)/PAR(1)/DT/2.)+1
      N1= (NPTS-N2)/2
      SIGT= SIGT-DT*FLOAT(N1)
      A = 6.283185*PAR(1)
      B = 4.*PAR(1)/PAR(5)
      DO 61 I=1,N1+1
        S(1,I)= 0.
   61 CONTINUE
      T = 0.
      F = 0.
      N1= N1+1
      DO 62 I=N1+1,NPTS
        T=T+DT
        S(1,I)=0.
        TRED=B*T
        C=1./TRED-2.+TRED
        C=C*3.141593*PAR(5)/PAR(2)
        IF(C.LT.70.) S(1,I)= SIN(A*T+PAR(3))*EXP(-C)
        IF(PAR(4).NE.0.) S(1,I)=S(1,I)*PAR(4)
   62 CONTINUE
      RETURN
C
      END
C
C=======================================================================
C
C     
C
      SUBROUTINE PLSIG(KPLOT,NUMS,NUM,MPTS,NPTS,TL,DT,S,AMP,N1,N2,FILPS)
      INTEGER KPLOT,NUMS,NUM,MPTS,NPTS,N1,N2
      REAL S(2,NPTS),TL,DT,AMP
      CHARACTER*(*) FILPS
C
      EXTERNAL ERROR,PLTIM,PLOTN,PLOT,NUMBER
      INTEGER LU6,I,N3,N4,NUM1
      REAL EPS,SMALL,T1,T2,T3,T4,A,X,Y
      PARAMETER (LU6=3)
C
      IF(MPTS.GT.NPTS) THEN
C       SS-11
        WRITE(*,'(2(A,I6))') ' MPTS=',MPTS,' NPTS=',NPTS
        CALL ERROR('SS-11: MPTS greater than NPTS')
      END IF
C
C     Maximum amplitude:
      AMP= 0.
      DO 1 I=1,NPTS
    1 AMP= AMAX1(AMP,ABS(S(1,I)))
      CALL RSEP3R('SMALL',SMALL,0.002)
      EPS= SMALL*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+DT*FLOAT(N1)
      T2= TL+DT*FLOAT(NPTS-N2-1)
      T3= TL+DT*FLOAT(N3-1)
      T4= TL+DT*FLOAT(N4-1)
      A = T2-T1
      WRITE(LU6,'(I4,5F11.3,E11.3)') NUM,T3,T1,T2,T4,A,AMP
      IF(FILPS.EQ.' '.OR.KPLOT.LE.0.) RETURN
C
C     Plotting the signal:
      NUM1= MOD(NUMS-1,14)+1
      CALL PLOTN(FILPS,(NUMS-1)/14)
      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
      IF(AMP.EQ.0.) THEN
        Y = 0.
      ELSE
        Y = S(1,I)/AMP
      END IF
      CALL PLOT(X,Y,2)
   11 X = X+A
      RETURN
      END
C
C=======================================================================
C
C     
C
      SUBROUTINE PLOPN
C
      EXTERNAL PLOTS,PLOT
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
      REAL T,T3,T4,B
      EXTERNAL PLOT,NUMBER
      REAL A
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
      REAL D,SN,SH,FNW,AA,W1,W2,Q1,Q2
      INTEGER INU,I,IL,K,NKK,NCK,LCK,L2K,LX,LA,LS,NW,ICK,J,J1,JH,JH1,ID,
     *        JJ
      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 'error.for'
C     error.for
      INCLUDE 'sep.for'
C     sep.for
      INCLUDE 'gse.for'
C     gse.for
      INCLUDE 'forms.for'
C     forms.for
      INCLUDE 'length.for'
C     length.for
      INCLUDE 'calcops.for'
C     calcops.for
C
C=======================================================================
C