SciELO - Scientific Electronic Library Online

 
vol.17 issue4Multiband FSS Analysis and Synthesis Based on Parallel Non Coupled Metallic Strips Using WCIP Method author indexsubject indexarticles search
Home Pagealphabetic serial listing  

Services on Demand

Journal

Article

Indicators

Related links

Share


Journal of Microwaves, Optoelectronics and Electromagnetic Applications

On-line version ISSN 2179-1074

J. Microw. Optoelectron. Electromagn. Appl. vol.17 no.4 São Caetano do Sul Oct./Dec. 2018

http://dx.doi.org/10.1590/2179-10742018v17i41219 

Article

Comparison of FDTD and TDWP Methods for Simulating Electromagnetic Wave Propagation above Terrains

Mahmoud Fallah1 

Amir. H Nazeri1 

1School of Electrical Engineering, Iran University of Science and Technology, Tehran, Iran

Abstract

Researchers have always been in face of path loss calculation in different media for applications such as telecommunication link design. Wave propagation calculations in large spaces using the FDTD method is time-consuming and imposes a great computational burden. For this reason, to replace the classical FDTD method for wave propagation simulation and path loss calculation in large spaces, optimized methods, namely TDWP, have been provided. In this paper, the use of the TDWP method for wave propagation simulation and path loss calculation above a terrain is investigated. Longitudinal components of ground-waves are taken into account (including direct waves, ground reflections, and surface waves). Propagation space is longitudinally divided into smaller FDTD windows with finite length. The electromagnetic pulse travels through these windows from left to right to the desired point. But despite its capability in reducing computational burden and increasing processing speed, TDWP has lower precision in instantaneous field simulation and calculation of propagation coefficients, so that results obtained from FDTD and TDWP are clearly different. In this paper, some efficient methods are proposed, which yield an increase in method accuracy.

Index Terms Electromagnetic wave propagation; Groundwaves; Surface wave propagation; Terrain modeling; Consecutive windows; Finite-difference time-domain (FDTD); Time-domain wave propagation (TDWP)

I. INTRODUCTION

Modern communication and radar systems need effective tools to foresee propagation. In many cellular communication and HF and VHF radar sites, high precision propagation simulation is vital. The results are necessary for choosing proper places for ground stations and optimum height for transmitters and receivers, and also reducing the costs of implementing communication systems. One of the primary aims of propagation simulation is to forecast path loss (LP) with high precision. Path loss is defined as:

Where PT,GT,GR,PR,LP, and LS denote transmitted power, transmitter antenna gain, receiver antenna gain, propagation loss, and system loss, respectively. In (1), antenna gains are measured with respect to an isotropic reference antenna, and LS represents all system losses, including thermal loss, signal loss in cables, mismatch loss, etc. For example, this equation is used in cellular communication system design to foresee each antenna coverage and dead zones before they are installed. Propagation loss denotes all loss caused by factors such as propagation medium refraction index, multipath interference, and reflections from natural and man-made obstacles. So the simulator must take all such parameters into account.

LP=PTGTGRPRLS (1)

In the past two decades, researchers have been in favor of path loss calculation, with numerous efforts made in this direction. In [3], virtual tools for ground-wave propagation modeling and numerical simulation strategies are reviewed. The propagation loss of 71 floors in 17 different buildings, in an urban area, is measured in 4 frequencies ranging from 800 MHz to 8 GHz in [4]. In [5], a study of path loss models in the 2 – 73.5 GHz range in an Urban-Macro Environment is performed. In [6], electromagnetic (EM) wave measurement in 2.4 GHz in underwater environment is conducted. Also, it is shown that the propagation velocity and absorption coefficient in fresh water are frequency independent.

So far, several numerical methods have been introduced for the analysis of electromagnetic problems. Most of these methods are based on full-wave analysis in the time or frequency domain; and utilize one or two difference equations, plus boundary conditions (depending on the specific case) to solve the problem [7 - 10]. Each method is used to solve a specific class of problems. Although these methods are powerful and can be used to solve a wide variety of problems, due to improper usage and intrinsic limitations, sometimes the derived answer is erroneous. Unless errors such as overflow or underflow have occurred, computers always produce answers. The case is whether this answer correctly represents the physics of the problem in the real world.

