C
C Program CRTRAY converting the unformatted output of program CRT into C formatted files with rays suitable for plotting. C C Version: 5.80 C Date: 2003, August 14 C C Coded by: Ludek Klimes C Department of Geophysics, Charles University Prague, C Ke Karlovu 3, 121 16 Praha 2, Czech Republic, C http://sw3d.cz/staff/klimes.htm 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 LINES C designed to store 3-D curves, e.g., for the purposes of plotting. 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-R', and 'CRT-I' are C read by CRTRAY. C For general description of file CRTOUT refer to file C C writ.for. C Default: CRTOUT=' ' which means 'CRT-R'='r01.out', C 'CRT-I'='r01i.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 two-point rays. C Otherwise, the two-point rays are denoted by the receiver C index. C Description of file C 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 rays. C Description of file C SRC. C Default: SRC=' ' C RAYS='string' ... Name of the output formatted file with rays. C This file is designed, for instance, as an input file for C plotting rays. C Description of file C RAYS. C Default: RAYS='rays.out' C Parameters describing the form of output file RAYS: C NQ=integer ... Number of reals in each line of output file RAYS. 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=4 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 KSRC=integer... Modification of the initial point of the ray. C 0: No modification of the initial point of the ray. C 1: If the sourcer file is specified, the coordinates of C the initial point of a ray are replaced by the source C coordinates and the travel time is linearly C interpolated from the initial point to the source C point. C Default: KSRC=0 C KTWO=integer... Converting initial-value to two-point travel time. C 0: No modification of the initial-value travel time. C 1: The initial-value travel time is replaced by the C two-point travel time. C The two-point travel time is the initial-value travel C time minus the travel time at the initial point of the C ray. C Default: KTWO=0 C Optional parameters specifying the form of the real quantities C written in the output formatted files: C MINDIG,MAXDIG=positive integers ... See the description in file C C forms.for. C C C Output formatted file RAYS (format LIN): C (1) / (a slash). C (2) For each ray (2.1), (2.2), and (2.3): C (2.1) 'RAYTXT',/ C 'RAYTXT'... String in apostrophes describing the ray. See the C description of input parameter KALL. C /... An obligatory slash after the string, in place of the C coordinates of the reference point. C (2.2) For each point of the ray (2.1.1): C (2.2.1) (OUT(I),I=1,NQ),/ 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=4: C (2.2.1) X1,X2,X3,TT,/ C X1,X2,X3... Coordinates of the point of the ray. C TT... Arrival time at the point. C (2.3) / (a slash). C (3) / (a slash). C Description of format LIN 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,FORM2 C CRTOUT..File 'crtout.for'. C AP00... File 'ap.for' (called by CRTOUT). C FORM2...File 'forms.for'. C FORM1...File 'forms.for' (called by FORM2). C C----------------------------------------------------------------------- C CHARACTER*80 FILSEP,FCRT,CH INTEGER LU0 PARAMETER (LU0=1) C C Auxiliary storage locations: INTEGER LU1,LU2,LU3,MQ,MPTS,MOUT PARAMETER (LU1=1,LU2=2,LU3=3) PARAMETER (MQ=11,MPTS=2,MOUT=MQ*MPTS) C 1:X1, 2:X2, 3:X3, 4:TT C Optional: 5:P1, 6:P2, 7:P3, 8:AR1, 9:AI1, 10:AR2, 11:AI2 CHARACTER*80 FILE1,FILE2,FILE3,FILSRC,FILREC CHARACTER*(2+8*MQ) FORMAT CHARACTER*80 RAYTXT INTEGER NQ,KALL,NPTS,LENTXT,IQ REAL OUT(MOUT),OUTMIN(MQ),OUTMAX(MQ) C C----------------------------------------------------------------------- C C Reading name of SEP file with input data: WRITE(*,'(A)') '+CRTRAY: Enter input filename: ' FILSEP=' ' READ(*,*) FILSEP WRITE(*,'(A)') '+CRTRAY: Working ... ' C C Reading all data from the SEP file into the memory: IF (FILSEP.NE.' ') THEN CALL RSEP1(LU0,FILSEP) ELSE C CRTRAY-01 CALL ERROR('CRTRAY-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='r01.out' FILE2='r01i.out' IF (FCRT.NE.' ') THEN OPEN (LU1,FILE=FCRT,FORM='FORMATTED',STATUS='OLD') READ (LU1,*) FILE1,CH,FILE2,CH CLOSE (LU1) ENDIF CALL RSEP3T('REC',FILREC,' ') CALL RSEP3T('SRC',FILSRC,' ') CALL RSEP3T('RAYS',FILE3,'rays.out') C C Number of output quantities: CALL RSEP3I('NQ',NQ,4) C Switch between all rays and only two-point rays: CALL RSEP3I('KALL',KALL,0) C NQ=MIN0(NQ,MQ) CALL TXT1(LU3,FILSRC,FILREC) OPEN(LU1,FILE=FILE1,FORM='UNFORMATTED',STATUS='OLD') OPEN(LU2,FILE=FILE2,FORM='UNFORMATTED',STATUS='OLD') OPEN(LU3,FILE=FILE3) C C....................................................................... C C Loop for the points of rays 10 CONTINUE C Reading the results of the complete ray tracing NPTS=0 CALL CRTOUT(LU1,LU2,KALL,0,1,MQ,MPTS,NPTS,OUT,OUTMIN,OUTMAX) IF(IWAVE.LT.1)THEN C End of rays GO TO 80 END IF C Writing the results of the complete ray tracing FORMAT(1:1)='(' CALL FORM2(NQ,OUTMIN,OUTMAX,FORMAT(2:2+8*NQ)) IF(IPT.LE.1)THEN C New ray - recording the initial point CALL TXT2(KALL,0,ISRC,IWAVE,IRAY,IREC,LENTXT,RAYTXT) WRITE(LU3,'(A)') '/' WRITE(LU3,'(2A)') RAYTXT(1:LENTXT),' /' WRITE(LU3,FORMAT) OUT(1),(' ',OUT(IQ),IQ=2,NQ),' /' END IF WRITE(LU3,FORMAT) OUT(1+MQ*(NPTS-1)), * (' ',OUT(IQ),IQ=2+MQ*(NPTS-1),MQ*(NPTS-1)+NQ),' /' GO TO 10 C 80 CONTINUE WRITE(LU3,'(A)') '/' WRITE(LU3,'(A)') '/' CLOSE(LU3) CLOSE(LU2) CLOSE(LU1) WRITE(*,'(A)') '+CRTRAY: 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 C C======================================================================= C