## Abstract

A fast calculation method for computer generation of spherical holograms in proposed. This method is based on wave propagation defined in spectral domain and in spherical coordinates. The spherical wave spectrum and transfer function were derived from boundary value solutions to the scalar wave equation. It is a spectral propagation formula analogous to angular spectrum formula in cartesian coordinates. A numerical method to evaluate the derived formula is suggested, which uses only *N*(log*N*)^{2} operations for calculations on *N* sampling points. Simulation results are presented to verify the correctness of the proposed method. A spherical hologram for a spherical object was generated and reconstructed successfully using the proposed method.

© 2013 OSA

## 1. Introduction

Holograms are capable of recording and reproducing all three dimensional information like motion parallax, accomodation, occlusion, convergence and so on. For three-dimensional information storage and retrieval(of an object) to be complete, recording and reconstruction is to be done on all sides(directions) of an object, *i.e*, for 360° on both azimuth and polar directions. This is achieved by recording and reconstructing on a spherical surface surrounding the object, which is called as spherical holography. However available optical techniques and numerical methods have restricted holography to be done on plane and cylindrical surfaces only. The restriction due to optical techniques arises from achieving coherent illumination of object and reference wave on a spherical surface. One solution is to model the object on the computer and numerically generate the hologram. But still the following issues prevail, i) The unavailability of fast calculation methods or availability of only approximated fast calculation methods for spherical CGH and ii) the need to illuminate the whole spherical surface using a single reference wave front for optical reconstruction. However the availability of optic fibers and developments in multidimensional lasers [1, 2] that include microball lasers or lasing microspheres have provided hope in realizing spherical holography. Motivated by these appealing facts it was decided to make some developments to numerical methods for spherical holography. Accrodingly, this research work is an attempt to devise a new fast calculation method for computer generation of spherical holograms.

Wave propagation in the paraxial regime is defined by the Rayleigh Sommerfeld diffraction integral [3]. This formula enables us to calculate the optical field(complex amplitude) at an arbitrary point or surface(hologram) knowing the radiating source point or surface(object). This method also known as the direct integration method and is slow due to large sampling points in object and hologram. To speed up computation the Rayleigh Sommerfeld integral is expressed as convolution integral or as angular spectrum of plane waves and then evaluated using fast Fourier transform(FFT). This speeds up the calculation considerably, however demands the object surface and hologram surface to satisfy shift invariance. Accordingly, to generate a plane hologram using FFT, the object should also be a parallel plane surface or set of parallel plane surfaces. If the hologram surface is non planar or not parallel to the object surface then the direct integration method has to be used, which consumes *O*(*N*^{4}) or *O*(*N*^{3}) operations for *N* sampling points. But later on smart ways to use FFT for calculations even when the hologram surface is non planar or parallel were devised. Tommasi *et al.*, [4] developed a fast computation method where the hologram is a tilted plane with respect to the object surface. Sakamoto *et al.*[5] took advantage of the shift invariance in rotation and developed a method to compute cylindrical hologram for a plane object using FFT. Later the method was extended to a volume object by slicing it into planar segments [6]. Yamaguchi *et al.*[7] proposed a fast computation method for cylindrical holograms by segmenting the object and hologram surface into elemental plane and parallel surfaces. Sando *et al.*[8] came up with a brilliant idea by choosing the object to be a concentric cylindrical surface with the hologram surface which makes the system completely shift invariant. They derived the point spread function(PSF) for this configuration and developed the convolution integral to calculate wave progagation from one cylindrical surface(object) to other(hologram).The convolution method was evaluated using FFT, making it a fast computation method. Later Jackin *et al.*[9], devised a new and fast computation method for the same problem by deriving the Transfer Function and defining wave propagation in spectral domain. The advantage of this method is, it uses one less FFT operation for its computation compared to the convolution method. Then the method was extended to a volume object by slicing in to cylindrical surfaces [10].

