## Abstract

A simulator is designed in MATLAB code which gives the propagation characteristics of a general-type beam in turbulent atmosphere. When the required source and medium parameters are entered, the simulator yields the average intensity profile along the propagation axis in a video format. In our simulator, the user can choose the option of a “user defined beam” in which the source and medium parameters are selected as requested by the user by entering numerical values in the relevant menu boxes. Alternatively, the user can proceed with the option of “pre-defined beam” in which the average intensity profiles of beams such as annular, cos-Gaussian, sine-Gaussian, cosh-Gaussian, sinh-Gaussian, their higher-order counterparts and flat-topped can be observed as they propagate in a turbulent atmosphere. Some samples of the simulator output are presented.

© 2006 Optical Society of America

## 1. Introduction

The variation of the received average intensity profile along a turbulent path depends on the type of the beam used as incidence. This dependency is investigated by many researchers for various types of beams, such as fundamental mode [1,2], Laguerre Gaussian [3], higher order mode [4], cosh-Gaussian and cos-Gaussian [5,6], Hermite-sine-Gaussian and Hermite-sinh- Gaussian [7], elliptical Gaussian [8], dark hollow [9], flat-topped [10], and higher-order annular Gaussian [11]. In all of these studies, the formulation is developed starting with the intended incidence. In this paper, we have managed to combine the formulation of the average received intensity in turbulence for most of the mentioned beams in one compact form and designed a simulator in MATLAB code. For this purpose, we have formulated the average received intensity by using general-type beam incidence [12] and transformed this formulation into a form of simulator. The simulator input can be either “user defined beam” in which the general beam source and medium parameters are entered by the user in the relevant menu boxes or “pre-defined beam” in which the known beam source and medium parameters can be entered. As the output of the simulator, one can observe in video format the average received intensity distribution of the “user defined beam” or “pre-defined beam” as the beam propagates along a link having atmospheric turbulence. Here, the video format refers to the progress of the average intensity profiles along the propagation axis, not the timely variation of the intensity. Our simulator is based on horizontal links where the structure constant of turbulence is constant, and it also handles the free space (no turbulence) propagation characteristics of general-type beams as special cases.

## 2. Formulation

A general beam, excluding the multimode composition, emitted from a source co-centric with its own plane will be given as [12]

$$\times {H}_{{m}_{\ell}}\left({a}_{y\ell}{s}_{y}+{b}_{x\ell}\right)\mathrm{exp}\left[-\left(0.5k{\alpha}_{y\ell}{s}_{y}^{2}+j{V}_{yt}{s}_{y}\right)\right],$$

where, *s*_{x}
and *s*_{y}
are the *x* and *y* components of the source plane vector **s**, such that **s**=(*s*_{x}
, *s*_{y}
). All *ℓ* subscripted terms establish the specific parameters of the individual beams comprising the general beam through summation. In this manner, *N* denotes the number of beams, *A*_{ℓ}
is the amplitude factor, *θ*_{ℓ}
is the phase, ${H}_{{n}_{\ell}}\left({a}_{x\ell}{s}_{x}+{b}_{x\ell}\right)$ and ${H}_{{m}_{\ell}}\left({a}_{y\ell}{s}_{y}+{b}_{x\ell}\right)$ are Hermite polynomials governing the beam variations for *s*_{x}
and *s*_{y}
directions, where *n*_{ℓ}
and *m*_{ℓ}
are the order, *a*_{xℓ}
and *a*_{yℓ}
characterize the width, *b*_{xℓ}
and *b*_{yℓ}
are the complex displacement parameters, *V*_{xℓ}
and *V*_{yℓ}
are the complex parameters used to create physical location displacement and phase rotation or a combination of both, named as the displacement parameters, *α*_{xℓ}
and *α*_{yℓ}
are related to the source sizes, *α*_{sxℓ}
, *α*_{syℓ}
and source focusing parameters, *F*_{xℓ}
, *F*_{yℓ}
along *s*_{x}
and *s*_{y}
directions via

