Integral equations formulation of plasmonic transmission lines

In this paper, a comprehensive integral equation formulation of plasmonic transmission lines is presented for the first time. Such lines are made up of a number of metallic strips with arbitrary shapes and dimensions immersed within a stack of planar dielectric or metallic layers. These lines support a number of propagating modes. Each mode has its own phase constant, attenuation constant, and field distribution. The presented integral equation formulation is solved using the Method of Moments (MoM). It provides all the propagation characteristics of the modes. The new formulation is applied to a number of plasmonic transmission lines, such as: single rectangular strip, horizontally coupled strips, vertically coupled strips, triangular strip, and circular strip. The numerical study is performed in the frequency (wavelength) range of 150-450 THz (0.66-2.0 μm). The results of the proposed technique are compared with those obtained using Lumerical mode solution, and CST. Very good agreement has been observed. The main advantage of the MoM is its intrinsic speed for this type of problem compared to general purpose solvers. ©2014 Optical Society of America OCIS codes: (250.5403) Plasmonics; (240.7370) Surface plasmons; (350.5500) Propagation; (230.7370) Waveguides; (060.4510) Optical communications. References and links 1. P. R. West, S. Ishii, G. V. Naik, N. K. Emani, V. M. Shalaev, and A. Boltasseva, “Searching for better plasmonic materials,” Laser Photonics Rev. 4(6), 795–808 (2010). 2. Q. Li and M. Qiu, “Plasmonic wave propagation in silver nanowires: guiding modes or not?” Opt. Express 21(7), 8587–8595 (2013). 3. E. Ozbay, “Plasmonics: merging photonics and electronics at nanoscale dimensions,” Science 311(5758), 189– 193 (2006). 4. G. Veronis, S. E. Kocabas, D. A. B. Miller, and S. Fan, “Modeling of plasmonic waveguide components and networks,” J. Computational Theoretical Nanoscience 6(8), 1808–1826 (2009). 5. A. Hosseini, H. Nejati, and Y. Massoud, “Design of a maximally flat optical low pass filter using plasmonic nanostrip waveguides,” Opt. Express 15(23), 15280–15286 (2007). 6. S. A. Maier, Plasmonics: Fundamentals and Applications (Springer, 2007). 7. H. Wei and H. Xu, “Nanowire-based plasmonic waveguides and devices for integrated nanophotonic circuits,” Nanophotonics 1(2), 155–169 (2012). 8. T. Srivastava and A. Kumar, “Propagation characteristics of channel plasmon polaritons supported by a dielectric filled trench in a real metal,” J. Appl. Phys. 106(4), 043104 (2009). 9. J. Park, K.-Y. Kim, I.-M. Lee, H. Na, S.-Y. Lee, and B. Lee, “Trapping light in plasmonic waveguides,” Opt. Express 18(2), 598–623 (2010). 10. I. E. Hashem, N. H. Rafat, and E. A. Soliman, “Dipole nantennas terminated by traveling wave rectifiers for ambient thermal energy harvesting,” IEEE Trans. Nanotechnol. 13(4), 767–778 (2014). 11. T. Holmgaard and S. I. Bozhevolnyi, “Theoretical analysis of dielectric-loaded surface plasmon-polariton waveguides,” Phys. Rev. B 75(24), 245405 (2007). 12. K. Leosson, T. Nikolajsen, A. Boltasseva, and S. I. Bozhevolnyi, “Long-range surface plasmon polariton nanowire waveguides for device applications,” Opt. Express 14(1), 314–319 (2006). 13. S. I. Bozhevolnyi, “Plasmonic waveguides: challenges and opportunities,” in Optical Fiber Communication Conference and Exposition and the National Fiber Optic Engineers Conference, Los Angeles, CA, USA, 4–8 March 2012. #214024 $15.00 USD Received 12 Jun 2014; revised 26 Jul 2014; accepted 26 Aug 2014; published 9 Sep 2014 (C) 2014 OSA 22 September 2014 | Vol. 22, No. 19 | DOI:10.1364/OE.22.022388 | OPTICS EXPRESS 22388 14. X. Sun, M. Z. Alam, J. S. Aitchison, and M. Mojahedi, “Comparison of confinement and loss of plasmonic waveguides,” IEEE Photonics Conference (IPC), pp. 618–619, Burlingame, CA, USA, 23–27 Sept. 2012. 15. A. R. Camara, P. M. P. Gouvêa, A. C. M. S. Dias, A. M. B. Braga, R. F. Dutra, R. E. de Araujo, and I. C. S. Carvalho, “Dengue immunoassay with an LSPR fiber optic sensor,” Opt. Express 21(22), 27023–27031 (2013). 16. S. Park and S. H. Song, “Polymeric variable optical attenuator based on long range surface plasmon polaritons,” Electron. Lett. 42(7), 402–404 (2006). 17. F. Liu, Y. Rao, X. Tang, R. Wan, Y. Huang, W. Zhang, and J. Peng, “Hybrid three-arm coupler with long range surface plasmon polariton and dielectric waveguides,” Appl. Phys. Lett. 90(24), 241120 (2007). 18. J. Lee, F. Lu, and M. A. Belkin, “Broadly wavelength tunable bandpass filters based on long-range surface plasmon polaritons,” Opt. Lett. 36(19), 3744–3746 (2011). 19. T. Nikolajsen, K. Leosson, and S. I. Bozhevolnyi, “Surface plasmon polariton based modulator and switches operating at telecom wavelengths,” Appl. Phys. Lett. 85(24), 5833 (2004). 20. D. Dai and S. He, “A silicon-based hybrid plasmonic waveguide with a metal cap for a nano-scale light confinement,” Opt. Express 17(19), 16646–16653 (2009). 21. G. Veronis and S. Fan, “Modes of subwavelength plasmonic slot waveguides,” J. Lightwave Technol. 25(9), 2511–2521 (2007). 22. D. F. P. Pile, D. K. Gramotnev, M. Haraguchi, T. Okamoto, and M. Fukui, “Numerical analysis of coupled wedge plasmons in a structure of two metal wedges separated by a gap,” J. Appl. Phys. 100(1), 013101 (2006). 23. T. T. Minh, K. Tanaka, and M. Tanaka, “Complex propagation constants of surface plasmon polariton rectangular waveguide by method of lines,” Opt. Express 16(13), 9378–9390 (2008). 24. H. Giefers, C. Plessl, and J. Forstner, “Accelerating finite difference time domain simulations with reconfigurable dataflow computers,” in Proceedings of 4th International Symposium on Highly-Efficient Accelerators and Reconfigurable Technologies, pp. 33–38, Edinburgh, Scotland, UK, 13–14 June 2013. 25. J. Hoffmann, C. Hafner, P. Leidenberger, J. Hesselbarth, and S. Burger, “Comparison of electromagnetic field solvers for the 3D analysis of plasmonic nano antennas,” Proc. SPIE 7390, 73900J (2009). 26. J. M. Taboada, J. Rivero, F. Obelleiro, M. G. Araújo, and L. Landesa, “Method-of-moments formulation for the analysis of plasmonic nano-optical antennas,” J. Opt. Soc. Am. A 28(7), 1341–1348 (2011). 27. A. M. Kern and O. J. F. Martin, “Modeling near-field properties of plasmonic nanoparticles: a surface integral approach,” Proc. SPIE 7395, 739518, 739518-13 (2009). 28. E. A. Soliman, G. A. E. Vandenbosch, E. Beyne, and R. P. Mertens, “Full-wave analysis of multiconductor multislot planar guiding structures in layered media,” IEEE Trans. Microw. Theory Tech. 51(3), 874–886 (2003). 29. A. Gholipour, R. Faraji-Dana, G. A. E. Vandenbosch, and S. Safavi-Naeini, “Surface impedance modeling of plasmonic circuits at optical communication wavelengths,” J. Lightwave Technol. 31(20), 3315–3322 (2013). 30. R. Araneo, G. Lovat, and P. Burghignoli, “Graphene nanostrip lines: dispersion and attenuation analysis,” IEEE 16th Workshop on Signal and Power Integrity, pp. 75–78, Sorrento, Italy, 13–16 May 2012. 31. M. Vrancken and G. A. E. Vandenbosch, “Semantics of dyadic and mixed potential field representation for 3-D current distributions in planar stratified media,” IEEE Trans. Antenn. Propag. 51(10), 2778–2787 (2003). 32. E. A. Soliman and G. A. E. Vandenbosch, “Green's functions of filament sources embedded in stratified dielectric media,” Prog. Electromagnetics Res. 62, 21–40 (2006). 33. P. B. Johnson and R. W. Christy, “Optical constants of the noble metals,” Phys. Rev. B 6(12), 4370–4379 (1972). 34. E. D. Palik, Handbook of Optical Constants of Solids (Academic Press, 1985). 35. S. D. Conte and C. de Boor, Elementary Numerical Analysis: An algorithmic Approach (McGraw-Hill, 1990). 36. MATLAB, version 7.9.0. Natick, Massachusetts: The MathWorks Inc., 2009. 37. W. C. Chew, Wave and Fields in Inhomogeneous Media (IEEE Press, 1995). 38. Lumerical Solutions, Inc. http://www.lumerical.com/tcad-products/mode/ 39. CST Microwave Studio, 2012. www.cst.com 40. F. Wooten, Optical Properties of Solids (Academic Press, 1972), Chap.3.