Fast computation solutions for spherical computer generated hologram emplyoing PSF(convolution method) was proposed by Tachiki *et al.*[11]. Though the object was assumed to be a concentric spherical surface with the hologram surface, shift invariance does not exist due to unequal sampling points on a spherial grid(ie., the grid points are more crowded at the poles). To facilitate fast calculation using FFT, an approximation to the convolution integaral was proposed, which forced the PSF to be spatially invariant. This research work is an attempt to achieve the same but by deriving the transfer function and defining wave propagation in spectral domain. In other words, it is to find a spectral propagation formula (for spherical surfaces in spherical coordinates) analogous to the angular spectrum formula(for plane surfaces in cartesian coordinates). The expected advantages compared to the earlier method are faster calculation and no approximations. Section 2 explains the theoretical development of the calculation method. The numerical procedure for evaluating the formula is described in section 3. Simulation results to verify the proposed method are presented in Section 4. Section 5 summarises the work with possible applications and future improvements to this work.

## 2. Theoretical background

We will describe points on the sphere by their latitude *θ* ∈ [−*π*/2, *π*/2] and longitude *ϕ* ∈ [0, 2*π*] as shown in Fig. 1. Throughout this work we will be dealing with square-integrable functions that span a space *L*^{2}(*S*) on the unit sphere *S*. Fast computation algorithms take advantage of the cyclic and peroidic properties of the transformation kernel for fast calculations. The cyclic and periodic properties exists only if the system is shift invariant between the transformation planes for that operation. Which means that, the object and hologram surface should be shift invariant in order to devise a fast computation scheme. To achieve this the object and hologram surface are chosen to be concentric spherical surfaces, so that they can remain shift invariant in rotation along *ϕ* and *θ* directions. But these surfaces are defined on a spherical grid, where the sampling points are more dense at the poles than at the equator. Hence shift invariance is not satisfied. However in the hologram generation process, the hologram and object are band-limited functions on *L*^{2}(*S*) space. The band-limited functions on *L*^{2}(*S*) space has a very useful and important property that *any rotated version of a band-limited function is also a band-limited function with the same bandwidth*[12, 13]. Thus, they are refered to as having *uniform resolution*, at all points on the sphere, meaning that they are shift-invariant. The triangular truncation and Gauss-legendre quardature method that occurs in the transformation operation is responsible for this property(explained in the next section). Hence the system shown in Fig. 1 does have shift invariance relationship between object and hologram surfaces and approves the possiblity of fast computation formula. With this assurance we proceed to develop the fast calculation method for computer generated spherical holography starting from the basic equations of electromagnetism.

An electromagnetic field is defined by Maxwell’s equations and its propagation by the Helmholtz wave equation. Hence for any particular system the complex amplitude of a propagating wave at any instance of time and any where in space, can be found by solving the wave equation, applying its constraints and conditions. Accordingly for the system shown in Fig. 1 the solution can be derived starting from the wave equation as follows, The time dependent vector wave equation **u**(*r*,*θ*, *ϕ*) is expressed as

*r*is the radius of the spherical surface of interest and

*θ*,

*ϕ*represent the azimuthal and polar angles in the surface. The Laplacian operator ∇

^{2}in spherical co-ordinates is defined as

**u**” can be discarded. The scalar wave equation given in spherical coordinates becomes

*m*must be an integer so that there is continuity and periodicity in Φ(

*ϕ*). Φ

_{1}and Φ

_{2}are constants.

The solution to polar Eq. (6) is

_{1}and Θ

_{2}are constants. ${Q}_{n}^{m}$ is not finite at the poles where

*cos*(

*θ*) = ±1, so this solution is discarded(Θ

_{2}= 0).

For the radial differential equation Eq. (7) the solutions are

where, ${h}_{n}^{(1)}$ and ${h}_{n}^{(2)}$ are the Spherical Hankel functions of the first and second kind respectively. Since we are interested only in outgoing wave, we can neglect the second term (*R*

_{2}= 0)

The solution to Eq.(8) is

Again here we assume*T*

_{2}= 0 due to the convention of time

The angle functions in Eq. (9) and Eq. (10) are conveniently combined into a single function called as spherical harmonic ${Y}_{n}^{m}$[14, 15] defined by

*is called as the Condon-Shortley phase. Hence the spherical harmonics can also be represented in a short form as shown below(neglecting the Condon-Shortley phase) Combining all the above equations, the travelling wave solution to Eq. (3) can be represented as*

^{m}The radiated field is completely defined when the coefficient *A _{mn}* is determined. This is achieved by using the orthonormal property of the spherical harmonics. Assume that the wave-field

*u*(

*r*,

*θ*,

*ϕ*) is known on a sphere of radius

*r*=