where ere *k*=2*π*/*λ* is the wave number with *λ* being the wavelength and *j*=(-1)^{0.5}.

The average intensity 〈*I*(**p**,*L*)〉, of a general beam on a receiver plane located at *L* distance away from the source can be found from extended Huygens-Fresnel integral in the following way [13]

$$\times \mathrm{exp}\{jk\left[{\left(\mathbf{p}-{\mathbf{s}}_{1}\right)}^{2}-{\left(\mathbf{p}-{\mathbf{s}}_{2}\right)}^{2}\right]\u20442L\}\u3008\mathrm{exp}\left[\psi ({\mathbf{s}}_{1},\mathbf{p})+{\psi}^{*}({\mathbf{s}}_{2},\mathbf{p})\right]\u3009,$$

where the ensemble average term within the integrand is [2]

with *ψ* being the fluctuations of the complex amplitude, *p*_{x}
and *p*_{y}
are the *x* and *y* components of the receiver plane vector **p**, such that **p**=(*p*_{x}
,*p*_{y}
), *ρ*
_{0}=(0.545 ${C}_{n}^{2}$*k*
^{2}
*L*)^{-3/5} is the coherence length of a spherical wave propagating in the turbulent medium and ${C}_{n}^{2}$
, the structure constant. We note that Eq.(4) is derived under the quadratic approximation for the Rytov’s phase structure function.

The source field *u*(**s**
_{1}), when multiplied by its conjugate *u*
^{*}(**s**
_{2}), indicated by the sign *, will be

$$\times {H}_{{n}_{{\ell}_{1}}}\left({a}_{x{\ell}_{1}}{s}_{1x}+{b}_{x{\ell}_{1}}\right)\mathrm{exp}\left[-\left(0.5k{\alpha}_{x{\ell}_{1}}{s}_{1x}^{2}+j{V}_{x{\ell}_{1}}{s}_{1x}\right)\right]$$

$$\times {H}_{{m}_{{\ell}_{1}}}\left({a}_{y{\ell}_{1}}{s}_{1y}+{b}_{y{\ell}_{1}}\right)\mathrm{exp}\left[-\left(0.5k{\alpha}_{y{\ell}_{1}}{s}_{1y}^{2}+j{V}_{y{\ell}_{1}}{s}_{1y}\right)\right]$$

$$\times {H}_{{n}_{{\ell}_{2}}}^{*}\left({a}_{x{\ell}_{2}}{s}_{2x}+{b}_{x{\ell}_{2}}\right)\mathrm{exp}\left[-\left(0.5k{\alpha}_{x{\ell}_{2}}^{*}{s}_{2x}^{2}-j{V}_{x{\ell}_{2}}^{*}{s}_{2x}\right)\right]$$

$$\times {H}_{{m}_{{\ell}_{2}}}^{*}\left({a}_{y{\ell}_{2}}{s}_{2y}+{b}_{y{\ell}_{2}}\right)\mathrm{exp}\left[-\left(0.5k{\alpha}_{y{\ell}_{2}}^{*}{s}_{2y}^{2}-j{V}_{y{\ell}_{2}}^{*}{s}_{2y}\right)\right],$$

After substituting Eqs. (4) and (5) into Eq. (3) and solving the resulting integral by the repeated use of Eq. 3.462.2 of Ref. [14], which is ${\int}_{-\infty}^{\infty}\mathit{dx}{x}^{n}\mathrm{exp}\left(-{\mathit{px}}^{2}+2\mathit{qx}\right)=n!\mathrm{exp}({q}^{2}\u2044p){(\pi \u2044p)}^{0.5}{(q\u2044p)}^{n}\sum _{t=0}^{[n\u20442]}{p}^{t}\u2044\left[{\left({4q}^{2}\right)}^{t}\left(n-2t\right)!\left(t\right)!\right],$, we obtain the following expression for the average intensity at the receiver plane

where $b=k\u20442L,{a}_{sx{\ell}_{1}}=0.5k{\alpha}_{x{\ell}_{1}}+1\u2044{\rho}_{0}^{2},{a}_{sx{\ell}_{2}}=0.5k{\alpha}_{x{\ell}_{2}}+1\u2044{\rho}_{0}^{2}$, ${a}_{sy{\ell}_{1}}$ and ${a}_{sy{\ell}_{2}}$ are the y counterparts of ${a}_{sy{\ell}_{1}}$ and ${a}_{sy{\ell}_{2}}$,

*E*_{y}
is attained by changing all *x* subscripts to *y* in *E*_{x}
,

$$\times {\left(-j\right)}^{{\ell}_{nx\phantom{\rule{.2em}{0ex}}1}+{\ell}_{nx\phantom{\rule{.2em}{0ex}}2}-2{k}_{x\phantom{\rule{.2em}{0ex}}1}-2{k}_{x\phantom{\rule{.2em}{0ex}}2}}{\left(-1\right)}^{{\ell}_{x\phantom{\rule{.2em}{0ex}}1}+{\ell}_{x\phantom{\rule{.2em}{0ex}}2}}{2}^{{n}_{{\ell}_{1}}+{n}_{{\ell}_{2}}-{\ell}_{x\phantom{\rule{.2em}{0ex}}1}-{\ell}_{x\phantom{\rule{.2em}{0ex}}2}-{\ell}_{nx\phantom{\rule{.2em}{0ex}}1}-{\ell}_{nx\phantom{\rule{.2em}{0ex}}2}}$$

$$\times \left(\begin{array}{c}{n}_{{\ell}_{1}}\\ 2{\ell}_{x1}\end{array}\right)\left(\begin{array}{c}{n}_{{\ell}_{2}}\\ 2{\ell}_{x2}\end{array}\right)\left(\begin{array}{c}{{n}_{\ell}}_{1}-2{\ell}_{x1}\\ {\ell}_{nx1}\end{array}\right)\left(\begin{array}{c}{\ell}_{nx1}-2{k}_{x1}\\ {\ell}_{nx11}\end{array}\right)\left(\begin{array}{c}{{n}_{\ell}}_{2}-2{\ell}_{x2}\\ {\ell}_{nx2}\end{array}\right)\frac{\left({\ell}_{nx1}\right)!}{\left({\ell}_{nx1}-2{k}_{x1}\right)!\left({k}_{x1}\right)!}$$

$$\times \frac{\left({\ell}_{nx11}+{\ell}_{nx2}\right)!}{\left({\ell}_{nx11}+{\ell}_{nx2}-2{k}_{x2}\right)!\left({k}_{x2}\right)!}{\left({a}_{x{\ell}_{1}}\right)}^{{\ell}_{nx1}}{\left({a}_{x{\ell}_{2}}^{*}\right)}^{{\ell}_{nx2}}{\left({b}_{x{\ell}_{1}}\right)}^{{n}_{{\ell}_{1}}-2{\ell}_{x1}-{\ell}_{nx1}}{\left({b}_{x{\ell}_{2}}^{*}\right)}^{{n}_{{\ell}_{2}}-2{\ell}_{z2}-{\ell}_{nx2}}{\left({\rho}_{0}^{2}\right)}^{{\ell}_{nx2}}$$

$$\times {\left[{\rho}_{0}^{4}\left({a}_{sx{\ell}_{1}}-jb\right)\left({a}_{sx{\ell}_{2}}^{*}+jb\right)-1\right]}^{-{\ell}_{nx11}-{\ell}_{nx2}+{\ell}_{x2}}{\left({a}_{sx{\ell}_{1}}-jb\right)}^{-{\ell}_{nx1}+{\ell}_{x1}+{\ell}_{x2}}{\left({V}_{x{\ell}_{1}}+2b{p}_{x}\right)}^{{\ell}_{nx1}-2{k}_{x\ell}-{\ell}_{nx11}}$$

$$\times {\left[-{\rho}_{0}^{2}\left({V}_{x{\ell}_{2}}^{*}+2b{p}_{x}\right)\left({a}_{sx{\ell}_{1}}-jb\right)+{V}_{x{\ell}_{1}}+2b{p}_{x}\right]}^{{\ell}_{nx11}+{\ell}_{nx2}-2{k}_{x2}},$$

*S*_{y}
is constructed by changing all *x* subscripts to *y* in *S*_{x}
given by Eq. (8), *T*_{ℓx}
_{1}=1×3×....(2*ℓ*_{xi}
-1) for *ℓ*_{x}
_{1}≠0 where *i*=1,2, all appearances in the form of $\left(\begin{array}{c}{B}_{1}\\ {B}_{2}\end{array}\right)$ represent binomial coefficients, hence $\left(\begin{array}{c}{B}_{1}\\ {B}_{2}\end{array}\right)={B}_{1}!\u2044\left[\left({B}_{1}-{B}_{2}\right)!{B}_{2}!\right]$, ! is the factorial notation, the square brackets, [], placed in the upper limits of some summations mean that the integral part of the expression within the square brackets is to be taken.

## 3. Results and discussion

We have designed a simulator that will generate the average receiver intensity of a general beam, when the necessary source and medium parameters, namely *N*, *A*_{ℓ}
, *θ*_{ℓ}
, *α*_{sxℓ}
, *α*_{syℓ}
, *F*_{xℓ}
, *F*_{yℓ}
, *n*_{ℓ}
, *m*_{ℓ}
, *V*_{xℓ}
, *V*_{yℓ}
, *λ*,${C}_{n}^{2}$
,*L* are supplied to Eq. (6). Presently, the remaining parameters are internally set as *α*_{xℓ}
=1/*α*_{sxℓ}
, *a*_{yℓ}
=1/*α*_{syℓ}
, *b*_{xℓ}
=*b*_{yℓ}
=0.

Our simulator offers the choice of “user defined beam” and “pre-defined beam”. We emphasize that Eq. (6) together with Eqs. (7), (8) along with the associated definitions given in the text unify the expressions for propagation of a wide variety of basic beam types into a single expression facilitating numerical investigation. This unification is realized by designing a simulator that implements the propagation of general-type beams as provided in Eq. (6). A selection of such beams, also named as “pre-defined” types when running the simulator, are tabulated below with the relevant parameter usage.

For all beams, there exists *x*-*y* symmetry. In this context, *n*, *m* are positive integers, *α*_{s}
, *F* and *V* are positive numeric values conforming to the following equalities for cos-Gaussian, cosh-Gaussian, sine-Gaussian, sinh-Gaussian, Hermite-cosh-Gaussian, Hermite-sine-Gaussian, Hermite-sinh-Gaussian, Hermite-cos-Gaussian beams, *n*=*n*
_{1}=*n*
_{2}, *m*=*m*
_{1}=*m*
_{2}, *α*_{s}
=*α*_{sx}
_{1}=*α*_{sx}
_{2}=*α*_{sy}
_{1}=*α*_{sy}
_{2}, *F*=*F*_{x}
_{1}=*F*_{x}
_{2}=*F*_{y}
_{1}=*F*_{y}
_{2}, *V*=*V*_{x}
_{1}=*V*_{x}
_{2}=*V*_{y}
_{1}=*V*_{y}
_{2}. For higherorder annular beam *α*_{s}
_{1}>*α*_{s}
_{2}.

Upon launching, the screen layouts after having made the choice of “user defined beam” or “pre-defined beam” are shown in Figs. 1(a) and 1(b), respectively. User defined beams are those whose parameters are freely selected by the user by inserting the numerical value in the relevant menu boxes displayed to the upper right hand side of the simulator’s page. This way, any type of beam can be constructed.

Within the context of this paper, pre-defined beams are confined to fundamental annular, cos-Gaussian, sine-Gaussian, cosh-Gaussian and sinh-Gaussian, flat-topped beams and their higher orders. In all these cases, except the flat-topped case, *N*=2, phase factors are *θ*_{ℓ}
=[0 0]. The amplitude factors alternate as follows; *A*_{ℓ}
[0.5 0.5] for cos-Gaussian and cosh-Gaussian beams, while *A*_{ℓ}
[0.5–0.5] for sinh-Gaussian and annular beams and *A*_{ℓ}
[0.5*j*-0.5*j*] for sine-Gaussian beam. The displacement parameters, and *V*_{xℓ}
*V*_{yℓ}
are purely real for cos and sine beams, and purely imaginary for cosh and sinh beams, they are also *x*-*y* symmetric and equal in magnitude but with opposite signs for the first and second beams of the summation, except for the annular beam, where *V*_{xℓ}
=*V*_{yℓ}
=0. In the flat-topped case, beams with descending source sizes and binomial coefficient related amplitude factors are summed up to the given flatness order, *N*. Hence the corresponding settings are ${A}_{\ell}=\frac{{\left(-1\right)}^{\ell}}{N}\left(\begin{array}{c}N\\ \ell \end{array}\right)$, α* _{sxℓ}* = α

*= α*

_{syℓ}*/ √ℓ = α*

