ESSA TR ERL 165-ITS 106 A UNITED STATES DEPARTMENT OF COMMERCE PUBLICATION * f 4rts o» / ESSA Technical Report ERL 165-ITS 106 U.S. DEPARTMENT OF COMMERCE Environmental Science Services Administration Research Laboratories Theoretical LF, VLF Field Calculations With Spherical Wave Functions of Integer Order J. RALPH JOHLER CARLENE MELLECKER BOULDER, COLO. APRIL 1970 ESSA RESEARCH LABORATORIES. The mission of the Research Laboratories is to study the oceans, inland waters, the lower and upper atmosphere, the space environment, and the earth, in search of the under- standing needed to provide more useful services in improving man's prospects for survival as influenced by the physical environment. Laboratories contributing to these studies are: Earth Sciences Laboratories: Geomagnetism, seismology, geodesy, and related earth sciences; earthquake processes, internal structure and accurate figure of the Earth, and distribution of the Earth's mass. Atlantic Oceanographic and Meteorological Laboratories: Oceanography, with emphasis on the geology and geophysics of ocean basins, oceanic processes, sea-air interactions, hurricane research, and weather modification (Miami, Florida). Pacific Oceanographic Laboratories: Oceanography; geology and geophysics of the Pacific Basin and margins; oceanic processes and dynamics; tsunami generation, propaga- tion, modification, detection, and monitoring (Seattle, Washington). Atmospheric Physics and Chemistry Laboratory: Cloud physics and precipitation; chem- ical composition and nucleating substances in the lower atmosphere; and laboratory and field experiments toward developing feasible methods of weather modification. Air Resources Laboratories: Diffusion, transport, and dissipation of atmospheric con- taminants; development of methods for prediction and control of atmospheric pollution (Silver Spring, Maryland). Geophysical Fluid Dynamics Laboratory: Dynamics and physics of geophysical fluid systems; development of a theoretical basis, through mathematical modeling and computer simulation, for the behavior and properties of the atmosphere and the oceans (Princeton, New Jersey). National Severe Storms Laboratory: Tornadoes, squall lines, thunderstorms, and other severe local convective phenomena toward achieving improved methods of forecasting, detecting, and providing advance warnings (Norman, Oklahoma). Space Disturbances Laboratory: Nature, behavior, and mechanisms of space disturb- ances; development and use of techniques for continuous monitoring and early detection and reporting of important disturbances. Aeronomy Laboratory: Theoretical, laboratory, rocket, and satellite studies of the physical and chemical processes controlling the ionosphere and exosphere of the earth and other planets. Wave Propagation Laboratory: Development of new methods for remote sensing of the geophysical environment; special emphasis on propagation of sound waves, and electro- magnetic waves at millimeter, infrared, and optical frequencies. Institute for Telecommunication Sciences: Central federal agency for research and services in propagation of radio waves, radio properties of the earth and its atmosphere, nature of radio noise and interference, information transmission and antennas, and meth- ods for the more effective use of the radio spectrum for telecommunications. Research Flight Facility: Outfits and operates aircraft specially instrumented for re- search; and meets needs of ESSA and other groups for environmental measurements for aircraft (Miami, Florida). ENVIRONMENTAL SCIENCE SERVICES ADMINISTRATION BOULDER, COLORADO 80302 U.S. DEPARTMENT OF COMMERCE Maurice H. Stans, Secretary ENVIRONMENTAL SCIENCE SERVICES ADMINISTRATION Robert M. White, Administrator RESEARCH LABORATORIES Wilmot N. Hess, Director ESSA TECHNICAL REPORT ERL 165-ITS 106 Theoretical LF, VLF Field Calculations With Spherical Wave Functions of Integer Order J. RALPH JOHLER CARLENE MELLECKER CO - i INSTITUTE FOR TELECOMMUNICATION SCIENCES BOULDER, COLORADO April 1970 For sale by the Superintendent of Documents, U. S. Government Printing Office, Washington, D. C. 20402 Price 75 cents. FOREWORD This work was completed on the DIDS project of the U. S. Army Strategic Communications Command, Fort Huachuca, Arizona 85613 as part of the work of interagency work order DAHC -20-69-C -Oil 6. The computer programs were developed for the Defense Atomic Support Agency under interagency work order 805-69o Test computa- tions made in connection with this development are also included. Permission to quote data points plotted in figures 30, 31, 32 and 33 (Morgan, 1966) was granted by the Department of the Army, Office of the Secretary of the Army, Office of Civil Defense. The FORTRAN listing of the computer program described in Section 7 is available from the authors upon request. TABLE OF CONTENTS Page FOREWORD ii ABSTRACT iv 1. INTRODUCTION 1 2. COMPUTATION AND DISCUSSION 2 3. THEORY OF n-WAVES 19 4. CONCLUSIONS 24 5. REFERENCES 26 6. GLOSSARY 29 FIGURES 34 7. COMPUTER PROGRAMS 77 in ABSTRACT The propagation of LF and VLF terrestrial radio waves is simulated on the computer with a precise theory of spherical wave functions of integer order. Anisotropic reflections from an electron and ion density profile of the ionosphere are accommodated by the analysis. Variations in the profile and the magnetic field 'with distance along the propagation path can also be introduced ad hoc. The analysis method is incorporated into a computer program which is written in a flexible manner so that it can be applied to a variety of scientific studies of low-frequency ionospheric wave propagation. Calculations of current interest are presented to illustrate use of the technique. Amplitude perturbations of approximate- ly 1 5 decibels are described at 60 kHz as a result of an electron density profile perturbation. It is also found that the 26 kHz signal is compara- tively insensitive to a variation in magnetic dip angle when propagation occurs in the daytime near the magnetic equator. Thus, using values of the magnetic field vector at the center of the propagation path is a good approximation in certain circumstances. KeyWords: MF, LF, VLF; low frequency wave propagation; propagation; terrestrial radio waves; low frequency radio waves; radio waves; spherical waves; spherical wave functions. THEORETICAL LF, VLF FIELD CALCULATIONS WITH SPHERICAL WAVE FUNCTIONS OF INTEGER ORDER J. Ralph Johler and Carlene Mellecker 1. INTRODUCTION The formulation of a theory of propagation for MF, LF, and VLF ionospheric radio waves using spherical wave functions of integer order as building blocks was presented in a previous work (Johler, 1970). Ordinarily, the spherical wave functions of integer order are found to be useful only at ELF or for propagation problems in which ka is small, where a is the radius of the sphere and k is the wave number (Johler and Lewis, 1969). In the approach presented in this paper, this restriction is no longer necessary because of the properties of the j -series summation used. Indeed, numerical testing of the method at medium frequencies (MF) indicates that the upper frequency limit of usefulness of spherical wave functions had not as yet been found. The theory has been trans- lated into a FORTRAN computer program for use in scientific investi- gations of the ionospheric radio waves. The philosophy in the computer program design emphasized flexibility in the application of the analysis method to a variety of scientific studies concerned with MF, LF and VLF wave propagation. Obviously, extensive calculation of the field, such as that employed in system studies, would require optimization of the computer program for a specific system application to insure economy of computer time. This paper demonstrates the application of the method of spheri- cal wave functions of integer order to LF and VLF. Calculations of interest are presented and compared with experimental data. The FORTRAN computer programs necessary for the calculations pre- sented are also described. 2. COMPUTATION AND DISCUSSION The theoretical approach used in this paper is the geometric series propagation formula (Johler, 1970) CO E = E„_ + ) E . . (1) 'r ?s ° Zj t> ^ 3=1 where E r is the ground wave (Johler et al. , 1956; Johler and Morgenstern, 1965; Johler, 1969a, b) and E i are the terms of the geometric series, j = 1, 2, 3 • * • , depicted as optical rays in figure 1. It should be emphasized that there is no set of optical rays of the homogeneous, plane wave type that will give the exact solution pre- sented in this analysis except as a very approximate estimate. How- ever, a rigorous set of n-waves, discussed in detail in section 3 will be referenced in the physical discussion to follow. The n-waves are almost equivalent to homogeneous plane waves with respect to the ionosphere boundaries. If the terrestrial waveguide is assumed to be sharply bounded at the top with a reflection coefficient, T, # = -1, various ionospheric waves, E rj j , j = 1, 2, 3 * * " , can be calculated (Johler, 1970) from t the formula , tQuantities B,G,R, ,p etc. , are defined in section 6. -2- co -4(-l) J G(9,n) pi Ri" 1 E r j = B ) • U) C -1 a C -1 a |_ -^ n C -1 a + ~^~ ln C - "3 a J ft=o Each term of this summation (2) is a spherical wave of integer order or an n-wave (Johler, 1970). The theoretical basis of such n-waves is discussed in section 3. The spherical reflection coefficient for a TM-wave (vertically polarized electric vector) on the ground, R, , is a function of the ground conductivity, o, mhos/m, as is the wave number, k_ 2 . If o = 5 mhos/m (sea water), the amplitude, |E r i \, is illustrated in figure Z, j =1 through 10. The height of the waveguide is h = 65 km. The frequency, f = 10 kHz. The curves presented are therefore independent of the complicated ionosphere reflection coeffi- cient and (except at very short distances) illustrate the basic character of the various terms of the j -series (or hop series). Each term, E r ,j, at short distance, d/ j, is a monotonic curve. Undulations then appear as a function of distance as d/ j is increased. In the region of the geometric -optical horizon of the ray theory, and beyond, as depicted in figure 1, the wave theory gives a curve that shows an exponential type of attenuation characteristic of diffracted waves. The standing wave as a function of distance in the intermediate region is evident at a certain critical angle of incidence, cp i , on the ground. This sug- gests lateral wave displacement (Brekhovskikh, I960), since a critical domain of the complex angle of incidence on the ground or the iono- sphere is necessary to excite such waves (Johler, 1970). Since, T«, = -1» these are lateral waves along the ground. At great distance, the waves diffract around the curvature of the ground. In such a region a considerable number of n-waves or terms of the summation over n, discussed in section 3, are important in the series (2). As a consequence of this fact, the curves are not 3- precisely exponential (linear on a semilog graph). In figure 3, the ground wave, E r , is also shown for comparison. A higher frequency (2 6 kHz) is illustrated in figure 4, assuming a ground conductivity = 0. 005 mhos/ m. At short distances the j -series terms ordinarily have small amplitude. With T,. = -1, the amplitudes illustrated in figures 2, 3 and 4 are abnormally large. Thus, except at very short distances, the basic features shown in the figures 2, 3, 4 dominate the terms of the j -series. More complicated effects can be introduced by using more realistic reflection coefficients for the ionosphere. An exponential model of the ionosphere and a two layer model ionosphere developed by Thrane* in 1969 are illustrated in figure 5„ An alternatie "two -layer model" profile developed by Bain and May (1967) is illustrated in figure 6. In a "two -layer model" the electron density tends to form a clear maximum in the D-region near 60 km. The electron density then decreases before increasing again. Both electron density and d-c conductivity for the electrons, o~ # , Ne 2 ,~a a = (3) e m e V are shown in the figures as a function of altitude above the ground level, km. A theoretical daytime noon model developed by W. S. Knapp** in 1967 based on computer programs of the type given by D. C. Hyovalti*** in 1965 , is shown in figure 7. Both the electron density and the positive ion density are given. The negative ion density can be obtained from conservation of change. N + = N e + N- . (4) * E. V. Thrane, Private Communication, ITS, ESSA, Boulder, Colo. *# W. S. Knapp, Private Communication, G.E. Tempo, Santa Barbara, Calif. #**D. C. Hyovalti, Unpublished Report, ITS, ESSA, Boulder, Colo. _4_ Then, the total d-c conductivity is given by 0. , = H^L + 2 2k* , (5) where v ' = v Els (a = e,+) a '° a '° m a + m ' v and where the mass m, of oxygen and the mass m. of the electron can be used. Of course, the effect of ions is quite small in the normal ionosphere since m e < < m+ . (6) Collision frequencies used by Johler and Berry (1965) were used in all calculations involving v #J . Positive ion collision frequencies, v + , and negative ion collision frequencies are also given by Hirsch (1967). Under disturbed conditions such as polar cap events or nuclear events (Johler and Berry, 1965), it may be necessary to consider conductivity caused by ions. A daytime profile similar to the profile in figure 7 is given in figure 8. However, the electron densities at the lower altitudes between 40 and 70 km become negligible more abruptly with decreasing altitude. There is some uncertainty concerning the precise form of the lower boundary. This can affect the propagation of terrestrial radio waves in the waveguide below the ionosphere. Belrose and Segal (1967) give comparative profiles of electron number density at the lower boundary of the ionosphere. The model shown in figure 8 is compared in figure 9 with a nighttime model developed from similar computer calculations. Also, the corresponding collision frequency is given. Of course, other pro- files are available in the literature that can be used to model the ionosphere (Belrose, 1968). To demonstrate the use of spherical wave •5- function methods we shall use the models given in figures 5 through 9 in the computer program discussed in section 6. The models given in figures 7 and 8 are regarded by us as more representative of daytime noon conditions of the ionosphere. The lat- ter is a more sharply bounded model. The models given in figures 5 and 6 are special models that may be used to explain special pheno- mena, such as the effect of a tendency to form a lower layer of electrons. In such discussions the exponential model given in figure 5 is often referenced. The available data on the nighttime ionosphere is not as abundant as the daytime noon model. The nighttime model given in figure 9 may or may not be typical. In (2), let p, = (-1)* pi Rl" 1 ; (?) then, an isotropic reflection coefficient can be introduced by replacing -1 with T e-J Cj = (T e J* ^ Rf 1 , (8) where T 09 is the reflection coefficient for a wave incident on the ionosphere with the electric vector of the incident and reflected wave in the plane of incidence. This is called a TM-wave (transverse magnetic) reflection. If the ionosphere is anisotropic, C 3 is more complicated. Thus, (Johler, 1961), Ci - T 9e >g = XV e J. ee i £\. m X 6m J. me C 2 - R e T ee + R ra T- m T, 2 rp 3 ee C 3 - 2R e Rn T 9e T em T me + R T, (9) "* ^m. -"-mm ■'■am me 6- A simple system for working out the first few Cj is depicted in figure 10 using a geometric notion of rays reflected to and from between the two waveguide boundaries. The rigorous matrix form of this solution is given in other papers, Johler (1970) and Johler (1966). Recently, Berry, et ah (1969) have worked out this matrix expansion so that ionosphere reflection coefficient effects can be separated from ground reflection coefficient effects as was originally postulated in previous works (Berry, 1964; Berry and Chrisman, 1965). It is quite apparent from figure 10 that the expressions for the Cj ' s can become quite complicated and uninteresting. However, the first few terms of the geometric series are instructive from a physical point of view. Thus, considering the first term, j = 1, of the geometric series, the incident vertical electric polarization gives rise only to TM -wave propagation between the transmitter and the receiver. While it is true that an abnormal component is generated at the ionosphere by the conversion coefficient, T em , the assumed vertical polarization of the source and observer dipoles does not involve this component in the waveguide propagation mechanism shown. For two reflections, j = 2, the reflected TE -wave (dashed ray) is reflected at the ground and then converted back into a TM-wave at the second reflection by the operation T em T ne . Of course, the conversion coefficient, T em , indicates conver- sion occurs in a region [2, 1] on the second term j = 2, first reflection k = 1 of the set [j,k], while T n0 indicates conversion occurs in a region \_2, 2] at a greater distance from the transmitter. Obviously the process can be continued ad infinitum, as indicated by j = 3 and 4 in figure 10. Thus, in general a region labeled [j,k] is the k th reflec- tion (beginning at the transmitter and numbering the regions with the k = 1 , 2, 3 • • ■ integers) of E r> j . When j = 3 or greater, the TE - wave reflection coefficient, T n n , becomes significant to propagation in which both transmitter and receiver are vertical dipoles. This is a consequence of the conversion effect of the terrestrial magnetic field at the ionosphere since the source does not give rise to TE -waves. Thus, TM and TE propagation modes become coupled at the ionosphere. As the terrestrial magnetic field is theoretically forced to vanish, T en = T B e = 0, and propagation is controlled by T e , (TM-waves) exclu- sively. Note also that the ground reflection coefficient, R m , for TE - waves is not used in the propagation mechanism if the ionosphere is isotropic. It is clear from the indexing [j,k] that regions of the ionosphere have been identified along the propagation path between transmitter and receiver. These regions can be tagged with the aid of points represent- ing the center of the reflecting region [j,k] by drawing geometric - optical rays depicted in figure 1. The waveguide height has been postu- lated as h. This is a physical parameter. In actuality, the geometric series expansion is ordinarily performed in a thick bottom layer repre- senting the earth-ionosphere waveguide as is illustrated in figure 11. In a rigorous solution, the ionosphere is then represented by concentric spherical shells to great height such that decreasing the shell thickness or increasing the number of such shells or decreasing the height to which such shells are emplaced produces negligible effect on the field propagated from S to O. The precise thickness of the waveguide shell in which the j -series expansion is performed is then primarily a matter of phase reference for the reflection coefficient calculation, although a choice can often be found in which the phase of the reflection coefficient, T c0 varies slowly in the summation process (2). This slowly varying phase reference can be used to improve the efficiency of the calculation, since the number of values of reflection coefficient calculated can be decreased. Values of T e , = T, e (n) are found for each n by interpolation in a calculated set of T e# (cp 1 ) as a function of cp t , related to n by (see sec. 3), ■8- r ( 2 )' sincp i = [ _i %r] ' (10) V n(n+ 1) sincpi ^ ^^ i , (n) K o g where g = a + h, k = — , and a is the radius of the terrestrial sphere, a~ 63 60 km. Although the computer program can accommodate complex angles of incidence, cp 4 , it is evident from the approximation (11) that the angle is almost real. In the presence of magneto -ionic coupling of the wave into the terrestrial magnetic field lines, the phase of the reflection coefficient does not usually display such a minimum phase variation with angle of incidence, cp t . Under such conditions a quite high density of values of reflection coefficient as a function of cp t or n must be calculated for high accuracy. Ordinarily, the bottom of the lowest spherical shell (figs. 1 through 11) at a radial distance r = a + h, is used as the reference point for the phase of the reflection coefficient. If a criterion of reflec- tion height is used in such a way that the phase is slowly varying as a function of cp ± , the phase of the reflection coefficient should be re- referenced to the height used (Budden 1961). Thus if « means "is replaced by" , T, t « T.. exp [2ik coscp i (h - h )] , (12) where h is the bottom of the ionosphere where the reflection coeffi- cient is calculated, and h is the height to which the reflection coeffi- cient is to be re -referenced. It is quite apparant that the phase of the reflected wave is retarded if h > h , and advanced if h < h . Equa- tion (12) is restricted in usage to plane wave reflection from a plane ionosphere. The use of spherical shells would require more compli- cated modification (Johler, 1970). In many calculations the use of avail- able reflection coefficients that are based on plane wave, plane iono- sphere models is quite accurate and convenient. A computer program is described in section 7 for computing the reflection coefficient of an anisotropic ionosphere using an exact method for the reflection of a plane -wave from a plane -ionosphere with continuous stratification of the profile with altitude. It is used with calculations performed with the spherical wave function programs which are also given in section 7. Some further discussion of the geometric series expansion in the thick bottom isotropic shell, (fig. 11) containing the transmitter S and the receiver O (the source and observer) is appropriate to a further understanding of the nature of the reflection coefficient. After the height of the bottom shell, h, has been established, the jth term of the series (1) can be calculated by the summation (2). Each term of the summa- tion (2) is a spherical wave of integer order (Johler, 1970) or an n-wave. The theoretical basis of such n-waves is discussed in section 3. One such n-wave is depicted in figure 11, for j = 1 and j = 2. The upgoing wave is shown with a solid line. The downgoing wave is shown with a large dashed line. Since a particular n-wave has been chosen, the angle of incidence with the ionosphere is given by (10) or approximately by (11). A short-dashed line shows the angle of incidence of the geometric - optical ray. Note that the angle of incidence of the geometric -optical ray and the n-wave do not necessarily correspond. One can make the conjecture that important contributions to the sum of the n-waves (2) will occur when the angle of incidence of the n-waves approaches the geometric -optical wave. This, indeed, has been found to be precisely the case when the series is summed for models presented in this paper. In fact, the n-waves with very steep angles of incidence, cp t , on the ionosphere were found to be severely attenuated as were the almost -10 grazing incidence waves. In the latter case, the angle of incidence became complex (-Imsincpi > > 0) . Suppose the n-wave depicted in figure 1 1 is the lower limit n-wave below which smaller n values given negligible contribution. Then as equation (2) is summed, it can be imagined that the regions \.\, l~\, [2,1], [.2,2] ■ • • , are the important regions for the terms of the j -series. Then, instead of an illuminated point used in geometric -optical theory, we have an illuminated point set extending over a finite region of the propagation path. It is also clear that the n-waves of small order n must exhibit considerable lateral displacement along the ionosphere. One must at the same time admit lateral displacement along the ground. This suggests physical modifi- cations of the analysis which can take into account local regions of the ionosphere [1,1], [2,1], [2,2] *•• and local regions of the ground [1,1], [1,2], [2,1], [2,2], [2,3] ••• , with different reflection coefficients, T,,(j,k), R,(j,k) ••• . The reflection coefficients T,,(j,k) may be a consequence of profile or magnetic field variations along the propagation path. The reflection coefficient R, (j,k) may be a consequence of ground conductivity or ground impedance changes along the propagation path. Suppose for example, the ionosphere is disturbed over a region of radius R, as depicted in figure 12. If the region R were of suffi- cient size, it could be represented by a spherical reflection coefficient T,, and a reference height that may be the reflection height, h . Also, a disturbance at some shorter distance (fig. 13) may have no effect on the term E r fl of the geometric series and the normal reflec- tion coefficient would be appropriate. But the E r j2 term of the series would require a special reflection coefficient and reflection height. The E r 3 wave may be unaffected. But the total field, E r , of equation (2) will be affected by this disturbance. The case depicted in figure 12 may also affect the E r 3 term of the geometric series, as 11- shown in figure 14 to the exclusion of the E r 2 term. Finally, we note the E r 4 term could also be disturbed by a disturbance in the region [4,2], i.e., the fourth hop, second reflection (fig. 15). The flexibility described above has been incorporated into the computer program described in section 7„ With some imagination, it is not difficult to envision further research on the nature of propagation via the D -region using these computer programs as a tool. It should be emphasized that a sufficient number of n-waves should be included to assure convergence of the n-series. Also, a sufficient number of j -waves or terms of the j series must be included for the desired accuracy. Usually 3 or 4 j -waves will suffice. If a large number are required, the above described additional flexibility feature becomes difficult to apply for obvious reasons. Finally, the lower limit of the n-series summation should be carefully determined with appropriate tests to be sure that the steeply incident waves are completely negligible. A rigorous method for calculating plane -wave, plane -ionosphere reflections in the presence of arbitrary magnetic induction and arbitrary variation of ion or electron number density and collision frequency with altitude was given by Johler and Harper (19 62). The use of this method was reviewed more recently by Johler (1967). Thus, in addition to the profile at the reflecting region [j,k], the magnetic azimuth, cp a , the magnetic dip angle, I, and the magnetic intensity, \ti\ , are required to calculate a reflection coefficient matrix T, T _ r T ee T em -1 ^ (13) T T me mm where, T= T(j,k) . (14) 12- Quite frequently, a particular investigation of low frequency- propagation will require only average parameters over the entire propagation path. Thus, the magnetic field parameters at the center of a long propagation path in the daytime at VLF usually suffice to describe the propagated field. An example of a 26 kHz field propagated at a magnetic azimuth, cp a ~ 147* from Hawaii will require magnetic para- meters given in table 1 for each member of the set Lj>k3, for a distance of 6180 km between transmitter and observer. Table 1. Magnetic Parameters for Figures 16 and 17 [j,k] 1, 1 2, 1 2, 2 3, 1 3, 2 3, 3 4, 1 4, 2 4, 3 4, 4 U|A/m cp a degrees I, degrees 25.8 147.8 - 3.1 26.0 147.8 21.2 28.4 144.9 -26.9 26. 6 147.0 28.2 25.8 147.8 - 3.1 29.7 147. 2 -33.8 27.0 146. 5 31.4 25. 6 148.2 9.5 2 6. 8 146.7 -15. 5 30.4 142. 3 -37.0 The amplitude of the field |E r | and the fields JE r> j | , j = 1, 2, 3, * • • , were calculated as a function of distance (fig. 16). The normal daytime noon profile (fig. 8) was employed. The magnetic parameters, [1,1], table 1, at the center of the propagation path for d = 3090 km were used for the entire path. In figure 17, the magnetic parameters for each reflecting region, [j,k], given in table 1, were inserted into the calculation so that a separate reflection coefficient for each reflecting region [j,k] was calculated. As the distance changes it was assumed that the reflecting regions [j,k] each remained the same, -13- It is noted that only minor changes in the field occur as a consequence of a variable magnetic field at each reflecting region. It should be noted that this propagation crosses the magnetic equator in such a way that the dip angle, I, table 1, changes from positive to negative along the propagation path. We therefore conclude that the value of magnetic field at the center of the propagation path may be a good assumption for daytime propagation over such a path. Thus, instead of calculating 10 sets of reflection coefficients, one set for the center of the propaga- tion path would suffice if accuracy of the order of a decibel were sufficient. The detailed nature of the terms of the geometric series, E r> j , j = to 7 illustrate typical behaviour for such terms as a function of distance. Thus, at great distance, d/ j, the amplitude decays almost exponentially (linear semilog graph). At shorter distances, d/j, i l small undulations exist in the curves. The total field, ^~ , does not L $=° reflect much of this detail but does show undulations as a function of distance resulting from phase interference between terms of the geometric series. The exponential and two -layer profiles, Thrane (1969) (fig. 5) and the Bain and May (1967) two -layer profile (fig. 6) were used to calculate the amplitude and the phase correction as a function of distance at 20 kHz, see figures 18, 19, 20, 21, 22, 23, 24, 25, and 2 6. The phase correction, cp c , is defined cp Q = -Arg E r - kid . (15) The amplitude is normalized to unity dipole current moment. Thus, the radiated power, P r , is given by: P r = 0.4(10- ls ) w 2 (I o is illustrated in figure 24 for this two layer model. Also, the corres- ponding phase difference, ArgE r>1 -ArgE r> r, o » is given in figure 2 6. The corresponding curves for the exponential profile in figure 5 are shown in figures 25 and 27. These curves were calculated using sea water ground constants, 0" = 5, £% - 80. However, the difference between the land, 0=0. 005, e 2 = 1 5 curves and the sea water curves is small. This can be ascertained by com- paring figures 18 with 19 and 20 with 21. A standing wave, as a function of distance at short distances, is shown in figure 24 which is apparently a consequence of a wave trapped by the two-layer model. It is also clear that the phenomenon completely disappears in the exponential model but also exists in the Thrane (1969) two -layer model of figure 5. It would be possible to observe such a wave if an experimental method to eliminate -15- the ground wave could be devised. Thus, the presence of the undulation as a function of distance, (i. e. , a standing wave), would be an indication of multiple D -region layers. The model of the D -region illustrated in figure 7 is typical of the D -region in the United States. The 60 kHz field as a function of distance was calculated for the Ft. Collins, Colorado, NBS transmitter WWV. The magnetic azimuth was assumed to be cp a =90° and 270° . The dip angle was assumed to be I = 70° . The magnetic intensity was assumed to be, \ 34 \ = 40 A/ m. The ground constants a = 0. 005 mhos/ m and e 2 = 15 were assumed. The two theoretical curves, cp a = 90° and 270°, are given as solid lines in figures 30, 3l, 32, and 33. Measured data of Morgan (1966) are plotted for various radials from the transmitter. The absolute field for 13 kw radiated power is given on the vertical scale in decibels above one volt per meter. Figure 30 gives the mea- sured data on a radial toward Palm Beach, Florida. Figure 31 gives the measured data on a radial toward Cape Fear, North Carolina. Figure 32 gives the measured data on a radial toward Brownsville, Texas. Finally, figure 33 gives the measured data on a radial toward Nantucket, Massachusetts. No attempt was made to determine ground constants more accurately or adjust profile so that data would fit more accurately. The agreement with theory is considered to be quite excellent. These data are quite remarkable, since data of amplitude vs distance is rare and can be obtained only after considerable effort. One could of course modify the profile and other parameters to obtain a more precise fit of the data for the various propagation directions, but determining a para- meter or parameters to change to obtain the desired result would require more research. To gain some notion of the extent of the amplitude perturbation that may occur as the ground constants of a propagation path are changed or as the ionosphere changes with time, a series of 60 kHz curves of 16- amplitude as a function of distance are presented in figures 34 through 43. These curves are normalized to unity dipole current moment. Figures 34, 35, 37, 39, 40, and 42 employ the model given in figure 7. Figures 3 6, 38, 41, and 43 employ the model given in figure 8, i.e. , the alternate normal daytime noon model with somewhat less ionization between 40 and 70 km. This latter model is considered to be more sharply bounded than the former (fig. 7). A single nighttime model given in figure 9 is used on the nighttime curve on all the figures 34 through 43. Figures 34 and 35 illustrate propagation northward, cp a = 0° and eastward, cp a = 90° assuming a ground conductivity 0" = 0. 005 mhos/m. Although both daytime and nighttime curves of amplitude as a function of distance are given, the principal change as a result of the magnetic azimuth change occurs at night. The change comprises a shifting of the undulations such that at a fixed distance both an increase and a decrease in amplitude are possible with a change in cp a . The numerical details of the computation indicate a strong coupling of energy into the whistler mode. This causes the phase of the reflection coeffi- cient to vary rapidly with angle of incidence on the ionosphere, i. e. , with distance. This phenomenon, together with the normal interference between the ground wave and the ionospheric waves generates a compli- cated curve as a function of distance. Comparison of figure 35 with figure 3 6 shows only a change in the daytime curve. In figure 3 6 we have introduced the alternate, more sharply bounded daytime model, figure 8, with the lower electron number density between 40 and 70 km. Although this model does not fit the data shown in figure 30 as well, it is, non the less a physically possible model. This gives us some quantitative intuition concerning the magnitude of perturbations possible at this frequency as the lower D-region changes from high to lower electron -17 density. Figure 37 should be compared with figure 34 since the only- change between these figures is an increase in the ground conductivity from a = 0.005 to o = 0. 01 mhos/m. Similarily, figure 38 should be compared with figure 36 (cp a = 90°, alternate daytime noon ionosphere model fig. 8). Figure 28 should be compared with figure 35. Again, if the conductivity is lowered from = 0. 01 to o = 0. 001 mhos/m, compare figure 37 with figure 40 and figure 3 6 with figures 43 (o" = 0. 005 to o = 0. 001). It is apparent that one cannot guarantee that a particular observer may not experience as much as 1 5 dB decrease or increase in signal at a particular locality as a result of ionosphere temporal changes. Thus, if a disturbance of the ionosphere occurs in the region R in figure 12, the first ionospheric wave, E r>1 , will be affected. Also (figure 14) the third ionospheric wave, E rj2 , will also be affected. It is then probable, if the disturbance is local that the second ionospheric wave, E rj2 , (figure l3)and the fourth ionospheric wave, E rj4 , will be unaffected. This disturbance may be a change in the profile or a lower- ing of the reflecting region. Then, wave interference between the dis- turbed and the undisturbed ionospheric waves could also produce fading similar to that described above at 60 kHz in which the entire propagation path was disturbed. The computer program described in section 7 will accommodate separate reflection coefficients for each of the reflecting regions, [j,k]. Also, a large number of such disturbances Lj»k] can be accommodated depending upon the dimension size of the j,k matrices. Under such circumstances, the principal interest is amplitude and phase as a function of time rather than distance. Numerical studies of this type are reserved for future work. -18- 3. THEORY OF n-WAVES The method of summation of the series (2) although surprisingly- simple tends to mask the physical significance of the n-wave s. This is a consequence of the fact that the factor, G(8,n), contains the Legendre function, P n (cos8), which is a real number with alternating sign. It would appear therefore that each n-wave propagates only in the radial direction. It would also appear that each n-wave is a standing wave in the direction of increasing distance, d = a8. However such a standing wave must be composed of traveling waves in the 8 -direction super- imposed upon each other, but traveling in different if not opposite directions. Using results obtained by Hob son [1896] and Barnes [1907], P n (-cos 9) = P n (cos8) exp (T inn) Q n (cos 8) sinnn . (17) This can be written for integral n only as P n (cos8) = (-l) ft P n (-cos8) = (-l) R P n [cos(tt -8)]. (18) It is found from this that P n (cos8) ~ ( " 1)n {(cosa + isin^) + (cosQ -isinQ)} (19) V 2rr(n+l) sin 8 where If,* TT n = (n+ J) (tt-9) -ii . (20) d = g e , (2i) *See section 6 for definitions of symbols. -19- the lateral displacement of the wave along the ionosphere, 1 and n+f s — k og K c (22) (23) then, using (19), P fi (cos0)~ V i(-l)M-l) ^_ ( exp (_ ikoDs) + iexp(ik Ds)}, (24) where s] 2TT(n+l) sin9 D 7 = D - 2ng , the negative or reverse net lateral displacement (distance) at the ionosphere around the world the long way. Hence, (2) can be written as _ _ M-o c Io*> V n(n+l)(2n+l) cfi Cw (1 + RJ 3 j r E - 3 -"sTT F^ ^ / 7TT f +n • 5 Po ' -i ^=0 V 2n (n+1) sin 9 i (-1)* +1 [exp( -ik Ds) + iexp(ik Ds) } . (25) The factor involving the exponential functions is the sum of two plane waves for either outgoing £i a or ingoing £i a type waves or rays tan> gent to the harmonic sphere of radius r n . Actually two sets of ortho- gonal rays tangent to a harmonic sphere, r fi , exist. Thus, in the meridian plane of the terrestrial sphere, B = £±1 = %. , (26) or n+ 2 _ k g g 9 n+ 2 k • r a = -^ . (27) -20- According to (10), s « sincpj . Thus, in general, there exist two ,(1,2) f P a (-cos9) iQ fl ( - cos 9) ^, orthogonal waves C, I ± — w_j.)» or four type waves (corresponding to superscripts l,2,and+, -signs) if Q n is the second solution ( Abramowitz and Stegun, 1964) of Legendre' s differen- tial equation. As n increases, r n becomes greater than a + h, cp x and sincp i become complex, and Re sincpj becomes greater than unity, i.e., the waves become inhomogeneous. It is interesting to note that (25) can be written as 00 00 E rj j = A Y B(n) exp (-ik Ds) + iA V B(n) exp(ik Ds) . (28) n=° ft=o More generally (25) or (2) can be written as F - C S CM l 1\» T P ^- C ° S9 ) * Qn(-CQS9) "1 E r , I - C I G(n) (-1)- |_ —^ i - j + C^G(n)(-l)» [ P -'- c 2 ° s9 ' +i Q,(-coB6) "j _ (29) a =o The quantities A, B(n), C, and G(n) are defined by reference to (25) and (2). In equations (25) and (29), the second terms, involving a summation over n, represent the wave traveling around the terres- trial sphere in the reverse (-9) direction. This field, is, of course, negligible, except near the antipode of the transmitter at low frequencies. Equation (2) is therefore valid near the antipode and should produce a standing wave as a consequence of the wave in the reverse direction. •21- At short distances, as a consequence primarily of the great lateral dis- placement of the n-wave depicted in figure 11 and the resultant severe attenuation, the second summation term can usually be neglected. Multiple forward (+8) and reverse (-9) waves can be readily demonstrated with the aid of Hobson' s ( 1896) convergent series for TT 5TT T < < — r ; here, 6 5 r(n) cos{(n+|)9-^ } i cos l(n+f)8-lg } P » (cose) = n F(ZTi) L - . " 1 + 2 • 2n+ 3 ,. . „,3/2 (2 sinB) (2 sin 8) 1 cosl(n + |e -gj ^ + 2- 4- 2n+ 3- 2n+ 5 ,„ . ..3/2 + # "J (30) (2 sin 9) and «, «■ r- r(,)r -" nt(n+i)6 -i 1 ! ^t(n + |)e-^} Q n (cosB)=Vn __L___ — z^teiT ., . „,3/2 2' (2sm8) 2 (2sm8) sin {(n+ |)8 -^} - - ••' , (31) 2-4- 2n+3- 2n+ 5 7, . ..3/2 J (2 sm8) which can be combined as , / ax rw ax r/ \ ' - exp [ Ti(n+i)(TT-9) ± i j ] ^(-cosS) . Q R (-cos8) _ ( . a r(n) j 4 2 X tt l "*' Hn+i) I r(n+ ^ L (2TTsin0)^ j exp[ =F i(n+ |)(TT-9) ± i~^] + 2- 2n+3 ~ . 173/2 (2n sm9) expl T i(n + |)(TT-9) ±~ + _ - — + 2 • 4- 2n + 3 • 2n + 5 .„ . AX 5/2 (2n sin 9) -22- }. '32) where T(n) 1 1 It is apparent that the equation (32) splits the series (29) into waves traveling in the positive ( + 9) or negative (-0) direction; the +iQ fl refers to waves in the -0 direction and the - iQ n refers to waves in the + 6 direction. The higher order terms of the series (32) can usually be neglected if n is large or n~ kia; the remaining terms comprise the important terms for the radio wave propagation problem under discussion. Equation (32) can be written, P n (-cos9) . Q n (-cos9) = 2 a TT ( " 1} HnTTT l r , . fi -,J exp L " lk ° Sl D J v ' [2n sin9 J ^ + i r D "I + exp I -ik s 2 j (2.2n+3)[2rrsin9] J/ ^ + i 2 5/2 (2. 4. 2n+3. 2n+5) [2n sin9] exp ik s 3 £ j + ••• } , (33) -23 where , n+ z -I S2 = — i k g S3 = — i Thus, the value given for s approximately in (22) is a real number. The angle of incidence of a single plane wave can be replaced more precisely with a sum of plane waves with angles of incidence si , s 2 , s 3 • • • . In short, s = sincpj is a complex angle of incidence. This is consistent with the general theory complex angle for the n-wave (10). 4. CONCLUSIONS Flexibility, computation accuracy and precision control and a quite general propagation theory for VLF, LF and MF terrestrial radio wave can be achieved with the use of spherical wave functions of integer order -when such functions are used as building blocks for computer programs. The analysis technique suggests a wide variety of scientific applications to the study of VLF, LF and MF ionospheric wave propagation, only a very few of which have been illustrated in this paper. It has been noted that wave interference and terrestrial magnetic field phenomenon and profile electron density changes can affect the amplitude of the 60 kHz signal by as much as 15 db at a fixed distance, ■24- especially at night. It has also been found that the 26 kHz signal is comparatively insensitive in the daytime to the variations in the magnetic dip angle when propagation occurs across the magnetic equator. Thus, the value of the magnetic field vector at the center of the path is a good approximation. -25- 5. REFERENCES Abramowitz, M. , and I. A. Stegun (1964), "Handbook of mathematical functions", National Bureau of Standards, AMS55, (U„ S. Government Printing Office, Washington, D. C. 20402). Bain, W. C. , and B. R„ May (19 67), D -region electron-density distributions from propagation data, Proc. Inst. El. Eng. , 114, No. 11, 1593-1597 Barnes, E. W. (1907), On generalized Legendre functions, Quart. J. Math. Q2/36, 179-184. Belrose, J. S. , and B. Segal (1967), Some comments on the deter- mination of D -region electron density distributions from the reflection of long waves, Ground Based Radio Wave Propagation Studies of the Lower Ionosphere, Conference Proceedings 11-15 April, 1966; Ottawa, Canada, (Defense Research Tele- communications Establishment, Radio Physics Laboratory, Shirley Bay, Ottawa, Ontario, Canada). Belrose, J. S. , (19 68), Low and very low frequency radio wave propagation, RADIOWAVE PROPAGATION, AGARD Lecture 1 Series XXDC, IV-1 to 115 (Technical Editing and Reproduction Ltd. Hanford House, 7-9 Charlotte St., London, W. 1.). Berry, L. A. , (1964), Wave hop theory of long distance propagation of LF radio waves, Radio Sci. J. Res. NBS 68D, 1275-1289. Berry, L. A., and M. E. Chrisman (1965), The path integrals of LF/VLF wave hop theory, Radio Sci. J. Res. NBS, '69D, 1469-1480. Berry, L. A., G. Gonzales and J. Lloyd (1969), The wave hop series for an anisotropic ionosphere, RADIO SCIENCE, Vol. 4, No. 11, 1025-1027 -26- Brekhovskikh, Leonid M. (i960), "Waves in layered media", (Academic Press, New York, London). Budden, K„ G. , (1961), Radio Waves In The Ionosphere (Cambridge at the University Press). Goldstein, M. , and R. M. Thaler (1959), Recurrence techniques for the calculation of Bessel functions, MTAC XIII (13), 102-108. Hirsh, N. M. (1967), VLF studies in simulated sub-D region (G. C. Dewey Corp. , 331 East 38th St. , New York, N. Y.). Hobson, E. W. (1896), On a type of spherical harmonics of unrestricted degree, order and argument, Phil. Trans. Roy. Soc. 187 (A), 443-531. Johler, J. R. , W. J. Kellar and L. C. Walters (1956), Phase of the low radiofrequency ground wave, NBS Circular 573 (U. S. Gov. Print. Off., Washington, D. C. 20402). Johler, J. Ralph (1961), On the analysis of LF ionospheric radio propagation phenomena, Jour, of Res. of NBS, D (Radio Prop.) 65D, No. 5 (Sept-Oct 1961), 507-528. Johler, J. Ralph (1966), Zonal harmonics in low frequency terrestrial radio wave propagation, NBS Tech. Note 335 (U. S. Govern- ment Print. Off., Washington, D. C. 20402). Johler, J. R. , (1967), Theory of propagation of low frequency terres- trial radio waves -- mathematical techniques for the inter- pretation of D -region propagation studies, Ground Based Radio Wave Propagation Studies of the Lower Ionosphere, Conference Proceedings 11 -15 April, 1966; Ottawa, Canada (Defense Research Telecommunications Establishment, Radio Physics Laboratory, Shirley Bay, Ottawa, Ontario, Canada). Johler, J. R. , (19 69a), Ground wave propagation in a normal and an ionized atmosphere, ESSA Technical Rept. ERL 121 -ITS 85 (U. S. Gov. Print. Off., Washington, D. C. 20402). -27- Johler, J. R. , (1969b), Loran C,D phase corrections over irregular, inhomogeneous terrain -- simplified computer methods, ESSA Tech. Rept. ERL 116-ITS 83 (U. S. Government Print. Off., Washington, D. C. 20402). Johler, J. R. , (1970), Spherical wave theory for MF, LF and VLF propagation (to be published). Johler, J. R. and L. A. Berry (1962), Propagation of terrestrial radio waves of long wavelength --theory of zonal harmonics with improved summation techniques, Jour, of Res. of NBS 66D (Radio Prop.) No. 6, 737-773. Johler, J. R. , and L. A. Berry (1965), On the effect of heavy ions on LF propagation, with special reference to a nuclear environ- ment, NBS Tech Note 313 (U. S. Government Print. Office, Washington, D. C. 20402). Johler, J. R. , and J. D. Harper, Jr. (1962), Reflection and trans- mission of radio waves at a continuously stratified plasma with arbitrary magnetic induction, Jour, of Res. of NBS, D (Radio Prop.) 66, No. 1 (Jan-Feb), 81-99. Johler, J. Ralph and Lewis, R. , (1969), Extra low-frequency terrestrial radio -wave field calculations with the zonal harmonics series, Jour, of Geo. Res. Vol. 74, No. 10 May 15, 1969, 2459-2469. Johler, J. R. , and J. C. Morgenstern (19 65), Propagation of the ground wave electromagnetic signal, with particular reference to a pulse of nuclear origin, Proc. IEEE, 53 , No. 12. Morgan, G. I. (1966), Measured daytime field intensities in the United States at VLF, LF and MF (available from Defense Documenta- tion Center, Washington, D. C. 20402). -28- 6. GLOSSARY E r = vertical electric field in spherical system (r,9,cp) r = radial direction. E r = ground wave field. E r . = ionospheric wave field, j = 1, 2, 3 • • • B = constant. B "= \ig !£-, [± Q = permeability of space, p. = 4tt(10" 7 ) Henry/ m. I Q i> - source dipole current moment, ampere -meters. k - = wave number of air, (U ID k , = -Tix ~ - = k (Tix ~ 1.0001 to 1.0003). - 1 c c c = speed of light, c - 2. 997925(10 8 )m/ s f = frequency, Hertz, f = W)/ 2tt . a = radius of the earth, a ~ 6. 3 6739(1 s )m. k_ 2 = wave number of the ground, i m • ° -2 c N e uo e 2 = dielectric constant of the ground. 0" = ground conductivity, mhos/m e = permittivity of space, e ~ 8. 8541 9(10" 12 ) G(0,n) = n(n+ l)(2n + 1) P n (cos9) d = distance along a great circle, d = a6 (n) = the set of integers, n = 0, 1, 2, 3 • • • . P n (z) = the Legendre function, which obeys the recurrence relation: P M-i(*) = -^71 zP n (z) -JL- P B .i(z) , -29- where P,(a) = 1 Pi(z) = z. tn k_ R e = - la k_ ^ln' _2a Po r (l) r r ( 2 ) r ( s ) J L r (i) J C^a = d X) (k_ x a) £ ( -h =ci a) (k- x a) tog - ^» (k g) -2 a ili a (k. 2 a) ~ I C £> a Ci 1,s) (z)=7=f H 2 ) ^n C-ia = logarithmic derivative of £-i a b-l a (i,2)'_ a (1,2) r a r (i,2), . i c - ia - ai c - ia - L~ z c * (z) J, N + = positive ion number density. N_ = negative ion number density. N e = electron number density. m+ = mass of the ion. m e = mass of the electron. e = electric charge. V = collision frequency, sec" . In' \|i-i a - t-=^- In ^- 2a R„ = --In C-i a + ; ln Y-s Cj = effective reflection coefficient defined by the -31- matrix equation: where t= r Tee Te,a "] , T T -"-Bis mm G _ r R 9 on r 1 On 9 " L -1 J ' u * '" Po L -R. J T eo = reflection coefficient (Johler and Harper, 1962) for vertical electric polarization of the incident and reflected ionosphere wave, T mm = reflection coefficient (Johler and Harper, 1962) for vertical magnetic polarization of the incident and reflected ionosphere wave, T 0m = conversion coefficient (Johler and Harper, 1962) from vertical electric incident to vertical magnetic reflected ionosphere wave, T m8 = conversion coefficient (Johler and Harper, 1962) from vertical magnetic incident to vertical electric reflected ionosphere wave, cp t = angle of incidence on the ionosphere, sincp, =[-i C fI! (k '*> l, Ci' (Kg) g '= radius to lower ionosphere boundary, g = a + h, h = height of lower ionosphere boundary above the ground, P r = radiated power, P r = 0.4(10- 12 )uo 2 (l l) 2 /z , •32 cp a = magnetic azimuth, measured clock-wise from magnetic north, degrees or radians, |3i| = magnitude of the terrestrial magnetic field vector W , I = magnetic dip or inclination, measured downward from the horizontal. , 1 n+ 2 s k g r n+ 2 1 n k k„ r>-0 c d = g e Q n (z) = the second solution of Legendre* s differential equation (Abramowitz and Stegun, 19 64). 33- CO 3 co CM CO CM II X CQ w o nj • H £ o M •H bO •rH ft CJ CQ a) o r— 1 CQ ' >N CO Rj rt s W CM CtJ rH M CVJ •» GO CJ i-t • H -^ II •rH Q O X bo •H -34.- t~\ — r o IT O _ E hi — 1 J£ h- II II ^ ,, ZH <3 UJ X UJ 3 (/) o _L_I I L ro CM II +-> i—i II 0) •i-i J •* CU N CO cu > 2 3 flj tuO o < 0) I— 1 > cr o II UJ u > VH lll cu 5 o ft CD f ) 6 _J CO CO •H ^ co UJ o 0) TJ 3 < 1- 3 un (S) <+-i a II o Q • 1-1 O TJ o • - 5 -t i—i I i-H m II ft rrt fi CO o <; rt H CD U OX) •H 83i3w/sinoA ni aanindwv -35- i r FREQUENCY = 100 kHz Tee = -' H = 65 km TYPICAL LAND LU UJ o > UJ Q => CL < DISTANCE IN KILOMETERS Figure 3. Amplitude of the ionospheric waves, j = 1, 2, 3 • • • , and the ground wave, j = 0, as a function of distance; waveguide height = 65 km; T e0 = -1, 0" = 0. 005 mhos/m; f = 10 kHz. -36- TEE = - I f = 26.1 kHz TYPICAL LAND 1,000 2,000 3,000 4,000 5,000 6,000 7,000 OISTANCE IN KILOMETERS 8,000 9,000 10,000 Figure 4. Amplitude of the ionospheric waves, j = 1, 2, 3 ••• and the ground wave, j = 0, as a function of distance; waveguide height = 65 km; T ee = -1; a - 0. 005mhos/ m; e 2 =15; f = 26 kHz. -37- Oli _D jo ul 6o| E a> X o> • H CO U so ej ^ •■- | Pi o a. to o u 0) ^ a) rH CXI I o I- (0 u CO rt ft d v 2 cj « .2 O t3 C3 §3 if} a> i-t r-H ■)-> ' ? § C n ^r o .?« w d 2 ft CO o o a ,d *l d g d u d o ■H CQ P o d o o -p OI>) 4 1H9I3H -42- .s i-H a, o o Ed u +J to 2 r4 d •* fl o ■rH CO ti ri a • X iJ H — 3 O bJD -— Q) h > o tf ^-^ 1> s o •n X CD CO CD • 1-1 u CD s o bO CM II u CD +-> (0 Cti -t-» CQ 0) •H o CO i— < CD a o o CO CO ci ii u CD CO (1) CD CO CO o CO o CO o S Q IT) CO bfi -48- i!ii — i — i r aaiaw/snoA ni aannidwv ■49- I II I — I 1 1 TT CD f£ o| r S o g » I- o uj cc 2 Q- o H w < i- o X 1I5 1 1 1 i i i i i_ ri +j bo t3 ri U •r* co I 8 ft-s CO S3 s ° o m CO to o o S ri co O ,3 LTl O CD O rS O N vO f-l ri a =270° 0.0 500.0 1000.0 1500.0 2000.0 D1ST. KM 2500.0 3000.0 Figure 30. Theoretical amplitude of the total field as a function of distance for the NBS Ft. Collins, Colorado transmitter (normalized to 13 kw radiated power) using the normal daytime noon profile model, figure 7; £ = 60 kHz a = 0. 005 mhos/m; e 2 = 1 5; cp a = 90° ; and 270°; 1=70°; |tt|=40A/m. Data (Morgan, 1966) measured on a radial from the transmitter toward Palm Beach, Florida are also shown. -63- 100 CD CZi I a. 60- 30+ 20 J - 10«- Datum Point 0.0 500.0 1000.0 1500.0 2000.0 DIST. KM 2500.0 3000.0 Figure 31. Theoretical amplitude of the total field as a function of distance for the NBS Ft. Collins, Colorado transmitter (normalized to 13 kw radiated power) using the normal daytime noon profile model, figure 7; f = 60 kHz; a = 0. 005 mhos/m; e 2 = 1 5; cp a = 90 u and 270* 1= 70°; |tf| = 40 A /m. Data (Morgan, 1966) measured on a radial from the transmitter toward Cape Fear, N. C. are also shown. 64- 100?- GQ ? 50+ 40- 30- 20" 10- Theoretica <£ a = 270° 0.0 500.0 1000.0 Figure 32. 1500.0 2000.0 D!ST. KM 2500.0 3000.0 Theoretical amplitude of the total field as a function of distance for the NBS Ft. Collins, Colorado transmitter (normalized to 13 kw radiated power) using the normal daytime noon profile model, figure 7; f = 60 kHz; a = 0.005 mhos/m; e 2 = 15; cp a = 90° and 270 s I = 70°; \U\ = 40 A /m. Data (Morgan, 1966) measured on a radial from the transmitter toward Brownsville, Texas are also shown. -65- 100 60- en ? 50 a. z: 40- 30- 20'- 10- Datum Point 0.0 500.0 1000.0 1500.0 2000.0 D1ST. KM 2500.0 3000.0 Figure 33. Theoretical amplitude of the total field as a function of distance for the NBS Ft. Collins, Colorado transmitter (normalized to 13 kw radiated power) using the normal daytime noon profile model, figure 7; f = 60 kHz; a = 0.005 mhos/m; e 2 = 15; cp a = 90° and 270" I = 70°; |tt| = 40 A/ m. Data (Morgan, 1966) -measured on a radial from the transmitter toward Nantucket, R. I. are also shown. -66- -100 — 110- -200 time: DAY ? NIGHT 0.0 500.0 1000.0 1500.0 2000.0 DIST. KM 2500.0 3000.0 Figure 34. Theoretical amplitude of the total field as a function of distance using the normal daytime noon profile model figure 7; and night- time model, figure 9; f = 60 kHz; a = 0. 005 mhos/m; e 2 =15; cp a = 0* ; I = 70°; \U\ = 40 A/m. -67 -100 -110- TIME $ DAY 7 NIGHT 500.0 1000.0 1500.0 2000.0 2500.0 3000.0 DIST. KM Figure 3 5. Theoretical amplitude of the total field as a function of distance using the normal daytime noon profile model figure 7; and the nighttime model, figure 9; f = 60 kHz; a = 0.005 mhos/m; e 2 = 15; cp a = 90°; I = 70°; \u\ = 40 A/ m. ■68- -100 -110" -190- -200 TIME $ DAY V NIGHT 0.0 500.0 1000.0 1500.0 2000.0 DIST. KM 2500.0 3000.0 Figure 36. Theoretical amplitude of the total field, as a function of distance using the alternate normal daytime noon profile model, figures 8, 9 and the nighttime model, figure*?; f = 60 kHz; O = 0. 005 mhos/m; e 2 = 15; cp a = 90°; I = 70° ; \U\ = 40 A/m. -69- -100 TIME DAY 7 NIGHT 500.0 1000.0 1500.0 2000.0 2500.0 3000.0 DIST. KM Figure 37. Theoretical amplitude of the total field, as a function of distance using the normal daytime noon profile model, figure 7; and the nighttime model, figure 9; f = 60 kHz; a = 0. 01 mhos/m; e 2 = 15; cp a = 0°; I = 70°; \u\ = 40 A/ m. -70- -100 - 1 1 0- - 'J .'J TIME $ DAY 7 NIGHT 500.0 1000.0 1500.0 2000.0 2500.0 DIST. KM 3000.0 Figure 38. Theoretical amplitude of the total field as a function of distance using the alternate normal daytime noon profile model, figures 8, 9 and the nighttime model, figure 9; f = 60 kHz; a = 0. 01 mhos / m; e 2 = 1 5; cp a = 0° ; I = 70°; |w| = 40 A/m. -71- -100 -110- -190- -200 0.0 500.0 1000.0 1500.0 2000.0 2500.0 3000.0 DIST. KM Figure 39. Theoretical amplitude of the total field as a function of distance using the normal daytime noon profile model, figure 7 and the nighttime model, figure 9; f = 60 kHz; 0" = 0. 01 mhos/m; e 2 = 15; cp a = 90°; I = 70°; |m| = 40 A/m. TIME $ DAY 7 NIGHT -72- -110 -120<- -130-- -lao- en ^-150+ a. r: ■160" -170- -180- -190- -200 time: $ DAY 7 NIGHT 0.0 500.0 1000.0 1500.0 2000.0 2500.0 3000.0 D1ST. KM Figure 40. Theoretical amplitude of the total field as a function of distance using the normal daytime noon profile model, figure 7 and the nighttime model, figure 9; f = 60 kHz; a = 0.001 mhos/m; e 2 =15; cp a = 0° ; 1 = 70°; \ W | = 40 A / m. -73- -100 - 1 10 time: £ DAY 7 NIGHT 500.0 1000.0 1500.0 2000.0 2500.0 DIST. KM 3000.0 Figure 41. Theoretical amplitude of the total field as a function of distance using the alternate normal daytime noon profile model, figures 8, 9 and the nighttime model, figure 9; f = 60 kHz; a = 0. 001 mhos/m; e 2 =15; cp t = 0° ; I = 70°; \u\ = 40 A/m. -74- -100 0.0 500.0 1000.0 1500.0 2000.0 2500.0 5000.0 PIST. KM Figure 42. Theoretical amplitude of the total field as a function of distance using the normal daytime noon profile model, figure 7 and the nighttime model, figure 9; f = 60 kHz; o = 0.001 mhos/m; e 2 = 15; cp a = 90°; I = 70°; |g| = 40 A/m. TIME $ PAY V NIGHT -75- -100 -110- -200 TIME $ DAY 7 NIGHT 0.0 500.0 1000.0 1500.0 2000.0 DIST. KM 2500.0 3000.0 Figure 43. Theoretical amplitude of the total field as a function of distance using the alternate normal daytime noon profile model, figures 8, 9 and the nighttime model, figure 9; f = 60 kHz; a =0.0 a = 0.001 mhos/m; e 2 = 1 5; cp a = 90°; I = 70°; \U\ = 40 A/m. -76- 7. COMPUTER PROGRAMS Computer programs RAYGUIDE and RAGIDOO are shown in the computation plan, figure 44. A program DIONS necessary for the cal- culation of reflection coefficients from profiles and SORTDION, which smooths the phase of the reflection coefficients and references the phase to the reflection height are also presented. Thus, DIONS generates data cards for SORTDION, and SORTDION generates data cards for RAYGUIDE or RAGIDOO. RAYGUIDE calculates four ionospheric waves or hops simultaneously and generates data cards for RAGIDOO, which can add, if necessary, a large number of additional hops or ionospheric waves. RAGIDOO can start with the zero order, j = 0, and calculate hops sequentially. But it is not as fast as RAYGUIDE, which calculates four hops simultaneously. For the majority of applications the first four hops suffice. The computer program SORTDION was written by Wieder* (19 67) and the program DIONS was written by Berry** (19 65). The latter was discussed in some detail recently by Johler (1967). The program RAYGUIDE was written by the authors. RAGIDOO is a modification of RAYGUIDE to calculate a large number of terms of the j -series, which was originally written by Wilson*** (19 67) and modified by the authors. The Legendre function, P n (z) is calculated exactly in RAYGUIDE for integral n, from the recursion formula, P n (z) = Un ; 1)z P B -i(z) -^ P n - 2 (z) , (33) PoU) = 1 (34) n n where * B. Wieder (1967), Private Communication, ITS, ESSA, Boulder, Colo. ** L. A. Berry (1965), Private Communication, ITS, ESSA, Boulder, Colo. ***J. H. Wilson (1967), Private Communication, ITS, ESSA, Boulder, Colo. -77- and Pi(z) = z . (35) Both integer and non-integer values of n can be used in a formula based on the saddle point approximation of the integral of LaPlace (letting z = cos 8) as 2 cos f(n+J) 6 -in P n (cos8) ~ — . (36) J 2n (n+ 1) sin9 This formula becomes inaccurate for sin 8 ~ (short distances, d = a 8) and n small. For n large the alternative (3 6) is available in the computer program. The spherical wave functions are also calculated from the recursion formulas, r (i,3) . 2n + 1 (12).. r (i,2) . . U + i (z) = — - — £ n (z) -C B -i (z) , (37) z starting with the values, and d 1,3) (z) = *iexp[ ±iz] (38) G ( _V 2) U) = exp [ ± iz ] . (39) Using a computation technique of Goldstein and Thayler [1959], the computer calculates the spherical wave functions (37), (38), and (39), using subroutines RHOSBN and ZETAIM. The ratio of the derivative of a spherical wave function to the function can also be determined by- recursion formulas ( Johler and Berry, 1962), which yield, ■78- c?'">'w . 1 n d a ' 2) (z) n z Cn-l (z) Gr-1 (z) z (40) where r (i, 2 ) y , v r (l, 2 ) / X Co (z) For large complex argument, z = k 2 a, the Debye approximation can be employed. This yields ^ d l} (he a) _ / n(n + 1) d l} (kea) ,/ (k 2 a) 2 - 1 • (41) (Im k 2 a < 0) This very simple approximation is quite excellent for n(n + 1) ^ (k 2 a) 2 But the series of zonal harmonics converges to great accuracy before n(n + 1) grows as large as |(k 2 a) 2 | . In fact, the series converges rapidly just beyond n~kia. Alternately, the function n . 2 . can (1) (1)' * R (ksa) be readily used in (40) by replacing the C n and C fi function with \|i R and ty R functions. Thus, C, (k 2 a) C (k 2 a) d l} (k 2 a) ^ (k 2 a) (42) /TTz where ^) n (z) = / — J R+ l(z) and J B+ jjz) is Bessel 1 s function of order v Z 2 2 n + \ and argument, z (Abramowitz and Stegun, 1964). Also, -79- — ill (z) = \li«'(z). These functions (40, 41 and 42) are either calculated dz T ft \ fa directly in RAYGUIDE or function PSIR is used. The latter function uses real argument only. -80- TERRESTRIAL MAGNETIC FIELD DATA AVAILABLE AS COMPUTER PROGRAM FOR [j,k] INPUT ANO CONTROL OATA IONOSPHERE PROFILE OATA FOR EACH REFLECTING REGION [),k] FIG's 1,12,13,14,15 ii PROVIOE COMPUTATION WITH ELECTRON AND ION COLLISION FREQUENCIES AS A SUBROUTINE PROGRAM DIONS COMPUTE AMPLITUDE AND PHASE OF REFLECTION COEFFICIENT MATRIX AS A FUNCTION OF REAL OR COMPLEX ANGLE OF INCIDENCE , USING i JOHLER ANO HARPER (1962) METHOD o o Q o c 5 ' PROGRAM RAYGUIDE ASSIGN REFLECTION COEFFICIENTS OF PROPER REFLECTING REGION [j,k,] AND CALCULATE EQ's (I) ond (2) USING REFLECTION COEFFICIENT MATRIX (9), FIGURE 10; INTERPOLATES ON SORTDION DATA FOR EACH TERM OF n - SERIES SUBROUTINES PROVIDE RECURSIVE CALCULATION OF SPHERICAL WAVE FUNCTIONS i 1 SUM FOUR IONOSPHERIC WAVES USING EXPLICIT MATRIX FORMULAS PLUS THE GROUND WAVE 4 I 1-0 SUBROUTINE PROVIDES GROUNO WAVE CALCULATION ' r * PROGRAM RAGIDOO CALCULATE SEQUENTIALLY ANT NUMBER OF ADDITIONAL TERM OF j- SERIES CO IF REQUIRED I i-o NUMERICAL SOLU COEFFICIE USING IMPLICIT ION OF REFLECTION NT MATRIX Figure 44. Conceptual plan of Program RAYGUIDE. GPO 858 - 620 II- PENN STATE UNIVERSITY LIBRARIES ADDDD720aDDia