One of the most ubiquitous methods is the finite-difference time domain, or FDTD method for short. This method is performed in the time domain and has been used for decades. Although real wave propagation takes place in three-dimensional space, a two-dimensional FDTD approximation can be used, if the problem is symmetric relative to one of the coordinate axes.

In this paper, the propagation factors derived using the TDWP and FDTD methods are compared. In order to find differences in the results obtained from the two methods and their causes, spontaneous fields in different time periods might be studied. In the first part, the TDWP method and its application in solving problems are explained. In the next part, a simulation in the whole space is performed using both TWDP and FDTD, and the results are compared. Finally, the simulation process is modified to reduce these differences.

II. TDWP METHOD

The wave propagation problem may be set in city streets with numerous buildings, in a vast desert, or in a mountainous area. The classical FDTD method faces shortcomings when simulating wave propagation along long paths. Suppose we would like to study wave front propagation with spatial Gaussian distribution (along the height axis) in a 10 km path and a height of 50 m. If we divide this space into a mesh of 10 cm cubic cells, there would be 5000000 cells in this space, not considering PML blocks located at the left, right, and top of the computational space. The propagated wave consists of three components, regardless of its polarization. Every time these three components are updated, the new values should be stored in matrices. So we need three matrices with 5 milli on elements each, all with 15 million elements updated in every time step. Considering the fact that the utilized program uses some of the computer's RAM and CPU cycles, it can be concluded that processing such huge amounts of data with personal computers is impossible. Even if sufficient RAM capacity is obtained, processing the data takes considerable time.

To overcome this problem, instead of updating all components for every computational cell at every moment, a rectangular window around the pulse can be presumed, containing most of the propagated pulse energy. In this manner, the computational space can be divided into smaller windows, as shown in fig. 1 and 2.

Fig 1 Dividing natural terrain into smaller computational regions 

Fig 2 Dividing urban terrain into smaller computational regions 

These smaller rectangular windows must be surrounded by PML from the left, right, and top. PML blocks surrounding the computational space prevent reflection and refraction of the waves along the boundaries. The aim is to compute propagation loss along the wave's path and also along the height axis. To this end, we trigger the computational space with a short duration pulse. The trigger source is located at the beginning of the first window. When it begins to stimulate field components in the adjacent cells, the wave front begins to move along the length of the window to its border. When the wave front first reaches the last column of cells in this window, the magnitude of field components in this column are stored in some part of the computer memory, until the effective pulse energy reaches the end of the current window. After the resulting data is stored, the first window is closed and a new computational space is defined as the second window, containing the physical properties of the next compartment. In the new window, all cell field components are initially zero, which means that none of the cells are stimulated. The field components of cells in the last column of the preceding window are immediately copied as field components of cells of the first column of the current window (fig.3). Similar to the first window, the pulse travels to the second window, and the same procedure is repeated. By continuing in the same manner for numerous windows, pulse flight can be simulated for paths as long as tens of km.

Fig 3 Storing and restoring field components in two consequent windows 

In fact, in this manner, a rectangular computational space with very great length is divided into rectangular computational spaces with smaller length, but the same height as the huge rectangle, so that the simulation can be readily done with a personal computer.

III. Problem formulation

A. Wave Propagation Formulation

Although the physics of propagation is taking place in a three-dimensional space, for a symmetrical problem with respect to one of the Cartesian coordinates, we can use a two-dimensional FDTD approximation. Maxwell's equations and the constitutive parameters, in the form of vector differentials, in an isotropic space containing magnetic and dielectric materials, are as follows:

Where B and D are associated to the H and E vectors. σ And σ* are the electric and magnetic conductivity of the material, and ρ and ρ* are the electric and magnetic charge, respectively. By assuming TMy (Hy,Ex,Ez) polarization for the propagated wave and considering a zero value for the electric and magnetic conductivity coefficient, equation 2 and 3 can be simplified as:

×E=Btσ*H×H=Dt+σE.D=ρ.B=ρ* (2)
B=μH=μrμoHD=εE=εrεoE (3)