_{sx}*/ √ℓ,*

_{sy}*F*=

_{xℓ}*F*= ∞,

_{yℓ}*n*=

_{ℓ}*m*= 0,

_{ℓ}*V*=

_{xℓ}*V*= 0. Upon selecting the pre-defined beam option, the default values for source and medium are assigned as shown in the menu boxes. The user can modify the source and medium parameters manually. We note that Figs. 2–7 are plotted for

_{yℓ}*λ*=1.55 µm, ${C}_{n}^{2}$ =1×10

^{-15}m

^{-2/3}.

Below, sample outputs from the simulator are provided. Figure 2 illustrates three pictures of a general beam during the program run captured at *L*=1 km with the settings given in the Fig. Here the upper left picture is the intensity plot of the beam viewed perpendicular to the propagation axis, where more brightness means higher intensity spots. The second picture to the right in the upper row refers to the contour plot of the beam. Finally the right one of the second row is the 3D picture of the beam. Figure 3 displays the same for a pre-defined Hermite-cosh-Gaussian beam at *L*=2 km, where all the settings are again given in the Fig. In Figs. 2 and 3, 〈*I*_{r}*N*〉 means that for the receiver intensity, the following normalization is applied

where the denominator of the right hand side corresponds to the maximum value of the source plane intensity. We note that Fig. 3 matches exactly the second intensity profile picture in Fig. 6 of Ref. [7].

