PROGRAM   R A Y P L O T
***********************

  PROGRAM RAYPLOT IS DESIGNED FOR PLOTTING OF RAY DIAGRAMS,
TRAVEL TIMES AND AMPLITUDES OF SEISMIC BODY WAVES FROM ANY FILE
WHICH HAS A FORM OF THE FILE LU GENERATED FOR THIS PURPOSE BY
PROGRAM SEIS.
  ALTERNATIVELY, PROGRAM RAYPLOT CAN BE USED FOR THE INSPECTION
OF THE CONTENT OF THE FILE LU, WHICH IS WRITTEN IN BINARY FORM.
THE FILE LU CAN BE GENERATED E.G. BY PROGRAM SEIS.

***************************************************************

A SHORT DESCRIPTION OF THE PROGRAM.
**********************************

  RAY DIAGRAMS, TRAVEL TIMES AND AMPLITUDES ARE PLOTTED IN THREE
DIFFERENT FRAMES. FOR EACH ELEMENTARY WAVE, FOR WHICH COMPUTATIONS
WERE PERFORMED, IT IS POSSIBLE TO PLOT ALL THE PICTURES OR ONLY
SOME OF THEM. SEVERAL AMPLITUDE PICTURES MAY BE PLOTTED FOR ONE
ELEMENTARY WAVE, CORRESPONDING TO DIFFERENT COMPONENTS OR, IF
SLIGHT ABSORPTION IS CONSIDERED, TO DIFFERENT FREQUENCIES. THE
HORIZONTAL COORDINATE IN THE RAY DIAGRAM CORRESPONDS TO THE RANGE
IN THE MODEL, AND IT INCREASES FROM THE LEFT TO THE RIGHT. THE
VERTICAL COORDINATE IN THE RAY DIAGRAM PLOTS CORRESPONDS TO THE
DEPTH. IN TRAVEL TIME AND AMPLITUDE PLOTS, THE HORIZONTAL AXIS
CORRESPONDS TO THE RANGE WHEN THE RECEIVERS ARE SITUATED ALONG THE
EARTH'S SURFACE OR ALONG AN INTERFACE. FOR RECEIVERS SITUATED
ALONG A VERTICAL PROFILE, THE HORIZONTAL AXIS OF TRAVEL TIME AND
AMPLITUDE PLOTS CORRESPONDS TO THE DEPTH IN THE MODEL. THE VERTICAL
COORDINATE CORRESPONDS TO THE TRAVEL TIME OR REDUCED TRAVEL TIME
PLOTS AND DECADIC LOGARITHM OF AN AMPLITUDE FOR AMPLITUDE PLOTS.
ALL THE RAYS STORED IN LU OR ONLY SELECTED RAYS FROM IT CAN BE
PLOTTED. THE SELECTION CAN BE DONE IN TWO WAYS:
   A) GLOBAL SELECTION, ACCORDING TO THE POSITION OF THE
TERMINATION POINT OF THE RAY. THE INDEX "IND" (SEE THE DESCRIPTION
OF PROGRAM SEIS) IS USED FOR THIS SELECTION. IN THIS WAY, IT
IS POSSIBLE TO PLOT THE RAYS WITH THE TERMINATION POINTS ALONG
THE EARTH'S SURFACE, ALONG THE LEFT VERTICAL BORDER, RIGHT VERTICAL
BORDER, OR ALONG SOME INTERFACE. IT IS EVEN POSSIBLE TO SELECT AND
PLOT THE RAYS WHICH DO NOT TERMINATE PROPERLY IN CORRESPONDENCE
WITH THE NUMERICAL CODE OF THE WAVE (IND=8, IND=9, ETC.). GLOBAL
SELECTION IS NOT APPLICABLE TO THE FILE LU GENERATED BY PROGRAM
SEIS. THE FILE GENERATED BY PROGRAM SEIS CONTAINS RAYS WHICH
ARE ALREADY SELECTED. THEY HAVE TERMINATION POINTS ALONG THE
EARTH'S SURFACE (IND=3), AT A SPECIFIED SYSTEM OF RECEIVERS.
   B) MANUAL SELECTION, USING SUCCESSIVE NUMBERS OF RAYS, SEE
