C
C======================================================================= C PROGRAM MTT C to calculate Multivalued Travel Times on a 3-D grid of points C C---------------------------------------------------------------------- C C The post processing program to the CRT program performs interpolation C of travel times to the gridpoints of arbitrary regular rectangular C 3-D grid of points. The interpolation inside ray tubes formed by C vertices of homogeneous triangles created during 3-D two-point C ray tracing is realized using the controlled initial-value ray tracing C algorithm. C If the I-th and (I+1)-th points of the three rays forming the ray tube C are situated inside the same complex block, the program performes C interpolation within ray cells whose bottom is formed by the I-th C points and top by (I+1)-th points of the rays. The best results can C thus be obtained if the travel time is the independent variable for C numerical integration in CRT program (default). C C This version of the program interpolates travel times within ray cells C trilinearly (in the coordinate system connected with each cell), C causing large travel-time errors inside wide ray cells. This should C be accounted for when interpolating travel times in the areas of large C geometrical spreading. C C The boundaries between ray cells are determined numerically, with C different rounding errors in different cells. The rounding errors C produce narrow strips along the boundaries between the cells, where C the travel times of the corresponding history are randomly missing C or doubled. If a gridpoint falls, by a chance, into the strip, C the corresponding arrival at the gridpoint may be missing or doubled. C A small percentage of gridpoints with missing or doubled arrivals is C thus inevitable with the present interpolation algorithm. C C The interpolation is not performed in demarcation belts between C different ray histories. The maximum width of the belts is controlled C by input parameter AERR of the CRT C program. The typical width of the belts, measured in take-off angles, C is 0.0001 radians in order of magnitude. The belts may become wide C in areas of large geometrical spreading. The division of rays into C different histories is, by default, very detailed and is controlled C by input parameter PRM0(3) of the CRT C program. Such a behaviour is useful for two-point ray tracing but has C some awkward consequences on the interpolation. For example, the C demarcation belt between the rays incident at different sides of the C computational volume for ray tracing extends across the whole model, C causing an artificial gap within the continuous travel times. C If the ray history does not account for the termination of a ray at C different sides of the computational volume, the gaps corresponding C to the edges of the computational volume are removed, but the corners C along the edges are not filled with the ray cells any more. Then the C computational volume for ray tracing should sufficiently exceed the C extent of target volume covered by the grid for interpolation. C Input parameter PRM0(3) of the CRT C program may thus considerably influence the results of the C interpolation. However, gaps due to the demarcation belts C corresponding to structural interfaces cannot be removed within C the present interpolation algorithm. C C....................................................................... C C C Main input data read from external interactive device (*): C The data consist of character strings, read by list C directed (free format) input. The strings have thus to be C enclosed in apostrophes. C 'SEP','INPUTFILE',/ C 'SEP'... Name of the file describing the configuration of C the target grid, where the output quantities are C to be interpolated. C Description of file SEP C 'INPUTFILE' ... Name of file with input data. C Defaults: 'SEP'='grd.h', 'INPUTFILE'='mtt.dat' C C Structure of the input data file 'INPUTFILE': C The data consist of following lines of C character strings, read by list directed (free format) input: C (1) 'CRT-R','CRT-I','CRT-T',/ C 'CRT-R'...Name of the input unformatted file with the quantities C stored along rays (see C.R.T.5.5.1). C Description of file CRT-R C 'CRT-I'...Name of the input unformatted file with the quantities C at the initial points of rays, corresponding to file C CRT-R (see C.R.T.6.1). C Description of file CRT-I C 'CRT-T'...Name of the input formatted file with the indices of C rays forming the homogeneous triangles in the C ray-parameter domain. C Description of file CRT-T C Defaults: 'CRT-R'='r01.out', 'CRT-I'='r01i.out', C 'CRT-T'='t01.out'. C (2) 'NUM',/ C 'NUM'... Name of the output formatted file with N1*N2*N3 array of C integer values, containing the number of arrivals at C each gridpoint of the target grid. C Description of file NUM C Default: 'NUM'='num.out' C (3) Zero to MOUT times following line: C 'FOUT(i)','OUT(i)',/ C 'FOUT(i)'.Name of the output formatted file with the array of C values, containing for each gridpoint the quantities C of all determined arrivals. The number of values C thus equals the sum of the integers in file 'NUM'. C The kind of computed quantities written into this C file is given by value of 'OUT(i)'. C Description of files 'FOUT(i)' C 'OUT(i)'..Character string, describing which quantities are to be C written into corresponding output file 'FOUT(i)'. C Following is a list of possible values of 'OUT(i)' and C corresponding quantities written into output files: C 'tt' ....... Travel times C 'hi' ....... Ray histories C Defaults: 'FOUT(1)' ='mtt.out', 'OUT(1)' ='tt', C 'FOUT(2)' =' ' , 'OUT(2)' =' ' , C . . C . . C 'FOUT(MOUT)'=' ' , 'OUT(MOUT)=' ' . C C C C Input formatted input SEP parameter file: C The file with coordinates of the origin of the target grid C O1, O2, O3, numbers of gridpoints N1, N2, N3 and spatial step C between gridpoints D1, D2, D3. The file has the form of SEP file, C for the description of the SEP format refer C to the file 'sep.for'. C N1=positive integer... Number of gridpoints along the X1 axis. C Default: N1=1 C N2=positive integer... Number of gridpoints along the X2 axis. C Default: N2=1 C N3=positive integer... Number of gridpoints along the X3 axis. C Default: N3=1 C O1=real, O2=real, O3=real... Coordinates of the first point of the C grid. Default: O1=0, O2=0, O3=0. C D1=real... Grid interval along the X1 axis. Default: D1=1. C D2=real... Grid interval along the X2 axis. Default: D2=1. C D3=real... Grid interval along the X3 axis. Default: D3=1. C Example of data SEP C C C Output formatted file 'NUM': C N1*N2*N3 integers: For each gridpoint, the number of travel C times determined at the gridpoint. The innermost loop corresponds C to N1, the outermost loop corresponds to N3. C C C Output formatted file 'FOUT(i)': C N numbers, where N is the sum of integers in file NUM. C For each gridpoint, all quantities determined at the gridpoint. C The kind of quantities is determined by corresponding C parameter OUT(i), for example if OUT(1)='TT', then travel times C are written into output file named 'FOUT(1)'. C C....................................................................... C C Date: 1997, October 25 C C Coded by Petr Bulant C bulant@seis.karlov.mff.cuni.cz C C----------------------------------------------------------------------- C C Common block /RAM/ to store the information about triangles and rays: INCLUDE 'mtt.inc' C ........................... C Common block /CIGRD/ to store the parameters of the C interpolation grid. REAL O1,O2,O3,D1,D2,D3 INTEGER N1,N2,N3 COMMON/CIGRD/O1,O2,O3,D1,D2,D3,N1,N2,N3 SAVE /CIGRD/ C O1,O2,O3 ... Coordinates of the origin of the grid. C D1,D2,D3 ... Steps per gridpoint. C N1,N2,N3 ... Numbers of gridpoints. C....................................................................... C External functions required: C C Auxiliary storage locations: INTEGER IRAY1,IRAY2,IRAY3 INTEGER I1,I2,I3,I4,I5,INDX,IT1 INTEGER I,J,K,L,M,N INTEGER LU0,LU1,LU2,LU3,LU4,LU5 PARAMETER (LU0=10,LU1=1,LU2=2,LU3=3,LU4=4,LU5=5) CHARACTER*80 FILE0,FILE1,FILE2,FILE3,FILE4,FILE5 C C----------------------------------------------------------------------- C WRITE(*,'(''+'',79('' ''))') WRITE(*,'(A)') '+Reading input data.' C Reading in a name of the file with the data for interpolation C grid, and a name of the main input data file: FILE4='grd.h' FILE0='mtt.dat' WRITE(*,'(A)') * ' Enter the names of the input data files (SEP,INPUTFILE): ' READ(*,*) FILE4,FILE0 OPEN(LU0,FILE=FILE0,FORM='FORMATTED',STATUS='OLD') C Reading in filenames of the files with C computed triangles and rays: FILE1='r01.out' FILE2='r01i.out' FILE3='t01.out' READ(LU0,*) FILE1,FILE2,FILE3 OPEN(LU1,FILE=FILE1,FORM='UNFORMATTED',STATUS='OLD') OPEN(LU2,FILE=FILE2,FORM='UNFORMATTED',STATUS='OLD') OPEN(LU3,FILE=FILE3,FORM='FORMATTED',STATUS='OLD') C Reading in filenames of the output files together C with quantities describing what is to be written into the files: OUT(1)='tt' OUT(2)=' ' OUT(3)=' ' OUT(4)=' ' OUT(5)=' ' FILE5='num.out' FOUT(1)='mtt.out' FOUT(2)=' ' FOUT(3)=' ' FOUT(4)=' ' FOUT(5)=' ' READ(LU0,*) FILE5 I1=1 2 CONTINUE READ(LU0,*,END=4) FOUT(I1),OUT(I1) I1=I1+1 IF (I1.GT.MOUT) THEN C MTT-19 PAUSE 'Error MTT-19: Small array FOUT for the output files.' STOP C This error may be caused by wrongly prescribed C input data FOUT and OUT, or by C small parameter MOUT defined at the top of this file. ENDIF GOTO 2 4 CONTINUE C Checking the requests for output files: NOUT=0 DO 6, I1=1,MOUT CALL LOWER(OUT(I1)) IF ((OUT(I1)(1:2).NE.' ').AND.(OUT(I1)(1:2).NE.'tt').AND. * (OUT(I1)(1:2).NE.'hi')) THEN C MTT-20 PAUSE 'Error MTT-20: Wrong value of OUT.' STOP C This error may be caused by wrongly prescribed value of C input data OUT. ELSEIF (((OUT(I1)(1:2).EQ.' ').AND.(FOUT(I1).NE.' ')) .OR. * ((OUT(I1)(1:2).NE.' ').AND.(FOUT(I1).EQ.' '))) THEN C MTT-21 PAUSE 'Error MTT-21: Inconsistency between FOUT and OUT.' STOP C Either both of FOUT(i) and OUT(i) or none of them must be C prescribed. See the description of C input data FOUT and OUT. ELSEIF ((OUT(I1)(1:2).NE.' ').AND.(FOUT(I1).NE.' ')) THEN NOUT=NOUT+1 FOUT(NOUT)=FOUT(I1) OUT(NOUT)=OUT(I1)(1:2) ENDIF 6 CONTINUE NQ=5+NOUT C C Reading in the file with parameters of target grid: CALL CIGRID(LU4,FILE4) C C C Reading in file with computed triangles, C sorting the rays in triangles: NT=0 NRAMP=0 IRAYMI=999999 IRAYMA=0 10 CONTINUE READ(LU3,*,END=20) IRAY1,IRAY2,IRAY3 IF (IRAY1.GT.IRAYMA) IRAYMA=IRAY1 IF (IRAY2.GT.IRAYMA) IRAYMA=IRAY2 IF (IRAY3.GT.IRAYMA) IRAYMA=IRAY3 IF (IRAY1.LT.IRAYMI) IRAYMI=IRAY1 IF (IRAY2.LT.IRAYMI) IRAYMI=IRAY2 IF (IRAY3.LT.IRAYMI) IRAYMI=IRAY3 IF (IRAY1.LT.IRAY2) THEN I=IRAY1 IRAY1=IRAY2 IRAY2=I ENDIF IF (IRAY2.LT.IRAY3) THEN I=IRAY2 IRAY2=IRAY3 IRAY3=I ENDIF IF (IRAY1.LT.IRAY2) THEN I=IRAY1 IRAY1=IRAY2 IRAY2=I ENDIF NRAMP=NRAMP+1 IF (NRAMP.GT.MRAMP) CALL CIEROR(1) IRAM(NRAMP)=IRAY1 NRAMP=NRAMP+1 IF (NRAMP.GT.MRAMP) CALL CIEROR(1) IRAM(NRAMP)=IRAY2 NRAMP=NRAMP+1 IF (NRAMP.GT.MRAMP) CALL CIEROR(1) IRAM(NRAMP)=IRAY3 NT=NT+1 GOTO 10 20 CONTINUE CLOSE(LU3) NR=IRAYMA-IRAYMI+1 C C Sorting the triangles: DO 40, I1=NRAMP-5,1,-3 DO 30, I2=1,I1,3 L=I2+3 IF (IRAM(I2).GT.IRAM(L)) THEN J=I2+1 K=I2+2 M=I2+4 N=I2+5 I =IRAM(I2) IRAM(I2)=IRAM(L) IRAM(L) =I I =IRAM(J) IRAM(J) =IRAM(M) IRAM(M) =I I =IRAM(K) IRAM(K) =IRAM(N) IRAM(N) =I ENDIF 30 CONTINUE 40 CONTINUE C C C Forming an auxiliary array with information about when can be rays C erased from the memory ("deleting array"): IF (NRAMP+NR.GT.MRAMP) CALL CIEROR(1) DO 50, I1=NRAMP+1,NRAMP+NR IRAM(I1)=0 50 CONTINUE NRAMP=NRAMP+NR ORAYE=IRAYMI-1-3*NT DO 60, I1=1,3*NT,3 IRAM(IRAM(I1 )-ORAYE)=IRAM(I1) IRAM(IRAM(I1+1)-ORAYE)=IRAM(I1) IRAM(IRAM(I1+2)-ORAYE)=IRAM(I1) 60 CONTINUE C C C Forming an auxiliary array with information about addresses C of the ends of records for rays in array RAM ("addressing array"): C "Ray" IRAYMI-1: NRAMP=NRAMP+1 IF (NRAMP.GT.MRAMP) CALL CIEROR(1) IRAM(NRAMP)=NRAMP+NR C All other rays: IF (NRAMP+NR.GT.MRAMP) CALL CIEROR(1) DO 70, I1=NRAMP+1,NRAMP+NR IRAM(I1)=0 70 CONTINUE NRAMP=NRAMP+NR ORAYA=IRAYMI-1-3*NT-NR-1 C C C Loop for all the triangles: IT1=1 I1=-1 100 CONTINUE I2=NINT(100*IT1/3/NT) IF (I2.NE.I1) THEN WRITE(*,'(''+'',79('' ''))') WRITE(*,'(A,1I4,A)') '+Interpolating. ',I2,'% completed.' I1=I2 ENDIF C C If necessary, reading in new rays: IF * (((IRAM(IRAM(IT1)-ORAYA+1).EQ.0).AND.(IRAM(IT1).NE.IRAYMA)) .OR. * ((IRAM(IRAM(IT1)-ORAYA ).EQ.0).AND.(IRAM(IT1).EQ.IRAYMA))) * CALL CIREAD(LU1,LU2,IT1) C C Empting the array RAM to enable writing of possible interpolated C quantities: IF ((MRAMP-NRAMP).LT.(NINT(MRAM/10))) CALL CIERAS C C Interpolation along the rays of the triangle: CALL CITUBE(IT1) C IT1=IT1+3 IF (IT1.LT.3*NT) GOTO 100 C End of the loop for all the triangles. CLOSE(LU1) CLOSE(LU2) C C C Sorting the results according to the travel times: DO 340, I4=1,NOUT IF (OUT(I4)(1:2).EQ.'tt') THEN DO 330, I3=0,N3-1 DO 320, I2=0,N2-1 DO 310, I1=0,N1-1 300 CONTINUE INDX=I3*N1*N2+I2*N1+I1 INDX=MRAM-INDX INDX1=IRAM(INDX) IF (INDX1.EQ.0) GOTO 306 302 CONTINUE INDX2=IRAM(INDX1) IF (INDX2.EQ.0) GOTO 306 IF (RAM(INDX2+I4).LT.RAM(INDX1+I4)) THEN IRAM(INDX1)=IRAM(INDX2) IRAM(INDX2)=INDX1 IRAM(INDX)=INDX2 GOTO 300 ENDIF INDX=INDX1 INDX1=INDX2 GOTO 302 306 CONTINUE 310 CONTINUE 320 CONTINUE 330 CONTINUE ENDIF 340 CONTINUE C C Writing the results of interpolation: C Writing numbers of computed arrivals: IF (N1*N2*N3.GT.MRAMP) CALL CIEROR(1) NRAMP=0 I5=0 DO 130, I3=0,N3-1 DO 120, I2=0,N2-1 DO 110, I1=0,N1-1 NRAMP=NRAMP+1 I4=0 INDX=I3*N1*N2+I2*N1+I1 INDX=MRAM-INDX 105 CONTINUE C Check for the consistentcy of the results IF (INDX.LE.MRAMP.OR.INDX.GT.MRAM) THEN C MTT-14 PAUSE 'Error MTT-14: Disorder in array RAM' STOP C This error should not appear. C Please contact the author. ENDIF INDX=IRAM(INDX) IF (INDX.EQ.0) THEN IRAM(NRAMP)=I4 I5=I5+I4 GOTO 110 ENDIF I4=I4+1 GOTO 105 110 CONTINUE 120 CONTINUE 130 CONTINUE IF (NRAMP.NE.N1*N2*N3) THEN C MTT-15 PAUSE 'Error MTT-15: Disorder in array RAM' STOP C This error should not appear. C Please contact the author. ENDIF CALL WARRAI(LU5,FILE5,'FORMATTED',.FALSE.,0,.FALSE.,0,N1*N2*N3, * IRAM) C C C Writing interpolated quantities: IF (I5.GT.MRAMP) CALL CIEROR(1) DO 240, I4=1,NOUT NRAMP=0 DO 230, I3=0,N3-1 DO 220, I2=0,N2-1 DO 210, I1=0,N1-1 INDX=I3*N1*N2+I2*N1+I1 INDX=MRAM-INDX IF (IRAM(INDX).EQ.0) GOTO 206 INDX=IRAM(INDX) 205 CONTINUE C Check for the consistentcy of the results IF (INDX.LE.MRAMP.OR.INDX.GT.MRAM) THEN C MTT-16 PAUSE 'Error MTT-16: Disorder in array RAM' STOP C This error should not appear. C Please contact the author. ENDIF NRAMP=NRAMP+1 IF (OUT(I4)(1:2).EQ.'tt') THEN C Travel time: RAM(NRAMP)=RAM(INDX+I4) ELSEIF (OUT(I4)(1:2).EQ.'hi') THEN C Ray history: IRAM(NRAMP)=IRAM(INDX+I4) ENDIF INDX=IRAM(INDX) IF (INDX.NE.0) GOTO 205 206 CONTINUE 210 CONTINUE 220 CONTINUE 230 CONTINUE IF (NRAMP.NE.I5) THEN C MTT-17 PAUSE 'Error MTT-17: Disorder in array RAM' STOP C This error should not appear. C Please contact the author. ENDIF IF (OUT(I4)(1:2).EQ.'tt') THEN CALL WARRAY(LU5,FOUT(I4),'FORMATTED',.FALSE.,0.,.FALSE.,0.,I5, * RAM) ELSE CALL WARRAI(LU5,FOUT(I4),'FORMATTED',.FALSE.,0,.FALSE.,0,I5, * IRAM) ENDIF 240 CONTINUE C STOP END C C C======================================================================= C SUBROUTINE CIREAD(LU1,LU2,IT1) C C---------------------------------------------------------------------- C Subroutine to read the unformatted output of program CRT and C to write it into array (I)RAM used in program MTT. C Reading the output files is completed by a simple invocation of C subroutine AP00 of file 'ap.for'. C INTEGER LU1,LU2,IT1 C Input: C LU1 ... Number of logical unit corresponding to the file with C the quantities stored along rays (see C.R.T.5.5.1). C LU2 ... Number of logical unit corresponding to the file with C the quantities at the initial points of rays, C corresponding to file 'RAYPOINTS' (see C.R.T.6.1). C IT1 ... Position of the first ray of the actually processed C triangle in the array IRAM. C No output. C C Coded by Petr Bulant C 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 Common block /RAM/ to store the information about triangles and rays: INCLUDE 'mtt.inc' C....................................................................... C C Auxiliary storage locations: INTEGER IRAY0,I1 C IRAY0.. Index of the last ray read in into the array RAM. C C----------------------------------------------------------------------- C Loop for the points of rays: 10 CONTINUE IF ((NRAMP+2*NQ).GT.MRAMP) THEN C Empting the memory: CALL CIERAS IF ((NRAMP+2*NQ).GT.MRAMP) CALL CIEROR(1) ENDIF C Reading the results of the complete ray tracing: IRAY0=IRAY CALL AP00(0,LU1,LU2) IF (IWAVE.LT.1) THEN C End of rays: IF (IRAY0.LT.IRAM(IT1)) THEN C MTT-02 PAUSE * 'Error MTT-02: end of file with rays, wanted ray not found' STOP C This error may be caused by K2P C not equal to zero, then only two-point rays are stored in C output files of CRT. ENDIF RETURN ENDIF IF (IRAY.LT.IRAYMI) GOTO 10 IF (IRAM(IRAY-ORAYE).NE.0) THEN C Writing the results of the complete ray tracing - recording C new point on a ray: IF (IPT.LE.1) THEN IF (ISHEET.EQ.0) THEN C MTT-18 PAUSE * 'Error MTT-18: a ray with history equal to 0 in CIREAD.' STOP C All the rays are of the same history, the interpolation C in such a case makes no sense. Check the value of the C parameter IPOINT in file C rpar.for. ENDIF C New ray - recording the initial point: C Coordinates: RAM(NRAMP+1)=YI(3) RAM(NRAMP+2)=YI(4) RAM(NRAMP+3)=YI(5) C Index of surface: IRAM(NRAMP+4)=0 C Sequential index of point: IRAM(NRAMP+5)=0 C Quantities to be interpolated: DO 20, I1=1,NOUT IF (OUT(I1)(1:2).EQ.'tt') THEN C Travel time: RAM(NRAMP+(5+I1))=YI(1) ELSEIF (OUT(I1)(1:2).EQ.'hi') THEN C Ray history: IRAM(NRAMP+(5+I1))=ISHEET ENDIF 20 CONTINUE NRAMP=NRAMP+NQ ENDIF C Recording the new point: C Coordinates: RAM(NRAMP+1)=Y(3) RAM(NRAMP+2)=Y(4) RAM(NRAMP+3)=Y(5) C Index of surface: IRAM(NRAMP+4)=ISRF C Sequential index of point: IRAM(NRAMP+5)=IPT C Quantities to be interpolated: DO 30, I1=1,NOUT IF (OUT(I1)(1:2).EQ.'tt') THEN C Travel time: RAM(NRAMP+(5+I1))=Y(1) ELSEIF (OUT(I1)(1:2).EQ.'hi') THEN C Ray history: IRAM(NRAMP+(5+I1))=ISHEET ENDIF 30 CONTINUE NRAMP=NRAMP+NQ ENDIF IRAM(IRAY-ORAYA)=NRAMP IF ((IPT.LE.1).AND.(IRAY.GT.IRAM(IT1)).AND.(IRAY.NE.IRAYMA)) * RETURN GOTO 10 C END C C C======================================================================= C SUBROUTINE CIERAS C C---------------------------------------------------------------------- C Subroutine for empting the array (I)RAM. All the parameters C of all the rays, which will no more be used, are erased. C C No input. C No output. C C Coded by Petr Bulant C C ........................... C Common block /POINTC/ to store the results of complete ray tracing: INCLUDE 'pointc.inc' C IRAY .. Index of the ray being actually read in by CIREAD. C This procedure supposes, that any ray with higher C index than IRAY was not read in. C None of the storage locations of the common block are altered. C ........................... C Common block /RAM/ to store the information about triangles and rays: INCLUDE 'mtt.inc' C....................................................................... C Auxiliary storage locations: INTEGER I1,I2,J1 INTEGER IADDRP C I1 ... Controls main loop over rays. C I2 ... Controls the loop over parameters of ray I1. C J1 ... address of the last used record of array RAM. C C----------------------------------------------------------------------- J1=IRAM(IRAYMI-1-ORAYA) IADDRP=J1 C Loop for the rays: DO 20, I1=IRAYMI,IRAY IF (IRAM(I1-ORAYE).GE.(IRAY-1)) THEN C This ray is not to be erased: DO 10, I2=IADDRP+1,IRAM(I1-ORAYA) J1=J1+1 IRAM(J1)=IRAM(I2) 10 CONTINUE IADDRP=IRAM(I1-ORAYA) IRAM(I1-ORAYA)=J1 ELSE C This ray is to be erased: IRAM(I1-ORAYE)=0 IADDRP=IRAM(I1-ORAYA) IRAM(I1-ORAYA)=J1 ENDIF 20 CONTINUE NRAMP=J1 RETURN END C======================================================================= C SUBROUTINE CITUBE(IT1) C C---------------------------------------------------------------------- C Subroutine for interpolation within ray tube formed by the rays C IRAM(IT1), IRAM(IT1+1), IRAM(IT1+2). C INTEGER IT1 C Input: C IT1 ... The address of the index of the first ray of the ray tube, C in which the interpolation is to be performed. C No output. C C Coded by Petr Bulant C C ........................... C Common block /RAM/ to store the information about triangles and rays: INCLUDE 'mtt.inc' C....................................................................... INTEGER I1A,I2A,I3A,I1B,I2B,I3B INTEGER I1MA,I2MA,I3MA,J1,J2 C I1A,I2A,I3A,I1B,I2B,I3B ... Integers controling main loop over C points on the rays (along ray tube). Numbers 1,2,3 denote C the first, second and third ray, character A denotes bottom of C the ray cell and B denotes top of the ray cell. C Each integer contains the address just before the parameters C of the corresponding point: C the first parameter of the first point: RAM(I1A+1) C J1 ... Counts the number of identical points of the ray cell C ( J1=0,1,2,3 ). C J2 ... Displaies maximum sequential index of a point allowed C for actual ray cell. C C----------------------------------------------------------------------- IF ((IRAM(IRAM(IT1 )-ORAYA).EQ.0).OR. * (IRAM(IRAM(IT1+1)-ORAYA).EQ.0).OR. * (IRAM(IRAM(IT1+2)-ORAYA).EQ.0)) THEN C MTT-03 PAUSE * 'Error MTT-03: parameters of a ray not found in memory in CITUBE' STOP C This error may be caused by K2P C not equal to zero, then only two-point rays are stored in C output files of CRT. ENDIF C I1MA=IRAM(IRAM(IT1 )-ORAYA)-NQ I2MA=IRAM(IRAM(IT1+1)-ORAYA)-NQ I3MA=IRAM(IRAM(IT1+2)-ORAYA)-NQ C Skipping the first ray cell: I1A=IRAM(IRAM(IT1 )-ORAYA-1)+NQ I2A=IRAM(IRAM(IT1+1)-ORAYA-1)+NQ I3A=IRAM(IRAM(IT1+2)-ORAYA-1)+NQ C Loop for points on the rays (loop for ray cells): 10 CONTINUE J1=3 J2=MAX0(IRAM(I1A+5),IRAM(I2A+5),IRAM(I3A+5)) IF ((IRAM(I1A+5).EQ.IRAM(I2A+5)).AND. * (IRAM(I1A+5).EQ.IRAM(I3A+5))) J2=J2+1 IF ((IRAM(I1A+4).EQ.0).AND.(IRAM(I1A+5).LT.J2)) THEN I1B=I1A+NQ J1=J1-1 ELSE I1B=I1A ENDIF IF ((IRAM(I2A+4).EQ.0).AND.(IRAM(I2A+5).LT.J2)) THEN I2B=I2A+NQ J1=J1-1 ELSE I2B=I2A ENDIF IF ((IRAM(I3A+4).EQ.0).AND.(IRAM(I3A+5).LT.J2)) THEN I3B=I3A+NQ J1=J1-1 ELSE I3B=I3A ENDIF IF (J1.EQ.3) THEN IF ((IRAM(I1A+4).EQ.IRAM(I2A+4)).AND. * (IRAM(I1A+4).EQ.IRAM(I3A+4))) THEN C Crossing the interface: I1A=I1A+NQ I2A=I2A+NQ I3A=I3A+NQ J1=0 IF ((I1A.GT.I1MA).OR.(I2A.GT.I2MA).OR.(I3A.GT.I3MA)) * RETURN IRAM(I1A+4)=0 IRAM(I2A+4)=0 IRAM(I3A+4)=0 GOTO 10 ELSE IF (IRAM(I1A+4).EQ.0) THEN I1B=I1A+NQ J1=J1-1 ELSE I1B=I1A ENDIF IF (IRAM(I2A+4).EQ.0) THEN I2B=I2A+NQ J1=J1-1 ELSE I2B=I2A ENDIF IF (IRAM(I3A+4).EQ.0) THEN I3B=I3A+NQ J1=J1-1 ELSE I3B=I3A ENDIF IF (J1.EQ.3) THEN IF ((I1A.GE.I1MA).AND.(I2A.GE.I2MA).AND.(I3A.GE.I3MA)) * RETURN C MTT-04 PAUSE * 'Error MTT-04: the points reached different interfaces.' STOP C This error should not appear. C Please contact the author. ENDIF ENDIF ENDIF IF ((I1B.GT.I1MA).OR.(I2B.GT.I2MA).OR.(I3B.GT.I3MA)) THEN C MTT-05 PAUSE 'Error MTT-05: point exceded the ray.' STOP C This error should not appear. C Please contact the author. ENDIF C C The ray cell formed by points, whose parameters can be found C just behind the addresses I1A,I2A,I3A,I1B,I2B,I3B, C is prepared for interpolation: IF (((RAM(I1A+1).NE.RAM(I2A+1)).OR.(RAM(I1A+2).NE.RAM(I2A+2)) * .OR.(RAM(I1A+3).NE.RAM(I2A+3))) .AND. * ((RAM(I1A+1).NE.RAM(I3A+1)).OR.(RAM(I1A+2).NE.RAM(I3A+2)) * .OR.(RAM(I1A+3).NE.RAM(I3A+3))) .AND. * ((RAM(I2A+1).NE.RAM(I3A+1)).OR.(RAM(I2A+2).NE.RAM(I3A+2)) * .OR.(RAM(I2A+3).NE.RAM(I3A+3))) .AND. * ((RAM(I1B+1).NE.RAM(I2B+1)).OR.(RAM(I1B+2).NE.RAM(I2B+2)) * .OR.(RAM(I1B+3).NE.RAM(I2B+3))) .AND. * ((RAM(I1B+1).NE.RAM(I3B+1)).OR.(RAM(I1B+2).NE.RAM(I3B+2)) * .OR.(RAM(I1B+3).NE.RAM(I3B+3))) .AND. * ((RAM(I2B+1).NE.RAM(I3B+1)).OR.(RAM(I2B+2).NE.RAM(I3B+2)) * .OR.(RAM(I2B+3).NE.RAM(I3B+3)))) THEN C Standard ray cell formed by 4, 5 or 6 points. The bottom C and the top of the ray cell are triangles (i.e. they are C formed by different points), the three sides of the ray cell C are formed by lines, triangles or tetragons. CALL CIINTP(I1A,I2A,I3A,I1B,I2B,I3B) ELSE C TADY TO JESTE NENI HOTOVE !!!!!!! CONTINUE ENDIF I1A=I1B I2A=I2B I3A=I3B GOTO 10 C End of the loop along the ray tube. END C======================================================================= C SUBROUTINE CIINTP(I1A,I2A,I3A,I1B,I2B,I3B) C C---------------------------------------------------------------------- C Subroutine for interpolation within standard ray cell formed by the C points I1A,I2A,I3A,I1B,I2B,I3B. C INTEGER I1A,I2A,I3A,I1B,I2B,I3B C Input: C I1A,I2A,I3A,I1B,I2B,I3B ... Integers specifying the six C points on the rays, whch define the ray cell. C Numbers 1,2,3 denote the first, second and third ray, C character A denotes bottom of the ray cell and B C denotes top of the ray cell. C Each integer contains the address just before the parameters C of the corresponding point: C the first parameter of the first point: RAM(I1A+1) C No output. C C Subroutines and external functions required: EXTERNAL CIPPP,CILSIG,CICUB,CIKVA,CICUBR,CIKVAR,CIDET3,CISUBD REAL CIPPP,CICUB,CIKVA,CIDET3,CISUBD LOGICAL CILSIG C C Coded by Petr Bulant C C ........................... C Common block /RAM/ to store the information about triangles and rays: INCLUDE 'mtt.inc' C ........................... C Common block /CIGRD/ to store the parameters of the C interpolation grid. REAL O1,O2,O3,D1,D2,D3 INTEGER N1,N2,N3 COMMON/CIGRD/O1,O2,O3,D1,D2,D3,N1,N2,N3 SAVE /CIGRD/ C O1,O2,O3 ... Coordinates of the origin of the grid. C D1,D2,D3 ... Grid steps. C N1,N2,N3 ... Numbers of gridpoints. C....................................................................... REAL ZERO PARAMETER (ZERO=0.000001) INTEGER K1MI,K1MA,K2MI,K2MA,K3MI,K3MA INTEGER I1,I2,I3,I4,I5 REAL X1,X2,X3 REAL Y1,Y2,Y3 REAL XX(3,3),ZZ(3,3) REAL AAA,BBB,CCC,DDD,EEE,FFF INTEGER NRR REAL RR(4),CRR(4),ROOT(3) INTEGER IR REAL R1,R2 INTEGER INDX REAL ERRMAX REAL WB,WC,W1,W2,W3 REAL A11,A21,A31,A12,A22,A32,A13,A23,A33,B11,B21,B31,B12,B22,B32 REAL A1,A2,UU11,UU21,UU12,UU22,VV1,VV2,DETUU,TT C----------------------------------------------------------------------- C Choosing gridpoints, which may be contained in the ray cell: K1MI=MAX0(0,NINT(( * ( AMIN1(RAM(I1A+1),RAM(I2A+1),RAM(I3A+1), * RAM(I1B+1),RAM(I2B+1),RAM(I3B+1)) - O1 ) / D1)+0.5) ) K1MA=MIN0(N1-1,NINT(( * ( AMAX1(RAM(I1A+1),RAM(I2A+1),RAM(I3A+1), * RAM(I1B+1),RAM(I2B+1),RAM(I3B+1)) - O1 ) / D1)-0.5) ) IF (K1MA.LT.K1MI) RETURN K2MI=MAX0(0,NINT(( * ( AMIN1(RAM(I1A+2),RAM(I2A+2),RAM(I3A+2), * RAM(I1B+2),RAM(I2B+2),RAM(I3B+2)) - O2 ) / D2)+0.5) ) K2MA=MIN0(N2-1,NINT(( * ( AMAX1(RAM(I1A+2),RAM(I2A+2),RAM(I3A+2), * RAM(I1B+2),RAM(I2B+2),RAM(I3B+2)) - O2 ) / D2)-0.5) ) IF (K2MA.LT.K2MI) RETURN K3MI=MAX0(0,NINT(( * ( AMIN1(RAM(I1A+3),RAM(I2A+3),RAM(I3A+3), * RAM(I1B+3),RAM(I2B+3),RAM(I3B+3)) - O3 ) / D3)+0.5) ) K3MA=MIN0(N3-1,NINT(( * ( AMAX1(RAM(I1A+3),RAM(I2A+3),RAM(I3A+3), * RAM(I1B+3),RAM(I2B+3),RAM(I3B+3)) - O3 ) / D3)-0.5) ) IF (K3MA.LT.K3MI) RETURN C C Checking whether all the rays of the top (bottom) of the ray cell C are on the same side of the bottom (top) of the ray cell: IF (I1A.NE.I1B.AND.I2A.NE.I2B.AND.I3A.NE.I3B) THEN IF (.NOT.CILSIG( * CIPPP(I1A,I2A,I1A,I3A,I1A,RAM(I1B+1),RAM(I1B+2),RAM(I1B+3)), * CIPPP(I1A,I2A,I1A,I3A,I1A,RAM(I2B+1),RAM(I2B+2),RAM(I2B+3)), * CIPPP(I1A,I2A,I1A,I3A,I1A,RAM(I3B+1),RAM(I3B+2),RAM(I3B+3)), * 0.)) RETURN IF (.NOT.CILSIG( * CIPPP(I1B,I2B,I1B,I3B,I1B,RAM(I1A+1),RAM(I1A+2),RAM(I1A+3)), * CIPPP(I1B,I2B,I1B,I3B,I1B,RAM(I2A+1),RAM(I2A+2),RAM(I2A+3)), * CIPPP(I1B,I2B,I1B,I3B,I1B,RAM(I3A+1),RAM(I3A+2),RAM(I3A+3)), * 0.)) RETURN ENDIF C ERRMAX=(RAM(I1A+1)**2+RAM(I2A+1)**2+RAM(I3A+1)**2 * +RAM(I1B+1)**2+RAM(I2B+1)**2+RAM(I3B+1)**2 * +RAM(I1A+2)**2+RAM(I2A+2)**2+RAM(I3A+2)**2 * +RAM(I1B+2)**2+RAM(I2B+2)**2+RAM(I3B+2)**2 * +RAM(I1A+3)**2+RAM(I2A+3)**2+RAM(I3A+3)**2 * +RAM(I1B+3)**2+RAM(I2B+3)**2+RAM(I3B+3)**2)/6.0E12 C C Loop for all the selected gridpoints: DO 100, I1=K1MI,K1MA DO 90, I2=K2MI,K2MA DO 80, I3=K3MI,K3MA X1=O1+D1*FLOAT(I1) X2=O2+D2*FLOAT(I2) X3=O3+D3*FLOAT(I3) INDX=I3*N1*N2+I2*N1+I1+1 C C C Checking the position of the gridpoint with respect to C the planes bounding the ray cell: C IF (I1A.NE.I1B.AND.I2A.NE.I2B.AND.I3A.NE.I3B) THEN IF (.NOT.CILSIG( * CIPPP(I1B,I2B,I1B,I3B,I1B,RAM(I1A+1),RAM(I1A+2),RAM(I1A+3)), * CIPPP(I1B,I2B,I1B,I3B,I1B,X1,X2,X3) , 0. , 0. )) * GOTO 79 C IF (.NOT.CILSIG( * CIPPP(I1A,I2A,I1A,I3A,I1A,RAM(I1B+1),RAM(I1B+2),RAM(I1B+3)), * CIPPP(I1A,I2A,I1A,I3A,I1A,X1,X2,X3) , 0. , 0. )) * GOTO 79 C IF (CILSIG( * CIPPP(I1A,I2B,I2A,I1B,I1A,RAM(I2A+1),RAM(I2A+2),RAM(I2A+3)), * CIPPP(I1A,I2B,I2A,I1B,I1A,RAM(I3A+1),RAM(I3A+2),RAM(I3A+3)), * CIPPP(I1A,I2B,I2A,I1B,I1A,RAM(I3B+1),RAM(I3B+2),RAM(I3B+3)), * 0.)) THEN IF (.NOT.CILSIG( * CIPPP(I1A,I2B,I2A,I1B,I1A,RAM(I2A+1),RAM(I2A+2),RAM(I2A+3)), * CIPPP(I1A,I2B,I2A,I1B,I1A,X1,X2,X3),0.,0.)) GOTO 79 ELSEIF (CILSIG( * CIPPP(I1A,I2B,I2A,I1B,I2A,RAM(I1A+1),RAM(I1A+2),RAM(I1A+3)), * CIPPP(I1A,I2B,I2A,I1B,I2A,RAM(I3A+1),RAM(I3A+2),RAM(I3A+3)), * CIPPP(I1A,I2B,I2A,I1B,I2A,RAM(I3B+1),RAM(I3B+2),RAM(I3B+3)), * 0.)) THEN IF (.NOT.CILSIG( * CIPPP(I1A,I2B,I2A,I1B,I2A,RAM(I1A+1),RAM(I1A+2),RAM(I1A+3)), * CIPPP(I1A,I2B,I2A,I1B,I2A,X1,X2,X3),0.,0.)) GOTO 79 ENDIF C IF (CILSIG( * CIPPP(I2A,I3B,I3A,I2B,I2A,RAM(I3A+1),RAM(I3A+2),RAM(I3A+3)), * CIPPP(I2A,I3B,I3A,I2B,I2A,RAM(I1A+1),RAM(I1A+2),RAM(I1A+3)), * CIPPP(I2A,I3B,I3A,I2B,I2A,RAM(I1B+1),RAM(I1B+2),RAM(I1B+3)), * 0.)) THEN IF (.NOT.CILSIG( * CIPPP(I2A,I3B,I3A,I2B,I2A,RAM(I3A+1),RAM(I3A+2),RAM(I3A+3)), * CIPPP(I2A,I3B,I3A,I2B,I2A,X1,X2,X3),0.,0.)) GOTO 79 ELSEIF (CILSIG( * CIPPP(I2A,I3B,I3A,I2B,I3A,RAM(I2A+1),RAM(I2A+2),RAM(I2A+3)), * CIPPP(I2A,I3B,I3A,I2B,I3A,RAM(I1A+1),RAM(I1A+2),RAM(I1A+3)), * CIPPP(I2A,I3B,I3A,I2B,I3A,RAM(I1B+1),RAM(I1B+2),RAM(I1B+3)), * 0.)) THEN IF (.NOT.CILSIG( * CIPPP(I2A,I3B,I3A,I2B,I3A,RAM(I2A+1),RAM(I2A+2),RAM(I2A+3)), * CIPPP(I2A,I3B,I3A,I2B,I3A,X1,X2,X3),0.,0.)) GOTO 79 ENDIF C IF (CILSIG( * CIPPP(I3A,I1B,I1A,I3B,I1A,RAM(I3A+1),RAM(I3A+2),RAM(I3A+3)), * CIPPP(I3A,I1B,I1A,I3B,I1A,RAM(I2A+1),RAM(I2A+2),RAM(I2A+3)), * CIPPP(I3A,I1B,I1A,I3B,I1A,RAM(I2B+1),RAM(I2B+2),RAM(I2B+3)), * 0.)) THEN IF (.NOT.CILSIG( * CIPPP(I3A,I1B,I1A,I3B,I1A,RAM(I3A+1),RAM(I3A+2),RAM(I3A+3)), * CIPPP(I3A,I1B,I1A,I3B,I1A,X1,X2,X3),0.,0.)) GOTO 79 ELSEIF (CILSIG( * CIPPP(I3A,I1B,I1A,I3B,I3A,RAM(I1A+1),RAM(I1A+2),RAM(I1A+3)), * CIPPP(I3A,I1B,I1A,I3B,I3A,RAM(I2A+1),RAM(I2A+2),RAM(I2A+3)), * CIPPP(I3A,I1B,I1A,I3B,I3A,RAM(I2B+1),RAM(I2B+2),RAM(I2B+3)), * 0.)) THEN IF (.NOT.CILSIG( * CIPPP(I3A,I1B,I1A,I3B,I3A,RAM(I1A+1),RAM(I1A+2),RAM(I1A+3)), * CIPPP(I3A,I1B,I1A,I3B,I3A,X1,X2,X3),0.,0.)) GOTO 79 ENDIF ENDIF C C C Computation of the coefficients AAA,BBB,CCC,DDD C of the cubic equation: C XX(1,1)=RAM(I1A+1)-RAM(I1B+1) XX(1,2)=RAM(I1A+2)-RAM(I1B+2) XX(1,3)=RAM(I1A+3)-RAM(I1B+3) XX(2,1)=RAM(I2A+1)-RAM(I2B+1) XX(2,2)=RAM(I2A+2)-RAM(I2B+2) XX(2,3)=RAM(I2A+3)-RAM(I2B+3) XX(3,1)=RAM(I3A+1)-RAM(I3B+1) XX(3,2)=RAM(I3A+2)-RAM(I3B+2) XX(3,3)=RAM(I3A+3)-RAM(I3B+3) C ZZ(1,1)=RAM(I1B+1)-X1 ZZ(1,2)=RAM(I1B+2)-X2 ZZ(1,3)=RAM(I1B+3)-X3 ZZ(2,1)=RAM(I2B+1)-X1 ZZ(2,2)=RAM(I2B+2)-X2 ZZ(2,3)=RAM(I2B+3)-X3 ZZ(3,1)=RAM(I3B+1)-X1 ZZ(3,2)=RAM(I3B+2)-X2 ZZ(3,3)=RAM(I3B+3)-X3 C AAA=CIDET3(XX) BBB=CISUBD(ZZ,XX) CCC=CISUBD(XX,ZZ) DDD=CIDET3(ZZ) C C C Solving the equation according to their coefficients: IF (AAA.EQ.0.) THEN IF (BBB.EQ.0.) THEN C Linear equation: IF (CCC.EQ.0.) THEN C MTT-06 PAUSE * 'Error MTT-06: wrong coefficients of equation in CIINTP' STOP C This error should not appear. C Please contact the author. ENDIF R1=-DDD/CCC IF (R1.GE.0..AND.R1.LE.1.) THEN IR=1 ROOT(1)=R1 ELSE IR=0 ENDIF ELSE C Kvadratic equation: CALL CIKVAR(BBB,CCC,DDD,IR,R1,R2) IF (IR.GE.1) THEN IF ((R1.LT.0.).OR.(R1.GT.1.)) THEN C MTT-11 PAUSE 'Error MTT-11: First root of quadratic eq. wrong' STOP C This error should not appear. C Please contact the author. ENDIF ROOT(1)=R1 ENDIF IF (IR.GE.2) THEN IF ((R2.LT.0.).OR.(R2.GT.1.)) THEN C MTT-12 PAUSE 'Error MTT-12: Second root of quadratic eq. wrong' STOP C This error should not appear. C Please contact the author. ENDIF ROOT(2)=R2 ENDIF ENDIF ELSE C Cubic equation: C Computation of the coefficients EEE,FFF,CCC of the derivative C of the cubic equation: EEE=3.*AAA FFF=2.*BBB C C Solving the equation for the derivative: CALL CIKVAR(EEE,FFF,CCC,IR,R1,R2) C Deciding, whether the cubic equation is to be solved: RR(1)=0. IF (IR.GE.1) THEN RR(2)=R1 NRR=3 ELSE NRR=2 ENDIF IF (IR.GE.2) THEN RR(NRR)=R2 NRR=NRR+1 ENDIF RR(NRR)=1. DO 10, I4=1,NRR CRR(I4)=CICUB(AAA,BBB,CCC,DDD,RR(I4)) 10 CONTINUE IR=0 I4=1 20 CONTINUE IF ((CRR(I4)*CRR(I4+1)).LE.0.) THEN IR=IR+1 I4=I4+1 ELSE DO 30, I5=I4,NRR-1 CRR(I5)=CRR(I5+1) RR(I5)=RR(I5+1) 30 CONTINUE NRR=NRR-1 ENDIF IF (I4.LT.NRR) GOTO 20 IF (IR.EQ.0) GOTO 79 C C C Solving the cubic equation: CALL CICUBR(AAA,BBB,CCC,DDD,RR,IR,ROOT) ENDIF C DO 70, I4=1,IR WB=ROOT(I4) IF (WB.EQ.0.) GOTO 79 WC=1.-WB IF ((WB.LT.0.).OR.(WC.LT.0.)) THEN C MTT-13 PAUSE 'Error MTT-13: Root outside the cell' STOP C This error should not appear. C Please contact the author. ENDIF A11=WB*RAM(I1A+1) + WC*RAM(I1B+1) A21=WB*RAM(I1A+2) + WC*RAM(I1B+2) A31=WB*RAM(I1A+3) + WC*RAM(I1B+3) A12=WB*RAM(I2A+1) + WC*RAM(I2B+1) A22=WB*RAM(I2A+2) + WC*RAM(I2B+2) A32=WB*RAM(I2A+3) + WC*RAM(I2B+3) A13=WB*RAM(I3A+1) + WC*RAM(I3B+1) A23=WB*RAM(I3A+2) + WC*RAM(I3B+2) A33=WB*RAM(I3A+3) + WC*RAM(I3B+3) A11=A11-A13 A21=A21-A23 A31=A31-A33 A12=A12-A13 A22=A22-A23 A32=A32-A33 A13=X1-A13 A23=X2-A23 A33=X3-A33 A1=SQRT(A11*A11+A21*A21+A31*A31) A2=SQRT(A12*A12+A22*A22+A32*A32) B11=A11*A2+A12*A1 B21=A21*A2+A22*A1 B31=A31*A2+A32*A1 B12=A11*A2-A12*A1 B22=A21*A2-A22*A1 B32=A31*A2-A32*A1 UU11=B11*A11+B21*A21+B31*A31 UU21=B12*A11+B22*A21+B32*A31 UU12=B11*A12+B21*A22+B31*A32 UU22=B12*A12+B22*A22+B32*A32 VV1 =B11*A13+B21*A23+B31*A33 VV2 =B12*A13+B22*A23+B32*A33 DETUU=UU11*UU22-UU12*UU21 C Determinant eq zero in case of infinitely thin ray cell, C in such a case gridpoint cannot lie inside the cell. IF (DETUU.NE.0.) THEN W1=(UU22*VV1-UU12*VV2)/DETUU W2=(UU11*VV2-UU21*VV1)/DETUU W3=1.-W1-W2 IF ((W1.GE.0.).AND.(W1.LE.1.).AND. * (W2.GE.0.).AND.(W2.LE.1.).AND. * (W3.GE.0.).AND.(W3.LE.1.)) THEN C The gridpoint is situated inside the ray cell. C C Checking accuracy of interpolation coefficients: Y1=WB*(W1*RAM(I1A+1)+W2*RAM(I2A+1)+W3*RAM(I3A+1))+ * WC*(W1*RAM(I1B+1)+W2*RAM(I2B+1)+W3*RAM(I3B+1)) Y2=WB*(W1*RAM(I1A+2)+W2*RAM(I2A+2)+W3*RAM(I3A+2))+ * WC*(W1*RAM(I1B+2)+W2*RAM(I2B+2)+W3*RAM(I3B+2)) Y3=WB*(W1*RAM(I1A+3)+W2*RAM(I2A+3)+W3*RAM(I3A+3))+ * WC*(W1*RAM(I1B+3)+W2*RAM(I2B+3)+W3*RAM(I3B+3)) IF (((Y1-X1)**2 + (Y2-X2)**2 + (Y3-X3)**2).GT.ERRMAX) THEN C MTT-23 PAUSE 'Error MTT-23: Inexact interpolation coefficients' STOP C Interpolation coefficients are not exact enough to C provide exact coordinates of gridpoints. Thus the error C of the same order will appear in interpolation of C travel times and other quantities. C In principle this error should not appear. C On the other hand, if you wish to continue in C interpolation even with inaccurate interpolation C coefficients, put the letter 'C' to the beginning of the C previous line with 'STOP'. ENDIF C C Interpolation of output quantities: INDX=MRAM-INDX+1 40 CONTINUE IF (IRAM(INDX).NE.0) THEN INDX=IRAM(INDX) GOTO 40 ENDIF NRAMT=NRAMT-(NOUT+1) IF (NRAMT.LE.NRAMP) CALL CIEROR(1) IRAM(INDX)=NRAMT IRAM(NRAMT)=0 MRAMP=NRAMT-1 DO 50, I5=1,NOUT IF (OUT(I5)(1:2).EQ.'tt') THEN C Travel time: TT=WB*W1*RAM(I1A+(5+I5))+WB*W2*RAM(I2A+(5+I5))+ * WB*W3*RAM(I3A+(5+I5))+ * WC*W1*RAM(I1B+(5+I5))+WC*W2*RAM(I2B+(5+I5))+ * WC*W3*RAM(I3B+(5+I5)) RAM(NRAMT+I5)=TT ELSEIF (OUT(I5)(1:2).EQ.'hi') THEN C Ray history: IRAM(NRAMT+I5)=IRAM(I1A+(5+I5)) ENDIF 50 CONTINUE GOTO 79 ENDIF ENDIF 70 CONTINUE 79 CONTINUE 80 CONTINUE 90 CONTINUE 100 CONTINUE END C C======================================================================= C REAL FUNCTION CIPPP(IU1,IU2,IV1,IV2,IA,X1,X2,X3) C C---------------------------------------------------------------------- INTEGER IU1,IU2,IV1,IV2,IA REAL X1,X2,X3 C Subroutine designed to indicate the position of C a point X: (X1,X2,X3) with respect to C a plane defined by two vectors u: IU1 --> IU2, v: IV1 --> IV2, C and by a point IA. C The subroutine computes the scalar product a.w of vector C a: IA --> X with vector w: w = u x v being normal to the plane. C The sign of this product then indicates, on which side of the C plane is the point X located. C All the computation is performed in 3-D Cartesian space. C C Input: C IU1,IU2 .. Two points defining vector u. C IV1,IV2 .. Two points defining vector v. C IA ... Point of the plane. C Each integer contains the address just before the parameters C of the corresponding point: C the first parameter of the point IA: RAM(IA+1) C X1,X2,X3 ... Coordinates of the examined point. C Output: C CIPPP ... The value of the scalar product. C C Coded by Petr Bulant C INCLUDE 'mtt.inc' C C...................................................................... REAL U1,U2,U3,V1,V2,V3,W1,W2,W3,A1,A2,A3 C----------------------------------------------------------------------- U1=RAM(IU2+1)-RAM(IU1+1) U2=RAM(IU2+2)-RAM(IU1+2) U3=RAM(IU2+3)-RAM(IU1+3) V1=RAM(IV2+1)-RAM(IV1+1) V2=RAM(IV2+2)-RAM(IV1+2) V3=RAM(IV2+3)-RAM(IV1+3) W1=U2*V3-U3*V2 W2=U3*V1-U1*V3 W3=U1*V2-U2*V1 A1=X1-RAM(IA+1) A2=X2-RAM(IA+2) A3=X3-RAM(IA+3) CIPPP=W1*A1+W2*A2+W3*A3 END C C======================================================================= C REAL FUNCTION CICUB(AAA,BBB,CCC,DDD,XX) C C---------------------------------------------------------------------- REAL AAA,BBB,CCC,DDD,XX C Subroutine designed to calculate the value of cubic function C given by coefficients AAA,BBB,CCC,DDD in point XX. C...................................................................... CICUB=AAA*XX*XX*XX+BBB*XX*XX+CCC*XX+DDD END C C======================================================================= C REAL FUNCTION CIKVA(AAA,BBB,CCC,XX) C C---------------------------------------------------------------------- REAL AAA,BBB,CCC,XX C Subroutine designed to calculate the value of kvadratic function C given by coefficients AAA,BBB,CCC in point XX. C...................................................................... CIKVA=AAA*XX*XX+BBB*XX+CCC END C C======================================================================= C REAL FUNCTION CIDET3(A) C C---------------------------------------------------------------------- REAL A(3,3) C Subroutine designed to calculate the determinant of 3x3 matrix AA. C Input: C A ... 3x3 real matrix. C Output: C CIDET3 . Determinant. C C Coded by Petr Bulant C...................................................................... CIDET3=A(1,1)*(A(2,2)*A(3,3)-A(3,2)*A(2,3))+ * A(2,1)*(A(3,2)*A(1,3)-A(1,2)*A(3,3))+ * A(3,1)*(A(1,2)*A(2,3)-A(2,2)*A(1,3)) END C C======================================================================= C REAL FUNCTION CISUBD(A,B) C C---------------------------------------------------------------------- REAL A(3,3),B(3,3) C Subroutine designed to calculate the sum C C SUM ( e * A * B * B ) C ijk ijk 1i 2j 3k C C where i,j,k=1..3 and e is a Levi-Civita symbol. C ijk C Input: C A ... 3x3 real matrix. C B ... 3x3 real matrix. C Output: C CISUBD . Value of the sum mentioned above. C C Coded by Petr Bulant C...................................................................... CISUBD=A(1,1)*(B(2,2)*B(3,3)-B(3,2)*B(2,3))+ * A(2,1)*(B(3,2)*B(1,3)-B(1,2)*B(3,3))+ * A(3,1)*(B(1,2)*B(2,3)-B(2,2)*B(1,3))+ * A(1,2)*(B(2,3)*B(3,1)-B(3,3)*B(2,1))+ * A(2,2)*(B(3,3)*B(1,1)-B(1,3)*B(3,1))+ * A(3,2)*(B(1,3)*B(2,1)-B(2,3)*B(1,1))+ * A(1,3)*(B(2,1)*B(3,2)-B(3,1)*B(2,2))+ * A(2,3)*(B(3,1)*B(1,2)-B(1,1)*B(3,2))+ * A(3,3)*(B(1,1)*B(2,2)-B(2,1)*B(1,2)) END C C C======================================================================= C SUBROUTINE CIKVAR(AAA,BBB,CCC,IR,R1,R2) C C---------------------------------------------------------------------- C Subroutine to solve the kvadratic equation on the interval <0,1>. INTEGER IR REAL AAA,BBB,CCC,R1,R2 C Input: C AAA,BBB,CCC ... Coefficients of the equation. C Output: C IR ... Number of real roots of the equation from <0,1>. C R1,R2. Real roots of the equation from <0,1>, R1. INTEGER IR REAL AA0,BB0,CC0,DD0,RR(4),ROOT(3) C Input: C AA0,BB0,CC0,DD0 ... Coefficients of the equation. C IR ... Expected number of real roots from the interval <0,1>. C RR ... Array of points, which divide the interval <0,1> into C subintervals, each subinterval containing just one C root of the cubic equation. C RR(1) ... 0. C RR(NRR) ... 1. C for IR=1 NRR=2 C for IR=2 NRR=3, RR(2) is a point in which C the cubic equation has zero derivative C for IR=3 NRR=4, RR(2) and RR(3) are points in which C the cubic equation has zero derivative. C Output: C ROOT ... Real roots of the equation from <0,1> or undefined. C C Subroutines and external functions required: EXTERNAL CICUB,CIKVA REAL CICUB,CIKVA C C Coded by Petr Bulant C C....................................................................... REAL ZERO PARAMETER (ZERO=0.000001) REAL TRIAA0,DVEBB0,POM1,POM2,POM3,DER,RA,RB INTEGER I1,I2 C----------------------------------------------------------------------- TRIAA0=3.*AA0 DVEBB0=2.*BB0 DO 20, I2=1,IR RA=RR(I2) RB=RR(I2+1) I1=0 10 CONTINUE I1=I1+1 IF (I1.GT.100) THEN IF (ABS(RA-RB).LE.ZERO) RETURN C MTT-08 PAUSE 'Error MTT-08: infinite loop in CICUBR' c STOP C This error should not appear. C Please contact the author. ENDIF C POM1=CICUB(AA0,BB0,CC0,DD0,RB) POM2=CICUB(AA0,BB0,CC0,DD0,RA) DER=CIKVA(TRIAA0,DVEBB0,CC0,RA) IF (DER.NE.0.) THEN ROOT(I2)=RA-POM2/DER IF((ROOT(I2).LE.RA).OR.(ROOT(I2).GE.RB))ROOT(I2)=(RA+RB)/2. ELSE ROOT(I2)=(RA+RB)/2. ENDIF POM3=CICUB(AA0,BB0,CC0,DD0,ROOT(I2)) IF (POM1*POM3.LT.0.) THEN RA=ROOT(I2) ELSE RB=ROOT(I2) ENDIF C POM1=CICUB(AA0,BB0,CC0,DD0,RB) POM2=CICUB(AA0,BB0,CC0,DD0,RA) POM3=POM1-POM2 IF (POM3.EQ.0.) THEN C MTT-09 PAUSE 'Error MTT-09: division by zero in CICUBR.' STOP C This error should not appear. C Please contact the author. ENDIF ROOT(I2)=RA-(RB-RA)*POM2/POM3 POM3=CICUB(AA0,BB0,CC0,DD0,ROOT(I2)) IF (POM1*POM3.LT.0.) THEN RA=ROOT(I2) POM2=POM3 ELSE RB=ROOT(I2) POM1=POM3 ENDIF IF (ABS(RA-RB).GT.ZERO) GOTO 10 IF (RA.NE.RB) THEN IF (POM1.EQ.0.) THEN ROOT(I2)=RB ELSE IF (POM2.EQ.0.) THEN ROOT(I2)=RA ELSE IF (POM1-POM2.NE.0.) THEN C Linear interpolation of the root: RA=RA-ROOT(I2) RB=RB-ROOT(I2) POM3=(RA*POM1-RB*POM2)/(POM1-POM2) ROOT(I2)=ROOT(I2)+POM3 ELSE C MTT-22 PAUSE 'Error MTT-22: Equal functional values in CICUBR.' STOP C This error should not appear. C Please contact the author. STOP END IF END IF 20 CONTINUE C END C C======================================================================= C LOGICAL FUNCTION CILSIG(R1,R2,R3,R4) C C---------------------------------------------------------------------- REAL R1,R2,R3,R4 C Subroutine to compare signs of real quantities. All the quantities C with absolute value less then parameter ZERO are considered to C be equal 0. C C Input: R1,R2,R3,R4 ... Real quantities. C Output: CILSIG ... .TRUE. for all the quantities nonnegative or C all th quantities nonpositive. C .FALSE. in other cases. C....................................................................... REAL ZERO PARAMETER (ZERO=.000001) C ZERO ... Constant used to decide whether C the real quantiti .EQ. zero. C----------------------------------------------------------------------- IF (ABS(R1).LT.ZERO) R1=0. IF (ABS(R2).LT.ZERO) R2=0. IF (ABS(R3).LT.ZERO) R3=0. IF (ABS(R4).LT.ZERO) R4=0. IF (((R1.LE.0.).AND.(R2.LE.0.).AND.(R3.LE.0.).AND.(R4.LE.0.)).OR. * ((R1.GE.0.).AND.(R2.GE.0.).AND.(R3.GE.0.).AND.(R4.GE.0.))) * THEN CILSIG=.TRUE. ELSE CILSIG=.FALSE. ENDIF END C C C======================================================================= C SUBROUTINE CIGRID(LU4,FILE4) C C---------------------------------------------------------------------- C Subroutine to read in the file with parameters of an interpolation C grid. The file is assumed to be in form of a SEP-like parameter or C header file. The read quantities are then written into common block C CIGRD. INTEGER LU4 CHARACTER*80 FILE4 C Input: C LU4 ... Number of logical unit corresponding to the file with C the input quantities for interpolation grid. C FILE4 . Name of the file. C No output. C C Coded by Petr Bulant C C ........................... C Common block /RAM/ to store the information about triangles and rays: INCLUDE 'mtt.inc' C ........................... C Common block /CIGRD/ to store the parameters of the C interpolation grid. REAL O1,O2,O3,D1,D2,D3 INTEGER N1,N2,N3 COMMON/CIGRD/O1,O2,O3,D1,D2,D3,N1,N2,N3 SAVE /CIGRD/ C O1,O2,O3 ... Coordinates of the origin of the grid. C D1,D2,D3 ... Grid steps. C N1,N2,N3 ... Numbers of gridpoints. C....................................................................... INTEGER I1 C----------------------------------------------------------------------- C Reading all the data from file FILE4 to the memory C (SEP parameter file form): CALL RSEP1(LU4,FILE4) C C Recalling the data specifying grid dimensions C (arguments: Name of value in input data, Variable, Default): CALL RSEP3R('O1',O1,0.) CALL RSEP3R('O2',O2,0.) CALL RSEP3R('O3',O3,0.) CALL RSEP3R('D1',D1,1.) CALL RSEP3R('D2',D2,1.) CALL RSEP3R('D3',D3,1.) CALL RSEP3I('N1',N1,1) CALL RSEP3I('N2',N2,1) CALL RSEP3I('N3',N3,1) C IF ((D1.LT.0.).OR.(D2.LT.0.).OR.(D3.LT.0.)) GOTO 10 IF ((N1.LE.0).OR.(N2.LE.0).OR.(N3.LE.0)) GOTO 10 IF (D1.EQ.0.) THEN IF (N1.EQ.1) THEN D1=1. ELSE GOTO 10 ENDIF ENDIF IF (D2.EQ.0.) THEN IF (N2.EQ.1) THEN D2=1. ELSE GOTO 10 ENDIF ENDIF IF (D3.EQ.0.) THEN IF (N3.EQ.1) THEN D3=1. ELSE GOTO 10 ENDIF ENDIF MRAMP=MRAM-N1*N2*N3 IF (MRAMP.LE.1) CALL CIEROR(1) NRAMT=MRAMP+1 DO 5, I1=NRAMT,MRAM IRAM(I1)=0 5 CONTINUE RETURN C 10 CONTINUE C MTT-10 PAUSE 'Error MTT-10: This specification of the interpolation * grid may cause problems. Please, specify * D1,D2,D3 and N1,N2,N3 greater than 0. * Di may equal 0 in case that corresponding Ni equals 1. ' STOP C END C C C======================================================================= C SUBROUTINE CIEROR(IERR) C C----------------------------------------------------------------------- INTEGER IERR C C Subroutine designed to print error messages of different C CI* subroutines using command 'PAUSE'. C C Input: C IERR ... Index of the error. C No output. C Coded by Petr Bulant C----------------------------------------------------------------------- C IF (IERR.EQ.001) THEN C MTT-01 PAUSE 'Error MTT-01: array (I)RAM too small.' STOP C This error may be caused by too small dimension of array C RAM. Try to enlarge the parameter MRAM in common block C RAM. ENDIF END C C======================================================================= C INCLUDE 'ap.for' C ap.for INCLUDE 'forms.for' C forms.for INCLUDE 'length.for' C length.for INCLUDE 'sep.for' C sep.for C C======================================================================= C