C
C Subroutine file 'mttq.for' to deal with the user-defined quantities C to be interpolated by program MTT. C The user is welcome to edit this file according to the instructions C provided. He will be able to read his input parameters, and to C define quantities in the points of the rays. The quantities, which C are assumed to be real-valued, are then bilinearly interpolated by C program MTT. C C Version: 5.40 C Date: 2000, May 31 C C Coded by Petr Bulant and Ludek Klimes C Department of Geophysics, Charles University Prague, C Ke Karlovu 3, 121 16 Praha 2, Czech Republic, C E-mails: bulant@seis.karlov.mff.cuni.cz, C klimes@seis.karlov.mff.cuni.cz C C....................................................................... C C A list of the quantities which may already be interpolated by MTT can C be found in the file mtt.for. C C If the user wishes to interpolate other quantities, he may add them by C editting this file. Any user-defined quantity may be interpolated. C The quantities stored in common block POINTC C are at user's disposal for the calculation of the user-defined C quantity at a point of ray. The user may define the value of his C quantity at the initial point of a ray in subprogram entry MTTQRI. C Similarly, the user may define the value of his quantity at the other C points of a ray in subprogram entry MTTQRA. C To let program 'mtt.for' know about the new quantity for C interpolation, the quantity should be registered in subroutine MTTQ. C The user should supplement subroutine MTTQ with reading the filename C of the output file with the interpolated quantity. At the same place, C the user may also wish to supplement reading some other input C parameters to be used in entries MTTQRI and MTTQRA for the calculation C of the new quantity. The supplemented input parameters must be C declared in the SAVE statement in the beginning of subroutine MTTQ. C The filename (and possible other input parameters) should be read by C individual invocations of the subroutines RSEP3T, RSEP3I and RSEP3R of C file sep.for. C C Places to update when defining a new quantity for interpolation: C SAVE statement for possible input data C Reading possible input data C Registering a new quantity C Calculating the value of a new quantity at the C initial point of a ray C Calculating the value of a new quantity at the C other points of a ray C You may also wish to add your description of a new quantity to the C list of the names of output C files with interpolated quantities in program MTT. C Please, send your updated files to the author C to include your improvements in the next version of 'mttq.for'. C C======================================================================= C SUBROUTINE MTTQ C C Subroutine designed to read the input data describing the names C of the output files with the user-defined interpolated quantities, C and also to read any other input parameters, which will be used by the C user. All the parameters to be read are assumed to be specified in C the SEP format in the C MTT input file. C C....................................................................... C C Dummy arguments of entries: CHARACTER*80 QNAME REAL QVALUE C C Subroutines and external functions required: EXTERNAL CIEROR,ERROR,RSEP3T,RSEP3R,RSEP3I C CIEROR ... File mtt.for. C ERROR ... File error.for. C RSEP3T,RSEP3R,RSEP3I ... C File sep.for. C C Common block /POINTC/ to store the results of complete ray tracing: INCLUDE 'pointc.inc' C pointc.inc. C C Common block /MTTC/ to store the information about triangles and rays. INCLUDE 'mtt.inc' C mtt.inc C C----------------------------------------------------------------------- C C Input data: REAL GBE11,GBE21,GBE31,GBE12,GBE22,GBE32 REAL GBR11,GBR12,GBR22,GBY11,GBY22 REAL AMPMAX SAVE GBE11,GBE21,GBE31,GBE12,GBE22,GBE32 SAVE GBR11,GBR12,GBR22,GBY11,GBY22 SAVE AMPMAX C Input data read in this subroutine in order to be used in entries C MTTQRI and MTTQRA to calculate the quantities for interpolation C must be declared in the SAVE statement. C C Values calculated in entry MTTQRI and used in entry MTTQRA: REAL R11,R12,R22,C11,C21,C12,C22,CDET SAVE R11,R12,R22,C11,C21,C12,C22,CDET C The values must be declared in the SAVE statement. C C Temporary storage locations: CHARACTER*80 TXTERR INTEGER I REAL C31,C32,H11,H21,H31,H12,H22,H32,H13,H23,H33,UU,U,U1,U2,U3 REAL Q111,Q211,Q121,Q221,Q112,Q212,Q122,Q222 REAL GREEN(32) C C....................................................................... C Reading the input data: C CALL RSEP3R('GBE11',GBE11,1.) CALL RSEP3R('GBE21',GBE21,0.) CALL RSEP3R('GBE31',GBE31,0.) CALL RSEP3R('GBE12',GBE12,0.) CALL RSEP3R('GBE22',GBE22,1.) CALL RSEP3R('GBE32',GBE32,0.) CALL RSEP3R('GBR11',GBR11,0.) CALL RSEP3R('GBR12',GBR12,0.) CALL RSEP3R('GBR22',GBR22,0.) CALL RSEP3R('GBY11',GBY11,0.) CALL RSEP3R('GBY22',GBY22,0.) CALL RSEP3R('AMPMAX',AMPMAX,999999.) C C C User-defined quantity: C Integer value IVALUE of integer input SEP parameter NAMEI is read C by statement * CALL RSEP3I('NAMEI',IVALUE,0) C where 0 is the default value of the parameter. C C Real value VALUE of the integer or real-valued input SEP parameter C NAMER is read by statement * CALL RSEP3R('NAMER',VALUE,0.) C where 0. is the default value of the parameter. C C....................................................................... C C Registering quantities for bilinear interpolation within ray cells C and reading the filenames of the respective output files: C C Components of the 4*4 paraxial ray propagator matrix in RCCS: CALL RSEP3T('MPQ11',FOUT(NOUT+1),' ') IF (FOUT(NOUT+1).NE.' ') THEN NOUT=NOUT+1 IF (NOUT.GT.MOUT) CALL CIEROR(27,TXTERR) CHOUT(NOUT)='MPQ11' ENDIF CALL RSEP3T('MPQ21',FOUT(NOUT+1),' ') IF (FOUT(NOUT+1).NE.' ') THEN NOUT=NOUT+1 IF (NOUT.GT.MOUT) CALL CIEROR(27,TXTERR) CHOUT(NOUT)='MPQ21' ENDIF CALL RSEP3T('MPQ31',FOUT(NOUT+1),' ') IF (FOUT(NOUT+1).NE.' ') THEN NOUT=NOUT+1 IF (NOUT.GT.MOUT) CALL CIEROR(27,TXTERR) CHOUT(NOUT)='MPQ31' ENDIF CALL RSEP3T('MPQ41',FOUT(NOUT+1),' ') IF (FOUT(NOUT+1).NE.' ') THEN NOUT=NOUT+1 IF (NOUT.GT.MOUT) CALL CIEROR(27,TXTERR) CHOUT(NOUT)='MPQ41' ENDIF CALL RSEP3T('MPQ12',FOUT(NOUT+1),' ') IF (FOUT(NOUT+1).NE.' ') THEN NOUT=NOUT+1 IF (NOUT.GT.MOUT) CALL CIEROR(27,TXTERR) CHOUT(NOUT)='MPQ12' ENDIF CALL RSEP3T('MPQ22',FOUT(NOUT+1),' ') IF (FOUT(NOUT+1).NE.' ') THEN NOUT=NOUT+1 IF (NOUT.GT.MOUT) CALL CIEROR(27,TXTERR) CHOUT(NOUT)='MPQ22' ENDIF CALL RSEP3T('MPQ32',FOUT(NOUT+1),' ') IF (FOUT(NOUT+1).NE.' ') THEN NOUT=NOUT+1 IF (NOUT.GT.MOUT) CALL CIEROR(27,TXTERR) CHOUT(NOUT)='MPQ32' ENDIF CALL RSEP3T('MPQ42',FOUT(NOUT+1),' ') IF (FOUT(NOUT+1).NE.' ') THEN NOUT=NOUT+1 IF (NOUT.GT.MOUT) CALL CIEROR(27,TXTERR) CHOUT(NOUT)='MPQ42' ENDIF CALL RSEP3T('MPQ13',FOUT(NOUT+1),' ') IF (FOUT(NOUT+1).NE.' ') THEN NOUT=NOUT+1 IF (NOUT.GT.MOUT) CALL CIEROR(27,TXTERR) CHOUT(NOUT)='MPQ13' ENDIF CALL RSEP3T('MPQ23',FOUT(NOUT+1),' ') IF (FOUT(NOUT+1).NE.' ') THEN NOUT=NOUT+1 IF (NOUT.GT.MOUT) CALL CIEROR(27,TXTERR) CHOUT(NOUT)='MPQ23' ENDIF CALL RSEP3T('MPQ33',FOUT(NOUT+1),' ') IF (FOUT(NOUT+1).NE.' ') THEN NOUT=NOUT+1 IF (NOUT.GT.MOUT) CALL CIEROR(27,TXTERR) CHOUT(NOUT)='MPQ33' ENDIF CALL RSEP3T('MPQ43',FOUT(NOUT+1),' ') IF (FOUT(NOUT+1).NE.' ') THEN NOUT=NOUT+1 IF (NOUT.GT.MOUT) CALL CIEROR(27,TXTERR) CHOUT(NOUT)='MPQ43' ENDIF CALL RSEP3T('MPQ14',FOUT(NOUT+1),' ') IF (FOUT(NOUT+1).NE.' ') THEN NOUT=NOUT+1 IF (NOUT.GT.MOUT) CALL CIEROR(27,TXTERR) CHOUT(NOUT)='MPQ14' ENDIF CALL RSEP3T('MPQ24',FOUT(NOUT+1),' ') IF (FOUT(NOUT+1).NE.' ') THEN NOUT=NOUT+1 IF (NOUT.GT.MOUT) CALL CIEROR(27,TXTERR) CHOUT(NOUT)='MPQ24' ENDIF CALL RSEP3T('MPQ34',FOUT(NOUT+1),' ') IF (FOUT(NOUT+1).NE.' ') THEN NOUT=NOUT+1 IF (NOUT.GT.MOUT) CALL CIEROR(27,TXTERR) CHOUT(NOUT)='MPQ34' ENDIF CALL RSEP3T('MPQ44',FOUT(NOUT+1),' ') IF (FOUT(NOUT+1).NE.' ') THEN NOUT=NOUT+1 IF (NOUT.GT.MOUT) CALL CIEROR(27,TXTERR) CHOUT(NOUT)='MPQ44' ENDIF C C Widths of Gaussian beams: CALL RSEP3T('GBW11',FOUT(NOUT+1),' ') IF (FOUT(NOUT+1).NE.' ') THEN NOUT=NOUT+1 IF (NOUT.GT.MOUT) CALL CIEROR(27,TXTERR) CHOUT(NOUT)='GBW11' IF(GBY11.LE.0.) THEN C MTTQ-01 CALL ERROR('MTTQ-01: Input parameter GBY11 is not positive') C To calculte and interpolate the sum of squares of Gaussian C beam widths identified by input parameter GBW11, positive C real-valued parameter GBY11 must be specified in the input C SEP parameter or history file. END IF ENDIF CALL RSEP3T('GBW22',FOUT(NOUT+1),' ') IF (FOUT(NOUT+1).NE.' ') THEN NOUT=NOUT+1 IF (NOUT.GT.MOUT) CALL CIEROR(27,TXTERR) CHOUT(NOUT)='GBW22' IF(GBY22.LE.0.) THEN C MTTQ-02 CALL ERROR('MTTQ-02: Input parameter GBY22 is not positive') C To calculte and interpolate the sum of squares of Gaussian C beam widths identified by input parameter GBW22, positive C real-valued parameter GBY22 must be specified in the input C SEP parameter or history file. END IF ENDIF CALL RSEP3T('GBW1',FOUT(NOUT+1),' ') IF (FOUT(NOUT+1).NE.' ') THEN NOUT=NOUT+1 IF (NOUT.GT.MOUT) CALL CIEROR(27,TXTERR) CHOUT(NOUT)='GBW1' IF(GBY11.LE.0.) THEN C MTTQ-03 CALL ERROR('MTTQ-03: Input parameter GBY11 is not positive') C To calculte and interpolate the sum of squares of Gaussian C beam widths identified by input parameter GBW1, positive C real-valued parameter GBY11 must be specified in the input C SEP parameter or history file. END IF ENDIF CALL RSEP3T('GBW2',FOUT(NOUT+1),' ') IF (FOUT(NOUT+1).NE.' ') THEN NOUT=NOUT+1 IF (NOUT.GT.MOUT) CALL CIEROR(27,TXTERR) CHOUT(NOUT)='GBW2' IF(GBY22.LE.0.) THEN C MTTQ-04 CALL ERROR('MTTQ-04: Input parameter GBY22 is not positive') C To calculte and interpolate the sum of squares of Gaussian C beam widths identified by input parameter GBW2, positive C real-valued parameter GBY22 must be specified in the input C SEP parameter or history file. END IF ENDIF CALL RSEP3T('AMP',FOUT(NOUT+1),' ') IF (FOUT(NOUT+1).NE.' ') THEN NOUT=NOUT+1 IF (NOUT.GT.MOUT) CALL CIEROR(27,TXTERR) CHOUT(NOUT)='AMP' ENDIF C C User-defined quantity: C Select string NAME identifying a new quantity and, simultaneously, C an input SEP parameter describing the name of the output file with C the new quantity. Then uncomment and modify the following lines C designed to read the filename of the output file with the new C quantity and to register the quantity if the output filename is C not blank. * CALL RSEP3T('NAME',FOUT(NOUT+1),' ') * IF (FOUT(NOUT+1).NE.' ') THEN * NOUT=NOUT+1 * IF (NOUT.GT.MOUT) CALL CIEROR(27,TXTERR) * CHOUT(NOUT)='NAME' * insert lines to check the values of input data * ENDIF C Repeat the above lines for each new quantity. C C End of subroutine MTTQ: RETURN C C======================================================================= C ENTRY MTTQRI(QNAME,QVALUE) C C Entry designed to compute the user-defined quantity at an initial C point of a ray. C C Input: C QNAME...Character string, which characterizes the quantity. C C Output: C QVALUE..Value of the quantity at the initial point of a ray. C C All the quantities computed by CRT at the initial point of the ray C under consideration are now stored in the common block C pointc.inc. The user may use all the C quantities from the common block and all his input parameters to C compute the value VALUE of his quantity at the initial point of the C ray. This value must be then assigned to the output quantity QVALUE. C C----------------------------------------------------------------------- C C Components of the 4*4 paraxial ray propagator matrix in RCCS: IF(QNAME(1:3).EQ.'MPQ') THEN IF(QNAME.EQ.'MPQ11') THEN QVALUE=1. ELSE IF(QNAME.EQ.'MPQ21') THEN QVALUE=0. ELSE IF(QNAME.EQ.'MPQ31') THEN QVALUE=0. ELSE IF(QNAME.EQ.'MPQ41') THEN QVALUE=0. ELSE IF(QNAME.EQ.'MPQ12') THEN QVALUE=0. ELSE IF(QNAME.EQ.'MPQ22') THEN QVALUE=1. ELSE IF(QNAME.EQ.'MPQ32') THEN QVALUE=0. ELSE IF(QNAME.EQ.'MPQ42') THEN QVALUE=0. ELSE IF(QNAME.EQ.'MPQ13') THEN QVALUE=0. ELSE IF(QNAME.EQ.'MPQ23') THEN QVALUE=0. ELSE IF(QNAME.EQ.'MPQ33') THEN QVALUE=1. ELSE IF(QNAME.EQ.'MPQ43') THEN QVALUE=0. ELSE IF(QNAME.EQ.'MPQ14') THEN QVALUE=0. ELSE IF(QNAME.EQ.'MPQ24') THEN QVALUE=0. ELSE IF(QNAME.EQ.'MPQ34') THEN QVALUE=0. ELSE IF(QNAME.EQ.'MPQ44') THEN QVALUE=1. END IF RETURN END IF C C Widths of Gaussian beams: IF (QNAME(1:3).EQ.'GBW') THEN C Ray-centred coordinate system at the initial point UU=YI(6)**2+YI(7)**2+YI(8)**2 U=SQRT(UU) H13=YI( 6)/U H23=YI( 7)/U H33=YI( 8)/U H11=YI( 9) H21=YI(10) H31=YI(11) H12=H23*H31-H33*H21 H22=H33*H11-H13*H31 H32=H13*H21-H23*H11 C Minus slowness gradient in the given vectors and RCCS U1=(YLI(4)*GBE11+YLI(5)*GBE21+YLI(6)*GBE31)*UU U2=(YLI(4)*GBE12+YLI(5)*GBE22+YLI(6)*GBE32)*UU U3=(YLI(4)* H13+YLI(5)* H23+YLI(6)* H33)*UU C Transformation matrix from the given vectors to RCCS C11=H11*GBE11+H21*GBE21+H31*GBE31 C21=H12*GBE11+H22*GBE21+H32*GBE31 C31=H13*GBE11+H23*GBE21+H33*GBE31 C12=H11*GBE12+H21*GBE22+H31*GBE32 C22=H12*GBE12+H22*GBE22+H32*GBE32 C32=H13*GBE12+H23*GBE22+H33*GBE32 CDET=C11*C22-C12*C21 C Projection of the real part of matrix M to the given vectors R11=GBR11+U1*C31+C31*U1-C31*C31*U3 R12=GBR12+U1*C32+C31*U2-C31*C32*U3 R22=GBR22+U2*C32+C32*U2-C31*C31*U3 C Sum of squares of Gaussian beam widths corresponding to GBY11 IF(QNAME.EQ.'GBW11') THEN QVALUE=(C11*C11+C21*C21)/GBY11 RETURN END IF IF(QNAME.EQ.'GBW1') THEN QVALUE=SQRT((C11*C11+C21*C21)/GBY11) RETURN END IF C Sum of squares of Gaussian beam widths corresponding to GBY22 IF(QNAME.EQ.'GBW22') THEN QVALUE=(C12*C12+C22*C22)/GBY22 RETURN END IF IF(QNAME.EQ.'GBW2') THEN QVALUE=SQRT((C12*C12+C22*C22)/GBY22) RETURN END IF END IF C C Amplitudes: IF(QNAME.EQ.'AMP') THEN QVALUE=AMPMAX RETURN END IF C C User-defined quantity: C Uncomment and modify the following lines, where NAME is the name C of your quantity. * IF(QNAME.EQ.'NAME') THEN * insert lines to calculate the value of a new quantity * QVALUE=value of a new quantity * RETURN * ENDIF C Repeat the above lines for each new quantity. C C End of entry MTTQRI: RETURN C C======================================================================= C ENTRY MTTQRA(QNAME,QVALUE) C C Entry designed to compute the user-defined quantity at other than C initial point of a ray. C C Input: C QNAME...Character string, which characterizes the quantity. C C Output: C QVALUE..Value of the quantity at the point of a ray. C C All the quantities computed by CRT at the point of the ray C under consideration are now stored in the common block C pointc.inc together with the values at C the initial point of the ray. The user may use all the C quantities from the common block and all his input parameters to C compute the value VALUE of his quantity at the point under C consideration. This value must be then assigned to the output C quantity QVALUE. C C----------------------------------------------------------------------- C C Components of the 4*4 paraxial ray propagator matrix in RCCS: IF(QNAME(1:3).EQ.'MPQ') THEN IF(QNAME.EQ.'MPQ11') THEN QVALUE=Y(12) ELSEIF (QNAME.EQ.'MPQ21') THEN QVALUE=Y(13) ELSEIF (QNAME.EQ.'MPQ31') THEN QVALUE=Y(14) ELSEIF (QNAME.EQ.'MPQ41') THEN QVALUE=Y(15) ELSEIF (QNAME.EQ.'MPQ12') THEN QVALUE=Y(16) ELSEIF (QNAME.EQ.'MPQ22') THEN QVALUE=Y(17) ELSEIF (QNAME.EQ.'MPQ32') THEN QVALUE=Y(18) ELSEIF (QNAME.EQ.'MPQ42') THEN QVALUE=Y(19) ELSEIF (QNAME.EQ.'MPQ13') THEN QVALUE=Y(20) ELSEIF (QNAME.EQ.'MPQ23') THEN QVALUE=Y(21) ELSEIF (QNAME.EQ.'MPQ33') THEN QVALUE=Y(22) ELSEIF (QNAME.EQ.'MPQ43') THEN QVALUE=Y(23) ELSEIF (QNAME.EQ.'MPQ14') THEN QVALUE=Y(24) ELSEIF (QNAME.EQ.'MPQ24') THEN QVALUE=Y(25) ELSEIF (QNAME.EQ.'MPQ34') THEN QVALUE=Y(26) ELSEIF (QNAME.EQ.'MPQ44') THEN QVALUE=Y(27) END IF RETURN END IF C C Widths of Gaussian beams: IF (QNAME(1:3).EQ.'GBW') THEN C Transforming matrices of geometrical spreading to given vectors Q111= Y(12)*C11+Y(16)*C21 Q211= Y(13)*C11+Y(17)*C21 Q121= Y(12)*C12+Y(16)*C22 Q221= Y(13)*C12+Y(17)*C22 Q112=( Y(20)*C22-Y(24)*C12)/CDET Q212=( Y(21)*C22-Y(25)*C12)/CDET Q122=(-Y(20)*C21+Y(24)*C11)/CDET Q222=(-Y(21)*C21+Y(25)*C11)/CDET C Dependence of the widths on real part R11,R12,R22 of matrix M Q111=Q111+Q112*R11+Q122*R12 Q211=Q211+Q212*R11+Q222*R12 Q121=Q121+Q112*R12+Q122*R22 Q221=Q221+Q212*R12+Q222*R22 C Sum of squares of Gaussian beam widths corresponding to GBY11 IF(QNAME.EQ.'GBW11') THEN QVALUE=(Q111*Q111+Q211*Q211)/GBY11+(Q112*Q112+Q212*Q212)*GBY11 RETURN END IF IF(QNAME.EQ.'GBW1') THEN QVALUE=SQRT( * (Q111*Q111+Q211*Q211)/GBY11+(Q112*Q112+Q212*Q212)*GBY11) RETURN END IF C Sum of squares of Gaussian beam widths corresponding to GBY22 IF(QNAME.EQ.'GBW22') THEN QVALUE=(Q121*Q121+Q221*Q221)/GBY22+(Q122*Q122+Q222*Q222)*GBY22 RETURN END IF IF(QNAME.EQ.'GBW2') THEN QVALUE=SQRT( * (Q121*Q121+Q221*Q221)/GBY22+(Q122*Q122+Q222*Q222)*GBY22) RETURN END IF END IF C C Amplitudes: IF(QNAME.EQ.'AMP') THEN CALL AP21(GREEN) QVALUE=0. DO 10 I=15,32 QVALUE=QVALUE+GREEN(I)**2 10 CONTINUE QVALUE=SQRT(QVALUE) QVALUE=AMIN1(QVALUE,AMPMAX) RETURN END IF C C User-defined quantity: C Uncomment and modify the following lines, where NAME is the name C of your quantity. * IF(QNAME.EQ.'NAME') THEN * insert lines to calculate the value of a new quantity * QVALUE=value of a new quantity * RETURN * ENDIF C Repeat the above lines for each new quantity. C C End of entry MTTQRA: RETURN END C C======================================================================= C