INDEX NRAY IN INPUT DATA CARD NO.3.
  THE TRAVEL TIMES OBTAINED EXPERIMENTALLY OR IN SOME OTHER
WAY CAN BE PLOTTED TOGETHER WITH COMPUTED TRAVEL TIMES (FOR
APPLICATIONS IN MATHEMATICAL MODELLING). THE SAME APPLIES
TO AMPLITUDES. AMPLITUDES OF RADIAL, TRANSVERSE AND VERTICAL
COMPONENTS OF DISPLACEMENT VECTOR CAN BE PLOTTED.
  CALCOMP PLOTTING ROUTINES ARE REQUIRED.
  THE PROGRAM GENERATES A POSTSCRIPT FILE WITH THE DESIRED PLOT.

****************************************************************

INPUT DATA
**********

  INPUT DATA CONSIST PARTIALLY OF THE DATA GENERATED IN THE PROGRAM
SEIS (OR ANY OTHER PROGRAM WHICH GENERATES A FILE OF THE SAME
STRUCTURE AS SEIS), STORED IN THE BINARY FORM IN THE FILE LU,
AND PARTIALLY OF THE ADDITIONAL INPUT DATA CONTROLLING PLOTTING.
SOME OF THE ADDITIONAL INPUT DATA ARE READ IN FROM STANDARD INPUT
BY LIST-DIRECTED INPUT (FREE FORMAT). SOME OF THE DATA MAY BE ALSO
SPECIFIED IN THE FORM OF THE SEP PARAMETERS,
SEE ALSO THE SEP PARAMETERS CONTROLLING THE FORM OF THE OUTPUT
POSTSCRIPT FILE.


THE DATA STORED IN LU.
*********************
  FOR DETAILS, SEE THE DESCRIPTION OF THE CONTENTS OF THE FILE LU1
IN PROGRAM SEIS.
  THE DATA ARE STORED IN THE FOLLOWING ORDER:
1) ICONT,ITPR
2) NINT, (NPNT(I),I=1,NINT)
3) FOR I=1,2,...,NINT:
   (A1(J,I),B1(J,I),C1(J,I),D1(J,I),X1(J,I),III(J,I),J=1,NC),
   WHERE NC=NPNT(I)-1.
4) X0,Z0,ROS,VPS,VSS
5) N,IND
6) (X(J),Y(J),J=1,N)
   THE DATA 5,6 ARE REPEATED SUCCESSIVELY FOR ALL RAYS. AFTER
   THE LAST RAY, N=0.
7) NS
8) (INDI(I),X(I),T(I),TAS(I),ANG(I),AX(I),AY(I),AZ(I),PHX(I),PHY(I),
   PHZ(I),I=1,NS)
   AFTER THIS, DATA FOR A NEW RAY DIAGRAM ARE STORED, STARTING
   FROM 1). ICONT=0 IS AN INDICATION OF THE END OF THE FILE LU.

THE ADDITIONAL INPUT DATA
*************************
  THE ADDITIONAL INPUT DATA CAN BE INTRODUCED BY A CARD READER
OR BY A TERMINAL. THE PLACES, WHERE THE DATA FROM LU ARE READ IN,
ARE DENOTED BY **LU/1, **LU/2, ETC.

1)ONE CARD
   NEXT,ISHIFT,IPRINT
      NEXT ...  NEXT=1: CONTINUATION OF COMPUTATIONS
                NEXT=0: TERMINATION OF COMPUTATIONS.
      ISHIFT... SHIFT BETWEEN INDIVIDUAL PICTURES, IN CM.
                IF NOT SPECIFIED, ISHIFT=8.
      IPRINT... CONTROLS THE PRINTING OF SOME DATA AND RESULTS
                ON THE LINE PRINTER.
             IPRINT=0: ONLY ADDITIONAL INPUT DATA ARE REPRODUCED.
             IPRINT=1: IN ADDITION TO IPRINT=0, ALSO THE TABLES
                OF VALUES OF AMPLITUDES SHOWN IN THE PLOTS.
             IPRINT=2: IN ADDITION TO IPRINT=0, ALSO THE TABLES
                OF TRAVEL TIMES AND AMPLITUDES STORED IN THE FILE LU.
             IPRINT=3: IN ADDITION TO IPRINT=2, ALSO THE TABLES
                OF RAYS. SEE MORE DETAILS ABOUT THE OUTPUT IN
                THE END OF THE DESCRIPTION.