In Fig. 4, the configuration of the electric and magnetic field nodes are shown for TM polarization. In this figure, fields components within closed dash lines have the same index in the implementation of the algorithm. The nodes of the electric field, with lines to illustrate the direction of field components, are shown with solid circles; and the nodes of the magnetic field are shown with empty squares. Triangle shaped dashed lines represent a group of nodes that have identical array indices in the computer algorithm. For example, in the bottom left of the grid, all three nodes inside the triangle have (i = 1, k = 1) indices. The diagram on the bottom left of the grid shows nodes related to each field and their indices for a typical group. In the diagram on the right side of the grid, the same nodes and their array indices have been shown when used in a computer program. In this polarization, none of the nodes of the fields were located in places with the correct values of discretized space.

Fig 4 Configuration of the nodes of electric and magnetic fields for TM polarization used in simulation 

With the above explanations, and assuming a discrete time form of tn12=(n12)Δt , equations 4, 5 and 6 can be indexed:

μHyt=EzxExz (4)
εEzt=Hyx (5)
εExt=Hyz (6)

The components of the electric and magnetic field are calculated in the Δt, 2Δt, 3Δt, …, nΔt, and 12Δt,32Δt,52Δt,,(n+12)Δt moments, respectively.

If a node of the magnetic field is placed within the setting of the conductor, it is assumed that the entire cell containing that node is placed inside the conductor. Thus, all electric field components of that cell (which include that node of the magnetic field) are fixed to zero.

IV. Path loss calculation

In this problem, the propagation space is a semi-confined space, limited by z = 0 and z → ∞ vertically, and by x → ∞ and x → ∞ in the plane parallel to the ground. We would like to calculate the loss caused by the propagation medium, without considering the effect of the receiver and transmitter antenna. In other words, we would like to calculate the power lost in propagation, after the wave is transmitted by the transmitter antenna, before it is absorbed by the receiving antenna. Therefore, in (10), GT, GR, and Ls should be substituted with 1. In every given point, PR can be calculated in terms of field magnitude [2]:

Hyn+12(i,k)=Hyn12(i,k)+ΔtμΔx(Ezn(i+1,k)Ezn(i,k))ΔtμΔz(Exn(i,k+1)Exn(i,k)) (7)
Ezn(i,k)=Ezn1(i,k)+ΔtεΔx(Hyn12(i,k)Hyn12(i1,k)) (8)
Exn(i,k)=Exn1(i,k)ΔtεΔz(Hyn12(i,k)Hyn12(i,k1)) (9)
PR=E2Z0λ24π (10)

Therefore, to compute the propagation loss, we need to calculate the electric field magnitude at every given point in the propagation space.

To be able to compare the results for different media, we need to eliminate the effect of the reduction of wave power due to the distance from the source via the inverse square law, which is the only loss associated with a lossless medium without obstacles. The propagation factor at every given point is defined as |E/E0| (for any specific frequency). Where E and E0 are the measured and reference electric field at the desired point, respectively [2]. The reference field is the field propagated from the same source, not taking the terrain into account (fig.5). When the wave reaches the boundaries of the propagation space into free space, it must be ensured that it leaves the studied medium and will not be reflected back. To simulate the free space surrounding the propagation space, electromagnetic absorber layers must be used. These layers (fig.5) attenuate the wave, so that it doesn't return to the propagation medium after leaving it. one of the most appropriate types of such absorber layers, for simulations, is the “Perfectly Matched Layer” or PML for short [11]. So in this method, we first perform the simulation for the terrain and store the derived E (t) values for both a specific height and different distances, or also a specific distance and different heights. In this step, the propagation medium is surrounded by the PML from the top, left, and right. In the next step, the simulation is repeated with the terrain removed and the propagation space covered by PML from the bottom as well. In this step, the E0(t) values are stored in the same cells as the last simulation [12].

Fig 5 For every scenario, two simulation steps should be performed. In the first, the fields are calculated in the presence of the terrain. In the second, the fields are calculated without regarding the terrain. 

Finally, after both steps are done, by performing DFT analysis on the time-domain data available (in the desired bandwidth), |E/E0| for each desired frequency is derived. Usually, it would be favorable to calculate and plot this ratio in decibel units (2*10log|E/E0|) for a specific distance or height.

A. Source Stimulation

To stimulate the fields of the computational space, we need to model a physical source. Because finite time distributions, which also have limited frequency bandwidth, provide better responses in FDTD simulations, the pulse shape (in the time domain) is usually selected as Gaussian. However, here we prefer to use the derivative of the Gaussian function, in order to prevent DC effects that can cause shifting of the wave propagation frequency in some applications. However, we should note that the desired pulse for our source should cover the range of frequencies determined by the bandwidth. The Gaussian function derivative used in our algorithm is as follows:

f(t)=2/3α(tφ)e(α(tφ)2)α=2/74BW2 , φ=4α (11)

Where ϕ is the time shift of the function toward the positive direction of the time axis; which is equal to 4/α in our case. The pulse shape of the Gaussian function derivative in the time and frequency domains are given in fig. 6, respectively, for BW = 200 MHz.

Fig 6 a) The diagram of the Gaussian pulse derivative in terms of time for the bandwidth 200MHz, b) The diagram of the Fourier transform of the derivative Gaussian pulse. 

Now which i and j the simulation applies to depend on the type of stimulation. The distribution, e.g. pulse time distribution, can also be in different forms. One of the most popular spatial distributions along the vertical axis is the Gaussian distribution; where its centricity is shown with h _ source.

The statistical process of spatially distributing the pulse along the z axis is achieved by its variance.

G(z)=12πσvexp((zh_source)22σv) (12)

h_source is the height of the spatial distribution center for the source from the ground and σν is the spatial variance of the source. In fig.7, the spatial distribution of the source with centricity of h_source = 25m from the ground at a height of 50 m is shown for different variances.

Fig 7 Spatial distribution of the source with centricity of Z0 = 25 m from the ground at a height of 50 m for different variances 

Thus, the pulse function pulse(y, t) will be the multiplication of the two functions:

pulse(z,t)=G(z)F(t) (13)

In order to compare the FDTD and TDWP methods, a simulation is set up as shown in fig. 8, containing two conducting triangles in free space. The ground plane is also assumed to be a perfect electrical conductor. The dimensions of the simulated space are shown in fig.8.

Fig 8 Dimensions of the simulated space and structures used in simulation. Four different columns at four different distances, where the field components should be stored in order to compute the propagation factors 

The fields used in this simulation are Gaussian derivative functions of time with a central frequency of 50 MHz and a bandwidth of 100 MHz. The spatial distribution is a Gaussian derivative function of height as well. The center of the spatial distribution is taken to be 15 m high. The stimulation source is assumed to consist of a number of dipoles perpendicular to the ground. Therefore, the polarization is TM along the x direction. We would like to obtain the propagation factor as a function of height for four different distances, as shown in fig. 8. We compute the desired factors by two different procedures. In the first procedure, the simulation is performed for the whole space (using FDTD method) as shown in fig. 9 and the results are recorded.

Fig 9 Instantaneous images captured for comparing the FDTD and TDWP methods in the setup shown in fig. 8  

In the second procedure, the space depicted in fig.8 is cut into two halves, or in fact two windows (TDWP method), and the simulation is rerun. To achieve a better understanding of the propagation, the vertical component of the normalized electric field in the computational space is color-plotted. The instantaneous images of the wave front, yielded from the FDTD and TDWP methods, are shown in fig.9.

Propagation factor curves derived from these two methods are shown in fig.10. It can be seen that the two curves do not match perfectly. The reason can be deduced by analyzing the images in fig. 9.

Fig 10 Comparison between propagation losses for lines x1, x2, x3. and x4. obtained from FDTD and TWDP 

As it can be seen in frames 3 and 4 of fig.9 (in the TDWP method), when the propagated wave reaches the end of the first window, it is restored in the beginning of the second window. But, as seen in frames 6 of fig. 9, reflections from the right side triangle (situated in the second window) are neither propagated back into the first window, nor absorbed by the PML located at the left boundary of the second window. Instead, the reflection is re-reflected into the second window.

