Stochastic optimization for a tip-tilt adaptive correcting system

ARTICLE IN PRESS UNCORRECTED PROOF S0010-4655(02)00681-1/FLA AID:2511 Vol.•••(•••) P.1 (1-19) ELSGMLTM(COMPHY):m3 2002/09/24 Prn:1/10/2002; 16:39 cpc2511 by:PS p. 1 Computer Physics Communications ••• (••••)•••–••• www.elsevier.com/locate/cpc 149 250 351 452 553 654 755 856 957 10 58 11 59 12 60 13 61 14 62 15 63 16 64 17 65 18 66 19 67 20 68 21 69 22 70 23 71 24 72 25 73 26 74 27 75 28 76 29 77 30 78 31 79 32 80 33 81 34 82 35 83 36 84 37 85 38 86 39 87 40 88 41 89 42 90 43 91 44 92 45 93 46 94 47 95 48 96 Stochastic optimization for a tip-tilt adaptive correcting system✩ M.S. Zakynthinakia,∗,Y.G.Saridakisb aDepartament de Matemática Aplicada I, Universitat Politécnica de Catalunya, Diagonal 647, E-08028 Barcelona, Spain bApplied Mathematics and Computer Laboratory, Department of Sciences, Technical University of Crete, 73100, Chania, Greece Received 24 April 2002; received in revised form 26 June 2002; accepted 13 July 2002 Abstract We present computer simulations of a tip-tilt adaptive optics system, where stochastic optimization is applied to the problem of dynamic compensation of atmospheric turbulence. The system uses a simple measure of the light intensity that passes through a mask and is recorded on the image plane, to generate signals for the tip-tilt mirror. A feedback system rotates the mirror adaptively and in phase with the rapidly changing atmospheric conditions. Computer simulations and a series of numerical experiments investigate the implementation of the method in the presence of drifting atmosphere. In particular, the study examines the system’s sensitivity to the rate of change of the atmospheric conditions and investigates the optimal size of the mirror’s masking area and the algorithm’s optimal degree of stochasticity. 2002 Elsevier Science B.V. All rights reserved. PACS: 02.60.Pn; 07.05.Tp; 95.75.Qr; 95.75Pq Keywords: Stochastic optimization; ALOPEX algorithm; Adaptive wavefront correction; Tip-tilt mirror; Masking 1. Introduction Atmospheric layers act, due to turbulence, as a distorting optical system which causes parallel light rays to diverge. A typical short exposure image (exposure time <2–20 msec), formed by viewing a point source through turbulence, does not consist of a single diffraction pattern with a diameter fixed by the diffraction limit of the telescope, but rather of a number of su- ✩This work was supported by the Greek General Secretariat of Research and Technology under the grant PENED-107.527. *Corresponding author. Partially supported by the grant SB 2000-0396 from the Secretaría de Educación y Universidades and Fondo Social Europes. E-mail addresses: [email protected] (M.S. Zakynthinaki), [email protected] (Y.G. Saridakis). perimpozed speckles, microscopic images of the point source, which are distributed over a diameter determined by the severity of turbulence (see Section 2.1). This phenomenon is referred to as atmospheric seeing [1]. Due to seeing, the capabilities of even the larger ground based telescopes remain unexploited. Today it is technologically possible to smooth out the distorting effects of the earth’s atmosphere by the use of adaptive optics systems, which are able to adaptively cancel out, or at least minimize, atmospheric seeing. The basic idea of adaptive optics was formulated in 1953 by Babcock [2,3]. However, adaptive optics systems were developed only recently, due to the technological limitations of the implementation of such systems. The basic limitation of an adaptive optics system lies in the fact that a largenumberof highly 0010-4655/02/$ – see front matter 2002 Elsevier Science B.V. All rights reserved. PII:S0010-4655(02)00681-1 ARTICLE IN PRESS UNCORRECTED PROOF S0010-4655(02)00681-1/FLA AID:2511 Vol.•••(•••) P.2 (1-19) ELSGMLTM(COMPHY):m3 2002/09/24 Prn:1/10/2002; 16:39 cpc2511 by:PS p. 2 2M.S. Zakynthinaki, Y.G. Saridakis / Computer Physics Communications ••• (••••)•••–••• 149 250 351 452 553 654 755 856 957 10 58 11 59 12 60 13 61 14 62 15 63 16 64 17 65 18 66 19 67 20 68 21 69 22 70 23 71 24 72 25 73 26 74 27 75 28 76 29 77 30 78 31 79 32 80 33 81 34 82 35 83 36 84 37 85 38 86 39 87 40 88 41 89 42 90 43 91 44 92 45 93 46 94 47 95 48 96 expensive and complex elements are needed for the wavefront sensing. By use of such sensors, the distorted wavefront is repeatedly evaluated and continuously corrected, in real-time, througha rotating and/or a deformable mirror. In this way the spatial resolution of the image of a star can be effectively improved up to the theoretical diffraction limit of the telescope system under use. Adaptive optics systems are currently in use in the world’s largest astronomical observatories. The present study investigates the development of an alternative adaptive optics method based on stochastic optimization, in which the expensive wavefront sensors will be replaced by a simple and cost effective system. According to Noll [4], the so-called tip-tilt distortion accounts for 85% of the aberration induced upon a wavefront, which is passing though a turbulent atmosphere. The tip-tilt is referred to as the first order distortions which define, in the case of atmospheric seeing, the time varying wavefront gradients. Assuming that the centroid of an image is the centre of its intensity distribution, the tip-tilt distortion results in a displacement of thecentroid and hence a blurring of the image. The correcting element for the system under study is assumed to be a flat mirror, which can be tilted along two orthogonalplanes to correctthe motion of theimage. For the higherorder— but less significant—corrections, a deformable mirror [1,5] is necessary. However, the primary concern of any adaptive optics system is the first order tip-tilt correction, by which a significant part of the distortions caused by seeing may be eliminated. Tip-tilt correction mayalso accountforthe minimizationof anythermal effects within the telescope enclosure and above the surface of the mirror. In addition, any optical aberrations caused by micro-movements or other possible external factors can be smoothed out. In the simulated tip-tilt adaptive optics system presented here, the wavefrontsensor is to be replaced by a mask applied directlyto theimage,on theimageplane. The mask allows only a fraction of the light information to pass through its narrow aperture. Behind the mask there is a light detector, evaluating the total light intensity passing through the mask and being recorded on the central area of the image plane. This central image area behind the mask’s aperture will be referred to as the masking area. A computer, by means of the optimization algorithm, drives the tip-tilt mirror in such a way as to bring the centroid (i.e. the centre of the light intensity of the image) over the aperture of the mask. By means of this optimization method, the restoration of the centroid of the image is achieved, in real time, resulting in the improvement of the light distribution of the long exposure image recorded (exposure time >2–20 msec). The process is repeated at each time step. Fig. 1 illustrates the optimization process. For the simulation purposes, atmospheric turbulence is modeled following the principles of the relative Kolmogorov theory [6] and by using a sequence of phase screens. A phase screen is a two-dimensional distribution of phase fluctuations, introduced into the optical system to simulate the effects of the turbulent atmosphere. These phase fluctuations, are numerically generated [1,5,7] by the so-called random mid-point displacement algorithm, used for the same purposes by Lane et al. [7] (see Section 2.2). The phase fluctuations are transported across the telescope aperture according to local wind conditions. For the optimization of the image, ALOPEX stochastic optimization is used (see Section 3 and [1,5, 8–11]). ALOPEX optimization methods are driven by the parallel incoherent dithering of the control variables and the time dependence of the feedback. Their main advantage is that no knowledge of the dynamics of the system is required [8,9]. The method’s effectiveness and practicality depend on the ability to follow the motion of the turbulent wavefronts. ALOPEX stochastic optimization algorithm is chosen due to its speed of convergence and its simple implementation [1,5,10,11]. In the present study we apply the ALOPEX stochastic optimization method to the problem of real time sharpness restoration of the image of a point star. We investigate •the typical method’s behaviour during the process of image restoration (Section 4.2); •the optimal size of the masking area (Section 4.3); •the system’s sensitivity to the rate of change of atmospheric distortions (Section 4.4.1); •the relation of the algorithm’s optimal mean noise amplitude and the rate of change of atmospheric distortions, in the case of bad weather conditions (Section 4.4.2). ARTICLE IN PRESS UNCORRECTED PROOF S0010-4655(02)00681-1/FLA AID:2511 Vol.•••(•••) P.3 (1-19) ELSGMLTM(COMPHY):m3 2002/09/24 Prn:1/10/2002; 16:39 cpc2511 by:PS p. 3 M.S. Zakynthinaki, Y.G. Saridakis / Computer Physics Communications ••• (••••)•••–••• 3 149 250 351 452 553 654 755 856 957 10 58 11 59 12 60 13 61 14 62 15 63 16 64 17 65 18 66 19 67 20 68 21 69 22 70 23 71 24 72 25 73 26 74 27 75 28 76 29 77 30 78 31 79 32 80 33 81 34 82 35 83 36 84 37 85 38 86 39 87 40 88 41 89 42 90 43 91 44 92 45 93 46 94 47 95 48 96 Fig. 1. A time step of the simulated adaptive optics system. 2. Computer simulations and background 2.1. Atmospheric distortion and adaptive correction Let us assume a point source at infinity. Coherent light rays coming from the point source and passing through the turbulent atmospheric layers, deviate from planarity.The distortedinformationpassesthroughthe circular telescope aperture (ρ, θ) and forms an image of the point star on the image plane (r, ψ). The light distribution of any point on the image plane, in the presence of interfering atmospheric distortions, can be given by the following modified Fresnel–Kirchhoff integral [1,5,12], I(r,ψ)= 1  0 2π  0 eik[δ(ρ,θ)−rρ cos(θ−ψ)]ρdθdρ 2 ,(1) where kis the wave number of the incident wave. The real function δ(ρ,θ) in Eq. (1) includesthe optical effects of both the distorting atmosphere φ(ρ,θ) and the correcting mirror φc(ρ, θ). In the absence of any distortions δ(ρ,θ) =0, while in general δ(ρ,θ) =φ(ρ,θ)−φc(ρ, θ). (2) Function δ(ρ,θ) is of great interest in adaptive optics. However, it is not directlymeasurable,as far as optical wavelengths are concerned. The use of an alternative method for the estimation of the image quality and the corrections needed would therefore speed up the process of image optimization in any adaptive optics system. In the adaptive optics system presented here, the problem of minimization of δ(ρ,θ) (maximum possible elimination of atmospheric distortions) is solved via the maximization of the total light intensity passing through the mask. As a result of keeping the light intensity inside the masking area at maximum levels, less spreading of the light distribution is achieved. As- ARTICLE IN PRESS UNCORRECTED PROOF S0010-4655(02)00681-1/FLA AID:2511 Vol.•••(•••) P.4 (1-19) ELSGMLTM(COMPHY):m3 2002/09/24 Prn:1/10/2002; 16:39 cpc2511 by:PS p. 4 4M.S. Zakynthinaki, Y.G. Saridakis / Computer Physics Communications ••• (••••)•••–••• 149 250 351 452 553 654 755 856 957 10 58 11 59 12 60 13 61 14 62 15 63 16 64 17 65 18 66 19 67 20 68 21 69 22 70 23 71 24 72 25 73 26 74 27 75 28 76 29 77 30 78 31 79 32 80 33 81 34 82 35 83 36 84 37 85 38 86 39 87 40 88 41 89 42 90 43 91 44 92 45 93 46 94 47 95 48 96 suming a central masking area of radius R, the total light intensity inside this area can be represented by the integral Stip tilt ≡ R  −R 2π  0 I(r,ψ)rdψdr, (3) which can be used to evaluate the cost function Stip tilt that describes the simulated system. For a tip-tilt system, i.e. a system that only corrects the two lowest sources of aberrations, the correcting phase φc(ρ, θ) is a functionof theslope of the mirror’s surface and is for simplicity expressed in Cartesian coordinates as φc(x, y) =a1x+a2y. (4) The variables xand yrepresent the phase gradients in the x-andy-direction, while the tip and tilt components a1,a2are variables determiningthe slope of the mirror. Changes in the tip and tilt components lead to a movementof the centroid of the star’s image. Fig. 2 shows the image of a point star, generated by means of Eq. (1), under the assumption that no atmospheric distortion is present. Such an image consists of a number of diffraction rings around a central light disk referred to as the airy disk of the image. Theoreticallythe angular width of the airy disk can be calculated using the relation δθ ≈1.22 λ D≈0.25λ(µm) D(m)arcsec,(5) where λand Ddenote the wavelength of the incident light and the diameter of the telescope aperture, respectively. Obviously the quantity δθ gives a theoretical limit of the resolution of a telescope of a given diameter D, since two neighboring sources are distinguishable only when placed at a minimum distance δθ apart. In our simulation, the light was considered monochromatic with a wavelength of λ=0.5µm. The spatial resolution of ground based telescopes with diameter greater than D>30 cm is, however, considerably less than the diffraction limit given by Eq. (5), even under good weather conditions. For visible wavefronts the light distribution typically appears spread on an area 10–20 times larger than the one which is theoretically expected. If, for given seeing, (a) (b) Fig. 2. The image of a point star—no disturbing atmosphere is present, (a) as shown on the image plane and (b) irradiance distribution. the angular resolution of a star is δθ, then Eq. (5) can be re-written in the form δθ ≈1.22 λ r0,(6) wherer0istheFried’scoherencelengthofatmospheric turbulence[1,13]. For given λ, the parameter r0experimentally represents the maximum telescope diameter that leaves light information unaffected by turbulence. Theoretically r0gives the characteristic length above which turbulent atmospheric movement remain statistically uncorrelated. The value of r0is time-varying ARTICLE IN PRESS UNCORRECTED PROOF S0010-4655(02)00681-1/FLA AID:2511 Vol.•••(•••) P.5 (1-19) ELSGMLTM(COMPHY):m3 2002/09/24 Prn:1/10/2002; 16:39 cpc2511 by:PS p. 5 M.S. Zakynthinaki, Y.G. Saridakis / Computer Physics Communications ••• (••••)•••–••• 5 149 250 351 452 553 654 755 856 957 10 58 11 59 12 60 13 61 14 62 15 63 16 64 17 65 18 66 19 67 20 68 21 69 22 70 23 71 24 72 25 73 26 74 27 75 28 76 29 77 30 78 31 79 32 80 33 81 34 82 35 83 36 84 37 85 38 86 39 87 40 88 41 89 42 90 43 91 44 92 45 93 46 94 47 95 48 96 (a) (b) (c) (d) Fig. 3. Speckle samples. and depends on the severity of turbulence and the wavelength used. Evidently a typical telescope diameter Dis larger than the coherencelength r0. As a result, the telescope pupil can be thought of as a dense collection of Kmicroscopic lenses of diameter r0, which densely cover the pupil plane. This number Kcan be approximated by K≈D r02 . The Kmicroscopic tilted lenses result in a random interference of light rays and the formation of Kmicroscopic images on the image plane. Each microscopic image of the star is of diameter λ/D approximately, while the displacement of each image from the central image area is random. These microscopic images of the point source are referred to as speckles (see Fig. 3). The speckle pattern moves as a whole on the image plane. Due to its randommotion and changes in shape, the time integrated, long exposure image, appears spread on a large area on the image plane. It is worth noticing that the speckles observed in the images of Fig. 3 are the result of bad atmospheric seeing as, for the particular simulation, there was D/r0= 17.3. ARTICLE IN PRESS UNCORRECTED PROOF S0010-4655(02)00681-1/FLA AID:2511 Vol.•••(•••) P.6 (1-19) ELSGMLTM(COMPHY):m3 2002/09/24 Prn:1/10/2002; 16:39 cpc2511 by:PS p. 6 6M.S. Zakynthinaki, Y.G. Saridakis / Computer Physics Communications ••• (••••)•••–••• 149 250 351 452 553 654 755 856 957 10 58 11 59 12 60 13 61 14 62 15 63 16 64 17 65 18 66 19 67 20 68 21 69 22 70 23 71 24 72 25 73 26 74 27 75 28 76 29 77 30 78 31 79 32 80 33 81 34 82 35 83 36 84 37 85 38 86 39 87 40 88 41 89 42 90 43 91 44 92 45 93 46 94 47 95 48 96 The Strehl Ratio (SR)isdefinedasarestoration criterion,used to estimate the sharpnessof the restored images, SR =I(0) M(0),(7) with M(r) being the irradiance distribution of the undistorted image. The value of SR gives a measure of the central light intensity of a restored image, in comparison to the central intensity peak of the optimal undistorted image (see Fig. 2). Evidently the maximum value of the Strehl Ratio is SR =1. For the simulation of the point star (see Fig. 2) we normalized the peak of the irradiance distribution to unity. Thus, for an undistorted light distribution M(r),thereis SR =M(0)=1. Under this assumption, the value of the SR of any undistorted image can be calculated by the relation SR =I(0) (see Eq. (7)). 2.2. Atmospheric simulation The atmospherically distorted wavefront at the telescope pupil is modeled in terms of Kolmogorov theory [6], which predicts the statistical properties of the fluctuations of the refractive index and leads to equations which describe the turbulence induced by the thin layers within the atmosphere. For the simulation of the atmospheric distortion, a phase screen, i.e. a two-dimensional distribution of phase fluctuations satisfying a Kolmogorov spectrum, is introduced to the optical system at each time step of the optimizing algorithm (recall Fig. 1). Lane, Glindemann and Dainty [7]usedtherandom mid-point displacement method, previously applied to the computer generation of artificial landscapes, for a Kolmogorov turbulent screen generation. The same technique is applied in the present study as well (see also [1,5]), for the generation of a 65 ×65 Kolmogorov sampling grid. The algorithm is based on making a coarsely sampled approximation to the turbulent atmosphere and the subsequent refinement of the successive smaller localized regions. For more details see [1,5,7]. Let us assume the nth time step of the optimization process. The numerical generation of the time evolving turbulent phase screen φ(ρ,θ)(n) begins with four independent starting samples α(n),β(n),γ(n) and δ(n), generated by a combination of the Gaussian variables Rc(n) and Rd(n),ofvarianceσcand σd, respectively [7]: α(n) =Rc(n) α+0.5Rd(n) αδ , β(n) =Rc(n) β+0.5Rd(n) βγ , γ(n) =Rc(n) γ−0.5Rd(n) βγ , δ(n) =Rc(n) δ−0.5Rd(n) αδ (8) with σcand σdsatisfying the relations 2σ2 c+σ2 d 2=6.88(D/r0)5/3, 2σ2 c+σ2 d=6.88(√2D/r0)5/3. (9) The time dependence of the Gaussian variables Rc(n) and Rd(n) which provide the four starting samples α(n),β(n),γ(n) and δ(n) of the turbulent phase screen, as described by Eqs. (8), is directly related to local atmospheric conditions [1]: Assume the number x(n) which moves across the interval [0,1]at a rate which depends upon the rate of change of the atmospheric conditions. The Gaussian variables Ri(n) j,wherei=c,d and j=α, β, γ, δ, αδ, βγ ,asin Eqs. (8), are then calculated by the relation Ri(n) j=±σi−2.0·log(x(n)). For a more detailed description of the simulation of the time evolving phase screens, see [1]. Having calculated the four starting samples α(n), β(n),γ(n) and δ(n), a new point m(n) can be generated by linear interpolation and the addition of a displacement $(n): m(n) =α(n) +β(n) +γ(n) +δ(n) 4+$(n).(10) If h(n) is defined to be the distance between two adjacent samples on the sampling grid, the variance of the displacement $(n) can be required [7]tobe equal to σ2 $=0.6091h5/3. It is worth noticing that, as the sampling grid becomes finer, the variance of the displacementis also reduced.The procedurecontinues by another linear interpolation between the corner points and again by the addition of displacements. As in [7], the calculation of a new point by use of a point at the edge of the grid requires the addition of ARTICLE IN PRESS UNCORRECTED PROOF S0010-4655(02)00681-1/FLA AID:2511 Vol.•••(•••) P.7 (1-19) ELSGMLTM(COMPHY):m3 2002/09/24 Prn:1/10/2002; 16:39 cpc2511 by:PS p. 7 M.S. Zakynthinaki, Y.G. Saridakis / Computer Physics Communications ••• (••••)•••–••• 7 149 250 351 452 553 654 755 856 957 10 58 11 59 12 60 13 61 14 62 15 63 16 64 17 65 18 66 19 67 20 68 21 69 22 70 23 71 24 72 25 73 26 74 27 75 28 76 29 77 30 78 31 79 32 80 33 81 34 82 35 83 36 84 37 85 38 86 39 87 40 88 41 89 42 90 43 91 44 92 45 93 46 94 47 95 48 96 (a) (b) (a) (b) (a) (b) Fig. 4. Samples of turbulent phase screens (a) two-dimensional projection, (b) functional form. ARTICLE IN PRESS UNCORRECTED PROOF S0010-4655(02)00681-1/FLA AID:2511 Vol.•••(•••) P.8 (1-19) ELSGMLTM(COMPHY):m3 2002/09/24 Prn:1/10/2002; 16:39 cpc2511 by:PS p. 8 8M.S. Zakynthinaki, Y.G. Saridakis / Computer Physics Communications ••• (••••)•••–••• 149 250 351 452 553 654 755 856 957 10 58 11 59 12 60 13 61 14 62 15 63 16 64 17 65 18 66 19 67 20 68 21 69 22 70 23 71 24 72 25 73 26 74 27 75 28 76 29 77 30 78 31 79 32 80 33 81 34 82 35 83 36 84 37 85 38 86 39 87 40 88 41 89 42 90 43 91 44 92 45 93 46 94 47 95 48 96 a displacement η(n) with variance σ2 η=0.4471h5/3. On the left edge, for example, the new point m(n) αβ is formed according to the rule m(n) αβ =α(n)+β(n) 2+η(n). For a more detailed description anda discussion of this interpolation/displacement procedure see [1,5,7]. Fig. 4 shows three samples of simulated turbulent screens. For the purposes of the present study, phase screens simulating turbulence were calculated for a wide range of atmospheric conditions, starting from D/r0=3.25 to D/r0=43.3. 3. ALOPEX stochastic optimization For the purpose of optimization, let f≡f(x 1,x2, ...,x N;y1,y2,...,y M)denote the cost function,that is the function that describes the system. Function fdepends on a set of Nparameters x1,x2,...,x N, called the control variables, and a set of Mparameters y1,y2,...,y M, that are not under control, and may be internal or external parameters of the system. As the set of the Mparameters y1,y2,...,y Mare not under control, we can re-write the cost function fusing the control variables x1,x2,...,x Nonly, f≡f(x 1,x2,...,x N;y1,y2,...,y M) ≡f(x 1,x2,...,x N). The stochastic optimization algorithm used for the maximization of the cost function is the socalled ALOPEX (ALgorithm Of Pattern EXtraction) algorithm. It was originally devised in [8,9]forthe purposeofexperimentallydeterminingreceptivefields of individual neurons in the visual pathway. We modified the ALOPEX optimization in [5,10,11]and introducednew versions inaddition to the ones already known. The ALOPEX process operates as follows: •At each time step of the procedure, the variables that determine the cost function are changed simultaneously by small increments, and the cost function is re-evaluated. •The changesin a variabledependstochastically on the change of the cost function and the change in that variable over the previous two time steps. •Since the changes are cumulative, the value of each variable reflects at all times the dependence of the cost function to changes in that variable over all previous time steps. •The process is guided by two free parameters, which are the size of the increments and the degree of stochasticity, i.e. the amplitude of noise. For the purposes of the present simulation, the cost function fof the system equates to the function Stip tilt, defined by the integral shown in Eq. (3) and referring to the total light distribution measured on the masking area of the image. From Eqs. (3), (1), (2) and (4), it can be easily seen that the cost function of the system Stip tilt is a function of the tip and the tilt components of the mirror, a1and a2. The restoration of the centroidof the speckled imagesis to be achieved by optimization of the cost function Stip tilt(a1,a2),via the appropriate changes in the slope of the mirror, i.e. via changes in the control variables a1and a2. An optimization of Stip tilt(a1,a2)in real time will in turn improve the sharpness of the final, long exposure image. Let a(n) ibe the value of the ith control variable after the nth time step and let S(n)(a(n) 1,a(n) 2)be the value of the cost function. According to the version of the ALOPEX stochastic optimization algorithm presented in [1,11], given the (random) initial conditions a(0) i, a(1) iand a(2) i(i=1,2), the value of the ith control variable is evaluated, at each time step, according to the rule: a(n+1) i=a(n) i+c)a(n) i )S(n) |)S(n−1)|+g(n) i, i=1,2,n⩾3,(11) where )a(n) i=a(n) i−a(n−1) iand )S(n) =S(n) − S(n−1). The noise terms giare essential ingredients in the process, as they provide the agitation necessary to drive the process [1,5,8–11]. The dynamics of the process dependsstronglyon the mean amplitudeof the giterms. With use of the above version of the ALOPEX stochastic optimization algorithm, the value of parameter cremainsconstant duringthe optimizationprocess. For the results presented in the present paper, there is c=0.4. For a more detailed discussion on the version of ALOPEX stochastic optimization algorithm presented by Eq. (11) see [1,11]. ARTICLE IN PRESS UNCORRECTED PROOF S0010-4655(02)00681-1/FLA AID:2511 Vol.•••(•••) P.9 (1-19) ELSGMLTM(COMPHY):m3 2002/09/24 Prn:1/10/2002; 16:39 cpc2511 by:PS p. 9 M.S. Zakynthinaki, Y.G. Saridakis / Computer Physics Communications ••• (••••)•••–••• 9 149 250 351 452 553 654 755 856 957 10 58 11 59 12 60 13 61 14 62 15 63 16 64 17 65 18 66 19 67 20 68 21 69 22 70 23 71 24 72 25 73 26 74 27 75 28 76 29 77 30 78 31 79 32 80 33 81 34 82 35 83 36 84 37 85 38 86 39 87 40 88 41 89 42 90 43 91 44 92 45 93 46 94 47 95 48 96 Characteristics of the method, such as •effectiveness and speed of convergence in real time, •no required knowledge of the dynamics of the system or of the functional dependenceof the cost function on the control variables, •easy and cost effective implementation are the main factors that may lead to a successful implementation of the suggested method to real astronomical systems. 4. Numerical experimentations 4.1. Definitions We define the following parameters: •Parameter µ, µ=masking area diameter airy disk diameter investigates the performance of the method, as it pertains to the diameter of the masking area. Recall that the airy disk is the central disk of the irradiance distribution of the image of the point source (see Fig. 2). •Parameter µ1, µ1≡Drestored m Ddistorted m estimates the restoration of the light intensity spread around the central maximum. Drestored mis the diameter of the light distribution at half the intensity maximum of the restored long exposure image (note that this parameter can be found in the bibliography as the FWHM (Full Width Half Maximum). Ddistorted mis the diameter of the light distribution at half the intensity maximum of the distorted long exposure image. •Parameter µ2, µ2≡SRrestored image SRdistorted image compares the SR values of the restored long exposure images to that of the distorted long exposure image. A combination of parameters µ1and µ2is believed to give a sufficient measure of the restoration of the image sharpness. 4.2. The restoration process The restoration process is assumed to start with the tip-tilt mirror aligned parallel to the image and pupil plane (i.e. for n=1, φ(1) c(ρ, θ) =0). The nth time step (where n⩾3, as in Eq. (11)) starts with the calculation of the distorting phase φ(n)(ρ, θ) (i.e. the nth turbulent phase screen, see Section 2.2). As described by Eq. (2), the evaluation of the function δ(n)(ρ, θ) involves the subtraction of the correcting phase φ(n−1) c(ρ, θ), which corresponds to the mirror’s position at the previous time step n−1, from the distorting phase φ(n)(ρ, θ). The total light intensity on the masking area is then measured and the cost function S(n) tip tilt is evaluated (see Eq. (3)). This value of the cost function is given as feedback for ALOPEX, see Eq. (11), for the calculation of the new position of the mirror (as described in the previous sections, ALOPEX calculates the coefficients a(n) 1and a(n) 2by which the correcting phase φ(n) c(ρ, θ) is evaluated). We recall (see Eq. (11)) that ALOPEX optimization requiresinformationfrom the two previoustime steps, i.e. the values of S(n−1) tip tilt ,S(n−2) tip tilt ,a(n−1) 1and a(n−1) 2. Fig. 5 illustrates the optimization process for the three successive time steps n−2, n−1andn(n⩾3). The success of the method is demonstrated through Fig. 6, which shows random samples of speckles, taken on the image plane as the exposure time progressed. The images shown in (a) are the result of the distorting atmosphere (ALOPEX optimization turned off), while the images in (b) are the result of a replay of the same time sequence of the turbulent phase screens of (a) (i.e they are subjected to the same distortion), but with correction by help the tip-tilt mirror (ALOPEX optimization turned on). We note that the three pairs of images shown in Fig. 6 correspond to the short exposure images of time steps n=50, n=500 and n=1000 of a 1000-step optimization process. For these images, parameter µwas kept constant and equal to µ=1.875 (see section below) and the severity of turbulence was D/r0=8.6 (compare with Fig. 3). By simple observation of the images of Fig. 6, one may verify that, although turbulence tends ARTICLE IN PRESS UNCORRECTED PROOF S0010-4655(02)00681-1/FLA AID:2511 Vol.•••(•••) P.16 (1-19) ELSGMLTM(COMPHY):m3 2002/09/24 Prn:1/10/2002; 16:39 cpc2511 by:PS p. 16 16 M.S. Zakynthinaki, Y.G. Saridakis / Computer Physics Communications ••• (••••)•••–••• 149 250 351 452 553 654 755 856 957 10 58 11 59 12 60 13 61 14 62 15 63 16 64 17 65 18 66 19 67 20 68 21 69 22 70 23 71 24 72 25 73 26 74 27 75 28 76 29 77 30 78 31 79 32 80 33 81 34 82 35 83 36 84 37 85 38 86 39 87 40 88 41 89 42 90 43 91 44 92 45 93 46 94 47 95 48 96 of the centroid of the speckles. This mean velocity of speckles, ¯vspeckles, may be expressed in terms of the average diameter of the speckle pattern, Dspeckles,and in units of Dspeckles/(time step). 4.4.1. Keeping the mean noise amplitude constant The image restoration process was studied for various values of ¯vspeckles. For the results presented in Table 3, ¯vspeckles varies from 0.231 ·10−2to 4.709 · 10−2Dspeckles/(time step), while the mean noise amplitude is constant and equal to gi=1.964 · 10−2Dspeckles. Table 3 Restoration parameters versus the mean image velocity ¯vspeckles[Dspeckles/(time step)]µ1µ2 4.709 ·10−20.984 1.023 3.118 ·10−20.943 1.098 2.458 ·10−20.937 1.129 2.310 ·10−20.913 1.196 2.198 ·10−20.815 1.294 2.117 ·10−20.775 1.362 2.024 ·10−20.762 1.364 1.939 ·10−20.682 1.571 1.861 ·10−20.504 1.802 1.539 ·10−20.496 1.789 1.159 ·10−20.494 1.803 0.927 ·10−20.499 1.831 0.772 ·10−20.527 1.848 0.661 ·10−20.497 1.923 0.578 ·10−20.498 1.906 0.514 ·10−20.499 2.002 0.462 ·10−20.499 2.057 0.330 ·10−20.499 2.048 0.231 ·10−20.499 2.053 Figs. 11 and 12 illustrate the results presented in Table 3. One can observe that, for ¯vspeckles <1.939 ·10−2 Dspeckles/(time step)the mean amplitude of noise gi implemented can be kept optimally constant, even though its average numerical values are higher than the average values of the displacement of the centroid of the speckles at each time step. Thiscan be explained due to the intrinsic properties of ALOPEX stochastic optimization, see Eq. (11). For velocities ¯vspeckles > 1.939 ·10−2Dspeckles/(time step)the method fails to satisfactory restore the image. Note, however, that such velocities correspond to fast rate of change of the position of the centroid of the speckles and hence represent worse than average weather conditions. We believe that in such cases of strong atmospheric distortion, the problem of satisfactory image sharpness restoration will be overcome by appropriately changing the mean amplitude giof the noise. The section that follows justifies our hypothesis. 4.4.2. Changing the mean noise amplitude We investigate the effect of the implementation of increasing mean noise amplitudes giin the case where the mean image velocity ¯vspeckles exceeds the limit of 1.939 ·10−2Dspeckles/(time step).Asshown in Table 4, which summarizes the results of the simulation, there is an optimal mean noise amplitude gifor each ¯vspeckles, resulting in a satisfactory restoration of the image sharpness. In Figs. 13 and 14, which illustrate the data of Table 4, the results of implementing a mean noise amplitude gi, which increases as ¯vspeckles increases, are presentedwith solid lines.Dashed lines correspond Table 4 Restoration parameters versus mean image velocity and noise amplitude ¯vspeckles[Dspeckles/(time step)]gi[Dspeckles]µ1µ2 4.709 ·10−24.622 ·10−20.575 1.925 3.118 ·10−23.343 ·10−20.516 1.927 2.458 ·10−22.910 ·10−20.526 1.992 2.310 ·10−22.857 ·10−20.583 2.021 2.198 ·10−22.619 ·10−20.498 1.982 2.117 ·10−22.245 ·10−20.542 2.052 2.024 ·10−22.245 ·10−20.501 2.041 1.939 ·10−22.245 ·10−20.518 2.039 1.861 ·10−21.964 ·10−20.504 1.802 ARTICLE IN PRESS UNCORRECTED PROOF S0010-4655(02)00681-1/FLA AID:2511 Vol.•••(•••) P.17 (1-19) ELSGMLTM(COMPHY):m3 2002/09/24 Prn:1/10/2002; 16:39 cpc2511 by:PS p. 17 M.S. Zakynthinaki, Y.G. Saridakis / Computer Physics Communications ••• (••••)•••–••• 17 149 250 351 452 553 654 755 856 957 10 58 11 59 12 60 13 61 14 62 15 63 16 64 17 65 18 66 19 67 20 68 21 69 22 70 23 71 24 72 25 73 26 74 27 75 28 76 29 77 30 78 31 79 32 80 33 81 34 82 35 83 36 84 37 85 38 86 39 87 40 88 41 89 42 90 43 91 44 92 45 93 46 94 47 95 48 96 Fig. 13. Intensity spread restoration versus mean speckle velocity, in relation to the noise amplitude: ···keeping gi’s constant, — changing the gi’s. Fig. 14. SR restoration versus mean speckle velocity, in relation to the noise amplitude: ···keeping gi’s constant, — changing the gi’s. to the data of Figs. 11 and 12 (Table 3). The success of the implementation is apparent. Fig. 15 shows the necessary increase in the noise amplitude when the mean image velocities ¯vspeckles become largerthan 1.939·10−2Dspeckles/(time step),i.e. under weather conditions which are worse than average. As can be observed,the mean noise amplitudecan be assumedto be linearly dependedon themean veloc- ARTICLE IN PRESS UNCORRECTED PROOF S0010-4655(02)00681-1/FLA AID:2511 Vol.•••(•••) P.18 (1-19) ELSGMLTM(COMPHY):m3 2002/09/24 Prn:1/10/2002; 16:39 cpc2511 by:PS p. 18 18 M.S. Zakynthinaki, Y.G. Saridakis / Computer Physics Communications ••• (••••)•••–••• 149 250 351 452 553 654 755 856 957 10 58 11 59 12 60 13 61 14 62 15 63 16 64 17 65 18 66 19 67 20 68 21 69 22 70 23 71 24 72 25 73 26 74 27 75 28 76 29 77 30 78 31 79 32 80 33 81 34 82 35 83 36 84 37 85 38 86 39 87 40 88 41 89 42 90 43 91 44 92 45 93 46 94 47 95 48 96 Fig. 15. Optimal noise amplitude versus speckle velocity: ···keeping gi’s constant, — changing the gi’s. Fig. 16. Linear regression of optimal noise amplitude on the mean speckle velocity. ity of the speckles on the image plane, ¯vspeckles.More specifically, a linear regression of the mean noise amplitude gion the mean velocity of speckles, ¯vspeckles, can be calculated: gi=α·¯vspeckles +β, where α=0.883(time step)/Dspeckles and β=0.552· 10−2. The correlation coefficient for the above for- ARTICLE IN PRESS UNCORRECTED PROOF S0010-4655(02)00681-1/FLA AID:2511 Vol.•••(•••) P.19 (1-19) ELSGMLTM(COMPHY):m3 2002/09/24 Prn:1/10/2002; 16:39 cpc2511 by:PS p. 19 M.S. Zakynthinaki, Y.G. Saridakis / Computer Physics Communications ••• (••••)•••–••• 19 149 250 351 452 553 654 755 856 957 10 58 11 59 12 60 13 61 14 62 15 63 16 64 17 65 18 66 19 67 20 68 21 69 22 70 23 71 24 72 25 73 26 74 27 75 28 76 29 77 30 78 31 79 32 80 33 81 34 82 35 83 36 84 37 85 38 86 39 87 40 88 41 89 42 90 43 91 44 92 45 93 46 94 47 95 48 96 mula is r=0.97849. Fig. 16 shows the above linear regression. 5. Conclusions The present work investigates the optimization of an image of a point star using a computer simulated tip-tilt adaptive optics system. The optimization of the long exposure images is achieved by use of a stochastic optimizationalgorithm.The tip-tilt mirroris rotated in such a way that the position of the centroid of the short exposure images is restored, in real time, inside the central area of the image plane (the masking area). The optimal diameter of the masking area is found to be approximately twice the airy disk of the point star. We also found that for optimal image restoration and under conditions of good or medium seeing, the mean amplitude of the required noise in the algorithm can be kept constant during the optimization process. For worse than average weather conditions, the optimal mean noise amplitude and the mean velocity of the short exposure images are found to be linearly depended. By use of the optimization method presented here, the reduction of the distorting effects of the atmosphere are achieved in real time in a way which is far simpler, inexpensive and easy to implement than the currently used methods. Our optimizing system, being able to provide images of high quality and sharpness for all ground-based telescopes, is a highly competitive application to a very complicated and demanding problem. References [1] M.S. Zakynthinaki, Stochastic optimization for adaptive correction of atmospheric distortion in astronomical observation, PhD thesis, Technical University of Crete, 2001. [2] H.W. Babcock, PASP 65 (1953) 229. [3] H.W. Babcock, J. Optical Soc. Amer. 48 (1958) 500. [4] R.J. Noll, J. Optical Soc. Amer. 66 (1976). [5] Y.G. Saridakis, M.S. Zakynthinaki, T.E. Kalogeropoulos, Internat. J. Appl. Sci. Comput. 5 (3) (1999) 252. [6] A.N. Kolmogorov, in: S.K. Friedlander, L. Topper (Eds.), Turbulence, Interscience, New York, 1965. [7] R.G. Lane, A. Glindemann, J.C. Dainty, Waves in Random Media 2 (1992) 209. [8] E. Harth, E. Tzanakou, Vision Res. 14 (1974) 1475. [9] T. Tzanakou, R. Michalak, E. Harth, Biol. Cybernet. 35 (1979) 161. [10] T.E. Kalogeropoulos, Y.G. Saridakis, M.S. Zakynthinaki, Comput. Phys. Comm. 99 (1997) 255. [11] Y.G. Saridakis, M.S. Zakynthinaki, in: Proc. of the 3rd Hellenic–European Conference on Mathematics and Informatics, LEA, Athens, 1996, p. 251. [12] R.K. Tyson, Principles of Adaptive Optics, Academic Press, 1991. [13] D.L. Fried, J. Optical Soc. Amer. 55 (1965).