**LU/1

2) ONE CARD
   ITEXT
   ARBITRARY ALPHANUMERIC TEXT WHICH IS SHOWN BELOW THE PLOTS.

3) ONE CARD, VARIOUS SWITCHES.
   NTICX,NTICY,NTICT,NTICA,INDO,INDT,ICOMP,NRAY,
                  IBOUND,IRED,IRS,NDX,NDY
      NTICX.. NTICX.NE.0: THE NUMBER OF MARKED INTERVALS ON
                THE HORIZONTAL RANGE AXIS BETWEEN TWO NEIGHBOURHOOD
                TICS WITH CORRESPONDING COORDINATE VALUES.
              NTICX=0: THE PLOTTING TERMINATES, INDEPENDENTLY
                OF THE CONTENT OF LU.
      NTICY.. NTICY.NE.0: THE SAME AS NTICX, BUT FOR THE DEPTH AXIS.
              NTICY=0: NO RAY DIAGRAM WILL BE PLOTTED
      NTICT.. NTICT.NE.0: THE SAME AS NTICX, BUT FOR A TIME AXIS
              NTICT=0: NO PLOT OF TRAVEL TIMES.
      NTICA.. NTICA.NE.0: THE SAME AS NTICX, BUT FOR THE AMPLITUDE AXIS.
              NTICA=0: NO PLOT OF AMPLITUDES.
      INDO...   CONTROLS THE GLOBAL SELECTION OF RAYS IN THE RAY
                DIAGRAM. IN THE RAY DIAGRAM, ONLY THE RAYS WITH
                IND=INDO ARE PLOTTED.
                IF INDO=0, ALL RAYS STORED AT LU ARE PLOTTED, IF
                THEY ARE NOT EXCLUDED DURING THE MANUAL SELECTION
                OF RAYS, SEE NRAY.
      INDT...   CONTROLS THE GLOBAL SELECTION OF TRAVEL TIMES
                IN THE PLOT OF THE TRAVEL-TIMES. ONLY THE TRAVEL
                TIMES FOR IND=INDT ARE PLOTTED. FOR INDT=1,2,3,
                AMPLITUDES ARE ALSO PLOTTED. IF INDT NOT SPECIFIED,
                INDT=3.
      NRAY...   CONTROLS THE MANUAL SELECTION OF RAYS FOR PLOTTING.
             NRAY.NE.0: NRAY IS A NUMBER OF RAYS STORED IN LU
                WHICH WILL NOT APPEAR IN THE RAY DIAGRAM.
                THESE RAYS ARE SPECIFIED IN CARD NO.6.
              NRAY=0: ALL RAYS ARE PLOTTED (WHEN, OF COURSE,
                THEY ARE NOT EXCLUDED IN THE GLOBAL SELECTION).
      IBOUND... IABS(IBOUND)... NUMBER OF POINTS USED IN THE PLOT
                OF THE RAY DIAGRAMS TO APPROXIMATE INTERFACES.
              IBOUND.LT.0... TO DISTINGUISH INTERFACES FROM
                RAYS, EACH INTERFACE IS PLOTTED BY A HEAVY LINE
                (3 SLIGHTLY SHIFTED SINGLE LINES).
              IBOUND.GT.0... EACH INTERFACE IS PLOTTED BY A
                SINGLE LINE.
              IF NOT SPECIFIED, IBOUND=100 AUTOMATICALLY.
                NOTE: THE FICTITIOUS PARTS OF INTERFACES AND THE
                PARTS OF INTERFACES WHICH COINCIDE WITH THOSE
                DIRECTLY ABOVE THEM ARE NOT PLOTTED.
      IRED.. IRED=0: NO REDUCTION OF TRAVEL TIMES
             IRED=1: REDUCTION OF TRAVEL TIMES, THE REDUCTION
                VELOCITY IS GIVEN IN CARD NO.7.
                NOTE: FOR RECEIVERS SITUATED ALONG A VERTICAL
                PROFILE (ITPR=1), AUTOMATICALLY IRED=0.
      IRS...    IRS+1: NUMBER OF THE POINT (THE SAME ON ALL RAYS),
                FROM WHICH PLOT OF THE RAY SHOULD START.
                FOR IRS=0, RAYS ARE PLOTTED RIGHT FROM THE SOURCE.
      NDX,NDY...CONTROL THE PRECISION OF NUMBERS DESCRIBING THE
                COORDINATE AXES IN THE PLOTS.
             ND.GT.0: THE NUMBER OF DIGITS TO THE RIGHT OF THE
                DECIMAL POINT.
             ND.EQ.0: ONLY INTEGER PORTIONS OF THE NUMBERS WITH
                DECIMAL POINTS.
             ND.LT.0: INTEGERS.

