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.20 C Date: 1998, November 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 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 C Description of data files: C C Main input data read from external interactive device (*): C The data consist of character strings and integers, read by list C directed (free format) input. The strings have thus to be C enclosed in apostrophes. The interactive * external unit may be C redirected to the file containing the data. C (1) 'REC','SRC','POINTS','INITIALPOINTS','PTS',NQ,KALL,ISRC,KTT C 'REC'...If non-blank, the name of the file with the names of the C receiver points. The names are then used within the C strings describing the points of two-point rays. C Otherwise, the two-point rays are denoted by the receiver C index. C Description of file REC C 'SRC'...If non-blank, the name of the file with the name of the C source point. The name is then used within the strings C describing the rays. C Description of file SRC C 'POINTS'... File with the quantities stored at the points of C intersections of rays with a specified surface (see C C.R.T.5.5.2). C 'INITIALPOINTS'... File with the quantities at the initial points C of rays, corresponding to file POINTS (see C.R.T.6.1). C 'PTS'.. Name of the output formatted file with ray points. C Description of file PTS C NQ... Number of reals in each output line. 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 KALL... 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 ISRC.. -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 KTT... 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: 'REC'=' ', 'SRC'=' ', 'POINTS'='s01.out', C 'INITIALPOINTS'='s01i.out', 'PTS'='pts.out', C NQ=11, KALL=0, ISRC=0, KTT=0. C C Input unformatted file POINTS: C See the description within source code file 'writ.for'. C Description of file POINTS C C Input unformatted file INITIALPOINTS: C See the description within source code file 'writ.for'. C Description of file INITIALPOINTS 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 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 Opening input and output files: FILREC=' ' FILSRC=' ' FILE1='s01.out' FILE2='s01i.out' FILE3='pts.out' C Number of output quantities: NQ =MQ C Switch between all rays and only two-point rays: KALL=0 ISRC=0 KTT =0 WRITE(*,'(2A)') * ' Enter 5 filenames (REC,SRC,S01,S01I,PTS),', * ' and 4 integers (NQ,KALL,ISRC,KTT): ' READ(*,*) FILREC,FILSRC,FILE1,FILE2,FILE3,NQ,KALL,ISRC,KTT 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) 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