Several samples from the simulator’s recorded video outputs are shown in Figs. 4–7. Figure 4 contains the four plots from the progress of the general beam belonging to Fig. 2 along the propagation axis at distances of *L*=0, 1, 3 and 5 km. The particular case considered here is a general beam type, therefore it is difficult to make any conclusive remarks regarding its propagation characteristics. Figure 5 shows the progress at *L*=0, 1, 3 and 5 km, this time for a pre-defined Hermite-sine-Gaussian beam with the settings of *α*_{sxℓ}
=*α*_{syℓ}
=[2 2] cm, *F*_{yℓ}
=*F*_{yℓ}
=[∞ ∞], *n*_{ℓ}
[1 1], *m*_{ℓ}
[0 0], *V*_{xℓ}
=*V*_{yℓ}
=[120 -120]. It can be seen from the combined examination of Fig. 5 and the associated movie file, when propagating in turbulence, Hermite-sine-Gaussian beams keep their original shape for a certain link length, then with the appearance of additional lobes, the overall profile turns into sinh-positioned TEM beam [7]. Again it can be verified that pictures of Fig. 5 agree well with the results of Ref. [7]. Here, the comparisons are made only with our earlier results, since to our knowledge, such related results for Hermite-cosh-Gaussian and Hermite-sine-Gaussian beams propagating in turbulence are not reported elsewhere in the literature. Figure 6 is for a pre-defined higher-order annular beam with the settings of *α*_{sx}
_{1}=3 cm, *α*_{sx}
_{2}=1.5 cm, *F*_{xℓ}
=*F*_{yℓ}
=[∞ ∞], *n*_{ℓ}
=[1 1], *m*_{ℓ}
[0 0], *V*_{xℓ}
=*V*_{yℓ}
[0 0]. In however this is hardly visible in parts (c) and (d), since the intensity level of the outer part of the beam is quite low due to propagating a long distance. The movie obtained from Fig. 6 shows that, during propagation in atmospheric turbulence, the higher-order annular beam will first act to enlarge the initially smaller lobes near on-axis, then a primary beam like formation is observed, finally at excessive propagation lengths, the profile will turn into a pure Gaussian shape. Figure 7 presents the progress at *L*=0, 1, 6 and 10 km, for a pre-defined flat-topped beam with the settings of *N*=10, *α*_{sx}
=*α*_{sy}
=3 cm, *F*_{xℓ}
=*F*_{xℓ}
=∞, *n*_{ℓ}
=*m*_{ℓ}
=0 and *V*_{xℓ}
=*V*_{yℓ}
. Again in Fig. 7, as the beam propagates, outer rings are formed, however, they are not clearly seen in the plots due to their lower intensity levels when compared with the bright central part of the beam. Figure 7 and its associated movie file demonstrate that, flat-topped beam will initially develop a circular ring in the center. As the propagation advances, the circumference of this ring will become narrower, while downward peak will gradually emerge from the center of the beam. Eventually this peak will smooth out the initial ring, with the profile turning into a pure Gaussian shape.