4) ONE CARD, DESCRIPTION OF THE RANGE AXIS
   XMIN,XMAX,XLEN,DTICX,SC
      XMIN,XMAX... THE MINIMUM AND MAXIMUM VALUES ON THE RANGE
                AXIS (IN THE USERS LENGTH UNITS).
      XLEN...   LENGTH OF THE RANGE AXIS, IN CM.
      DTICX...  THE DISTANCE BETWEEN TWO NEIGHBOURING MARKS ON
                THE RANGE AXIS WHICH ARE DENOTED BY CORRESPONDING
                COORDINATE VALUES (IN THE USERS LENGTH UNITS).
             DTICX.GT.0.: TIC MARKS STARTING FROM XMIN AND
                APPEARING AT THE SUBSEQUENT POINTS XMIN+DTICX,
                XMIN+2.*DTICX,...
             DTICX.LT.0.: TIC MARKS START AND CONTINUED TO
                BE PLOTTED FROM THE FIRST INTEGER MULTIPLE OF
                ABS(DTICX) GREATER THAN XMIN.
      SC...     CONTROLS THE SCALES OF TICS AND ALPHANUMERICAL
                TEXTS. FOR SC=1.0, THE TICS ARE 0.15 CM LONG AND
                COORDINATES AND TEXTS DESCRIBING THE PLOTS ARE
                0.4 AND 0.45 CM HIGH, RESPECTIVELY. FOR SC LARGER
                (LESS) THAN 1.0, THE SCALES ARE PROPORTIONALLY
                INCREASED (DECREASED). DEFAULT VALUE, SC=1.0.

NOTE: THE ABOVE PARAMETERS SPECIFY THE HORIZONTAL AXIS IN
PLOTS OF RAY DIAGRAMS FOR ANY RECEIVER CONFIGURATION AND IN
PLOTS OF TRAVEL TIMES AND AMPLITUDES FOR RECEIVERS SITUATED
ALONG THE EARTH'S SURFACE OR ALONG AN INTERFACE. FOR RECEIVERS
SITUATED ALONG A VERTICAL PROFILE, PARAMETERS FROM INPUT
DATA CARD NO.5  ARE AUTOMATICALLY USED FOR THE SPECIFICATION
OF HORIZONTAL AXES OF PLOTS OF TRAVEL TIMES AND AMPLITUDES.

5) ONE CARD, DESCRIPTION OF THE DEPTH AXIS
   THIS CARD IS INCLUDED ONLY WHEN NTICY.NE.0.
   YMIN,YMAX,YLEN,DTICY
   THE SAME AS IN 4, BUT FOR THE DEPTH AXIS.

**LU/2
**LU/3
**LU/4

6) ONE CARD, MANUAL SELECTION OF RAYS
   INCLUDED ONLY WHEN NRAY.NE.0.
   (NR(I),I=1,NRAY)
      NR(I)...  SUCCESSIVE NUMBERS OF RAYS STORED IN LU WHICH
                ARE NOT TO BE PLOTTED IN THE RAY DIAGRAM.

