C
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 complex-valued amplitude C A11 or A13 or A31 or A33, depending on a P or S wave. C The first index of amplitude Aij corresponds to the C endpoint of the ray, the second index corresponds to C initial point of the ray. Indices i,j=1,2 correspond C to two S-wave polarizations, indices i,j=3 correspond C to a P-wave. C Ray amplitudes A11, A22 or A33 are normalized to 1 at C an initial surface or on a unit sphere around a point C source. C 9. Imaginary part of complex-valued amplitude C A11 or A31 or A13 or A33, depending on a P or S wave. C 10. Real part of complex-valued amplitude C A21 or A32 or A23, or 0, depending on a P or S wave. C 11. Imaginary part of complex-valued amplitude C A21 or A32 or A23, or 0, depending on a P or S wave. C 12. Real part of complex-valued amplitude A12, or 0, C depending on a P or S wave. C 13. Imaginary part of complex-valued amplitude A12, or 0, C depending on a P or S wave. C 14. Real part of complex-valued amplitude A22, or 0, C depending on a P or S wave. C 15. Imaginary part of complex-valued amplitude A22, or 0, C depending on a P or S wave. C OUTMIN,OUTMAX... Minimum and maximum values of corresponding C quantities stored in array OUT. C C SEP parameter: 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 source 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 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 AP00 C AP00... File 'ap.for'. C C Date: 2005, November 12 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) REAL TTINI SAVE TTINI INTEGER KSRC,KTWO,KPHASE SAVE KSRC,KTWO,KPHASE DATA KSRC/-100/ C C....................................................................... C IF(KSRC.EQ.-100) THEN CALL RSEP3I('KSRC',KSRC,0) CALL RSEP3I('KTWO',KTWO,0) CALL RSEP3I('KPHASE',KPHASE,0) END IF 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 Initial travel time IF(IPT.LE.1)THEN TTINI=YI(1) IF(KSRC.GE.1) THEN R(1)=YI(3) R(2)=YI(4) R(3)=YI(5) C The first point in the source file is selected: CALL SRC(1,R(1),R(2),R(3)) C Linear Taylor expansion of travel time: R(1)=R(1)-YI(3) R(2)=R(2)-YI(4) R(3)=R(3)-YI(5) TTINI=TTINI+YI(6)*R(1)+YI(7)*R(2)+YI(8)*R(3) END IF END IF 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 C CRTOUT-01 CALL ERROR * ('CRTOUT-01: Too many ray points to fit in memory') END IF NPTS=NPTS+1 C C Coordinates: R(1)=YI(3) R(2)=YI(4) R(3)=YI(5) IF(KSRC.GE.1) THEN C The first point in the source file is selected: CALL SRC(1,R(1),R(2),R(3)) END IF DO 31 IQ=1,MIN0(NQ,3) OUT(IQ,NPTS)=R(IQ) 31 CONTINUE C Travel time: DO 32 IQ=4,MIN0(NQ,4) OUT(IQ,NPTS)=TTINI IF(KTWO.GE.1) THEN C Converting initial-value to two-point travel time OUT(IQ,NPTS)=OUT(IQ,NPTS)-TTINI END IF 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),7 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 C CRTOUT-02 CALL ERROR * ('CRTOUT-02: Too many ray points to fit in memory') 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 IF(KTWO.GE.1) THEN C Converting initial-value to two-point travel time OUT(IQ,NPTS)=OUT(IQ,NPTS)-TTINI 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) IF(QDET.NE.0.) THEN AUX=SQRT(QDETI*VI*RHOI/(QDET*V*RHO)) ELSE AUX=0. END IF 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 DO 46 IQ=9,NQ,2 C Modifying amplitudes: AUX=SQRT(OUT(IQ-1,NPTS)**2+OUT(IQ,NPTS)**2) IF(KPHASE.EQ.0) THEN C Real and imaginary parts of the complex-valued amplitude IF(ABS(OUT(IQ-1,NPTS)).LT.0.000001*AUX) THEN OUT(IQ-1,NPTS)=0. END IF IF(ABS(OUT(IQ,NPTS)).LT.0.000001*AUX) THEN OUT(IQ,NPTS)=0. END IF ELSE C Absolute value and phase of the complex-valued amplitude IF(AUX.GT.0.) THEN OUT(IQ,NPTS)=ATAN2(OUT(IQ,NPTS),OUT(IQ-1,NPTS)) OUT(IQ-1,NPTS)=AUX END IF IF(ABS(OUT(IQ,NPTS)).LT.0.000001) THEN OUT(IQ,NPTS)=0. END IF END IF 46 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 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,ISRC,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,ISRC,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 output 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 ISRC... Index of the source. C IWAVE...Index of the elementary wave. C (not used for IABS(KALL).LT.2) C IRAY... Index of the ray within the elementary wave. C (not used for IABS(KALL).LT.1) C IREC... Index of the receiver for a two-point ray, determined in C subroutine C RPAR4. 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: 2008, April 7 C Coded by Ludek Klimes C C----------------------------------------------------------------------- C EXTERNAL UARRAY REAL UARRAY C INTEGER MPTS,NSRC,NREC,LENSRC,LENREC PARAMETER (MPTS=1000) CHARACTER*8 POINTS(MPTS) REAL UNDEF,COOR1(MPTS),COOR2(MPTS),COOR3(MPTS) SAVE NSRC,NREC,LENSRC,LENREC,POINTS,UNDEF,COOR1,COOR2,COOR3 C INTEGER I C UNDEF=UARRAY() 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 C CRTOUT-03 CALL ERROR * ('CRTOUT-03: Too many source names to fit in memory') 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 C CRTOUT-04 CALL ERROR * ('CRTOUT-04: Too many receiver names to fit in memory') 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,ISRC,IWAVE,IRAY,IREC,LENTXT,TXT) 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 C CRTOUT-05 CALL ERROR * ('CRTOUT-05: Source index exceeding number of sources') 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 C CRTOUT-06 CALL ERROR * ('CRTOUT-06: Receiver index exceeding number of receivers') 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