C
C======================================================================= C PROGRAM MTT C to calculate Multivalued Travel Times on a 3-D grid of points. C C Version: 5.30 C Date: 1999, June 11 C C---------------------------------------------------------------------- C C The post processing program to the CRT program performs interpolation C of travel times to the gridpoints of arbitrary regular rectangular C 3-D grid of points. The interpolation inside ray tubes formed by C vertices of homogeneous triangles created during 3-D two-point C ray tracing is realized using the controlled initial-value ray tracing C algorithm. C If the I-th and (I+1)-st points of the three rays forming the ray tube C are situated inside the same complex block, the program performes C interpolation within ray cells whose bottom is formed by the I-th C points and top by (I+1)-st points of the rays. The best results can C thus be obtained if the travel time is the independent variable for C numerical integration in CRT program (default). C C This version of the program performs bicubic interpolation of C travel times within ray cells (in the coordinate system connected C with each cell), other quantities are interpolated bilinearly. C C The interpolation is not performed in demarcation belts between C different ray histories. The maximum width of the belts is controlled C by input parameter AERR of the CRT C program. The typical width of the belts, measured in take-off angles, C is 0.0001 radians in order of magnitude. The belts may become wide C in areas of large geometrical spreading. The division of rays into C different histories is, by default, very detailed and is controlled C by input parameter PRM0(3) of the CRT C program. Such a behaviour is useful for two-point ray tracing but has C some awkward consequences on the interpolation. For example, the C demarcation belt between the rays incident at different sides of the C computational volume for ray tracing extends across the whole model, C causing an artificial gap within the continuous travel times. C If the ray history does not account for the termination of a ray at C different sides of the computational volume, the gaps corresponding C to the edges of the computational volume are removed, but the corners C along the edges are not filled with the ray cells any more. Then the C computational volume for ray tracing should sufficiently exceed the C extent of target volume covered by the grid for interpolation. C Input parameter PRM0(3) of the CRT C program may thus considerably influence the results of the C interpolation. However, gaps due to the demarcation belts C corresponding to structural interfaces cannot be removed within C the present interpolation algorithm. C C The interpolation is not performed within the cells whose bottom C or top is not formed by different points. This happens mainly C in the case of point source. To obtain the results of interpolation C also around the point source, CRT should be run with parameter C ADVANC. C....................................................................... C C C Main input data read from external interactive device (*): C The data consist of one character string, read by list C directed (free format) input. The string has thus to be C enclosed in apostrophes. C The interactive * external unit may be redirected to C the file containing the data. C (1) SEP / C Name of the formatted data file with input parameters. If blank, C default values of the corresponding data are considered. The file C has the form of SEP file, for the description of the SEP format C refer to the file 'sep.for'. C Only the parameters, which should differ from the default, must be C specified in file SEP. C Default: SEP=' '. C C C Input parameters in the input data file SEP: C Names of input files related to ray tracing (output of CRT): C CRTOUT='string'... File with the names of the output files of C program CRT. If blank, default names of the output files C of program CRT are considered. C Description of file CRTOUT. C Default: CRTOUT=' '. C Data specifying the parameters of the target grid: C O1=real, O2=real, O3=real ... Coordinates of the origin of the C grid. C Default: O1=0., O2=0., O3=0.. C D1=real... Grid interval along the X1 axis. C Default: D1=1.. C D2=real... Grid interval along the X2 axis. C Default: D2=1.. C D3=real... Grid interval along the X3 axis. C Default: D3=1.. C N1=positive integer... Number of gridpoints along the X1 axis. C Default: N1=1. C N2=positive integer... Number of gridpoints along the X2 axis. C Default: N2=1. C N3=positive integer... Number of gridpoints along the X3 axis. C Default: N3=1. C PROJ1=real, PROJ2=real, PROJ3=real ... Three components of the C projection vector which is used in 2-D computations to C project the coordinates of gridpoints to the plane defined C by rays. If not specified (default), a vector in the C direction of the lowest dimension of the interpolation C grid is used as the projection vector. C Default: PROJ1=0., PROJ2=0., PROJ3=0.. C Names of output formatted files with interpolated quantities: C NUM='string' ... Name of the output file with N1*N2*N3 array C of integer values, containing the number of arrivals at C each gridpoint of the target grid. C Description of file NUM. C Default: NUM='mtt-num.out'. C MTT='string' ... Name of the output file with the array of C values, containing for each gridpoint the travel times C of all determined arrivals. C Description of output files M*. C Default: MTT='mtt-tt.out'. C MHIST='string' ... Name of the output file with the array of C the ray histories of all determined arrivals. C Description of output files M*. C Default: MHIST=' ' (the file is not generated). C MP1='string' ... Name of the output file with the array of the C first component of the slowness vector. C MP2='string' ... Name of the output file with the array of the C second component of the slowness vector. C MP3='string' ... Name of the output file with the array of the C third component of the slowness vector. C Description of output files M*. C Defaults: MP1=' ', MP2=' ', MP3=' ' (no files C generated). C MPQ11='string' ... Name of the output file with the array of the C component 1,1 of the 4x4 ray propagator matrix. C MPQ21='string' ... Name of the output file with the array of the C component 2,1 of the 4x4 ray propagator matrix. C MPQ31='string' ... Name of the output file with the array of the C component 3,1 of the 4x4 ray propagator matrix. C MPQ41='string' ... Name of the output file with the array of the C component 4,1 of the 4x4 ray propagator matrix. C MPQ12='string' ... Name of the output file with the array of the C component 1,2 of the 4x4 ray propagator matrix. C MPQ22='string' ... Name of the output file with the array of the C component 2,2 of the 4x4 ray propagator matrix. C MPQ32='string' ... Name of the output file with the array of the C component 3,2 of the 4x4 ray propagator matrix. C MPQ42='string' ... Name of the output file with the array of the C component 4,2 of the 4x4 ray propagator matrix. C MPQ13='string' ... Name of the output file with the array of the C component 1,3 of the 4x4 ray propagator matrix. C MPQ23='string' ... Name of the output file with the array of the C component 2,3 of the 4x4 ray propagator matrix. C MPQ33='string' ... Name of the output file with the array of the C component 3,3 of the 4x4 ray propagator matrix. C MPQ43='string' ... Name of the output file with the array of the C component 4,3 of the 4x4 ray propagator matrix. C MPQ14='string' ... Name of the output file with the array of the C component 1,4 of the 4x4 ray propagator matrix. C MPQ24='string' ... Name of the output file with the array of the C component 2,4 of the 4x4 ray propagator matrix. C MPQ34='string' ... Name of the output file with the array of the C component 3,4 of the 4x4 ray propagator matrix. C MPQ44='string' ... Name of the output file with the array of the C component 4,4 of the 4x4 ray propagator matrix. C Description of output files M*. C Defaults: MPQ*=' ' (no files generated). C C C Input formatted file CRTOUT: C The data consist of character strings read by list directed (free C format) input. The strings have thus to be enclosed in apostrophes. C (1) CRT-R, CRT-S, CRT-I, CRT-T / C CRT-R .. Name of the input unformatted file with the quantities C stored along rays (see C.R.T.5.5.1). C Description of file CRT-R C CRT-S .. Not used here. C CRT-I .. Name of the file with the quantities at the initial C points of rays, (see C.R.T.6.1). C Description of file CRT-I. C CRT-T .. Name of the file with the indices of rays forming C the homogeneous triangles in the ray-parameter domain. C Description of file CRT-T. C Default: CRT-R='r01.out', CRT-I='r01i.out', CRT-T='t01.out'. C C C Output formatted file NUM: C N1*N2*N3 integers: For each gridpoint, the number of travel C times determined at the gridpoint. The innermost loop corresponds C to N1, the outermost loop corresponds to N3. C C C Output formatted files M*: C N numbers, where N is the sum of integers in file NUM. C For each gridpoint, all values of a single interpolated quantity C determined at the gridpoint. The values are stored in descending C order according to the travel time. The loop over gridpoints C is the same as in file NUM. C If M* is blank, the output file is not generated. C C....................................................................... C C C Coded by Petr Bulant C bulant@seis.karlov.mff.cuni.cz C C----------------------------------------------------------------------- C Subroutines and external functions required: EXTERNAL LDWARF LOGICAL LDWARF C C Common block /MTTC/ to store the information about triangles and rays: INCLUDE 'mtt.inc' C mtt.inc. C ........................... C Common block /CIGRD/ to store the parameters of the C interpolation grid. REAL O1,O2,O3,D1,D2,D3 REAL PROJ1,PROJ2,PROJ3 INTEGER N1,N2,N3 COMMON/CIGRD/O1,O2,O3,D1,D2,D3,N1,N2,N3,PROJ1,PROJ2,PROJ3 SAVE /CIGRD/ C O1,O2,O3 ... Coordinates of the origin of the grid. C D1,D2,D3 ... Steps per gridpoint. C N1,N2,N3 ... Numbers of gridpoints. C PROJ1,PROJ2,PROJ3 ... Three components of the projection vector, C which is used in 2-D computations to project the C coordinates of gridpoints to the plane defined by C rays. C....................................................................... C Auxiliary storage locations: INTEGER IRAY1,IRAY2,IRAY3 INTEGER I1,I2,I3,I4,I5,INDX,INDX1,INDX2,IT1,IIT1 INTEGER I INTEGER LU0,LU1,LU2,LU3,LU4,LU5 PARAMETER (LU0=10,LU1=1,LU2=2,LU3=3,LU4=4,LU5=5) REAL AUX CHARACTER*80 FILSEP,FILCRT,FILE1,FILE2,FILE3,FILE5,CH CHARACTER*52 TXTERR C 'FOUT(i)'.Name of the output formatted file with the array of C interpolated values. C 'OUT(i)'..Character string, describing which quantities are to be C written into corresponding output file 'FOUT(i)'. C Following is a list of possible values of 'OUT(i)' and C corresponding quantities written into output files: C 'tt' ... Travel time. C 'p1' ... First component of the slowness vector. C 'p2' ... Second component of the slowness vector. C 'p3' ... Third component of the slowness vector. C 'hi' ... Ray history. C 'pq11'...Component 1,1 of the 4x4 ray propagator matrix. C 'pq21'...Component 2,1 of the 4x4 ray propagator matrix. C 'pq31'...Component 3,1 of the 4x4 ray propagator matrix. C 'pq41'...Component 4,1 of the 4x4 ray propagator matrix. C 'pq12'...Component 1,2 of the 4x4 ray propagator matrix. C 'pq22'...Component 2,2 of the 4x4 ray propagator matrix. C 'pq32'...Component 3,2 of the 4x4 ray propagator matrix. C 'pq42'...Component 4,2 of the 4x4 ray propagator matrix. C 'pq13'...Component 1,3 of the 4x4 ray propagator matrix. C 'pq23'...Component 2,3 of the 4x4 ray propagator matrix. C 'pq33'...Component 3,3 of the 4x4 ray propagator matrix. C 'pq43'...Component 4,3 of the 4x4 ray propagator matrix. C 'pq14'...Component 1,4 of the 4x4 ray propagator matrix. C 'pq24'...Component 2,4 of the 4x4 ray propagator matrix. C 'pq34'...Component 3,4 of the 4x4 ray propagator matrix. C 'pq44'...Component 4,4 of the 4x4 ray propagator matrix. C C----------------------------------------------------------------------- C C Reading in a name of the file with the data for interpolation C grid, and a name of the main input data file: FILSEP=' ' WRITE(*,'(A)') ' Enter the name of the input data file SEP: ' READ(*,*) FILSEP WRITE(*,'(''+'',79('' ''))') WRITE(*,'(A)') '+MTT: Reading input data.' C C Reading all the data from the SEP file into the memory: CALL RSEP1(LU1,FILSEP) C C Reading in filenames of the files with computed rays C and triangles: CALL RSEP3T('CRTOUT',FILCRT,' ') FILE1='r01.out' FILE2='r01i.out' FILE3='t01.out' IF (FILCRT.NE.' ') THEN OPEN (LU2,FILE=FILCRT,FORM='FORMATTED',STATUS='OLD') READ (LU2,*) FILE1,CH,FILE2,FILE3 CLOSE (LU2) ENDIF C C Opening the CRT output files: OPEN(LU1,FILE=FILE1,FORM='UNFORMATTED',STATUS='OLD') OPEN(LU2,FILE=FILE2,FORM='UNFORMATTED',STATUS='OLD') OPEN(LU3,FILE=FILE3,FORM='FORMATTED',STATUS='OLD') C C Reading in filenames of the output files, recording C which quantities are to be written into the files: CALL RSEP3T('NUM',FILE5,'mtt-num.out') CALL RSEP3T('MTT',FOUT(1),'mtt-tt.out') IF (FOUT(1).NE.' ') THEN OUT(1)='tt' OUT(2)=' ' OUT(3)=' ' OUT(4)=' ' NOUT=4 CALL RSEP3T('MP1',FOUT(2),' ') IF (FOUT(2).NE.' ') OUT(2)='p1' CALL RSEP3T('MP2',FOUT(3),' ') IF (FOUT(3).NE.' ') OUT(3)='p2' CALL RSEP3T('MP3',FOUT(4),' ') IF (FOUT(4).NE.' ') OUT(4)='p3' ELSE NOUT=0 CALL RSEP3T('MP1',FOUT(NOUT+1),' ') IF (FOUT(NOUT+1).NE.' ') THEN OUT(NOUT+1)='p1' NOUT=NOUT+1 ENDIF CALL RSEP3T('MP2',FOUT(NOUT+1),' ') IF (FOUT(NOUT+1).NE.' ') THEN OUT(NOUT+1)='p2' NOUT=NOUT+1 ENDIF CALL RSEP3T('MP3',FOUT(NOUT+1),' ') IF (FOUT(NOUT+1).NE.' ') THEN OUT(NOUT+1)='p3' NOUT=NOUT+1 ENDIF ENDIF CALL RSEP3T('MHIST',FOUT(NOUT+1),' ') IF (FOUT(NOUT+1).NE.' ') THEN OUT(NOUT+1)= 'hi' NOUT=NOUT+1 ENDIF CALL RSEP3T('MPQ11',FOUT(NOUT+1),' ') IF (FOUT(NOUT+1).NE.' ') THEN OUT(NOUT+1)= 'pq11' NOUT=NOUT+1 ENDIF CALL RSEP3T('MPQ21',FOUT(NOUT+1),' ') IF (FOUT(NOUT+1).NE.' ') THEN OUT(NOUT+1)= 'pq21' NOUT=NOUT+1 ENDIF CALL RSEP3T('MPQ31',FOUT(NOUT+1),' ') IF (FOUT(NOUT+1).NE.' ') THEN OUT(NOUT+1)= 'pq31' NOUT=NOUT+1 ENDIF CALL RSEP3T('MPQ41',FOUT(NOUT+1),' ') IF (FOUT(NOUT+1).NE.' ') THEN OUT(NOUT+1)= 'pq41' NOUT=NOUT+1 ENDIF CALL RSEP3T('MPQ12',FOUT(NOUT+1),' ') IF (FOUT(NOUT+1).NE.' ') THEN OUT(NOUT+1)= 'pq12' NOUT=NOUT+1 ENDIF CALL RSEP3T('MPQ22',FOUT(NOUT+1),' ') IF (FOUT(NOUT+1).NE.' ') THEN OUT(NOUT+1)= 'pq22' NOUT=NOUT+1 ENDIF CALL RSEP3T('MPQ32',FOUT(NOUT+1),' ') IF (FOUT(NOUT+1).NE.' ') THEN OUT(NOUT+1)= 'pq32' NOUT=NOUT+1 ENDIF CALL RSEP3T('MPQ42',FOUT(NOUT+1),' ') IF (FOUT(NOUT+1).NE.' ') THEN OUT(NOUT+1)= 'pq42' NOUT=NOUT+1 ENDIF CALL RSEP3T('MPQ13',FOUT(NOUT+1),' ') IF (FOUT(NOUT+1).NE.' ') THEN OUT(NOUT+1)= 'pq13' NOUT=NOUT+1 ENDIF CALL RSEP3T('MPQ23',FOUT(NOUT+1),' ') IF (FOUT(NOUT+1).NE.' ') THEN OUT(NOUT+1)= 'pq23' NOUT=NOUT+1 ENDIF CALL RSEP3T('MPQ33',FOUT(NOUT+1),' ') IF (FOUT(NOUT+1).NE.' ') THEN OUT(NOUT+1)= 'pq33' NOUT=NOUT+1 ENDIF CALL RSEP3T('MPQ43',FOUT(NOUT+1),' ') IF (FOUT(NOUT+1).NE.' ') THEN OUT(NOUT+1)= 'pq43' NOUT=NOUT+1 ENDIF CALL RSEP3T('MPQ14',FOUT(NOUT+1),' ') IF (FOUT(NOUT+1).NE.' ') THEN OUT(NOUT+1)= 'pq14' NOUT=NOUT+1 ENDIF CALL RSEP3T('MPQ24',FOUT(NOUT+1),' ') IF (FOUT(NOUT+1).NE.' ') THEN OUT(NOUT+1)= 'pq24' NOUT=NOUT+1 ENDIF CALL RSEP3T('MPQ34',FOUT(NOUT+1),' ') IF (FOUT(NOUT+1).NE.' ') THEN OUT(NOUT+1)= 'pq34' NOUT=NOUT+1 ENDIF CALL RSEP3T('MPQ44',FOUT(NOUT+1),' ') IF (FOUT(NOUT+1).NE.' ') THEN OUT(NOUT+1)= 'pq44' NOUT=NOUT+1 ENDIF C NQ=5+NOUT C C Reading in the file with parameters of target grid: CALL CIGRID C C Computing parameter DWARF: DWARF=1. 7 CONTINUE DWARF=DWARF*0.1 AUX=1.+DWARF IF (LDWARF(AUX)) GOTO 7 DWARF=DWARF*100. C C Reading in file with computed triangles, C sorting the rays in triangles: WRITE(*,'(''+'',79('' ''))') WRITE(*,'(A)') '+MTT: Reading the file with triangles.' NT=0 NRAMP=0 IRAYMI=999999 IRAYMA=0 10 CONTINUE READ(LU3,*,END=20) IRAY1,IRAY2,IRAY3 IF (IRAY1.GT.IRAYMA) IRAYMA=IRAY1 IF (IRAY2.GT.IRAYMA) IRAYMA=IRAY2 IF (IRAY3.GT.IRAYMA) IRAYMA=IRAY3 IF (IRAY1.LT.IRAYMI) IRAYMI=IRAY1 IF (IRAY2.LT.IRAYMI) IRAYMI=IRAY2 IF ((IRAY3.NE.0).AND.(IRAY3.LT.IRAYMI)) IRAYMI=IRAY3 IF (IRAY1.LT.IRAY2) THEN I=IRAY1 IRAY1=IRAY2 IRAY2=I ENDIF IF (IRAY2.LT.IRAY3) THEN I=IRAY2 IRAY2=IRAY3 IRAY3=I ENDIF IF (IRAY1.LT.IRAY2) THEN I=IRAY1 IRAY1=IRAY2 IRAY2=I ENDIF IF (NRAMP+3.GT.MRAMP) CALL CIEROR(1,TXTERR) NRAMP=NRAMP+1 IRAM(NRAMP)=IRAY1 NRAMP=NRAMP+1 IRAM(NRAMP)=IRAY2 NRAMP=NRAMP+1 IRAM(NRAMP)=IRAY3 NT=NT+1 GOTO 10 20 CONTINUE CLOSE(LU3) NR=IRAYMA-IRAYMI+1 C C Check for the 'triangles' in case of 2D computation: IF (IRAM(3).EQ.0) THEN DO 25, I1=6,3*NT,3 IF (IRAM(I1).NE.0) THEN C MTT-32 CALL ERROR('MTT-32: Disorder in triangles.') C The first 'triangle' is formed by two rays with zero C instead of the third ray. This identifies 2-D calculation C and all the 'triangles' should be formed by two rays and C zero instead of the third ray. This is not the case of the C current triangle. Either file CRT-T