Two-point network shortest-path ray tracing within Fresnel volumes. Description of the algorithm employing programs NET and NETIND.

Given: The model, the source point, and the receiver point. Algorithm: (1) Determination of the Fresnel volume: (1.1) Create source point file 'SRC' containing the coordinates of the source point, and receiver point file 'REC' containing the coordinates of the receiver point. The file format is described within the source code 'net.for'. Choose the rectangular grid for the network ray tracing: N1*N2*N3 big bricks, each containing one grid point. We assume here that this original grid is not indexed. (1.2) Generate the grid velocities (file 'VEL') and, possibly, indices of geological blocks (file 'ICB'). If using the model specification software package 'MODEL', this task is performed by the program 'grid.for'. (1.3) Calculate travel-time field 'TT1' from the source, and the two-point travel time with its error. This step is performed by the program 'net.for'. (1.4) Calculate travel-time field 'TT2' from the receiver. This step is performed by the program 'net.for'. (1.5) Run 'netind.for' program with the travel-time fields 'TT1' and 'TT2', and with the calculated two-point travel time plus its estimated error (as the maximum sum of travel times limiting the Fresnel volume). The Fresnel volume is determined and the corresponding index file 'IND' is created. (2) Two-point ray tracing within the Fresnel volume: (2.1) Edit the grid: Divide each of N1*N2*N3 big bricks into L1*L2*L3 small bricks, each small brick having the gridpoint at its centre. If big bricks cannot be divided (e.g. because being short of the computer memory), the algorithm has to be terminated here, without reaching the required accuracy. This step is usually completed together with step (1.5) or (2.5) by the 'netind.for' program. (2.2) Generate the grid velocities (file 'VEL') and, possibly, indices of geological blocks (file 'ICB'). The same as (1.2) except for the index file 'IND', describing the Fresnel volume, specified. The new files 'VEL' and 'ICB' correspond to the index file 'IND' and to the new subdivision of big bricks into small bricks. (2.3) Calculate new, indexed, travel-time field 'TT1' from the source, and the two-point travel time with its error. If the two-point travel time is sufficiently accurate, the procedure may be terminated here and the file 'TT1' is not required. This step is performed by the program 'net.for'. (2.4) Calculate new, indexed, travel-time field 'TT2' from the receiver. This step is performed by the program 'net.for'. (2.5) Run 'netind.for' program with the travel-time fields 'TT1' and 'TT2', and with the calculated two-point travel time plus its estimated error (as the maximum sum of travel times limiting the Fresnel volume). Then the new Fresnel volume is determined and the new index file 'OUT' is created. Whereas the old index file 'IND' corresponds to N1*N2*N3 big bricks, the new index file 'OUT' corresponds to (N1*L1)*(N2*L2)*(N3*L3) big bricks. In this way, the new index file upgrades the old small bricks to big bricks. (3) Changing to the new, finer Fresnel volume: (3.1) Edit the grid: Upgrade old small bricks to big bricks, i.e. write (N1*L1),(N2*L2),(N3*L3) in place of N1,N2,N3, while L1,L2,L3 become undefined. This step is usually completed together with step (2.5) by the 'netind.for' program. (3.2) Replace the old index file 'IND' by the new index file 'OUT'. (3.3) Proceed to (2.1). Strategy of array dimensioning for network ray tracing within Fresnel volumes: In 'net.for', neglecting the template forward stars, the most of the memory is used by the following 7 arrays: IND(N123)... Indexing of Fresnel volumes, P(L1234)... Slownesses (always used), TT(L1234)... First-arrival times (always used), IPOSQ(L1234)... Queue (always used), IPRED(L1234)... Predecessors (used in Fresnel volumes), NFS(L1234)... Optimum sizes of forward stars (often used), ICB(L1234)... Indices of geological blocks (often used). Here the dimensions denote the number of storage locations required. If the index array IND(N123) is used, the minimum dimension N123 is the number of big bricks in the whole model volume. The minimum dimension L1234 of the other 4 to 6 arrays is the number of small bricks in the Fresnel volume. Let us denote: NDIM=2,3 ... For 2-D or 3-D calculation, respectively, NARR=4,5,6 ... Number of arrays dimensioned by L1234, MRAM... Total memory available for the above 7 arrays. The error of arrival time is proportional to the grid step: error=constant*step. The grid step in the last but one iteration is dependent on the total number of small bricks in the model volume in that iteration, that is the number N123 of all big bricks in the last iteration: step=constant*N123**(-1/NDIM). The ratio of the Fresnel volume to the model volume depends on the error in the last but one iteration: ratio=constant*error**((NDIM-1)/2). Finally: ratio=constant*N123**(-(NDIM-1)/(2*NDIM)). Then the number of small bricks within the Fresnel volume is: L1234=L1*L2*L3*N123*ratio as an objective, we choose to maximize the total number Y of small bricks within the whole model volume during the last iteration: Y=L1*L2*L3*N123 =L1234/RATIO =constant*L1234*N123**((NDIM-1)/(2*NDIM)). The memory available for all the above arrays is: MRAM=N123+NARR*L1234. The objective function Y=(constant/NARR)*(MRAM-N123)*N123**((NDIM-1)/(2*NDIM)) is extremal for N123=MRAM*(NDIM-1)/(3*NDIM-1), i.e. In 2-D: N123=MRAM/5, L1234=MRAM*4/(5*NARR), In 3-D: N123=MRAM/4, L1234=MRAM*3/(4*NARR). There is no considerable difference between this optimum choice and the smaller value of N123=MRAM/NARR. For example, if NARR=6, the difference of the grid step is 1/4 per cent in 2-D, and 1 per cent in 3-D. That is why in most cases the following two-iteration strategy will be sufficient: 1-st iteration: N123=0, L1234=MRAM*(NDIM-1)/(3*NDIM-1) 2-nd iteration: N123=MRAM*(NDIM-1)/(3*NDIM-1), L1234=(MRAM-N123)/NARR Input data for 'net.for' (first iteration): N1*N2*N3=MRAM*(NDIM-1)/(3*NDIM-1) Input data for 'netind.for': L1MAX=0, L2MAX=0, L3MAX=0 (defaults) In 3-D, the big bricks may be expected to be divided into 2*2*2 to 6*6*6 small bricks, depending on the complexity of the model. A finer division is unlike in inhomogeneous 3-D models. Although the number of network nodes in the first iteration is slightly greater than in the second iteration, the first iteration should be several times (2 to 6 times in 3-D) faster than the second iteration because of smaller optimized forward stars. On the other hand, if desirable to save the computation time, the following three-iteration strategy 1-st iteration: N123=0, L1234=MRAM*(NDIM-1)/(3*NDIM-1)/(2**NDIM), 2-nd iteration: N123=MRAM*(NDIM-1)/(3*NDIM-1)/(2**NDIM), L1*L2*L3=2**NDIM, 3-rd iteration: N123=MRAM*(NDIM-1)/(3*NDIM-1), L1234=(MRAM-N123)/NARR Input data for 'net.for' (first iteration): N1*N2*N3=MRAM*(NDIM-1)/(3*NDIM-1)/(2**NDIM) Input data for 'netind.for': L1MAX=2, L2MAX=2, L3MAX=2 (first run) L1MAX=0, L2MAX=0, L3MAX=0 (second run) may be little bit faster since nearly 2**(1.5*NDIM-0.5) times reducing the computation time of the first iteration. Let us note that the above strategies were suggested especially for 3-D, where the shortage of computer memory limiting the accuracy is a dominant feature. If, in 2-D, the computer memory does not limit the accuracy required, other strategies decreasing the computation time should be applied.