The reason the reflected wave from the second triangle is not reflected back into the first window is obvious. As consequent computational spaces are discrete, i.e. not continuously connected as in the FDTD method, when the wave front passes through a window for which calculations are over, no field will be regenerated in the previous windows. It is apparent that should we apply the method we use to store fields yielded from previous window and regenerate the field in current window in reverse direction, the result would be a heavy computational burden as FDTD. The waves reflected from the second triangle are absorbed by the PML blocks located at the right and top borders of the window; however, they are reflected by the left border. The reason is that the first column in the window is regenerating fields from the last window during the whole duration of the simulation. In other words, these cells do not take part in the wave propagation that happens in the right to left direction and solely represent the values stored from the previous window. It may seem that when the reflected wave from the second triangle reaches this column, the column's cells are not regenerating any field. But, in fact, the visible portion of the wave generated by these cells is only the tangible part, and these cells continue to regenerate fields, even though the resulting magnitudes are too small to be normally noticed. It is obvious then that as these cells do not take part in wave propagation, they practically act as a thin conducting wall and cause the reflected wave front from the second triangle to be re-reflected back into the second window. This is actually also the case for the cells generating the stimulation. As in this particular example it is possible to predict when the effective portion of the wave has passed through the previous window to the current window, the latter problem can be tackled by limiting the regeneration of field components in the first column. Moreover, considering the fact that the time domain Gaussian derivative function is known, the same trick can be used to limit the duration of the stimulation in the first window. Another trick which can be used for the stimulation is to set some threshold, so that when the source field magnitude falls below this threshold, the cells generating the source fields are set to zero. Fig 10, shows the comparison between the propagation loss results of two method (TDWP vs FDTD) in several lines.

B. Limiting field regeneration in the first column of the second window cells

In fig.11, instantaneous images taken from TDWP simulations with and without regeneration limit for the first column of the second window cells are compared.

Fig 11 Comparison of the instantaneous images taken, with field regeneration in the first column of the second window cells a) not limited b) limited. 

Frames 1 to 6 demonstrate snapshots taken from the 608th through 928th time steps of the simulation, with 64 time steps between each two consecutive frames.

If the propagation factor computed for lines x3 and x4 using the modified simulation method is compared with FDTD, it can be observed that the modified TDWP method provides results closer to the results derived from FDTD [5]. Fig 12, shows the comparison between the propagation loss results of two method, Optimized TDWP and FDTD, in Lines 3 and 4.

Fig 12 Comparison of the propagation loss calculated using FDTD and modified TWDP for lines x3 and x4. with field regeneration duration limited in the second window. 

As it can be sometimes impossible to predict the time when wave reflection from the terrain reaches the end of the previous window, field regeneration duration cannot be limited in general. This problem can be solved by applying the substitution F(x,y)=F(x,y)+f(t) . where F can be any field component Hx, Hy, or Hz, and f(t) represents the stored value of the component from the last stage. It might seem that the aforementioned substitution causes new values to sum up with old values; but it must be borne in mind that field regeneration takes place when the field component values are being updated, so that every time a value is updated, it does not contain the previous value anymore.

C. Analysis of Improving Factors of Optimized TDWP Method

In following some relations are defined which can help us to make a suitable analysis regarding the efficiency of proposed Technics, in order to improve traditional TDWP method.

If we denote the propagation loss factor of Traditional and our Optimized TDWP method as PFTDWP(op) and PFTDWP(re), respectively, the amount of generated Error in terms of y (axis) will be calculated for Traditional and Optimized methods as below, respectively. It should be emphasized that the one-window FDTD method is not applicable for long pathes. It is applied to this study just for a path of 100 m.

δop(y)=|PFTDWP(op)PFFDTD|δre(y)=|PFTDWP(re)PFFDTD| (14)

The ratio of δop(y) to δre(y) is defined as a Improving Factor for each height of y.

cf(y)=δop(y)δre(y) (15)

For investigating the improving factor of proposed method in one given line (like 0 < y < H, x = xo), we should first calculate the below Integrals for each line.

Δop=0Hδop(y)dyΔre=0Hδre(y)dy (16)

The ratio of Δre to Δop is defined as a Improving Factor for each vertical line x.

CF(x)=ΔopΔre (17)

In similar way, for calculating the propagation loss factor along a horizontal line (like 0 < x < L, y = yo) we will have :

δop(x)=|PFTDWP(op)PFFDTD|δre(x)=|PFTDWP(re)PFFDTD| (18)

The ratio of δop(x) to δre(x) can be defined as improving factor for each length of x.

cf(x)=δop(x)δre(x) (19)

For investigating the improving factor of proposed method in one given line (like 0 < y < H, y = yo), we should first calculate the below Integrals for each line.

Δop=0Lδop(x)dxΔre=0Lδre(x)dx (20)

The ratio of Δop to Δre is defined as a Improving Factor for each vertical line x.

