SUBROUTINE CRTOUT * (LU1,LU2,KALL,KREC,INI,NQ,MPTS,NPTS,OUT,OUTMIN,OUTMAX) C INTEGER LU1,LU2,KALL,KREC,INI,NQ,MPTS,NPTS REAL OUT(NQ,MPTS),OUTMIN(NQ),OUTMAX(NQ) C C Subroutine designed to prepare some output quantities of the CRT C program for printing. C C Input: C LU1,LU2... Logical unit numbers corresponding to files with ray C points and ray initial points. C KALL... KALL.LE.0: only two-point rays are considered, C KALL.GE.1: all rays are considered. C KREC... 0: No Taylor expansion of travel time. C 1: Linear Taylor expansion of travel time to the receiver. C 2: Quadratic Taylor expansion of travel time and linear C Taylor expansion of slowness vector to the receiver. C INI... INI.LE.0: initial points of rays are not considered, C INI.GE.1: initial points of rays are considered. C NQ... Number of reals in each output line. C MPTS... Maximum total number of ray points. C NPTS... Number of ray points already stored in array OUT. C OUT... Quantities at ray points already stored in the memory C during previous invocations. C OUTMIN,OUTMAX... Minimum and maximum values of corresponding C quantities stored in array OUT. C C Output: C NPTS... Number of ray points stored in array OUT. C Input value increased by 1 or 2 (if also the initial C point of a ray has been stored). C OUT... For each ray point, the first NQ quantities of the C following ones: 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, coresponding to P- or S1-polarization at C the initial point of the ray. C 9. Imaginary part of ray amplitude coresponding to P- or C S1-polarization at the initial point of the ray. C 10. Real part of ray amplitude coresponding to C S2-polarization at the initial point of the ray. C 11. Imaginary part of ray amplitude coresponding C S2-polarization at the initial point of the ray. C OUTMIN,OUTMAX... Minimum and maximum values of corresponding C quantities stored in array OUT. C C Common block /POINTC/ to store the results of complete ray tracing: INCLUDE 'pointc.inc' C None of the storage locations of the common block are altered. C C Subroutines and external functions required: EXTERNAL AP00 C AP00... File 'ap.for'. C C Date: 1996, April 17 C Coded by Ludek Klimes C C----------------------------------------------------------------------- C C Local storage locations: REAL GIANT PARAMETER (GIANT=1.E20) INTEGER IQ REAL QDETI,QDET,VI,V,RHOI,RHO,AUX REAL HI(18),H(18),HUI(9),RM(6),RN(6),R(3),P(3) C C....................................................................... C IF(NPTS.EQ.0) THEN DO 10 IQ=1,NQ OUTMIN(IQ)= GIANT OUTMAX(IQ)=-GIANT 10 CONTINUE END IF C 20 CONTINUE C Reading the results of the complete ray tracing CALL AP00(0,LU1,LU2) IF(IWAVE.LT.1)THEN C End of rays RETURN ELSE IF (KALL.GT.0.OR.IREC.GT.0) THEN C .............................................................. C Initial point of the ray: IF(INI.GT.0.AND.IPT.LE.1)THEN C New ray - recording the initial point IF(NPTS.GE.MPTS) THEN PAUSE 'ERROR: TOO MANY RAY POINTS TO FIT IN THE MEMORY' STOP END IF NPTS=NPTS+1 C C Coordinates: DO 31 IQ=1,MIN0(NQ,3) OUT(IQ,NPTS)=YI(2+IQ) 31 CONTINUE C Travel time: DO 32 IQ=4,MIN0(NQ,4) OUT(IQ,NPTS)=YI(1) 32 CONTINUE C Slowness vector: DO 33 IQ=5,MIN0(NQ,7) OUT(IQ,NPTS)=YI(1+IQ) 33 CONTINUE C Amplitudes: DO 34 IQ=8,NQ OUT(IQ,NPTS)=0. 34 CONTINUE DO 35 IQ=8,MIN0(NQ,NY-20),6 OUT(IQ,NPTS)=1. 35 CONTINUE C DO 39 IQ=1,NQ OUTMIN(IQ)=AMIN1(OUTMIN(IQ),OUT(IQ,NPTS)) OUTMAX(IQ)=AMAX1(OUTMAX(IQ),OUT(IQ,NPTS)) 39 CONTINUE END IF C .............................................................. C Other than initial ray point: IF(NPTS.GE.MPTS) THEN PAUSE 'ERROR: TOO MANY RAY POINTS TO FIT IN THE MEMORY' STOP END IF NPTS=NPTS+1 C C Coordinates: R(1)=Y(3) R(2)=Y(4) R(3)=Y(5) IF(KREC.GE.1) THEN CALL REC(IREC,R(1),R(2),R(3)) END IF DO 41 IQ=1,MIN0(NQ,3) OUT(IQ,NPTS)=R(IQ) 41 CONTINUE C Travel time: DO 42 IQ=4,MIN0(NQ,4) OUT(IQ,NPTS)=Y(1) IF(KREC.GE.1) THEN C Linear Taylor expansion of travel time: R(1)=R(1)-Y(3) R(2)=R(2)-Y(4) R(3)=R(3)-Y(5) OUT(IQ,NPTS)=OUT(IQ,NPTS)+Y(6)*R(1)+Y(7)*R(2)+Y(8)*R(3) IF(KREC.GE.2) THEN C Quadratic Taylor expansion of travel time: CALL AP03(0,HI,H,HUI) CALL AP08(0,H,HUI,RM,RN) P(1)=RN(1)*R(1)+RN(2)*R(2)+RN(4)*R(3) P(2)=RN(2)*R(1)+RN(3)*R(2)+RN(5)*R(3) P(3)=RN(4)*R(1)+RN(5)*R(2)+RN(6)*R(3) OUT(IQ,NPTS)=OUT(IQ,NPTS) * +0.5*(P(1)*R(1)+P(2)*R(2)+P(3)*R(3)) END IF END IF 42 CONTINUE C Slowness vector: DO 43 IQ=5,MIN0(NQ,7) OUT(IQ,NPTS)=Y(1+IQ) IF(KREC.GE.2) THEN C Linear Taylor expansion of slowness vector: IF(KREC.LT.2) THEN OUT(IQ,NPTS)=Y(1+IQ) ELSE OUT(IQ,NPTS)=Y(1+IQ)+P(IQ-4) END IF END IF 43 CONTINUE C Amplitudes: IF(NQ.GT.7) THEN C wrong! AUX=SQRT( YLI(1)**3*YLI(3)*ABS(YI(14)*YI(19)-YI(15)*YI(18))/ C old AUX=SQRT( YLI(1)**3*YLI(3)/ C old* (YL(1) *YL(3) *ABS( Y(20)*Y(25) - Y(21)*Y(24)))) CALL AP07(QDETI,QDET,VI,V,RHOI,RHO,IQ) AUX=SQRT(QDETI*VI*RHOI/(QDET*V*RHO)) DO 44 IQ=8,MIN0(NQ,NY-20) OUT(IQ,NPTS)=Y(20+IQ)*AUX 44 CONTINUE DO 45 IQ=NY-19,NQ OUT(IQ,NPTS)=0. 45 CONTINUE END IF C DO 49 IQ=1,NQ OUTMIN(IQ)=AMIN1(OUTMIN(IQ),OUT(IQ,NPTS)) OUTMAX(IQ)=AMAX1(OUTMAX(IQ),OUT(IQ,NPTS)) 49 CONTINUE C .............................................................. RETURN END IF GO TO 20 C END C C======================================================================= C SUBROUTINE TXT1(LU,SRCFIL,RECFIL) C C Subroutine designed to read the source and receiver names and prepare C them for entry TXT2. C C ENTRY TXT2(KALL,KTT,IWAVE,IRAY,IREC,LENTXT,TXT) C C Entry designed to generate the string describing the ray or its C endpoint. C C ENTRY REC(IREC,R1,R2,R3) C C ENTRY SRC(IREC,R1,R2,R3) C INTEGER LU,KALL,KTT,IWAVE,IRAY,IREC,LENTXT CHARACTER*(*) SRCFIL,RECFIL,TXT REAL R1,R2,R3 C C Input: C KALL... IABS(KALL)=0: only source (if the source file has been C submitted) and receiver (for two-point rays) information C is contained within the otput string. C IABS(KALL)=1: in addition, the string is prefixed with C the index of the ray. C IABS(KALL)=2: in addition, the string is prefixed also C with the index of the elementary wave. 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 IWAVE...Index of the elementary wave. C IRAY... Index of the ray within the elementary wave. C IREC... Index of the receiver for a two-point ray, C zero for other rays. C C Output: C LENTXT..Length of the string, including also apostrophes. C TXT... The string, beginning and terminating with apostrophes. C If KTT=1, the string is separated by C apostrophe,blank,apostrophe into the source and receiver C parts. C Examples: KTT IREC REC C --------- KALL SRC C ' ' 0 0 0 n C 'REC 13' 0 0 + n n C 'recnam' 0 0 + n y C 'srcnam' 0 0 0 y C 'srcnam, REC 13' 0 0 + y n C 'srcnam TO recnam' 0 0 + y y C 'RAY 112' 0 1 0 n C 'RAY 112, REC 13' 0 1 + n n C 'RAY 112 TO recnam' 0 1 + n y C 'RAY 112 FROM srcnam' 0 1 0 y C 'RAY 112 FROM srcnam, REC 13' 0 1 + y n C 'RAY 112 FROM srcnam to recnam' 0 1 + y y C 'WAVE 1, RAY 112' 0 2 0 n C 'WAVE 1, RAY 112, REC 13' 0 2 + n n C 'WAVE 1, RAY 112 TO recnam' 0 2 + n y C 'WAVE 1, RAY 112 FROM srcnam' 0 2 0 y C 'WAVE 1, RAY 112 FROM srcnam, rec 13' 0 2 + y n C 'WAVE 1, RAY 112 FROM srcnam TO recnam' 0 2 + y y C ' ' ' ' 1 0 0 n C ' ' 'REC 13' 1 0 + n n C ' ' 'recnam' 1 0 + n y C 'srcnam' ' ' 1 0 0 y C 'srcnam' 'REC 13' 1 0 + y n C 'srcnam' 'recnam' 1 0 + y y C 'RAY 112' ' ' 1 1 0 n C 'RAY 112' 'REC 13' 1 1 + n n C 'RAY 112' 'recnam' 1 1 + n y C 'RAY 112 FROM srcnam' ' ' 1 1 0 y C 'RAY 112 FROM srcnam' 'REC 13' 1 1 + y n C 'RAY 112 FROM srcnam' 'recnam' 1 1 + y y C 'WAVE 1, RAY 112' ' ' 1 2 0 n C 'WAVE 1, RAY 112' 'REC 13' 1 2 + n n C 'WAVE 1, RAY 112' 'recnam' 1 2 + n y C 'WAVE 1, RAY 112 FROM srcnam' ' ' 1 2 0 y C 'WAVE 1, RAY 112 FROM srcnam' 'REC 13' 1 2 + y n C 'WAVE 1, RAY 112 FROM srcnam' 'recnam' 1 2 + y y C C Subroutines and external functions required: INTEGER LENGTH EXTERNAL LENGTH C LENGTH..File 'length.for'. C C Date: 1996, August 30 C Coded by Ludek Klimes C C----------------------------------------------------------------------- C INTEGER MPTS,NSRC,NREC,LENSRC,LENREC PARAMETER (MPTS=1000) CHARACTER*8 POINTS(MPTS) REAL UNDEF,COOR1(MPTS),COOR2(MPTS),COOR3(MPTS) PARAMETER (UNDEF=-999999.) SAVE NSRC,NREC,LENSRC,LENREC,POINTS,COOR1,COOR2,COOR3 C INTEGER I,ISRC C C----------------------------------------------------------------------- C NSRC=0 LENSRC=0 IF(SRCFIL.NE.' ') THEN OPEN(LU,FILE=SRCFIL,STATUS='OLD') READ(LU,*) (POINTS(1),I=1,20) 11 CONTINUE IF(NSRC.GE.MPTS) THEN PAUSE 'ERROR: TOO MANY SOURCE NAMES TO FIT IN THE MEMORY' STOP END IF NSRC=NSRC+1 POINTS(NSRC)='$' COOR1(NSRC)=UNDEF COOR2(NSRC)=UNDEF COOR3(NSRC)=UNDEF READ(LU,*,END=12) * POINTS(NSRC),COOR1(NSRC),COOR2(NSRC),COOR3(NSRC) IF(POINTS(NSRC).EQ.'$') THEN GO TO 12 END IF LENSRC=MAX0(LENGTH(POINTS(NSRC)),LENSRC) GO TO 11 12 CONTINUE NSRC=NSRC-1 CLOSE(LU) END IF C NREC=NSRC LENREC=0 IF(RECFIL.NE.' ') THEN OPEN(LU,FILE=RECFIL,STATUS='OLD') READ(LU,*) (POINTS(NREC+1),I=1,20) 21 CONTINUE IF(NREC.GE.MPTS) THEN PAUSE 'ERROR: TOO MANY RECEIVER NAMES TO FIT IN THE MEMORY' STOP END IF NREC=NREC+1 POINTS(NREC)='$' COOR1(NREC)=UNDEF COOR2(NREC)=UNDEF COOR3(NREC)=UNDEF READ(LU,*) POINTS(NREC),COOR1(NREC),COOR2(NREC),COOR3(NREC) IF(POINTS(NREC).EQ.'$') THEN NREC=NREC-1 GO TO 22 END IF LENREC=MAX0(LENGTH(POINTS(NREC)),LENREC) GO TO 21 22 CONTINUE CLOSE(LU) END IF C RETURN C C----------------------------------------------------------------------- C ENTRY REC(IREC,R1,R2,R3) C IF(LENREC.GT.0.AND.IREC.GT.0.AND.IREC.LE.NREC-NSRC) THEN IF(COOR1(NSRC+IREC).NE.UNDEF) R1=COOR1(NSRC+IREC) IF(COOR2(NSRC+IREC).NE.UNDEF) R2=COOR2(NSRC+IREC) IF(COOR3(NSRC+IREC).NE.UNDEF) R3=COOR3(NSRC+IREC) END IF C RETURN C C----------------------------------------------------------------------- C ENTRY SRC(IREC,R1,R2,R3) C IF(LENSRC.GT.0.AND.IREC.GT.0.AND.IREC.LE.NSRC) THEN IF(COOR1(IREC).NE.UNDEF) R1=COOR1(IREC) IF(COOR2(IREC).NE.UNDEF) R2=COOR2(IREC) IF(COOR3(IREC).NE.UNDEF) R3=COOR3(IREC) END IF C RETURN C C----------------------------------------------------------------------- C ENTRY TXT2(KALL,KTT,IWAVE,IRAY,IREC,LENTXT,TXT) C C The name of the first point in the source file is selected: C ------ ISRC=1 C ------ C C Initial apostrophe: LENTXT=1 TXT='''' C C Index of the elementary wave: IF(IABS(KALL).GE.2) THEN TXT(LENTXT+1:LENTXT+10)='WAVE0000, ' WRITE(TXT(LENTXT+5:LENTXT+8),'(I4)') IWAVE LENTXT=LENTXT+10 END IF C C Index of the ray: IF(IABS(KALL).GE.1) THEN TXT(LENTXT+1:LENTXT+8)='RAY00000' WRITE(TXT(LENTXT+4:LENTXT+8),'(I5)') IRAY LENTXT=LENTXT+8 END IF C C Name of the source: IF(LENSRC.GT.0) THEN IF(ISRC.GT.0) THEN IF(ISRC.GT.NSRC) THEN PAUSE 'ERROR: SOURCE INDEX EXCEEDING NUMBER OF SOURCES' STOP END IF C Separator: IF(KALL.NE.0) THEN TXT(LENTXT+1:LENTXT+6)=' FROM ' LENTXT=LENTXT+6 END IF C Name: TXT(LENTXT+1:LENTXT+LENSRC)=POINTS(ISRC)(1:LENSRC) LENTXT=LENTXT+LENSRC END IF END IF C C Separation of source and receiver strings: IF(KTT.NE.0) THEN IF(LENTXT.LE.1) THEN LENTXT=2 END IF TXT(LENTXT+1:LENTXT+3)=''' ''' LENTXT=LENTXT+3 END IF C C Name of the receiver: IF(LENREC.GT.0) THEN C Separator: IF(KTT.EQ.0.AND.LENTXT.GT.1) THEN IF(IREC.GT.0) THEN TXT(LENTXT+1:LENTXT+4)=' TO ' END IF LENTXT=LENTXT+4 END IF C Name: IF(IREC.GT.0) THEN IF(IREC.GT.NREC-NSRC) THEN PAUSE 'ERROR: RECEIVER INDEX EXCEEDING NUMBER OF RECEIVERS' STOP END IF TXT(LENTXT+1:LENTXT+LENREC)=POINTS(NSRC+IREC)(1:LENREC) END IF LENTXT=LENTXT+LENREC END IF C C Index of the receiver: IF(LENREC.EQ.0) THEN C Separator: IF(KTT.EQ.0.AND.LENTXT.GT.1) THEN IF(IREC.GT.0) THEN TXT(LENTXT+1:LENTXT+2)=', ' END IF LENTXT=LENTXT+2 END IF C Index: IF(IREC.GT.0) THEN TXT(LENTXT+1:LENTXT+8)='REC00000' WRITE(TXT(LENTXT+4:LENTXT+8),'(I5)') IREC END IF LENTXT=LENTXT+8 END IF C C Terminating apostrophe: LENTXT=LENTXT+1 TXT(LENTXT:LENTXT)='''' C RETURN END C C======================================================================= C