```
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.
```