By launching the simulator program, it is possible to watch the entire movies of the beams belonging to Figs. 4–7 in the form of Figs. 2 and 3, as they travel along the propagation axis, starting at *L*=0 m up to *L*=5 km. In the simulator, the distance interval, at which an image of the beam is captured and displayed on the screen is 100 m. Finally we note that, by specifying ${C}_{n}^{2}$
=0 at the input menus of “user defined beam” and “pre-defined beam”, the user is able to acquire the propagation characteristics of such beams in free space, i.e., no turbulence versions.

## 4. Conclusion

A simulator is designed to provide the progress of the average received intensity profile of a general-type beam as it propagates in a turbulent horizontal link. This progress is shown on the screen at regular intervals of propagation length during the running of the program and the entire intensity plots are recorded in video format. Through the use of this simulator, after entering the appropriate source and medium parameters, the user is able to observe the average received intensity profiles in the above described formats in the form of “user defined beam” and “pre-defined beam” such as annular, cos-Gaussian, sine-Gaussian, cosh-Gaussian, sinh-Gaussian, their higher-order counterparts and flat-topped beams when propagating in turbulent atmosphere. As limiting case solutions, if the structure constant is entered as zero, our simulator also yields the intensity profiles of general-type beams in video format along the free space path.

## References and links