*a*. We also drop the time variable(which is not important) for simplicity.Now multiplying each side of Eq. (16) (evaluated at

*r*=

*a*) by ${Y}_{n}^{m}{(\theta ,\varphi )}^{*}$ and integrating over the sphere, gives

*θ*)d

*θ*d

*ϕ*′, is the solid angle for integrating on a sphere. Inserting the expression for

*A*back into Eq. (16), we get,

_{mn}Hence the wavefield at any spherical surface *u*(*r*,*θ*, *ϕ*) can be calculated knowing the wave-field at *u*(*a*,*θ*′, *ϕ*′).

For easy understanding and interpretation of Eq. (18), it is compared with the well known angular spectrum of plane waves equation [3] which is given by

*k*,

_{x}*k*). The propagation of the spectrum is defined by the quantity

_{y}*e*which is known as the transfer function. The propagated spectrum is recomposed into the wavefield at destination by the integral over

^{ikzz}*k*,

_{x}*k*. Thus this equation defines wave propagation from one plane surface

_{y}*u*(

*x*′,

*y*′,0) to another surface

*u*(

*x,y,z*). When comparing Eq. (18) with Eq. (19), the following can be deduced

- The decomposed wave components(spectra) are expressed by Eq. (13) which is composed of a travelling wave component in
*ϕ*and defined by*e*and a standing wave component in^{imϕ}*θ*given by ${P}_{n}^{m}(\text{cos}\theta )$ - The decomposed wave components in this system can be named as spherical wave components in analogous to plane wave components for planar system. Similarly Eq. (20) can be termed as Spherical wave spectrum as analogus to Angular spectrum of plane waves.
- The wavenumbers
*k*and_{x}*k*are imitated by the quantities_{y}*m/a*and*n/a*. Hence we can refer to the spherical wave spectrum as a k-space spectrum due to this anlogy. - Equation (20) can be viewed as a forward Fourier transform using ${Y}_{n}^{m}(\theta ,\varphi )$ as the basis function. In other words, the spectral space for spherical system is spanned by Spherical harmonic functions( ${Y}_{m}^{n}(\theta ,\varphi )$). This is also termed as Spherical Harmonic Transform (SHT)
- Hence the quantity
*h*(_{n}*kr*)/*h*(_{n}*ka*) can be refered to the Transfer function(TF) for spherical system, as opposed to the quantity*e*for planar system.^{ikzz}

The transfer function is crucial since it completely defines propagation and hence it is worth discussing some of its properties. During propagation amplitude and phase of the spectral components change with distance as defined by the transfer function. What is most important is the rate of change of phase of the transfer function which determines the sampling requirements. Accordingly the plot shown in Fig. 2 reveals that the phase change increases with increasing orders ’n’ of the transfer function (The plot was generated for wavelength 100*μm*, radius 10 and 0.5 cm, and upto 256 orders). Hence sampling requirements will be satisfied if highest order ’n’ of the transfer function is sampled according to Nyquist criteria. It can also be understood that, the rate of phase change becomes higher as the distance of propagation increases. This will demand very large sampling and also increase the numerical errors. It is worthy to note here that the spherical Hankel functions are asymptotic in nature. In the far field the spherical Hankel functions can be approximated by their asymptotic expressions which is as shown in Eq. 23.

Thus the devised formula which is analogus to the angular spectrum of plane waves(AS) formula defines wave propagation between spherical surfaces. The numerical procedure to implement the devised formula is discussed in the next section.

## 3. Numerical computation

The numerical computation of AS method heavily depends on the FFT operations for which a lot of tools and methods are available. But the numerical computation of the proposed method heavily depends on the SHT operations. Fast computation was guaranteed from the theory and from the geometry of the system. Now a numerical procedure that takes advantage of this is required. Fortunately lot of fast computation numerical methods have been reported for SHT and this research work could make use of it. Since FFT and numerical computation of AS method are well understood, this section intends to introduce SHT and numerical computation of proposed method in close anology to the former.

Continuing with the comparision from previous section, the numerical computation of wave propagation for planar and spherical systems according to Eq. (18) and Eq. (19) can be represented respectively as

where, FFT[...] and IFFT[...] denote the foward and inverse fast Fourier transform operations, while SHT[...] and ISHT[...] denote the forward and inverse spherical harmonic transform operations.*TF*and

_{s}*TF*are the transfer functions for spherical and planar wave propagations as defined by Eq. (18) and Eq. (19) respectively.Comparison reveals that the computation method (for spherical system) is analogous to the angular spetrum of plane waves method (for plane system) except the fact that, the Fourier Transform is replaced by the Spherical harmonic transform. Hence evaluation of SHT becomes the key, the others being only basic mathematical evalutations. Spherical harmonic transform have been studied extensively and fast computation algorithms and optimisation methods have been proposed. A method that requires only

_{c}*O*(

*N*

^{2}log

*N*) operations for

*N*grid points as opposed to the standard

*N*

^{3}operations was proposed by Chien

*et al.*[12]. They imposed truncations on the spectral components and used fast multipole method and fast Fourier transform for evaluation. Their method is called as “spectral truncation method” and the errors due to truncation were well within acceptable limits. Later on Healey

*et al.*[16] could achieve the same using

*O*(

*N*(log

*N*))

^{2}operations. They took advantage of the recursive properties of associated Legendre polynomial for fast computation. This method was found to be more efficient and hence was choosen for numerical evaluation the Spherical harmonic transforms in this work. A brief outline of the numerical evaluation is presented. For more details please refer to Healey [16] and Driscoll [13].

A 2D DFT is computed by separating the transformation kernel into its variables, and evaluating it as column transform followed by a row tranform. Similarly the spherical harmonic transform kernel
${Y}_{n}^{m}$ given by Eq. (15) is also variable separable and can be separated into *ϕ* component and *θ* component. Accordingly the transformation given by Eq. (20) can be represented as shown below

*U*(

^{m}*θ*) are evaluated for

*m*= −

*N*,...,

*N*as shown below where

*ϕ*= 2

_{i}*πi/I*for

*i*= 1,...,

*I*. The equispaced longitudes

*ϕ*enables the use of fast Fourier transform.

_{i}Second, the Legendre transform of the Fourier coefficients *U _{m}*(

*θ*) is to be evaluated for |

*m*| ≤

*n*≤

*N*. This is done using the Gaussian-Legendre quadrature as shown below,

*θ*and

_{j}*w*are respectively the Gauss nodes and weights and are calculated using the Fourier-Newton method as described by Swarztrauber [17]. The Gauss-Legendre quadrature replaces the integral by the sum. The fact that the summation runs only from |

_{j}*m*| to

*N*is refered to as the triangular truncation. The use of Gauss-Legendre quadrature method redistributes

*θ*into Gaussian nodes

*θ*. This along with the triangular truncation are responsible for uniform resolution on the latitudinal points. Again the triangular truncation along with the recurrance property of Legendre polynomial helps to achieve fast computation.

_{j}In the similar way the inverse spherical harmonic transform can be represented as

*O*(

*N*(log

*N*)

^{2}operations for the SHT computation.

There are a lot of software tools available on the internet to do the SHT operation. Most of them are tuned and dedicated for a geophysical process which requires only real SHT but holography requires complex SHT. However the package SHTools by Mark Wieczorek [18] could do a complex SHT operation in Fortran language. Though not tested by us this is the best one available that suits the work related to this paper. Hence we recommend using this if it is required to quickly reproduce the work in this paper. The next section describes the testing and verification of this numerical procedure using simulations.

## 4. Simulation results

The system considered for simulation experiments is shown in Fig. 1. The object(*O*(*a*,*θ*, *ϕ*)) is a spherical surface of radius 1*cm* and the hologram(*H*(*r*,*θ*, *ϕ*)) is another concentric spherical surface of radius 10*cm*. The reference is considered to be a virtual source emitting spherical waves from the center, ie., the wavefield due to reference has same phase and amplitude throughout the hologram plane. This is similar to using a plane reference wave with normal incidence in plane holography.

It is first required to verify the correctness of the method by expecting it to reproduce the already known diffraction results. Hence the proposed method is subjected to generate already reported diffraction patterns. For this the diffraction pattern reported by Tachiki et al [11] for spherical surfaces is used as reference. Accordingly, a simple object was chosen which is a spherial surface with two irradiating points at *ϕ* = −*π*/2 and *ϕ* = *π*/2, as shown in Fig. 3. The object and hologram are composed of 256 pixels in the longitudinal (north-south) direction and 512 pixels in the latitudinal(east-west) direction. The wavelength was chosen to be *λ* = 100*μm*, in order to reduce sampling requirements and visiblity of fringes. The procedure for numerical generation of hologram using the proposed method is expressed in an abstract form as shown below.