CF(y)=ΔopΔre (21)

At last, for attaining to a total improving factor in one window, relation 16 (or 20) should be employed over a given surface like Sδop(x,y)ds . As a case of study, for xo = 52.5m (Line3), the amounts of δορ(x) and δre(x) are 1.26 and 3.55 respectively in 19 m; Thus cf (19) = 0.35. by calculating the improving factor for xo = 52.5m in all heights we will have CF (x = 52.5) = 0.16, While this factor will be CF(x = 70) = 0.05 for xo = 70m (Line4). It shows that this technic make a better operation for Line4 in comparison to Line3.

V. Conclusion

In this paper, the TDWP method for wave propagation simulation and propagation loss computation was introduced. The benefits of this method in comparison to the FDTD method where discussed, and by comparing these two methods for a terrain with two simulation windows, the shortcomings of TDWP methods where demonstrated. After that, by proposing the approach of limiting field regeneration for first column cells of each window, the matching between the propagation factors curves derived from the two methods was ameliorated. The TDWP method enables wave propagation simulation with personal computers; but provides better results for plain surfaces rather than for terrains. Such as the simulation of wave propagation for islands located near coasts where surface impedance changes as a function of distance from the source. In order to achieve more accurate results for surfaces with terrains, modifications should be made, as discussed above.

REFERENCES

1. [1] Sevgi, L., Ponsford, A., Chan, H.C., “An Integrated Maritime Surveillance System Based on High-Frequency Surface-Wave Radars, Part1: Theoretical Background and Numerical Simulations”, IEEE Magazine, Vol. 43, No. 4, Agust, 2001. [ Links ]

2. [2] Akleman, F., Ozyalcin, M.O., Sevgi, L.,“Novel Time Domain Radiowave Propagators for Wireless Communication Systems”, Turk. J. Elec. Engin.(TUBITAK), Vol.10, No. 2, 2002. [ Links ]

3. [3] Sevgi, Levent. “Groundwave modeling and simulation strategies and path loss prediction virtual tools.” IEEE Transactions on Antennas and Propagation Vol.55, No.6, 2007. [ Links ]

4. [4] Okamoto, Hideaki, Koshiro Kitao, and Shinichi Ichitsubo. “Outdoor-to-indoor propagation loss prediction in 800- MHz to 8-GHz band for an urban area.” IEEE Transactions on vehicular Technology Vol.58, No.3, 2009. [ Links ]

5. [5] Thomas, Timothy A., et al. “A prediction study of path loss models from 2-73.5 GHz in an urban-macro environment.” Vehicular Technology Conference (VTC Spring), 2016 IEEE 83rd. IEEE, 2016. [ Links ]

6. [6] Sendra, Sandra, et al. “Underwater wireless communications in freshwater at 2.4 GHz.” IEEE Communications Letters Vol.17, No.9, 2013. [ Links ]

7. [7] Li, Yujia, and J-M. Jin. “A vector dual-primal finite element tearing and interconnecting method for solving 3-D large-scale electromagnetic problems.” IEEE Transactions on Antennas and Propagation Vol.54, No.10, 2006. [ Links ]

8. [8] Taflove, Allen, and Susan C. Hagness. Computational electrodynamics: the finite-difference time-domain method. Artech house, 2005. [ Links ]

9. [9] Peng, Zhen, Kheng-Hwee Lim, and Jin-Fa Lee. “Nonconformal domain decomposition methods for solving large multiscale electromagnetic scattering problems.” Proceedings of the IEEE Vol.101, No.2, 2013. [ Links ]

10. [10] Gibson, Walton C. The method of moments in electromagnetics. CRC press, 2014. [ Links ]

11. [11] Berenger, J. p., “Perfectly Matched layer for the FDTD Solution of Wave-Structure Interaction Problem”, IEEE Transactions on Antennas and Propagation, Vol. 44, No. 1, January 1996. [ Links ]

12. [12] Sevgi, L., Akleman, F., Felsen, L. B., “Groundwave Propagation Modeling: Problem-Matched Analytical Formulations and Direct Numerical Technics”, IEEE Magazine, Vol. 44, No. 1, February, 2002. [ Links ]

Received: January 26, 2018; Revised: February 09, 2018; Accepted: August 16, 2018

Creative Commons License This is an Open Access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.