C
C Program CRTPTS converting the unformatted output of program CRT into a C formatted file containing coordinates, travel times, slowness vectors, C and amplitudes at the endpoints of (usually two-point) rays. C C Version: 5.40 C Date: 2000, May 15 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 simple conversion program may serve as an example how to read and C process output files of the complete ray tracing program 'CRT'. C Reading the output files is completed by a simple invocation of C subroutine AP00 of file 'ap.for', called by means of subroutine CRTOUT C of file 'crtout.for'. C C The structure of the output file with rays is an extension of the C general file form POINTS C designed to store 3-D points. 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 Names of the input and output files: C CRTOUT='string'... File with the names of the output files of C program CRT. A single set of CRT output files is read C from CRTOUT. Only the files 'CRT-S', and 'CRT-I' are C read by CRTRAY. C For general description of file CRTOUT refer to C file writ.for. C Default: CRTOUT=' ' which means 'CRT-S'='s01.out', C 'CRT-I'='s01i.out' C REC='string' ... If non-blank, the name of the file with the C names of the receiver points. The names are then used C within the strings describing the points of two-point C rays. Otherwise, the two-point rays are denoted by the C receiver index. C Description of file REC C Default: REC=' ' C SRC='string' ... If non-blank, the name of the file with the C name of the source point. The name is then used within C the strings describing the points of the rays. C Description of file SRC C Default: SRC=' ' C PTS='string' .. Name of the output formatted file with ray points. C Description of file PTS C Default: PTS='pts.out' C Parameters describing the form of output file PTS: C NQ=integer ... Number of reals in each output line of file PTS. C If NQ exceeds parameter MQ (see the code below), it is set C to MQ. C The output reals represent: C 1. X1-coordinate. C 2. X2-coordinate. C 3. X3-coordinate. C 4. Travel time. C 5. P1 slowness-vector component. C 6. P2 slowness-vector component. C 7. P3 slowness-vector component. C 8. Real part of ray amplitude, normalized to 1 at an C initial surface or along on a unit sphere around a C point source, corresponding to P- or S1-polarization at C the initial point of the ray. C 9. Imaginary part of ray amplitude corresponding to P- or C S1-polarization at the initial point of the ray. C Printed only if greater than 0.000001 in abs value. C 10. Real part of ray amplitude corresponding to C S2-polarization at the initial point of the ray. C Printed only if greater than 0.000001 in abs value. C 11. Imaginary part of ray amplitude corresponding C S2-polarization at the initial point of the ray. C Printed only if greater than 0.000001 in abs value. C Default: NQ=11 C KALL=integer: C KALL.LE.0: only two-point rays are considered. C KALL.GE.1: all rays are considered. C Absolute value specifies the form of the strings C describing individual points. Here are some examples: C ABS(KALL)=0: 'REC 13' C 'recnam' C 'srcnam TO recnam' C ABS(KALL)=1: 'RAY 112' C 'RAY 112, REC 13' C 'RAY 112 TO recnam' C 'RAY 112 FROM srcnam' C 'RAY 112 FROM srcnam TO recnam' C ABS(KALL)=2: 'WAVE 1, RAY 112' C 'WAVE 1, RAY 112, REC 13' C 'WAVE 1, RAY 112 TO recnam' C 'WAVE 1, RAY 112 FROM srcnam' C 'WAVE 1, RAY 112 FROM srcnam TO recnam' C Values KALL=0 and KALL=1 specify the briefest strings. C Default: KALL=0 C ISRC=integer: C -1: Initial points of rays are written to the output file C instead of ray points situated on the storing surface. C 0: Ray points situated on the storing surface are written C to the output file. C 1: If the receiver file is specified, the coordinates C of ray endpoints are replaced by receiver coordinates C and the travel time is linearly interpolated to the C receivers. C 2: If the receiver file is specified, the coordinates C of ray endpoints are replaced by receiver coordinates C and the travel time is quadratically interpolated to C the receivers. C Default: ISRC=0 C KTT=integer: C 0: Creates a single string. C 1: Separates the string into 2 strings: the source and C receiver parts. This option is intended to generate C files with synthetic travel times. C Default: KTT=0 C C C Output formatted file PTS: C (1) / (a slash). C (2) For each ray endpoint (or initial point) (2.1): C (2.1) 'RAYTXT',(OUT(I),I=1,NQ),/ C 'RAYTXT'... One or two strings in apostrophes describing the ray. C See the description of input parameters KALL and KTT. C One string: output format PTS. C Two strings: output format FTT. C (OUT(I),I=1,NQ)... Output quantities at the ray point, see the C description of input parameter NQ. C /... An obligatory slash after at the end of line, in place C where the slowness vector components could be written. C For default NQ=11 and P-wave at the source one of: C (2.1) 'RAYTXT',X1,X2,X3,TT,P1,P2,P3,AR,AI,/ C (2.1) 'RAYTXT',X1,X2,X3,TT,P1,P2,P3,AR,/ C X1,X2,X3... Coordinates of the point of the ray. C TT... Arrival time at the point. C P1,P2,P3... Slowness vector. C AR... Real part of the complex-valued amplitude, normalized to C 1 on a unit sphere. C AI... Imaginary part of the amplitude if it is greater than C 0.000001. C /... An obligatory slash after at the end of line. C (3) / (a slash). C Description of format PTS C Description of format FTT C C----------------------------------------------------------------------- C C Common block /POINTC/ to store the results of complete ray tracing: INCLUDE 'pointc.inc' C pointc.inc C None of the storage locations of the common block are altered. C C Subroutines and external functions required: EXTERNAL CRTOUT,TXT1,TXT2,TTSORT,FORM2 C CRTOUT,TXT1,TXT2... File 'crtout.for'. C AP00... File 'ap.for' (called by CRTOUT). C LENGTH..File 'length.for' (called by CRTOUT). C TTSORT..File 'ttsort.for'. C INDEXX..File 'indexx.for' (called by TTSORT). C FORM2...File 'forms.for'. C FORM1...File 'forms.for' (called by FORM2). C C----------------------------------------------------------------------- C C Common block /RAMC/: INCLUDE 'ram.inc' C ram.inc C C Allocation of working arrays: INTEGER MPTS,MOUT PARAMETER (MPTS=MRAM/8,MOUT=MRAM-4*MPTS) INTEGER IRECS(MPTS),IWAVES(MPTS),IRAYS(MPTS),INDX(MPTS) REAL OUT(MOUT) EQUIVALENCE (IRECS ,RAM ) EQUIVALENCE (IWAVES,RAM( MPTS+1)) EQUIVALENCE (IRAYS ,RAM(2*MPTS+1)) EQUIVALENCE (INDX ,RAM(3*MPTS+1)) EQUIVALENCE (OUT ,RAM(4*MPTS+1)) REAL RECS(MPTS) EQUIVALENCE (IRECS,RECS) C C----------------------------------------------------------------------- C CHARACTER*80 FILSEP,FCRT,CH INTEGER LU0 PARAMETER (LU0=1) C C Auxiliary storage locations: INTEGER LU1,LU2,LU3,MQ PARAMETER (LU1=1,LU2=2,LU3=3) PARAMETER (MQ=11) C 1:X1, 2:X2, 3:X3, 4:TT, 5:P1, 6:P2, 7:P3, 8:AR1, 9:AI1, C 10:AR2, 11:AI2 CHARACTER*80 FILREC,FILSRC,FILE1,FILE2,FILE3 CHARACTER*(4+8*MQ) FORMAT CHARACTER*80 RAYTXT INTEGER NQ,KALL,ISRC,INI,KTT INTEGER MPTS1,NPTS,LENTXT,IQ,II,I,J REAL OUTMIN(MQ),OUTMAX(MQ) C C----------------------------------------------------------------------- C C Reading name of SEP file with input data: WRITE(*,'(A)') '+CRTPTS: Enter input filename: ' FILSEP=' ' READ(*,*) FILSEP WRITE(*,'(A)') '+CRTPTS: Working ... ' C C Reading all data from the SEP file into the memory: IF (FILSEP.NE.' ') THEN CALL RSEP1(LU0,FILSEP) ELSE C CRTPTS-01 CALL ERROR('CRTPTS-01: 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 C Reading input parameters from the SEP file: CALL RSEP3T('CRTOUT',FCRT,' ') FILE1='s01.out' FILE2='s01i.out' IF (FCRT.NE.' ') THEN OPEN (LU1,FILE=FCRT,FORM='FORMATTED',STATUS='OLD') READ (LU1,*) CH,FILE1,FILE2,CH CLOSE (LU1) ENDIF CALL RSEP3T('REC',FILREC,' ') CALL RSEP3T('SRC',FILSRC,' ') CALL RSEP3T('PTS',FILE3,'pts.out') C C Number of output quantities: CALL RSEP3I('NQ',NQ,11) C Switch between all rays and only two-point rays: CALL RSEP3I('KALL',KALL,0) CALL RSEP3I('ISRC',ISRC,0) CALL RSEP3I('KTT',KTT,0) C NQ=MIN0(NQ,MQ) MPTS1=MIN0(MPTS,MOUT/NQ) INI =MAX0(0,MIN0(-ISRC,1)) CALL TXT1(LU3,FILSRC,FILREC) FORMAT(1:1)='(' OPEN(LU1,FILE=FILE1,FORM='UNFORMATTED',STATUS='OLD') OPEN(LU2,FILE=FILE2,FORM='UNFORMATTED',STATUS='OLD') OPEN(LU3,FILE=FILE3) WRITE(LU3,'(A)') '/' C C....................................................................... C C Loop for the points of rays: NPTS=0 10 CONTINUE C Reading the results of the complete ray tracing CALL CRTOUT * (LU1,LU2,KALL,ISRC,INI,NQ,MPTS1,NPTS,OUT,OUTMIN,OUTMAX) IF(IWAVE.LT.1)THEN C End of rays GO TO 60 END IF NPTS=NPTS-INI IRECS(NPTS)=IREC IWAVES(NPTS)=IWAVE IRAYS(NPTS)=IRAY GO TO 10 60 CONTINUE C C....................................................................... C C Sorting two-point rays: IF(KALL.LE.0) THEN CALL TTSORT(NQ,NPTS,4,OUT,IRECS,RECS,INDX) ELSE DO 71 I=1,NPTS INDX(I)=I 71 CONTINUE END IF C C....................................................................... C C Writing ray points: C C Text describing the point: IF(ISRC.LT.0) THEN LENTXT=9 RAYTXT(1:LENTXT)='''000-000''' WRITE(RAYTXT(2:4),'(I3.3)') -ISRC END IF C C Writing: FORMAT(1:4)='(2A,' CALL FORM2(NQ,OUTMIN,OUTMAX,FORMAT(5:4+8*NQ)) DO 89 I=1,NPTS J=INDX(I) C C Text describing the point: IF(ISRC.LT.0) THEN WRITE(RAYTXT(LENTXT-3:LENTXT-1),'(I3.3)') IRECS(J) ELSE CALL TXT2(KALL,KTT,IWAVES(J),IRAYS(J),IRECS(J),LENTXT,RAYTXT) END IF C J=NQ*(J-1) DO 81 IQ=NQ,1,-1 IF(ABS(OUT(IQ+J)).GE.0.000001) THEN GO TO 82 END IF 81 CONTINUE 82 CONTINUE WRITE(LU3,FORMAT) * RAYTXT(1:LENTXT),(' ',OUT(II),II=1+J,IQ+J),' /' 89 CONTINUE C WRITE(LU3,'(A)') '/' CLOSE(LU3) CLOSE(LU2) CLOSE(LU1) WRITE(*,'(A)') '+CRTPTS: Done. ' STOP END C C======================================================================= C INCLUDE 'error.for' C error.for INCLUDE 'sep.for' C sep.for INCLUDE 'forms.for' C forms.for INCLUDE 'length.for' C length.for INCLUDE 'ap.for' C ap.for INCLUDE 'crtout.for' C crtout.for INCLUDE 'ttsort.for' C ttsort.for INCLUDE 'indexx.for' C indexx.for C C======================================================================= C