**LU/5
**LU/6
**LU/7
**LU/8

7) ONE CARD, DESCRIPTION OF THE TRAVEL-TIME AXIS
   INCLUDED ONLY WHEN NTICT.NE.0.
   TMIN,TMAX,TLEN,DTICT,VRED
   THE SAME MEANING AS IN 4, BUT FOR THE TRAVEL-TIME AXIS.
   VRED IS THE REDUCTION VELOCITY. FOR IRED=1, ALL THE QUANTITIES
   CORRESPOND TO THE REDUCED TRAVEL TIMES.
   NOTE: IN CASE OF ITPR=1, I.E. FOR RECEIVERS SITUATED ALONG
   A VERTICAL PROFILE), THE HORIZONTAL AXIS OF THE PLOT
   CORRESPONDS TO THE DEPTH, NOT TO THE RANGE. THE PARAMETERS
   FROM INPUT DATA CARD NO.5 ARE AUTOMATICALLY USED FOR THE
   DESCRIPTION OF THE HORIZONTAL AXIS OF THE PLOT

8) ONE CARD
   INCLUDED ONLY WHEN NTICY.NE.0.
   NEXP
   NUMBER OF POINTS OF THE ADDITIONAL TRAVEL-TIME CURVE WHICH
   IS TO BE PLOTTED IN ONE FRAME TOGETHER WITH THE COMPUTED
   TRAVEL-TIME CURVE.
   NEXP=0: NO ADDITIONAL TRAVEL-TIME POINTS ARE PLOTTED

9) ONE CARD
   INCLUDED ONLY WHEN NTICY.NE.0 AND NEXP.NE.0.
   (X(I),T(I),I=1,NEXP)
   RECEIVER COORDINATES AND CORRESPONDING TRAVEL TIMES IN AN
   ADDITIONAL TRAVEL-TIME CURVE. THE TIMES SHOULD BE IN
   A NON-REDUCED FORM. IN THE PLOT, THEY WILL BE REDUCED IN THE
   SAME WAY AS THE COMPUTED TIMES.

10) ONE CARD, DESCRIPTION OF THE AMPLITUDE AXIS
    INCLUDED ONLY WHEN NTICA.NE.0.
    AMIN,AMAX,ALEN,DTICA,FREQ,KABS,ICOMP,MSOUR
    AMIN,AMAX,ALEN,DTICA... THE SAME MEANING AS IN 4, BUT FOR
    THE AMPLITUDE AXIS (DECADIC LOGARITHMES OF AMPLITUDES).
    IF(ALEN.LT..00001), NO AMPLITUDE-DISTANCE CURVE WILL BE PLOTTED.
      FREQ..    FREQUENCY FOR WHICH AMPLITUDES ARE EVALUATED
                (IN CASE OF ABSORPTION, KABS.NE.0).
      KABS.. KABS=0... NO ABSORPTION IS CONSIDERED.
             KABS=1... NON-CAUSAL OR CAUSAL ABSORPTION
                (FUTTERMAN'S MODEL) IS CONSIDERED.
      ICOMP..ICOMP=0: AMPLITUDES OF VERTICAL COMPONENT.
             ICOMP=1: AMPLITUDES OF RADIAL COMPONENT.
             ICOMP=2: AMPLITUDES OF TRANSVERSE COMPONENT.
      MSOUR...  CONTROLS THE CHOICE OF THE TYPE OF THE POINT
                SOURCE.
             MSOUR=0: UNIT ISOTROPIC RADIATION PATTERN.
             MSOUR=1: SINGLE FORCE.
             MSOUR=2: DOUBLE COUPLE.
             MSOUR=3: EXPLOSIVE (IMPLOSIVE) SOURCE.

NOTE1: ABSORPTION EFFECTS CAN ONLY BE CONSIDERED IF GLOBAL
ABSORPTION FACTOR 'TAS', SEE DATA STORED IN LU/8, IS NON-ZERO.
OTHERWISE KABS HAS ONLY A FORMAL MEANING.

NOTE2: ABSORPTION AFFECTS ONLY AMPLITUDES STORED IN FILE LU,
NOT THE AMPLITUDES FORMING THE ADDITONAL AMPLITUDE CURVES
READ IN IN THE INPUT DATA CARD NO.13.

NOTE3: THE NOTE FROM THE DESCRIPTION OF THE INPUT DATA CARD
NO.7 HOLDS ALSO FOR AMPLITUDES.

11) ONE CARD, SPECIFICATION OF PARAMETERS OF THE SOURCE
    INCLUDED ONLY WHEN MSOUR.NE.0.
   IPAR(1),...,IPAR(4),PAR(1),...PAR(6)
     FOR THE MEANING OF THE ABOVE PARAMETERS SEE DESCRIPTION OF
     THE INPUT DATA CARD NO.7 IN PROGRAM SYNTPL IN THIS PACKAGE.

