====================================================================== U. S. Department of the Interior, U. S. Geological Survey xwell_radar: A Program for Calculating Crosswell Radar Traces for Layered Media Karl J. Ellefsen Open File Report 00-039 URL: http://greenwood.cr.usgs.gov/pub/open-file-reports/ofr-00-0039/ ====================================================================== 1. Overview Program 'xwell_radar' calculates radar traces that simulate crosswell radar data collected in horizontally layered media. The program is for computers with the Unix operating system and is written in the (ANSI) C and the Fortran (77) programming languages. This file 'README.TXT' describes (a) The contents and organization of the tape archive file (section 2). (b) Program 'xwell_radar', its installation, and its execution (section 3). (c) A data set that may be used to test whether the program is properly installed (section 4). If you find any errors or if you have any suggestions for improving the program, please contact: Karl J. Ellefsen U. S. Geological Survey MS 964, Box 25046 Denver, CO 80225 telephone: 303-236-7032 email: ellefsen@usgs.gov ====================================================================== 2. Contents and Organization of the Tape Archive File The name of the tape archive file (tar file) is 'xwell_radar.TAR', and its contents can be viewed with the Unix command: 'tar tvf xwell_radar.TAR'. In the top level of the tar file is 'README.TXT' --- a copy of this text file. The tar file also contains three directories. In directory 'xwell_radar' are 7 files, which are used to calculate the radar traces; these files are described in section 3.1. In directory 'bessel' are 32 files, which are used to calculate Bessel functions. In directory 'test_data' are 7 files, which are used to test the installation of the program; these files are described in section 4. ====================================================================== 3. Program 'xwell_radar' 3.1 Program Description The mathematical model, on which program 'xwell_radar' is based, is composed of horizontal layers. The top and the bottom layers are half-spaces. Each layer is homogeneous and isotropic; the constitutive relations are linear; and the electromagnetic properties are frequency-dependent. Locations are specified in circular cylindrical coordinates; the origin is at the bottom of the layer with the antennas. (Both antennas must be in the same layer and that layer cannot be the top layer or the bottom layer.) The transmitting antenna is on the z axis and is simulated by a vertical electric dipole with infinitesimal length. To simulate the receiving antenna, the vertical component of the electric field intensity is calculated. The simplest possible model consists of three layers: Layer 3 (half-space) --------------------------------------------------------- ^ |z axis | Layer 2 | # Receiving | Antenna | | Transmitting * Antenna | | -------------------------|==============================> Layer 1 (half-space) r axis Layers 1 and 3 are half-spaces; layer 2 contains the antennas. A more complex model might have, for example, 6 layers: Layer 6 (half-space) --------------------------------------------------------- Layer 5 --------------------------------------------------------- Layer 4 --------------------------------------------------------- ^ |z axis Layer 3 | Transmitting * # Receiving Antenna | Antenna -------------------------|==============================> r axis Layer 2 --------------------------------------------------------- Layer 1 (half-space) A mathematical solution for the electromagnetic fields in a three- layer model was derived by Wait (1995, p. 137), and this solution was modified to to handle models with many layers (Wait, 1995, p. 181- 185). The modified solution is used for program 'xwell_radar'. The solution accounts for both the near and the far fields, for the reflected and the refracted waves from all layers, and for the direct wave that propagates (directly) from the transmitting antenna to the receiving antenna. The solution is in the frequency domain; consequently, accounting for the frequency-dependence of the electromagnetic properties is easy. To make Wait's solution suitable for numerical computation, some minor modifications are necessary, and these are described in Ellefsen (1999). The equations used for numerical computation contain two integrals: One integral pertains to the radial wavenumber, and this integration is performed with the discrete wavenumber method (Bouchon, 1981). The other integral pertains to frequency, and this integration is performed with a fast Fourier transform. This last integral transforms Wait's solution into the time domain, yielding a radar trace. In the transmitting antenna, the current I is specified by a Ricker wavelet: I(t) = 0.25 / sqrt( pi ) * [ 1. - 0.5 * ( 2. * pi * f0 )^2 * ( t - t0 )^2 ] * exp[ - 0.25 * ( 2. * pi * f0 )^2 * ( t - t0 )^2 ] where t is time. ('*' means multiply; 'pi' means 3.14159; 'sqrt' means 'square root'; and '^2' means 'to the power 2'.) The center frequency is f0, and the time shift is t0. The time shift must be applied to make the wavelet causal; a very good value for the time shift can be calculated from this empirical formula: 0.85 / f0. For example, if the center frequency is 100 MHz, the time shift should be about 8.5 ns. The frequency-dependence in the dielectric permittivity is specified by the formula developed by Cole and Cole (1941): er(w) = efr + ( esr - efr ) / [ 1. + ( i * w * e_tau )^(1.-e_alpha) ] where er(w) is the relative dielectric permittivity at the angular frequency w, efr is the relative dielectric permittivity at infinite frequency, esr is the relative dielectric permittivity at zero frequency, i is the square root of -1, e_tau is the generalized relaxation time for the permittivity, and e_alpha accounts for multiple relaxation phenomena in the permittivity. The frequency- dependence in the magnetic permeability is specified with a similar formula: mr(w) = mfr + ( msr - mfr ) / [ 1. + ( i * w * m_tau )^(1.-m_alpha) ] where mr(w) is the relative magnetic permeability at the angular frequency w, mfr is the relative magnetic permeability at infinite frequency, msr is the relative magnetic permeability at zero frequency, m_tau is the generalized relaxation time for the permeability, and m_alpha accounts for multiple relaxation phenomena in the permeability. The program is written for the International System of Units; the following table lists the units for the various quantities that are either input to or output from the program: Where it Quantity appears Unit ------------------------ -------- ------------------ time input second distance input meter conductivity input Siemens/meter frequency input Hertz Hertz potential output Volt-meter electric field intensity output Volts/meter magnetic field intensity output Ampere-turns/meter The source code for program 'xwell_radar' is in 37 different files, of which 6 are in directory 'xwell_radar' and 31 are in directory 'bessel'. The following table lists the most significant files and gives a brief description of the operations performed by the code in those files: Brief Directory File Description ----------- ------------- ------------------------------------- xwell_radar xwell_radar.c This file contains the main part of the program. The code reads the input parameters, writes them, checks them, calculates the normalizations, applies the normalizations, calculates the traces, removes the normalizations from the traces, and writes the traces. xwell_radar calc_fld.f Calculates the electromagnetic field at one frequency. xwell_radar fft_d.f Performs a fast Fourier transform. xwell_radar prefft_d.f Calculates some quantities needed for a fast Fourier transform. xwell_radar file_open1.c Opens a file for reading or writing. xwell_radar mem_alloc1.c Dynamically allocates memory for arrays. bessel jn.f Acts as interface to subroutine zbesj. bessel zbesj.f Calculates, for a complex argument, Bessel functions of the first kind. The original versions of 'fft_d' and 'prefft_d' were published by Marple (1987) as 'fft' and 'prefft'; 'fft_d' and 'prefft_d' have been modified so that they perform double precision arithmetic. Subroutine 'zbesj' and its associated subroutines were written by D. E. Amos and are available also at 'http://www.netlib.org'. ---------------------------------------------------------------------- 3.2 Installation The installation requires three steps: 1. Copy the tar file to a suitable directory on your Unix computer. Extract the contents of the tar file using the command: 'tar xf xwell_radar.TAR'. During the extraction, three directories are created: 'xwell_radar', 'bessel', and 'test_data'. 2. Make the current directory 'bessel'. Use the file 'makefile' to compile the Fortran-language subroutines as object modules and then to put the modules into a library called 'bessel.a'. 3. Make the current directory 'xwell_radar'. Use the file 'makefile' to compile the c-language functions and Fortran- language subroutines as object modules and to load them as an the executuable file 'xwell_radar'. Program 'xwell_radar' was developed on a Hewlett-Packard computer (HP-9000) with the Unix operating system (version 10.20 of the HP-UX). Installation of program 'xwell_radar' on another computer may require some modification of the source code and the makefiles. ---------------------------------------------------------------------- 3.3 Input Data In directory 'test_data' is the sample input file 'in'. This input file is for a 3-layer model, for which a diagram is shown in section 3.1. Because the input parameters are fully described in this file, only a few remarks are necessary: 1. All input files must be similar to 'in': No text, except that related to the properties of the layers and the positions of the receiving antenna, may be changed --- only the numerical values may be changed. If another layer is added to the model, then more text and numerical values describing the properties of that layer must be added. The format for this additional text is identical to that in 'in'. Likewise, if more positions of the receiving antenna are added, then the format for the additional text is identical to that in 'in'. 2. The units for the input parameters are listed in section 3.1. 3. A suitable value for the normalization for time is roughly 1.e-9; a suitable value for the normalization for distance is roughly 1. If other values are chosen, they should be within a factor of 10 of these two values. 4. The parameters related to the dielectric permittivity and the magnetic permeability are for the formula developed by Cole and Cole (1941); see section 3.1. The sample input file 'in' shows parameters that make the dielectric permittivity and the magnetic permeability frequency-independent: For the dielectric permittivity, for example, the relative values at zero frequency and infinite frequency are identical, the generalized relaxation time is zero, and the exponent accounting for multiple relaxation phenomenon is zero. 5. The antennas cannot be in the top or the bottom layers, which are half-spaces. If the antennas must be in a half-space, there is a simple fix: Make an extra layer having the electromagnetic properties of the half-space, and put the antennas in it. 6. The thicknesses of the top and the bottom layers do not matter: Any values can be used. 7. The number of samples in a trace must be a power of 2: for example, 128, 256, or 512. 8. The selection of a suitable time shift for the Ricker wavelet is described in section 3.1. 9. The heights of the transmitting and the receiving antennas are measured from the bottom of the layer in which they are located. The radial distance from the transmitting antenna to the receiving antenna is measured parallel to the layer interface. (The radial distance is not the straight line distance!) Heights and radial distance are specified in this way to conform to the circular cylindrical coordinate system --- see section 3.1. 10. If the transmitting antenna or the receiving antenna or both are very close to the top or the bottom of a layer (that is, very close to an interface), numerical problems develop. 11. The source code is currently set for a maximum of 20 layers in the model and a maximum of 20 positions of the receiving antenna. These two values should be suitable for most situations. If these values must be increased, follow this procedure: Load file 'xwell_radar.c' into an editor and go to line 10. This line and the next are: #define MAX_POSITIONS 20 #define MAX_LAYERS 20 Change these c-language macros to the desired values. If MAX_LAYERS is changed then load file 'calc_fld.f' into the editor and go to line 27. This line is: parameter( stderr = 7, max_layers = 20, max_steps = 1000 ) Make the same change to parameter max_layers. Finally, recompile program 'xwell_radar'. ---------------------------------------------------------------------- 3.4 Execution To execute the program using the sample input file 'in',make the current directory 'test_data'. Type this command: '../xwell_radar/xwell_radar & list_test'. (For this input file, the total execution time for the program is about 10 seconds on an HP-9000 computer.) The output files are described in section 4. ====================================================================== 4. Test Data The 7 files in directory 'test_data' are described in sections 4.1 through 4.4. ---------------------------------------------------------------------- 4.1 File 'in' This is the input file, and it is described in section 3.3. This same input file was used by Ellefsen (1999) to demonstrate the accuracy of the traces computed with program 'xwell_radar'. To compute traces, follow the procedure in section 3.4. Five output files are created: 'list_test', 'trace_e_test', 'trace_h_test', 'trace_d_test', and 'trace_p_test'. ---------------------------------------------------------------------- 4.2 File 'list' This is the main output file, which lists the input parameters, the normalizations, and the traces that were computed. The file 'list_test' should be identical to 'list'. Both 'list_test' and 'list' have a warning about the parameters used to calculate the Ricker wavelet. This warning was generated because the time shift specified in file 'in' (11. ns) is slightly larger than the recommended value (0.85 ns). ---------------------------------------------------------------------- 4.3 Files 'trace_e', 'trace_h', 'trace_d', and 'trace_p' These four files contain traces that show components related to the the electromagnetic field: File Component ------- ---------------------------------------------------------- trace_e Vertical component of the electric field intensity trace_h Magnetic field intensity trace_d Vertical component of the electric field intensity for the direct wave trace_p Hertz potential Each file contains two traces corresponding to the two positions of the receiving antenna (see file 'in'). Each of the two traces consists of 512 samples, which are single-precision, floating point numbers (4 bytes). The numbers are stored in binary format, and they can be viewed with the command 'od -tfF trace_e | more', for example. The units associated with these numbers are listed in a table in section 3.1. For example, the unit for the numbers in file 'trace_e' is 'Volt/m'. The sample interval equals the duration of the trace divided by the number of the samples. For example, using the values in file 'in', the sample interval is 200. ns / 512, which equals= 0.3906 ns. The sample interval times the sample number, minus 1, equals the time for that sample. This table lists some times for the traces: Sample Sample Number Time Number minus 1 (ns) ------ ------- ------ 1 0 0. 2 1 0.3906 3 2 0.7813 ... ... ... 510 509 198.8 511 510 199.2 512 511 199.6 Notice that the time for the last sample is slightly less than 200 ns; this difference occurs because of the properties of the discrete Fourier transform. In practice, this difference is insignificant. Files 'trace_e_test', 'trace_h_test', 'trace_d_test', and 'trace_p_test' should be identical to 'trace_e', 'trace_h', 'trace_d', and 'trace_p', respectively. ---------------------------------------------------------------------- 4.4 File 'shell.plot' This file is a shell program that displays, in an X-window, the four files described in section 4.3. This shell program can only be used on those computers with the SU programs that are distributed by The Center for Wave Phenomenon at the Colorado School of Mines. To use this program, type 'shell.plot'. ====================================================================== 5. Disclaimer This open-file report was prepared by an agency of the United States Government. Neither the United States Government nor any agency thereof nor any of their employees makes any warranty, expressed or implied, or assumes any legal liability or responsibility for the accuracy, completeness, or usefulness of any information, apparatus, product, or process disclosed in this report or represents that its use would not infringe privately owned rights. Reference therein to any specific commercial product, process, or service by trade name, trademark, manufacturer, or otherwise does not constitute or imply its endorsement, recommendation, or favoring by the United States Government or any agency thereof. Although all data and software in this open-file report have been used by the USGS, no warranty, expressed or implied, is made by the USGS as to the accuracy of the data and related materials and (or) the functioning of the software. The act of distribution shall not constitute any such warranty, and no responsibility is assumed by the USGS in the use of these data, software, or related materials. ====================================================================== 6. Acknowledgements This work was supported by Mineral Resources Program of the U. S. Geological Survey. ====================================================================== 8. References Bouchon, M., 1981, A simple method to calculate Green's functions for elastic layered media: Bulletin of the Seismological Society of America, vol. 71, no. 4., p. 959-971. Cole, K. S., and Cole, R. S., 1941, Dispersion and adsorption in dielectrics, I., Alternating current characteristics: Journal of Chemical Physics, vol. 9, no. 4, p. 341-351. Ellefsen, K. J., 1999, Effects of layered sediments on the guided waves in crosswell radar data: Geophysics, vol. 6, no. 64, p. 1698- 1707. Marple, S. L., Jr., 1987, Digital spectral analysis with applications: Englewood Cliffs, New Jersey, 492 p. Wait, J. R., 1995, Electromagnetic Waves in Stratified Media: New York, MacMillan, 372 p. ======================================================================