Then the proposed method is tested to see whether it obeys the fundamental laws of diffraction and interference. In other words, it is to verify that it has the same qualities as angular spectrum method. For this two simulation experiments for qualitative analysis was performed. First it was intended to analyze the change in interference pattern with change in wavelength which is a fundamental law. Accordingly for the object shown in Fig. 3, the hologram was computed for wavelengths varying as a)150*μm*, b) 200 *μm*, c) 250 *μm*, d) 300 *μm*, e)350*μm* and f)400*μm* respectively. The corresponding holograms generated are shown in Fig. 5. As expected the fringe density decreases with increase in wavelength. Second another fundamental law which is the change in interference pattern with the change in distance between the coherent sources was verified. This also corresponds to the youngs double slit experiment. Accordingly, the hologram is computed for varying position of the pair of the point sources on the spherical surface. This setup is similar to the young double slit experiment. The positions of point sources were set to be a)(*ϕ* = −*π*/6, *ϕ* =*π*/6), b)(*ϕ* = −*π*/8, *ϕ* =*π*/8), c) (*ϕ* = −*π*/16, *ϕ* =*π*/16) and d) (*ϕ* = −*π*/32, *ϕ* = *π*/32) as shown in Fig. 6(a)–(d) respectively. The corresponding computed hologram pattern are shown in Fig. 6(e)–(h). As expected the fringe density decreases with decrease in distance between the point sources. Since the diffraction formula agrees well with the fundamental laws of diffraction it is confirmed that it behaves like the any other diffraction formula and can be used to simulate wave propagation. Moreover the above mentioned results also reveals that the wavefield on the spherical surface was computed correctly and as expected.