Introduction
With the ever increasing demand of high speed data communication, researchers are always seeking the development of miniaturized devices operating at high frequencies. Recent technology has enabled the fabrication of small devices down to the nano-scale. Nevertheless, researchers faced some limitations in designing classical devices operating above the tens of GHz range as a result of power dissipation and RC delays [1]. On the other hand, photonic devices are capable of providing a wide bandwidth. However, their dimensions, being in the order of micrometers, are incompatible with nano-sized electronic devices [1]. Plasmonic devices operating at optical frequencies offer a solution to this problem [2,3]. They are characterized by their small dimensions (sub-wavelength), bridging the size gap between optical and electronic devices and enabling the development of components essential for high capacity photonic and electronic circuits [3]. The operation of plasmonic devices is based on Surface Plasmon Polariton (SPP) phenomena where light interacts with free electrons of metals resulting in the excitation of an electromagnetic wave propagating at the metal/dielectric interface [4]. The excited wave is strongly confined to the interface and exponentially decays away from it with faster decay at the metal side [5]. In fact, the SPP phenomenon occurs in the optical frequency range where metal properties are different from their counterparts in the millimeter-wave range. At such high frequencies, the relative permittivity ε (ω) r of metals is frequency dependent. It consists in general of a negative real part, ε′ , representing the polarization strength, in addition to an imaginary part, ε′′ , representing (metallic) losses. The relative permittivity of metals can be expressed by the following Drude model [6]: where ε ∞ is the epsilon at infinity, ω p is the plasma frequency, and ω c is the collision frequency. Different plasmonic waveguide topologies have been presented in the literature [7][8][9][10][11][12]. Among these structures are metallic nanostrip [7], nano-trench and V-groove waveguides [8], metal-insulator-metal waveguides [9,10], dielectric-loaded SPP waveguides [11], and Long-Range SPP (LRSPP) [12]. Judging the performance of plasmonic waveguides is dependent on two factors: the mode confinement and the propagation distance of the excited mode [13,14]. The mode confinement is roughly defined by the ratio of the field penetrating the metal to that penetrating the dielectric [13]. Higher field penetration inside the metal implies better mode confinement. However, it also implies that the wave suffers from more attenuation due to losses encountered inside metals at optical frequencies. Waveguide structures like nano-strip, nano-trenches, and V-grooves are suitable for non-linear optical and bio-sensing applications where mode confinement and power localization are crucial [13]. On the other hand, LRSPP is a very good candidate for applications like optical sensors, attenuators, couplers, filters, and modulators [15][16][17][18][19] as it supports waves of long propagation length. Hybrid waveguides offer a compromise between mode confinement and long propagation length [20]. Various numerical techniques including the Finite-Difference Frequency-Domain (FDFD) [21], Finite-Difference Time-Domain (FDTD) [22], Finite Element Method (FEM) [11], the Method of Lines (MoL) [23], and the Effective Index Method (EIM) [11,23], were presented in the literature to analyze plasmonic waveguide topologies. In [21], the FDFD method is used to calculate the eigenmodes of various plasmonic transmission lines including the symmetric/asymmetric plasmonic slot, the modified plasmonic slot, and the plasmonic strip waveguide which consists of a metallic strip and metallic substrate. Accurate calculation using the FDFD solver requires a computational domain big enough to ensure that the fields are negligible at the boundaries [21]. In [22], the FDTD method is applied to study the dispersion and dissipation and the fields of asymmetric coupled nano-wedges at various angles. The main disadvantage of the FDTD method is the high need for computational resources, including computer memory and calculation time [24]. This is especially true for plasmonic structures where metals are highly dispersive [25].
In [11], the EIM and the FEM are used and compared for calculating the propagating modes in dielectric-loaded Surface Plasmon Polariton (SPP) waveguides, whereas in [23], the EIM and the MoL are used to calculate the propagating modes of rectangular hollow waveguides. In both papers, it has been demonstrated that the EIM is simple, but in some cases this method fails to give accurate results. For example, in [11], the accuracy of the EIM is limited to waveguides whose propagating modes are far from the cut-off frequency. The FDFD, FDTD, and FEM have the advantage of directly implementing the differential Maxwell equations [26]. However, the main drawback of these techniques is that the discretization is necessary not only for the plasmonic structures themselves, but also for the space surrounding them [27]. This results in an additional number of unknowns leading to huge memory requirements and calculation times.
The MoM offers a solution to this problem since it involves the discretization of the metallic object solely [27], resulting in a tremendous decrease in the number of unknowns. Furthermore, this technique is characterized by its high accuracy and stability [28] due to the fact that an integral formulation is adopted, and consequently most of the output parameters are expressed in an integral form. The MoM was used in [29], where the surface impedance model was presented to solve plasmonic circuits. This technique is suitable for LRSPP waveguides, which are constructed from nm-thick and μm-wide metal strips embedded in a dielectric medium, because only circumference currents are considered. In [30], the MoM has been applied to calculate the propagation characteristics of a nanostrip graphene transmission line. The nanostrip line is constructed from an infinitesimally thin strip of graphene which has a finite width and is mounted on a grounded substrate. Therefore, in this model only one component of the two transverse currents is considered.
In this paper, all currents physically flowing along the entire cross-sectional area of the strips are considered. The strips are assumed to be infinitely extended along one direction. As far as we can see, this method has not been adopted yet in the field of plasmonic transmission lines. Applying this formulation, the attenuation and phase constants of each propagating mode and the corresponding currents on the metal strips can be calculated. The paper is arranged as follows: in Section 2, the problem formulation is introduced. In Section 3, numerical results for different plasmonic transmission line topologies are presented and discussed. The propagating constants of the fundamental modes of the lines under investigation, as well as the current distribution along the metal strips are computed for each mode. Conclusions are drawn in Section 4. Figure 1 shows the topology of the general plasmonic transmission line under investigation. It consists of a number of metal strips embedded within a stack of flat layers. Each strip has an arbitrary cross-sectional shape and dimensions in the xz-plane, while it extends infinitely along the y-axis. Each flat background layer extends infinitely along the xy-plane, with a finite thickness along the z-axis. Each background layer can be a dielectric or a metallic layer. In the spectral domain, the electric field at the ith strip, i E  , due to current along the jth strip, j J  , can be expressed as follows [31]:

Problem formulation
where k x is the spectral counterpart of the spatial variable x, and k y is the unknown propagation constant to be obtained.
appropriate spectral domain Green's functions that link the current of the jth strip to the electric field at the ith strip [32]. The superscript of any Green's function consists of two parts: the first part is orientation of the electric field to be calculated, being either lateral E  , or z-directed z E . The second part is the type and orientation of the source, which is J  , J ′  , z J , or z J ′ indicating lateral current, derivative of lateral current, z-directed current, or derivative of z-directed current, respectively. Equation (2) can be mapped to the spatial domain as follows: respectively. The multiplication in the spectral domain by ( ) is replaced by differentiation with respect to x, ( ) x ∂ ∂ , in the spatial domain. In order to solve Eqs. (3)-(5) using the MoM, each metal strip is subdivided into segments. As shown in Fig. 2, rectangular segments with dimensions ∆ x and ∆ z along the x-and z-axis, respectively, are used. The unknown modal current components along each metal strip are represented by a superposition of known basis functions weighted by unknown coefficients: where N x , N y , and N z are the number of basis functions along x-, y-, and z-direction, respectively. A n , B n , and C n are the unknown weights of the nth basis function along x-, y-, and z-direction, respectively. ( )  It is clear from Eqs. (6)-(10) that J x and J z are represented by "roof-top" basis functions, while "rectangular prism" basis functions are used to represent J y . After meshing the metal strips and expanding the unknown currents along them, the next step is to satisfy the boundary conditions at these strips. Unlike for perfect conductors, an electric field exists inside the volume of plasmonic metal. The relationship between the electric field and the current along the ith metallic strip is governed by the following equation: where ε ri is the relative permittivity of the ith strip which can be obtained using the Drude model of Eq. (1), or using experimental data [33,34]. Substituting Eqs. (6) where the elements of the impedance matrix [ ] Z of Eq. (12), can be expressed as follows: , , , , where m and n are the orders of the test and basis functions, respectively. They also correspond to the row and column indices, respectively, of the associated sub-matrix. The operation , denotes inner product, while [ ] * denotes convolution, such that: where ( )