12) ONE CARD
    INCLUDED ONLY WHEN NTICA.NE.0.
   NEXP
    NUMBER OF POINTS IN THE ADDITIONAL AMPLITUDE CURVE WHICH
    IS TO BE PLOTTED IN THE SAME FRAME AS THE COMPUTED AMPLITUDE
    CURVE. FOR NEXP=0, NO ADDITIONAL AMPLITUDE CURVE IS PLOTTED.

13) ONE CARD
    INCLUDED ONLY WHEN NTICA.NE.0 AND NEXP.NE.0.
   (X(I),Y(I),I=1,NEXP)
     RECEIVER COORDINATES AND CORRESPONDING AMPLITUDES OF THE
     ADDITIONAL AMPLITUDE CURVE. REAL AMPLITUDES, NOT LOGARITHMES
     OF AMPLITUDES, SHOULD BE SPECIFIED.

****************************************************************

TERMINATION OF COMPUTATIONS
***************************

THE INPUT DATA CARDS NO.10-13 CAN BE REPEATED FOR THE GIVEN
ELEMENTARY WAVE ANY NUMBER TIMES (TO PERFORM PLOTTING OF AMPLITUDE-
DISTANCE CURVES OF THE ELEMENTARY WAVE E.G. FOR DIFFERENT FREQUENCIES,
IF ABSORPTION IS CONSIDERED, OR FOR DIFFERENT COMPONENTS, OR FOR
DIFFERENT SPECIFICATIONS OF THE SOURCE). THE LAST CARD SHOULD
BE A BLANK CARD (TO GUARANTEE ALEN.LT..00001). AFTER THE BLANK CARD,
A NEW SYSTEM OF CARDS 2-13 FOR A NEW ELEMENTARY WAVE STORED IN THE
FILE LU MAY FOLLOW.

AS SOON AS ICONT=0 IS FOUND IN LU, THE PLOTTING FROM THE FILE LU
TERMINATES. THE PLOTTING FROM THE FILE LU MAY BE ALSO TERMINATED
IF WE PUT NTICX=0 IN THE CARD NO.3 AND NEXT=0 in CARD NO.1, EVEN
IF LU STILL CONTAINS SOME OTHER DATA. IN BOTH CASES, THE FILE LU
IS CLOSED AND THE COMPUTATIONS TERMINATE.

******************************************************************

THE OUTPUT ON A LINE PRINTER
****************************

  THE OUTPUT IS CONTROLLED BY THE PARAMETER IPRINT, SEE INPUT DATA
CARD NO 1.
  FOR IPRINT=0, ALL THE ADDITIONAL DATA SHOWN ABOVE AND THE DATA
LU/1, LU/4, LU/5 AND LU/7 ARE PRINTED.
  FOR IPRINT=1, THE DATA FOR IPRINT=0 AND ALSO THE RECEIVER COORDINATES
AND VALUES OF CORRESPONDING AMPLITUDES WHICH APPEAR IN THE PLOTS ARE
PRINTED.
  FOR IPRINT=2, THE DATA FOR IPRINT=0 AND ALSO THE DATA LU/8 ARE
PRINTED.
  FOR IPRINT=3, THE DATA FOR IPRINT=2 AND ALSO ALL THE DATA FROM
LU ARE PRINTED, EVEN THE COORDINATES OF ALL RAYS.

 ****************************************************************