Since the theory is developed in context to computer generated holography, it is mandatory to verify its applicability to the same. For this it was decided to perform spherical hologram generation and then reconstruction from the hologram on the computer using the proposed method. The spherical object for which the hologram is to be made is shown in Fig. 7. The object and hologram were composed of 256 × 512 pixels. The wavelength for simulation was chosen to be *λ* = 30 *μm*. The generated hologram is shown in Fig. 8.

From this hologram the object was reconstructed back onto the original spherical surface using Eq. 18. This means that we are attempting to reconstruct a real image of the hologram and not the virtual one. For this the conjugate of the reference beam has to be used. Accordingly the numerical procedure of the for reconstruction is expressed in an abstract form as shown below.

*O*(

*N*(log

*N*))

^{2}operations for

*N*sampling points and hence it is a fast computation formula. The calculations were executed using a scripted language-python in a Dell Precision T7400 machine with 12 GB of RAM memory. A comparision of calculation time for the direct integration, convolution and spectral methods is shown in Table 1.

## 5. Concluding remarks

The solution to Helmholtz wave equation in spherical coordinates is derived using variable separable method. The spherical wave spectrum and transfer function were defined from the solution. A formula for computing the wave propagation from irradiating spherical surfaces is devised. A fast computation method that evaluates the wave propagation formula in *O*(*N*(log*N*))^{2} operations was suggested. The proposed method was tested for correctness using simulated experiments. Generation and reconstruction of a spherical hologram for a spherical object was also successful. Hence a new and fast computation method for computer genetated spherical holograms is realized.

This method can be extended by considering it as a farfield diffraction problem by utilizing its asymptotic form as shown in Eq. (23). This will be analogous to the Fresnel diffraction formula which will take this development more close to practical applications. More over other similar problems (in acoustics and geophysics) have shown that far field diffraction pattern is directional and not uniformly distributed over a spherical surface. Hence we do not need a spherical detector (which is not available now) or scan over the whole spherical surface in order to holographically image a spherical volume. This interesting fact takes it even more closer to realizing practical applications. The availability of optic fibers will also be a great aid in this process. When considering holographic display, printing a spherical hologram on a flexible film and reconstructing using a laser will be a challenging task. However with the availability of high precision lithographic machines and optic fibers this is not an impossible task. It is also worth refereing to Horvath [19] for more interesting facts on the importance of this research work to the future of lasers. Hence the days when computer generated spherical holography will have an application is not so far. We believe this research will be a small step in this regard and will aid the development of computer generated spherical holography.

## References and links

**1. **C. Gmachl, F. Capasso, E. E. Narimanov, J. U. Nockel, A. D. Stone, J. Faist, D. L. Sivco, and A. Y. Cho, “High power directional emission from microlasers with
chaotic resonators,” Science **280**, 1544–1545 (1998). [CrossRef]