**1. **Z. I. Feizulin and Y. Kravtsov, “Broadening of a laser beam in a turbulent medium,” Radiophys. Quantum Electron **10**, 33–35 (1967). [CrossRef]

**2. **S. C. H. Wang and M. A. Plonus, “Optical beam propagation for a partially coherent source in the turbulent atmosphere,” J. Opt. Soc. Am. **69**, 1297–1304 (1979). [CrossRef]

**3. **R. L. Phillips and L. C. Andrews, “Spot size and divergence for Laguerre Gaussian beams of any order,” Appl. Opt. **22**, 643–644 (1983). [CrossRef] [PubMed]

**4. **C. Y. Young, Y. V. Gilchrest, and B. R. Macon, “Turbulence induced beam spreading of higher order mode optical waves,” Opt. Eng. **41**, 1097–1103 (2002). [CrossRef]

**5. **H. T. Eyyuboğlu and Y. Baykal, “Average intensity and spreading of cosh-Gaussian laser beams in the turbulent atmosphere,” Appl. Opt. **44**, 976–983 (2005). [CrossRef] [PubMed]

**6. **H. T. Eyyuboğlu and Y. Baykal, “Analysis of reciprocity of cos-Gaussian and cosh-Gaussian laser beams in turbulent atmosphere,” Opt. Express **12**, 4659–4674 (2004). [CrossRef] [PubMed]

**7. **H. T. Eyyuboğlu and Y. Baykal, “Hermite-sine-Gaussian and Hermite-sinh-Gaussian laser beams in turbulent atmosphere,” J. Opt. Soc. Am. A **22**, 2709–2718 (2005). [CrossRef]

**8. **Y. Cai and S. He, “Average intensity and spreading of an elliptical Gaussian beam in a turbulent atmosphere,” Opt. Lett. **31**, 568–570 (2006). [CrossRef] [PubMed]

**9. **Y. Cai and S. He, “Propagation of various dark hollow beams in a turbulent atmosphere,” Opt. Express **14**, 1353–1367 (2006). [CrossRef] [PubMed]

**10. **H. T. Eyyuboğlu, Ç. Arpali, and Y. Baykal, “Flat topped beams and their characteristics in turbulent media,” Opt. Express **14**, 4196–4207 (2006). [CrossRef] [PubMed]

**11. **H. T. Eyyuboğlu, S. Altay, and Y. Baykal, “Propagation characteristics of higher-order annular Gaussian beams in atmospheric turbulence,” Opt. Commun. **264**25–34 (2006). [CrossRef]

**12. **Y. Baykal, “Formulation of correlations for general-type beams in atmospheric turbulence,” J. Opt. Soc. Am. A **23**, 889–893 (2006). [CrossRef]

**13. **L. C. Andrews and R. L. Phillips, *Laser Beam Propagation through Random Media*, (SPIE, Bellingham, Washington, 2005).

**14. **I. S. Gradshteyn and I. M. Ryzhik, *Tables of Integrals, Series and Products* (Academic Press, New York, 2000).