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.20
C Date: 1998, November 11
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 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 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) 'SEP',/
C     'SEP'...String with the name of the input SEP data file
C             containing the data describing the parameters of the
C             source time function and frequency-domain filter.
C             Description of file SEP
C     Default: 'SEP'='ss.h'.
C
C                                                     
C Input file 'SEP':
C     It has the form of SEP parameter file.
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 (NPTS-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=integer... 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 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 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: 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)*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                     of 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) 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) '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     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
      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
      PARAMETER (LU4=1,LU5=2,LU6=3,LU7=4)
C
      PARAMETER (UNDEF=-999999.)
      CHARACTER*80 TEXT1
      REAL PSGNL(10)
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*3  TCOM(NCOM)
      REAL         VCOM(NCOM)
      DATA TCOM/'XS1','XS2','XS3'/
C
C.......................................................................
C
C     Opening data files:
      FILDAT='ss.h'
      WRITE(*,'(A)') '+SS: Enter input filename: '
      READ(*,*) FILDAT
      WRITE(*,'(A)') '+SS: Working...            '
C
      CALL RSEP1(LU5,FILDAT)
      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('KEXP  ',KEXP  ,1)
      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('FHIGH',FR     ,FMAX)
      CALL RSEP3R('FMAX ',FMAX   ,0.5/DT)
      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
   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:
      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)
      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 NPTS 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
          WRITE(*,'(A)') '+SS: Nothing to do.                          '
        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'
      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.NE.0) THEN
            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
          END IF
          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)
      REAL S(2,NPTS),PAR(*)
C
      EXTERNAL ERROR
C
      GO TO (10,1,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     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)
      REAL S(2,NPTS)
      CHARACTER*(*) FILPS
C
      EXTERNAL ERROR,PLTIM,PLOTN,PLOT,NUMBER
      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.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+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
      Y = S(1,I)/AMP
      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
      EXTERNAL PLOT,NUMBER
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 'error.for'
C     error.for
      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