x e T z f R x e T z f m n i
Equation (12) can be rewritten in a compact form as follows: This equation is known as the characteristic equation, whose roots represent the unknown propagation constants, y k , of the fundamental modes supported by the multi-strip plasmonic transmission line. The solution of the characteristic equation is performed iteratively using the Müller method [35].

Numerical results
This section presents a study of different plasmonic transmission line topologies using the theory presented in the previous section. The topologies include single rectangular strip, horizontally coupled strips, vertically coupled strips, triangular strip, and circular strip, as shown in Fig. 3. All strips are made up of gold material due to its relatively low losses in the visible and near infrared range compared to other materials [1]. The cross-sectional dimensions of the transmission lines under investigation are shown symbolically in Fig. 3. The background medium is assumed to be free-space for the sake of concentrating on the direct coupling mechanism between the metal strips without interference with reflections from dielectric layers. In free-space, the spatial domain Green's functions are expressed as follows [37], where K 0 is the modified Bessel function of the second kind and zero order:

Single rectangular strip
Unlike microwave transmission lines, a single strip plasmonic transmission line can support a propagating mode. This is due to the fact that a single metal/dielectric interface can guide a SPP wave. The single strip line in Fig. 3(a) Comparing the results of the three simulators, it can be noticed that there is a very good agreement between them. Our solver shows better agreement with Lumerical with a maximum difference of 2.8%. It is also clear from Fig. 4 that as the frequency of operation increases, the effective refractive index and the attenuation constant increase. This can be explained from the fact that at optical frequency range, the skin depth δ is expressed as follows [40]: where λ 0 is the free-space wavelength, and Im(n) is the imaginary part of the complex refractive index of the metal, which can be calculated as follows: 2 2 Im( ) 0.5( ε ' ε '' ε ') n = + − (33) For gold, the skin depth versus frequency is plotted in Fig. 5. From this figure, it is clear that the skin depth increases with increasing the frequency of operation. This means more current penetration inside the lossy gold strip leading to higher effective refractive index and more attenuation for the propagating mode.  Using the calculated value of the propagation constant k y = 6.286 × 10 6 − j 2.275 × 10 5 rad/m at λ = 1.55 μm (frequency = 193.55 THz), the unknown modal current distribution along x-, y-, and z-axes is obtained and plotted in Fig. 6. The dominant current component is the longitudinal J y component along the direction of wave propagation, with maximum absolute value at the edges of the strip and minimum absolute value at its center. The J x component is maximum at the perpendicular right-and left-hand side edges of the strip with odd symmetry. On the other hand, the J z component can be considered as the rotated version of the J x component, showing maximum at the perpendicular upper and lower edges with odd symmetry. It is worth mentioning that the transverse currents J x , and J z are in phase, while the longitudinal current J y lags the transverse currents by 90°. Unlike the microwave transmission lines, the currents of this plasmonic line are no longer along the circumference of the strip, but they penetrate inside the surface of the metal strip. This is expected as the imaginary part of the permittivity of gold at optical frequencies is no longer much higher than the real part. In other words, metals at optical frequencies are not totally lossy and they allow field penetration within them to some extent. It is worth mentioning that the presented modal current distribution along the cross-section of the strip is the same as that calculated using CST. This comparison has been performed for all transmission lines investigated in this paper.