**2. **N. M. Lawandy and R. M. Balachandran “Random laser?,” Nature **373**, 204, (1995). [CrossRef]

**3. **J. W. Goodman, *Introduction to Fourier Optics* (McGraw-Hill, 1968).

**4. **T. Tommasi and B. Bianco, “Computer-generated holograms of tilted planes by a spatial frequency approach,” J. Opt. Soc. Am. A **10**, 299–305 (1993), http://www.opticsinfobase.org/josaa/abstract.cfm?URI=josaa-10-2-299. [CrossRef]

**5. **Y. Sakamoto and M. Tobise, “Computer generated cylindrical hologram,” in Practical Holography XIX: Materials and Applications Tung H. Jeong and Hans I. Bjelkhagen, eds., Proc. SPIE **5742**, 267–274 (2005).

**6. **A. Kashiwagi and Y. Sakamoto, “A fast calculation method of cylindrical computer-generated holograms which perform image-reconstruction of volume data,” in Adaptive Optics: Analysis and Methods/Computational Optical Sensing and Imaging/Information Photonics/Signal Recovery and Synthesis Topical Meetings on CD-ROM OSA Technical Digest (CD) (Optical Society of America, 2007), paper DWB7, http://www.opticsinfobase.org/abstract.cfm?URI=DH-2007-DWB7.

**7. **T. Yamaguchi, T. Fujii, and H. Yoshikawa, “Fast calculation method for computer-generated cylindrical holograms,” Appl. Opt. **47**, D63–D70 (2008), http://www.opticsinfobase.org/ao/abstract.cfm?URI=ao-47-19-D63. [CrossRef] [PubMed]

**8. **Y. Sando, M. Itoh, and T. Yatagai, “Fast calculation method for cylindrical computer-generated holograms,” Opt.Express **13**, 1418–1423 (2005), http://www.opticsinfobase.org/oe/abstract.cfm?URI=oe-13-5-1418. [CrossRef] [PubMed]

**9. **B. J. Jackin and T. Yatagai, “Fast calculation method for computer-generated cylindrical hologram based on wave propagation in spectral domain,” Opt.Express **18**, 25546–25555 (2010), http://www.opticsinfobase.org/oe/abstract.cfm?uri=oe-18-25-25546. [CrossRef] [PubMed]

**10. **B. J. Jackin and T. Yatagai, “360° reconstruction of a 3D object using cylindrical computer generated holography,” Appl. Opt. **50**, H147–H152 (2011), http://www.opticsinfobase.org/ao/abstract.cfm?URI=ao-50-34-H147. [CrossRef] [PubMed]

**11. **M. L. Tachiki, Y. Sando, M. Itoh, and T. Yatagai, “Fast calculation method for spherical computer-generated holograms,” Appl. Opt. **45**, 3527–3533 (2006), http://www.opticsinfobase.org/ao/abstract.cfm?URI=ao-45-15-3527. [CrossRef] [PubMed]

**12. **R. J. Chien and B. K. Alpert, “A fast spherical filter with uniform resolution,” J. Comput. Phys. **136**, 580–584 (1997). [CrossRef]

**13. **J. R. Driscoll and D. M Healy, “Computing Fourier transfroms and convolutions on the sphere,” Adv. Appl. Math. **15**, 201–250 (1994). [CrossRef]

**14. **N. N. Lebedev, *Special Functions and their Applications* (Prentice Hall, 1965).

**15. **G. B. Arfken and H. J. Weber, *Mathematical Method for Physicist* (Academic Press, 2001).

**16. **D. M. Healy Jr., D. Rockmore, P. J. Kostelec, and S. Moore, “FFTs for the 2-sphere - improvements and variations,” J. Fourier. Anal. Appl. **9**, No.4, 341–385(1998). [CrossRef]

**17. **P. N. Swarztrauber, “On computing the points and weights for Gauss-Legendre quadrature,” SIAM. J. Sci. Computing **24**. Issue. 3, 945–954(2002). [CrossRef]

**18. **M. Wieczorek, SHTools, http://shtools.ipgp.fr/.

**19. **Z. G. Horvath, “Beyond the beam: A history of multidimensional lasers,” Opt. Photonics News **23**. No. 7/8, 36–41(2012). [CrossRef]