Horizontally coupled strips
In this section, a transmission line consisting of two gold strips of width W = 50 nm and height H = 20 nm, are aligned horizontally with a spacing S = 20 nm between them, as shown in Fig. 3(b). This transmission line supports two propagating modes, which are referred to as the even and odd modes. The modal current distributions of these modes at f = 193.55 THz (λ = 1.55 μm), are shown in Figs. 7 and 8, respectively. According to Figs. 7(b) and 8(b), the dominant current component, J y , has even and odd symmetry, respectively. For the even mode, the currents are maximum at the external edges of the two strips. For the odd mode, the maximum current is at the inner edge of each strip. Similar to the single strip case, the longitudinal current J y lags the transverse currents J x and J z by 90° for both modes. Figure 9 shows the propagation characteristics of the even and odd modes of this transmission line versus frequency. Our solver shows a very good agreement with both Lumerical and CST for the two propagating modes. For the odd mode, it is clear that our solver gives results closer to Lumerical with a maximum difference of 4.88%. Both modes have effective refractive index and insertion loss increasing with the increase of the frequency of operation, which is similar to the single strip case. With respect to the even mode, the odd mode has higher n eff , which means higher percentage energy in the metal strips, leading to more insertion losses.

Vertically coupled strips
Now, the two strips with W = 50 nm and H = 20 nm are placed on top of each other with vertical spacing S = 20 nm, as shown in Fig. 3(c). Similar to the previous transmission line, this structure supports two propagating modes, whose modal current distribution at λ = 1.55 μm are plotted in Figs. 10 and 11. The longitudinal current components suggest that these are the even and odd modes, respectively, with 90° phase difference between the transverse and longitudinal currents. The effective refractive index as well as the insertion loss of the two propagating modes are plotted versus frequency in Fig. 12. Good agreement with Lumerical, and CST can be noticed for the two propagating modes, where the difference does not exceed 4.56%. Comparing the results of the horizontally and vertically coupled strips, the insertion loss is higher for the later case. This can be attributed to the fact that metallic surfaces facing each other are wider in the later case, which results in more losses within the metal volume.

Triangular and circular strips
In order to demonstrate the flexibility of the proposed formulation, strips with nonrectangular cross-section are considered. Specifically, the triangular and circular strips of Figs. 3(d) and 3(e) are investigated. The dimensions of these transmission lines are: W = 50 nm, H = 25 nm, and R = 50 nm. For such structures, the cross-section is approximated by a stair-case mesh. The cell size used for both examples has dimensions Δ x = Δ z = 2.5 nm. Figure  13 shows the effective refractive index and the insertion loss of these transmission lines as calculated by our solver, Lumerical, and CST. From the figure, it is clear that the effective refractive index and the insertion loss increase as the frequency increases similar to the previous examples of strips with rectangular cross-section. It is also clear that the agreement between the three solvers is good, which validates the obtained results and proves the versatility of the proposed formulation. The MoM is considered a powerful tool for the calculation of the propagation characteristics of transmission lines. The strength of the MoM arises from the fact that it requires meshing for only the surface of the metallic strips. On the other hand, Lumerical and CST, which are based on the finite difference method, require substantial additional meshing for the space surrounding the strips. This additional space should be made wide enough to ensure that the fields are not strongly disturbed by the truncation boundaries. Consequently, the number of unknowns is way larger than that of the MoM. Table 1 shows a comparison of the number of unknowns used in calculating the propagation characteristics of the single rectangular strip, horizontally coupled strips, vertically coupled strips, triangular strip, and circular strip using the three solvers at 193.55 THz (1.55 μm). It is clear that our solver requires much less number of unknowns compared to Lumerical and CST software, illustrating the intrinsic supremacy of this technique for this type of problem.

Conclusion
In this paper, a generalized formulation for plasmonic transmission lines is presented. The line consists of a number of metallic strips with arbitrary shapes and dimensions. The formulation is based on the integral equations solved using the Method of Moments (MoM). Such formulation is more accurate and faster than other numerical techniques. Moreover, it requires less computation resources. This comes from the fact that the MoM solves for unknown currents along the metallic strips only, instead of solving for field unknowns along the entire volume surrounding the strips. In addition, the MoM assumes open boundary with no need for applying truncation boundary conditions that result in errors. Finally, the fact that the MoM is based on integral equations rather differential equations makes it more stable and less sensitive to numerical noise. The proposed formulation is applied to simple test structures that include single strip of various cross-sections and two coupled strips. The obtained propagation characteristics are compared with those of commercial software and found to be in good agreement.