Lunar Scintillometer Array for Atmospheric Turbulence Measurement by Marie-Andree Dumais B.Eng., Ecole Polytechnique de Montreal, 2003 A THESIS S U B M I T T E D IN P A R T I A L F U L F I L M E N T O F T H E R E Q U I R E M E N T S F O R T H E D E G R E E O F Master of Science in The Faculty of Graduate Studies (Astronomy) The University Of British Columbia February 2006 © Marie-Andree Dumais 2006 11 Abstract The performance of large telescopes is directly related to the quality of the site. Of particular importance is the turbulence profile within the lower few hundred meters of the atmosphere. We have developed a new technique for measuring this profile, and the Fried parameter ro, by measuring spatial correlations in the scintillation of moonlight. Our lunar scintillometer array uses 12 photodiodes to measure the intensity fluctuation covariance function, and from this we obtain the index of refraction profile C2N to a height of about 200 m with a vertical resolution of order 2%. This paper describes the instrument, the technique, and results of initial tests. Contents Abstrac t 1 1 Contents u i Lis t of Tables v List of Figures • v i Acknowledgements v n 1 Introduct ion 1 2 T h e o r y 2 2.1 Optical path of fluctuations 4 2.2 Scintillation of light from an extended source 6 2.3 Fried parameter t 8 3 O t h e r technologies 9 3.1 Differential image motion monitor 9 3.2 MASS 10 3.3 SCIDAR 11 3.3.1 Principle 11 3.3.2 Generalized SCIDAR 12 3.3.3 S L O D A R 12 3.4 S O D A R 12 4 Des ign 17 4.1 Overview 17 4.2 Golomb Ruler 17 4.3 Photodiodes 22 4.4 Electrical circuit . ; 22 5 Observations 24 6 D a t a Analys i s " 26 6.1 Profile 26 6.2 Downhill simplex 28 6.3 Simulated Annealing 28 Contents 6.4 Monte Carlo 29 6.5 Lagrange multiplier 30 6.6 Simulations 30 6.7 Turbulence profiles from data 61 7 Conc lus ion 67 B ib l iography • • 68 V List of Tables 4.1 Photodiode positions 20 4.2 Separation Calculation 22 6.1 Profile Descriptions 30 vi List of Figures 2.1 Geometry of the cone of light 6 3.1 D I M M principle 10 3.2 M A S S principle . . 11 3.3 SCIDAR principle 12 3.4 Binary correlation .. . 14 3.5 SCIDAR Acquisition Plane '. . . . 15 3.6 generalized SCIDAR principle 15 3.7 S L O D A R principle 16 3.8 S O D A R principle 16 4.1 Instrument 17 4.2 Lunar scintillometer principle 18 4.3 Golomb 3 marks 19 4.4 Golomb 4 marks . 19 4.5 Picture of photodiodes 22 4.6 Block diagram 23 5.1 Labview Software 24 6.1 Data 27 6.2 Simplex downhill 32 6.3 Simulation models . . . 33 6.4 Simulations 60 6.5 Turbulence profiles recovered I 62 6.6 Turbulence profiles recovered II 63 6.7 Turbulence profiles recovered III 64 6.8 Turbulence profiles recovered IV 65 6.9 Fried parameter 66 V l l Acknowledgements The first people to thank are my supervisor Paul Hickson for bringing me in this project and this whole new field that is site testing and his wife Suzanne for keeping us awake during observations with tea and goodies. I would like to thank my parents and my brother for supporting me in my projects even if they don't exactly understand what I am doing. I also thank my roommates Lara and Lateef for all those moments shared in our crazy trips in the mountains and in the comfort of our apartment. While all is the center of the universe, Lateef exclaims, "that makes sense to me!" Thanks to Tom Felton for spending so much time fixing the electronics of the instrument. Finally, a special thank to Maryse and Cedrik, friends from my undergraduate days in Montreal, for all those conversations relating to physics or life challenges. 1 Chapter 1 Introduction The resolution of ground-based telescopes is limited by the turbulence in the atmosphere which causes image degradation. In order to reduce the effects of turbulence, sites are chosen to be at high altitude. Corrections are sometimes made by adaptive optics to compensate the atmospheric turbulence. Absorp-tion, scattering and scintillation are effects of the Earth's atmosphere. Scat-tering, caused by dust and molecules, will induce an intrinsic brightness in the sky background. The main source of scattered light is artificial, such as street lighting. Hence, a good site must be far from built-up areas. Light is also absorbed in the atmosphere due to molecular absorption by atmospheric gases. Visible light intensity is reduced by propagation through the atmosphere. Infrared radiation is absorbed by water vapour. Hence, a site < must be at high altitude to reduce air paths and to keep the water content as low as possible. Seeing is characterized by the change in the wavefront due to the variation of density through the atmosphere. Seeing results mainly from low level turbulence from the layer above the ground and high layers in the atmosphere. The ground-layer is affected by the ground structure and landscape. Thus, a good site requires a steady atmosphere to improve seeing. The choice of a site is critical for a large telescope as a lot of money is involved. It is also useful for adaptive optics to have a good idea of the at-mospheric turbulence profile for a certain site. It is then easier to build the adaptive optics system. Atmospheric turbulence can be characterized by three main parameters: the full width half maximum (FWHM) of the image of a star, the Fried parameter and the refractive index fluctuations [C2N) profile. The seeing is defined by the blurring of a source image. The Fried parameter ro gives an idea of the resolution of the telescope. It's the diameter of a circle over which the root mean square (rms) phase variation is 1 radian. This parameter describes the limit for a telescope diameter. For a diameter bigger than ro, the primary effect is the blurring of the instantaneous image. On the other hand, for a telescope diameter smaller than ro, the main effect is motion of the image that blurs the image during long exposures. Turbulence can be caused by the dome and the telescope as well the surface boundary layer and planetary boundary layer. The ground layer and high at-mosphere are the main concerns for site testing. The effects from the dome and the telescope can be minimized by careful design. 2 Chapter 2 Theory The turbulence that causes image degradation is divided into four categories: turbulence associated with the dome and the telescope, turbulence in the sur-face boundary layer or due to ground convection, turbulence in the planetary boundary layer or associated with orographic disturbances and turbulence in the tropopause or above. The first type is due to the telescope and the dome can easily be avoided as it is due to temperature gradients. Atmospheric seeing results from refractive index fluctuations. The refractive index of air is a function of the temperature, pressure and the concentration of water vapour [1]: 77 fi x 1 f l - 6 v n(X,P,T,u) - 1 = " - 0 T U (1 + 7.52 x 1 0 " 3 A - 2 ) ( P + 4810-). (2.1) P and T a r e in Pascal (Pa) and Kelvin (K) respectively and A is in nm. According to Barletti et al. [2], below 4 km, the C% profile is mostly affected by. orographic disturbances and above, the behaviour of the turbulence is inde-pendent of the location. A flow is characterized as turbulent when its Reynolds number (Re) is greater than a critical value. Re = VOLO/I/Q varies with the geometrical structure according to the characteristic velocity VQ, the character-istic size of the flow LQ and the kinematic viscosity of the fluid VQ. According to Kolomogorov [3], the kinetic energy of large scale motions is transferred to smaller and smaller scale motions for a fully developed turbulence flow, assum-ing motions on a small scale are isotropic. When the flow becomes laminar (Re decreasing to a small value), the kinetic energy dissipates into heat by viscous friction. The rate eo of viscous dissipation must equal the rate of production of turbulent energy for a stationary state. Therefore, the velocity V of motions at scale L depends only on L and eo for energy production and dissipation. Hence, by a dimensional reasoning, the relation V oc e ^ L 1 / 3 is obtained. Looking at the spectral analysis of the kinetic energy to express with the modulus K of the wave vector K, the energy E(n)dK between n and K + dK is proportional to V 2(re). Moreover, the characteristic length is inversely pro-portional to the wave vector modulus, L oc Hence, the Kolmogorov law is E(n)dK oc K ~ 2 / 3 or E(K) oc « T 5 / 3 . • (2.2) This law is valid only for the range LQ 1 C K C IQ1 , where Lrj describes the outer scale and IQ, the inner scale. The turbulence begins at the outer scale and the dissipation, at the inner scale. Chapter 2. Theory 3 The total energy at frequency K is related to the three-dimensional power spectrum <!>(«) = $(nx, K V , KZ). For the isotropic case considered here, E(K) = 4 T T K 2 $ ( K ) , (2.3) hence * ( K ) O C K T 1 1 / 3 . (2.4) On the other hand, temperature and humidity are also a function of height in atmosphere. Adding turbulence in the atmosphere creates inhomogeneities of temperature and humidity. The Kolmogorov law applies to the mixing of the air with water vapor or the mixing of cool air with warm air. Therefore, the power spectrum <T?T(«) of temperature fluctuations and the power spectrum $ c ( « 0 of humidity fluctuations are $ T ( K ) oc « T U / 3 and * c(/s) oc K ~ 1 1 / 3 . (2.5) By convention [4], the proportionality constant is written as $ T ( K ) = 0 . 0 3 3 C f .KT 1 1 / 3 (2.6) where Cf. is called the structure constant of temperature fluctuations. As seen above, the refractive index is a function of the temperature and the concentration of water vapour. The fluctuations for the refractive index are n = N-(N), (2.7) dN n = Q f @ (2-8) where 0 represents the temperature fluctuations Q — T — (T). In optical propa-gation for astronomical observations, terms related to water vapour concentra-tion are negligible [3]. Using the equation of the refractive index as a function of the temperature and the pressure, we obtain the value of our partial derivative. This value relates C%, the index structure constant to C\ by r2 -N — 7.9 x ICT 7 p C\. (2.9) Therefore, the power spectrum of the refractive index fluctuations is given by: $ J V ( K ) = 0 . 0 3 3 C ^ K - 1 1 / 3 . ( 2 . 1 0 ) We consider a monochromatic plane waves of wavelength A propagating downward from a star at zenith towards a ground-based observer described by <Sh(x) = \^h(x)\exp(id>h(x)) ( 2 . 1 1 ) where h is the altitude. We normalise to unity the unperturbed complex ampli-tude. The phase 4>h(x) is changed by the index fluctuations inside the layer. As Chapter 2. Theory 4 optical wavelengths observed are much smaller than the scale of the observed wavefront perturbations, the Fresnel approximation is used [3] where * denotes'convolution product. 2.1 O p t i c a l p a t h o f fluctuations Fluctuations of the phase (po{x) have an impact on the shape of the wavefront surface and are related to the angle of arrival of light. The radiation field is represented by $(x, z) where x is the plane perpendicular to the line of sight and z is the distance along the line of sight from the detector. This can be divided in two real function x i x , z) a n d <f>{x, z) z) = exp[x(a;, z) + i<j>{x, z)}. (2.13) The radiation field is constant at the top of the atmosphere. In order to nor-malize we take 0 = ^ = 0 [5]. Assuming a thin layer 5h at height h above the detector and defining ( as the zenith angle, the distance from the detector is z = / i s e c £ and 5z = Sh sec _. A phase fluctuation is introduced by a turbulent layer 5h. If one assumes the atmosphere is homogeneous everywhere except for a thin layer between h and h + Sh, it will introduce a phase shift 5<j)(x) 4>(x) =k n{x,z)dh. (2.14) Jh The quantity on the right hand side of equation 2.14 is infinitesimally small and proportional to Sh. Hence, the local change is 5\n[V{x,z)]=i6<t>(x,z). (2.15) At the detector, the propagation is described by Fresnel diffraction [3] 5V(x, 0) = 6{x, z) * exp I%-^-\ . (2.16) Xz \ Xz J With the derivative of the expression 2.15 and using the equation 2.16, we obtain ^ n ^ , 0 ) = ^ ^ , , ) ^ e x p ( ^ ) . (2.17) For small fluctuations ??, *(a;,0) ~ V{x,z). (2.18) so, 5\(x, 0) + i8<p(x, 0) = iS4>{x, z) * — exp ( -j- J (2.19) Chapter 2. Theory 5 or 5x(x,0)=54>(x,z)*j-zcos(^^ (2.20) and ^(a;,0) =J0(a: )z)* J - s i n f ^ - V (2.21) Az \ \z J In the plane of the detector [5], the spatial covariance functions are defined by Bx(r) =<x(x)x(x + r) >, (2-22) B+(r) =< 4>(x)4>(x + r) > . (2.23) Their Fourier transforms are the 2-d power spectra Wx(f) = J Bxexp(-2irrf)d2r, (2.24) W*'f) = J B^exp(-27rrf)d2r, (2.25) where / is the 2-d spatial frequency. We substitute the definition of \ a n d 4> m the previous expressions to get SWx(f) = 5W*(f,z)sm2(7r\zf2), (2.26) 6Wf{f) =6W*{f,z)coS2{7r\zf2), (2.27) where 6Wv(f,z) = j <(j){x,z)(l>{x + r,z)>e^{-2mrf)d2r. (2.28) On the other hand, according to Roddier [3], for Kolmogorov turbulence, the phase power spectrum has the form SWn,(f,z) = 0.38 sec (\~2 f-11/3C2N{h)5h. (2.29) Hence, substituting equation 2.29 in the Fourier transform for \ a n d 4> (equa-tions 2.26 and 2.27), the total power spectrum is: Wx{f) =0.38 sec C A " 2 / - 1 1 7 3 x J C2N{h) sin 2(7rA sec (f2)dh, (2.30) W+tf) = 0.38 s e c ( A - 2 / - 1 1 / 3 x j C2N(h) C O S 2 ( T T A secCj2)dh. (2.31) Chapter 2. Theory 6 Figure 2.1: Geometry of the cone of light received by a single photodiode. 2.2 S c i n t i l l a t i o n o f l i g h t f r o m a n e x t e n d e d s o u r c e In our experiment, light comes from a finite source, which is the moon, and falls on a finite size detector. First, we assume a circular active area for the detector of diameter a [5]. For a given cone, the ratio between y the aperture and z the height is a = | = 2 t a n ^ . (2.32) The extended area A(x,z), perpendicular to the line of sight at a height h is: A(x,z)-=u(^—\, (2.33) ,a + az Where H (y ) = { I H\. (2-34) The area of the beam is given by A(z) = (2.35) The normalized aperture function is defined by A(..»>-4$. (^ ) Since the intensity I is proportional to the amplitude squared x'' •> the intensity fluctuation covariance for two detectors separated by a distance r is four times the amplitude covariance fl,(r)=4flx(r). (2.37) Chapter 2. Theory 7 Substituting from the definition of Bx, Bi{r) = 0.38 j C2N(h)K(h,r)dh (2.38) where K{h,r) =&nsecC\-2 J sm2(irXhsecCf2)G(f,hsec()J0(2Trfr)r8/3df. (2.39) G{f,z) represents the window function for which the power spectra are mul-tiplied and integrated. It is the Fourier transform of the aperture function calculated above: G{f,z) = \ J A(x,z)exp{-2mx-f)d2x\2, (2.40) (2.41) •^(nia + az)/)12 [rc(a + az)f] The right hand side of equation 2.41 is integrated from 0 to x — (a + az)/2 in the case of the boundary limits of the function IT (equation 2.34). As the exponential changes between 1 and — 1, / varies between 0 and 1. Hence, the maximal value of the spatial frequency / is given by (a + az)'1. From equation 2.39, we would like to simplify the argument of the sine. Since f<(a + az)-1 (2.42) it follows that nXzf2 < f (2.43) (a + azy The maximal value of the expression of the right hand side can be found by equating to zero its derivative with respect to z ^ P ^ = 0. (2.44) (a + azy The maximum occurs when z = a/a. Substituting gives nXzf2 < (2.45) 4aa In our project, we have A = 500 nm, a = 5.8 mm, a = 2tan(#/2) and 6 = 0.5 deg so, wX/Aaa = 7.759 • 10~ 3 << 1. This quantity is sufficiently small that the sine may be approximated by its argument s i n ( 7 r A / i s e c £ / 2 ) ~ nXhsec(f2. Substituing this in equation 2.39 gives K(h,r) = 87r 3/i 2sec 3C J G{f,hsec()J0(2nfr)f4/3df (2.46) Chapter 2. Theory 8 or K(h,r) = 32irh2 sec3 C a + ahsecQf}2J0{2irfr)r2,Zdf. (2.47) (a + ah sec 2 We replace the integral with (2.48) with r (2.49) s = a + ah sec ( The function Q(s) is dependent of the shape of the moon as well as the detector size and frequency. [5] 2 . 3 F r i e d p a r a m e t e r One of the main aims of site testing- is to determine the degree of blurring produced by the atmosphere. This is quantified by Fried's coherence length ro: The resolving power is determined by the telescope if its diameter is smaller than the Fried's parameter. Otherwise, the atmosphere profile affects the re-solving power. The F W H M of a star image, in a large telescope, is directly related to the Fried parameter o o 1 - 3 / 5 r 2 (2.50) #i/2 = 0.976 — . ' r0 (2.51) The coherence time related the coherence length ro defined by to also determines the time scale over which the changes in turbulence becomes significant. 9 Chapter 3 Othe r technologies In the last few decades, site testing techniques have been developed to provide tools for the choice of the best location for a new telescope. Several technologies have been used to determined atmospheric structure. Here an overview of most of them is presented. 3 . 1 D i f f e r e n t i a l i m a g e m o t i o n m o n i t o r First used in the 60's, the Differential Image Motion Monitor (DIMM) [6] uses the variance of image motion technique to determine the seeing. The instru-ment is built to admit starlight through two circular holes at the entrance of a flat pupil to obtain two images as shown in figure 3;1. These subapertures of diameter D and separation d are created by a mask positioned at the entrance of a telescope. As the subaperture diameter D is smaller than the mean Fried parameter (ro), the images correspond to a blurred diffraction pattern. A wedge prism'separates light in the pupil plane into two images parallel to the aperture separation. Because of the small size of the aperture, scintillation due to turbulence affects the image. Because of the separation between apertures, instant intensity levels are not strongly correlated. Wavefront slope differences over the pupils result in differential motion of the two images. The differential motion will be evident as soon as the distance between the two apertures equals a few times their diameter. Hence a compact instrument can be built and still provide good sensitivity. The two images are related to two distinct paths of a single star. Hence, the wavefront perturbations caused by turbulence in the atmosphere cause the relative motion of the two images. On the other hand, the tracking errors and wind-shakes cause a common motion and affect both images in the same way. Two seeing measurements are made to improve the accuracy. A longitudinal measurement in the direction of the subaperture axis and a transverse perpen-dicular to the subaperture axis are averaged. One of the main advantages using this differential method is the low sensitivity to tracking errors. This instrument measures ro, and is not affected by wind during data acquisition. In order to study the turbulent profile, D I M M and MASS are used together. Chapter 3. Other technologies 10 Figure 3.1: D I M M principle: light hits the parabolic primary mirror M l and the hyperbolic secondary mirror M2 afterwards. At the entrance of the tele-scope, two holes (Wi and W2) allow light in. One of them is filled by a wedge prism to change the direction of the light, (http : //aanda.u — strasbg.fr :. 2002/articles/aas/full/2000/U/dsl885/node2.html) 3 . 2 M A S S The Multi-Aperture Scintillation Sensor (MASS) [7] determines the index of re-fraction profile from the scintillation of a single star. The starlight is collected by a telescope which redirects light into four concentric zones. Each zone has a spe-cific photodetector. This gives a vertical distribution of turbulence over several heights between 500 m and 16 km by measuring the differential scintillation between apertures. The wavefront is disturbed as it passes through different atmospheric layers. Each ring detects different intensity variations given the starlight scintillation. In order to avoid the effect of star brightness in the cal-culation of C%, the variance is computed of intensity normalized by the average intensity squared. As there are four rings, the differential scintillation index is calculated for six pairs. The diameters of the apertures are respectively 2, 3.7, 7.0 and 13 cm. This choice of configuration allows one to obtain integrated Cfj values for layers at 0.5, 1, 2, 4, 8 and 16 km. As seen in figure 3.2, light entering the telescope first arrives on mirror 2. It goes through a disk containing four different apertures. A beam splitter allows the long wavelength to exit the telescope and be used for other observations. The blue-green part of the spectrum is kept and directed by the segmentator 5 to four different photodetectors. The variance is computed in each ring and covariance for each pairs of rings. This instrument is often used with D I M M . Chapter 3. Other technologies 11 Figure 3.2: MASS principle: MASS receives light on a detector divided in four concentric sections. ( http : / / www.eso.org/gen — fac/pubs / astclim / par anal / asm / slodar / TheSTololoSLODAR-C ampaign.htm) 3.3 S C I D A R 3.3.1 Principle In 1973, Vernin and Roddier [8] invented the technique called Scintillation De-tection and Ranging to address the difficulty of measuring atmospheric turbu-lence at low altitude. Those layers typically cause 85% of the image degradation. SCIDAR works with an image acquisition system. Starlight is received by a tele-scope where it forms images of a virtual plane that is located below the pupil. This instrument measures the optical turbulence profile (C_) and horizontal ve-locity [v(h)) of turbulent layers. Using a binary star with an angular separation p will give an image of two similar patterns of intensity distribution separated by a distance d = ph for a unique turbulent layer at height h (figure 3.3). Assuming only one optically turbulent layer to simplify the understanding, this layer affects the pupil intensity at the ground on the detector. As seen on figure 3.4, on the image received, the rays passing through the same region of turbulence (X) are correlated. Thus, point B and C on this figure are correlated. The auto-correlation of several frames of the pupil image looks like figure 3.5 showing the expected main peak with two side-lobes. The motion of the principal peak of light intensity received on the detector (figure 3.4) gives the wind velocity by taking the displacement dx of the peak over the time dt. The main advantages are good spatial and temporal resolution and the access to wind profile. On the other hand, SCIDAR requires binary Chapter 3. Other technologies 12 Figure 3.3: S C I D A R uses a binary star and studies its separation at the pupil plane and determines the elevation of a turbulent layer, (http : //aanda.u — strasbg.fr : 2002/articles/aas/full/199S/10/ds708Q/node2.html) stars. 3.3.2 Generalized SCIDAR Detection of the ground level layer presents a challenge because the variance of scintillation is proportional to h5^6 for a layer at altitude h. The generalized SCIDAR, designed in 1998, has the analysis plane at a distance h g s below the pupil. Then the scintillation produced at an altitude h is proportional to (h — hgs)5/6 and the separation between the scintillation patterns is (h — hgs)p. 3.3.3 SLODAR SLOpe Detection And Ranging (SLODAR) [9] is similar to S C I D A R but uses a smaller aperture and allows one to scan the first kilometer. This instrument uses a Shack-Hartman wavefront sensor on a telescope and observes a double star. The resolution depends on the number of apertures and their sizes. Using measurements from close binary stars (5-7 arcsec), this instrument achieves a vertical resolution of about 1.5 km. The observation of wide binaries (55-60 arcsec) provides a resolution of about 150 m. 3 . 4 S O D A R In late 50's acoustic scattering from the atmosphere became a subject of inves-tigation. This instrument called Sonic Detection And Ranging uses propagation of sound to analyze turbulence in the atmosphere. S O D A R uses the same prin-ciple as sonar which detects objects under water except that the medium is air and the scattering is caused by temperature variations and wind flow. A short Chapter 3. Other technologies 13 pulse of sound is transmitted in the atmosphere and is scattered by small scale turbulence. The determination of the radial velocity is done by the measurement of the Doppler shift of the sound being scattered by temperature and humid-ity fluctuations. The range of turbulence is determined by the delay between the transmission of the acoustic pulse and the reception of the returned signal. Even if different scattering patterns exist, a part of the scattered acoustic en-ergy always reflects back towards the sensor. The measurements are done in two orthogonal directions to determine the wind direction. The intensity returned from reflections is proportional to the thermal structure of the atmosphere C\. S O D A R works as classical radar but it emits sound pulses (2 kHz). This emission is backscattered by temperature inhomogeneities present in air. The echoed sound pulses are then analyzed. The main advantage of this technology is the continuous operation. However, background sound and echoes can affect measurements. It is important to choose an area free of buildings. Moreover, higher frequencies are attenuated much more than lower frequencies. With increase of temperature or lower relative humidity, the attenuation of sound is increased. This technology offers good resolution within the first kilometre but the calibration of the instrument is challenging. Chapter 3. Other technologies 14 Figure 3.4: SCIDAR principle (http : //www.photonics.ic.ac.uk/scidar/scidarJn.html) Chapter 3. Other technologies 15 Figure 3.5: Simulated pupil cross correlation intensity (http I j www.photonics.ic.ac.uk / scidar / scidar Jn.html) M m W 6 Figure 3.6: Generalized SCIDAR allows one to study the turbu-lence profile at lower height. (http : //'www.astrosmo.unam.mx/ ~ r.avila/Scidar / Scidar Pinciple.html) Chapter 3. Other technologies 16 Figure 3.8: S O D A R principle: S O D A R uses the same prin-ciple as sonar to find inhomogeneities in atmosphere ( http : / j www.acoustics.sal f ord.ac.uk/acousticsjworld/svh/sodar.htm) 17 Chapter 4 Design 4 . 1 O v e r v i e w The instrument used in this thesis is designed to detect moonlight. Twelve pho-todiodes are mounted in an aluminium beam providing a maximum separation of 1.530 m. Tubes, 25.4 cm long, serve to collimate the light, restricting the solid angle seen by each photodiode to 2-degrees of sky. The photodiodes have an active area of 5.8 x 5.8 m m 2 . The instrument is set on a tripod to track the moon. A H G M - T i t a n Losmandy tripod was used to hold the instrument and track the moon. I Figure 4.1: Picture taken in the lab of the lunar scintillometer For any two photodiodes, the cone of light they receive overlaps at a certain height. The covariance of the intensity fluctuations seen by the two diodes is therefore sensitive to turbulence above this height. As seen in figure 4.2, the photodiodes have a particular distribution on the aluminium beam. As the separation between any two photodiodes determines the height at which they are sensitive to the Cjj profile, the choice of sensor po-sitions is important. In order to maximise the number of baselines, the Golomb configuration is used. 4 . 2 G o l o m b R u l e r The Golomb ruler [10] is a mathematical term referring to an assortment of non-negative integers such there are no two distinct pairs of numbers from this Chapter 4. Design 18 h Figure 4.2: Lunar scintillometer principle. The picture above shows how the signal of each pair of photodiodes are correlated. assortment having the same difference. A well-known example is a ruler built .such that the distance between any pairs is distinct. In our project, an optimal Golomb ruler has been used to determine the positions of photodiodes. This ruler refers to the shortest Golomb ruler possible for a given number of marks. For example, in a regular ruler of twelve marks, eleven different distances can be measured. The optimal Golomb ruler allows us to measure the same amount of distances with only five marks. Chapter 4. Design 19 0 1 3 Figure 4.3: Golomb ruler for 3 marks is illustrated above. 5 1 3 1 2 1 0 1 4 6 Figure 4.4: Golomb ruler for 4 marks is illustrated above. Chapter 4. Design 20 Photodiode Position (m) 1 0.000 2 0.036 3 ' 0.108 4 0.432 5 0.522 6 0.720 . 7 0.774 8 0.990 9 1.224 10 1.350 11 1.368 12 1.530 Table 4.1: Position of photodiodes on the array In our particular case, we choose the optimal ruler of 12 marks where the positions of the photodiodes on the beam are presented in table 4.1. This Golomb configuration brings us to a different baseline for each pair of photodiodes. Prom the photodiode positions, baselines can be calculated for each pair. Set # Photodiode 1 Photodiode 2 Baseline (m) 1 10 11 0.018 2 1 2 0.036 3 6 7 0.054 4 2 3 0.072 . 5 4 5 0.090 6 1 3 0.108 7 9 10 0.126 8 9 11 0.144 9 11 12 0.162 10 10 12 0.180 11 5 6 0.198 12 7 8 0.216 13 8 9 0.234 14 5 7 0.252 15 6 8 0.270 16 4 6 0.288 17 9 12 0.306 18 3 4 0.324 19 4 7 0.342 20 8 10 0.360 21 8 11 0.378 22 2 4 0.396 ' Table continue on next page Chapter 4. Design 21 Table continued ... 23 3 5 0.414 24 1 4 0.432 25 7 9 0.450 26 5 8 0.468 27 2 5 0.486 28 6 9 ' 0.504 29 1 5 0.522 30 8 12 0.540 31 4 8 0.558 32 7 10 0.576 33 7 11 0.594 34 3 6 0.612 35 6 10 0.630 36 6 11 0.648 37 3 7 0.666 38 2 6 0.684 39 5 9 0.702 40 1 6 0.720 41 2 7 0.738 42 7 12 0.756 43 1 7 0.774 44 4 9 0.792 45 6 12 0.810 46 5 10 0.828 47 5 11 0.846 48 3 8 0.882 49 4 10 0.918 50 4 11 0.936 51 • 2 8 . 0.954 52 1 8 0.990 53 5 12 1.008 54 4 12 1.098 55 3 9 1.116 56 2 9 1.188 57 1 9 1.224 58 3 10 1.242 59 3 11 1.260 60 2 10 1.314 61 2 11 1.332 62 1 10 1.350 63 1 11 1.368 64 3 12 1.422 65 2 12 1.494 Table continue on next page Chapter 4. Design 22 Table continued ... 66 1 12 1.530 Table 4.2: Calculation of separations between photodiodes As seen in table 4.2, 66 combinations of baselines, as oppose to only 11 on a regular ruler, are used for the calculation of the covariance for the turbulence profile. This increases the accuracy of our results. 4 . 3 P h o t o d i o d e s The photodectors shown in figure 4.5 used are Si photodiodes Hamamatsu S8746-01 which have a small package with a quartz window. Each photodi-ode has an internal 1 Gf l resistor with a cut-off frequency of 32 Hz and a 5 p F capacitor in parallel in the feedback. An external 100 Mil resistor has been connected for feedback. Hence, the total feedback resistance for the circuit is 90.9 Mil. This detector is sensitive to wavelengths between 190 and 1100 nm. Figure 4.5: 12 Si photodiodes S8746-01 are positioned on the beam 4 . 4 E l e c t r i c a l c i r c u i t The instrument contains 25 microcontrollers. These are used to set programmable gains (lx 2x 4x 8x), to control the analogue-to-digital converters, and to commu-nicate with the hub. These are connected in a star configuration with a central hub that connects to a laptop computer. A GPIB bus communicates between Chapter 4. Design 23 the laptop computer and the hub. A n RS-422 serial connection running at 690 K B / s allows communication between the hub and the remote photodiodes. The acquisition is done by Labview software that allows us to set D C and A C gain and trigger time via the GPIB bus. The intensity seen by each photodiode is recorded with 16-bit resolution at a maximum rate of 1000 samples (consisting of one reading for each of 12 photodiodes) per second. Data are acquired into dual 4080-bytes buffers. Each buffer contains 170 scans containing 2 bytes of A C signal for each of the twelve photodiodes. For every 255 A C signal scans, one D C signal scan is sent. Data are 2's compliment, 16-bit stored little endian. The signal from each photodiode is split into two paths as shown in figure 4.6. The D C signal has an inverted gain of one while the A C signal has a 5-pole low-band-pass filter with a gain of 64. Its cut-off frequency is 233 Hz. This filter is followed by an adjustable gain allowing setting a gain of one, two, four or eight. Data are acquired continuously until a stop command is sent to the instrument. The rate of acquisition can be chosen between 5 or 10 ms per sample by a function called trigger time. The instrument takes 170 ms to empty a buffer of 170 samples so the instrument is ready for another set of data. Photodiode DC gain of 64 AC ADC USB-GPIB HUB MICRO PROCESSOR LAPTOP COMPUTER Figure 4.6: Block diagram of the data acquisition system for the scintillometer array 24 Chapter 5 Observations During an observation, the first step is to acquire two scans of data. From this, the best gain for the measurements can be estimated by reading the level of the AC and DC signals. The second step is to set an observation time or a number of buffers acquired and the acquisition rate. The interface of the software is illustrated in figure 5.1. Figure 5.1: The Labview user interface software used to collect data from the instrument With the instrument set on the equatorial mount, we took measurements at the 6-m Liquid Mirror Observatory in Maple Ridge. It's located at 395 m above sea level and has an average seeing of about 2 arcsec. Observations can be taken in an interval of two days before and after the actual full moon and still be accurate. Measurements were taken during the night of Saturday August 20th to Sunday 21 s t . 10 observations were taken on that night for about 15 minutes each. During Chapter 5. Observations 25 the acquisition we were unable to track perfectly the moon because of problems with the lunar tracking mode of the tripod. The tripod was unable to track the moon. Later, we discovered that the setting was not properly done for its location and time zone. For the night of the observations, the moon was tracked manually by eyesight through the scope of the instrument. As this motion is not smooth, it may have induced noise in the data. Although a part of this effect can be correlated between the photodiodes. 26 Chapter 6 Data Analysis 6 . 1 P r o f i l e Once data have been acquired, the analysis is done by a C program. Every 256th scan, the instrument sends a D C value and the software updates the D C level. Data are separated into ten sets in order to have a better evaluation of the noise. The covariance _• is calculated DC-DC. where nc is the number of samples for each pair of photodiodes i and j; AC and DC are the signal read. The covariances calculated during observations for each set are shown on figure 6.1. The profile is calculated for a model consisting of a Dirac comb with an interval Sh between heights given by uniform logarithmic spacing. Starting with a constant value for the profile times Sh, we calculated the function Kij (equation 2.39) for all covariances and heights. The aim is to minimise the error or the difference between the model and the data. The first error is the x2 2 __ J_ [Bi ~ Bmodelji)}2 /g rib cr2 where a is the rms error and rib is the number of independent baselines. The second error is the Allan variance: °\ = ^ £ ^ + !) - l n ^ ) ] 2 • (6-3) where a is the coefficient of the profile and m the number of heights in the model. The Allan variance is used to smooth the profile assuming the turbulence does not change rapidly with height. The form of the function to minimize in our algorithm is E = <r£ + \cr\-a2 is defined by (x2 — l ) 2 - A is the Lagrange multiplier that is explained in section 6.5. The model vector is \og(C%) to force the profile to be positive. Chapter 6. Data Analysis 27 Chapter 6. Data Analysis 28 6.2 Downhill simplex The amoeba function is used to perform the minimization [11]. This function uses Powell's method [12] which starts at a point P and proceeds from there in a vector direction n in N-dimensional space. Any function f(P) is minimized along the line n by one-dimensional methods such line minimization. Hence, given inputs are P and n and we need to find F to minimize f{P + nT). The amoeba function is a downhill simplex method. This method looks for the minimum of a function of multiple independent variables. The simplex is represented by a geometric figure in N dimensions containing N + l points called vertices interconnected. For example, in two dimensions, a simplex is a triangle, in three, it is a tetrahedron. We first try a guess as a starting simplex including N + l points. The algorithm searches for a minimum. However, this minimum might be in fact a local minimum. Once the first points PQ is chosen, the N other points for the simplex can be Pi = P0 + Ta (6.4) where is N unit vectors and T is a constant that gives the problem's charac-teristic length scale. Y can be different for each vector direction. The algorithm goes on by choosing N + l points and then the value of the function for this set of points. It stops if it finds a good enough value. Otherwise, it will continue by choosing a new point closer the best points. The downhill simplex takes a series of steps where most of them are reflections. Each step tries to conserve the volume of the simplex until a "valley floor" is reached when no significant lower value of the function are found. A step can be a reflection, a reflection and expansion, a contraction or multiple contractions (in multiple directions) as seen on figure 6.2. The method stops when the vector distance moved in a step that is fraction-ally smaller in magnitude then tolerance criteria defined. When the algorithm found a minimum, the routine should be restarted at the point where it found the minimum. 6.3 Simulated Annealing This technique independently invented by Kirkpatrick, Gelatt, Vecchi in 1983 [13] and Cerny in 1985 [14], refers to the analogy of annealing solids. The idea is to start with an initial configuration and slowly change the value of parameters in order to reduce the temperature of the system (or minimize the energy of the system). Starting with a high temperature, all configurations are reachable which prevent us from becoming trapped in a local minimum. The temperature decreases as the simulation is running. This allows us to progressively find the minimum of the wells found by changing the configuration of the system and checking if this reduces or increases the energy. We keep the solution if it re-duces the energy. If the energy increases, we change to a new configuration with a probability exp(-dE/T), where dE is the quantity that represents the change Chapter 6. Data Analysis 29 between the new system and the old system. When energy stops decreasing, we decrease the temperature. When decreasing the temperature is inefficient, we stop the process. In our case, we want to reduce the sum of a\ and x 2-6 . 4 M o n t e C a r l o The Monte Carlo Markov chain was developed by Metropolis et al [15] in 1953 as they wanted to simulate the behavior of simple fluids. This method is used to estimate the most probable set of parameters for a non linear model. The posterior probability of an event X is p(X\D,I), (6.5) where X is proposition asserting the truth of a hypothesis, I is proposition representing our prior information and D is the proposition representing the data. X in this case is our model parameters. We are interested in the expectation value of a function f(X) of our parameters that can be calculated by integrating the function weight: </(*)> = / f(X)p(X\D,I)dX. , (6.6) The posterior probability is defined as the probability of X given D and I true and that is derived from or entailed by the specified value of D and I. A posterior probability distribution is the conditional probability distribution of the uncertain quantity given the data. To increase the efficiency of the simulation, the proposal density should have a similar shape to the posterior. The first step of the algorithm used is to choose the value to test XQ at t = 0. A function of motion q(Y, Xt) to find a new point to evaluate is chosen. A random acceptance criterion is chosen between 0 and 1, (7(0,1). The Metropolis ratio is defined as P(Y\D,I) q{Xt\Y) P(xt\D,i)q(Y\xty [-> If r, is larger than the acceptance criterion, Xt+i — Y with a probability r. Otherwise, Xt+i = Xt. On the following step, t is increased and we continue the algorithm until we have a satisfactory simulation. Usually, there is a burn-in period before the algorithm becomes efficient. q(Xt | Y), the proposal density, is chosen to be symmetric. Therefore, the ratio q(Xt \ Y)/q(Y \ Xt) is unity. Now, p is the probability for the system to be found. It is given by the Boltzman factor exp(—Ax/T). x 2 is the energy of the system, or the function to be minimized in our case. Chapter 6. Data Analysis 30 6 . 5 L a g r a n g e m u l t i p l i e r The Lagrange parameter is often used for optimisation problems. It allows one to optimize a function of several independent variables subject to one or more constraints. A new unknown scalar variable, called the Lagrange multiplier, is introduced for each constraint. A linear combination with the multipliers as coefficients is formed. For the purpose of the project, we have to minimize the function E = x2 + ^ A - (6-8) This parameter A is a tool to control the balance between the smoothness of the model and the goodness of fit. Once the Lagrange parameter is set, different levels of noise are applied to the simulation to understand how well the initial profile is recovered. 6 . 6 S i m u l a t i o n s In order to test the analysis program, we generated some simulated profiles in order to assess the accuracy with which the profile is recovered. We compare dif-ferent analysis methods including the downhill simplex, simulated annealing and Monte Carlo. Seven profiles were generated following a flat, step or exponential profiles. A delta function was added to represent turbulence at high altitude. As we employ a Lagrangian parameter, we first optimise this parameter in order to best recover the profile. A sum of delta function with 1 m spacing within the first 300 m of atmosphere and one at 10 km model the Cjy profile. The highest value of Cfj is 10~ 1 2 m 1 / 3 . The delta function at 10 km is 2 • 10~ 1 2 m 1 / 3 . Profile # Profile Type Cut off height (m) 0 flat step 10 1 flat step 30 2 flat step 100 3 only a 5 function at 10 km — 4 flat step with a 5 function at 10 km 10 5 exp(-/i/10) — 6 exp(—h/10) with a 5 function at 10 km — 7 /i/10exp(—h/50) with a S function at 10 km — Table 6.1: Description of the eight profiles simulated The first step is to calculate the covariance signal of the model without noise. A different set of A parameters are tried in order to find the smoothest solution and the closest one to the initial model. Below, a graph with some A parameters tested is shown for each model. The signal is normalized to have a variance of one part per million. Chapter 6. Data Analysis 31 The second step is to try different levels of noise in the data to see how well the profile could be recovered. The noise is added by choosing a percentage level, then multiplying by a random function, drawn from a Gaussian distribution, and adding this to the covariance profile calculated by the simulation program. The recovered Cjj profiles are shown below. For comparison, profiles recovered by the downhill simplex algorithm are also shown. Finally we compare the data with the original model. The term in parenthesis in the title refers to the parameter written in the legend. Chapter 6. Data Analysis Figure 6.2: Steps done by the algorithm of the simplex downhill Chapter 6. Data Analysis 33 c n 2 M o d e l prof i le 0 to 7 — , , , 1 , , , , 1 , , , r - y l og (He igh t ) Figure 6.3: Models sent to the optimisation program before adding noise. Height H;_C%[m->'*} Chapter 6. Data Analysis 34 Cl profile 0 (with noise.in the simulation) (DS) -12 -14 -16 1 Chapter 6. Data Analysis 35 profile 0 (Lagrangian parameter) 2 Chapter 6. Data Analysis 36 cj; profile 0 (with noise in the simulation) (SA) -12 -14 -16 -18 -20 3 Chapter 6. Data Analysis 37 cn2 profile (noise) 5000MC+SA -14 ~ -16 -18 -20 I I I I I I I I I I I I I I 1 j 1 1 1— ' I ° ffi • -e - # s * * * * • • -• °0.05 -_ 0 0 0.1 * 0.2 o * • X * 0.3 O * 0.5 * D ° A # • O A * ' • o * p -I I 1 1 I 1 , , , I , - , 4> 0 1 2 3 4 log(Height) 4 The profile 0 is a flat step of 10 m wide with a C2N of l O - 1 2 5 . In the DS plot, the program tends to smooth the step of the profile in one slope. The program assumes the profile does not fall drastically. In the low noise levels case (2-5%), the step is clear and decreases at 10 m. The SA profile is closer to the initial step for a A parameter of 1.5. Once the noise is added, the profile tends to be a slope. As the noise decreases, the slope is closer to a step. The last simulation is a Monte Carlo with 5000 iterations using the same A parameter. The results are similar to the SA plot but for high noise the profile increases around 10 m and drops. Chapter 6. Data Analysis 38 C§ p ro f i l e 1 (wi th no ise i n the s i m u l a t i o n ) (DS) -12 -16 . I I I I I I I I 8 A — ' 6 ' 1 1 . I " 1 8 -S I . O - " — O A i = i Amount of noise i - 0.01 • s 0.02 i _ 0 • 0.05 • i 2 O 0.1 -Q # & * • • O t Q * 0.2 0.25 0.3 I , , , , I , , 1 1 1 , 1 0.5 • , , , , I 0 1 2 3 4 log h (m) 5 Chapter 6. Data Analysis 39 Cj; profile 1 (Lagrangian parameter) -12 -14 -16 h 00 •2 -18 -20 -22 h 6 Chapter 6. Data Analysis 40 -15 -20 p ro f i l e 0 (with noise in the s i m u l a t i o n ) (SA) I i i i . I i i i " 1 • - - _ « i 4 i . i | i i , , | I s i t 8 * -Amount of noise - 0.01 ' 0.02 D 0.05 0.1 -° o • a • * o • i " 0 . 2 I * 0.25 0.3 0.5 -O , , , I , , . , T -0 1 - 2 3 4 log h (m) 7 The profile 1 is another flat step which show a cut-off at 30 m. The DS plot shows a clear step for noise levels between 2 and 5 %. For higher noise, the program smooth results to get a slope for the step. This slope is similar for all the noise levels used. The lagrangian parameter chosen is 1.5 in the SA simulations. A step is defined for a noise level of 1-2% otherwise, the profiles tend to be slopes. The slope is increasing with the noise level. Chapter 6. Data Analysis 41 Cl profile 2 (with noise in the simulation) (DS) -12 -13 QO O -14 -15 1 1 1 1 1 1 l i i i i 1 * * 5 a 0 -t * o w 0 • S 2 * I 1 * • i -o • * Amount of noise * * o - 0.01 o • 0.02 o 0 * ° 0.05 6 • • o 0.1 0.2 --• * 0.25 0.3 - 1 i I I 1 1 i t 1 o -0 1 . 2 3 4 logh(m) ° - 5 8 Chapter 6. Data Analysis 42 Cjj profile 2 (Lagrangian parameter) -12 -14 -16 Lagrangian Parameter - 0.005 0.05 0.5 ° 1 2 3 5 10 30 40 50 500 -18 -20 J I I I I I I , I I I I I I I I I I I I L 0 1 2 3 4 log h (m) 9 Chapter 6. Data Analysis 43 eg profile 2 (with noise in the simulation) (SA) -12 -14 i Amount of noise 0.01 0.02 ° 0.05 0.1 16 0.2 0.25 0.3 0.5 -18 -20 J I I 1_ I I I I I I I I I ! I I I I I I L 0 1 2 3 4 log h (m) 10 Chapter 6. Data Analysis 44 cn2 profile 2 5000MC+SA (noise) -12 -13 -15 -16 I „ I 1 1 1 1 1 1 1 1 1 1 ' 1 i i i i ' I o t — e - 0 g ini • o fl 0.1 _ • -fl • 0.3 0 - • o 0.5 fl o • - o _ o - fl • o > ° ° D I . ! -0 1 2 3 4 log(Heighl) 11 The fiat step that cuts off at 100 m in profile 2 simulated by DS gives a clear cut off at 100 m. The lagrangian parameter chosen is 5. In the case without noise, the curve on the simulation cut off around 100 m. When noise is added, the curve decreases around 10 m with a smooth slope. The 5 000 iterations Monte Carlo simulation gives similar result to the SA. For a profile without noise, it gives a step with a cut off at 100 m. With noise, the simulation gives a slope. Chapter 6. Data Analysis 45 -14.6 Cl profile 3 (with noise in the simulation) (DS) - i 1 1 1 i 1 1 1 r i i — i | i i i i r Amount of noise 0.01 -14.8 0.02 E -15 0.05 0.1 -15.2 -15.4 _L 1 2 log h (m) 12 Chapter 6. Data Analysis 46 Cg profile 3 (Lagrangian parameter) -16 £ -18 -20 i 1 1 1 1 1 1 1 1 1 _ • - Lagrangian Parameter _ - 0.005 0.05 • 0.5 o 1 - 2 o * 3 D O s -* • o s * • o a • * • 0 s 5 fi * • o 2 - K * • o fl _ 10 a ft 0 O s # X 0 • * o • fl 2 50 x e * • X 8 * • o 0 500 X 9 • o X * • - K *• -a 9 I 1 1 1 X 1 , , 1 0 1 2 3 4 log h (m) 13 Chapter 6. Data Analysis 47 profile 3 (with noise in the simulation) (SA) -15 r -16 -17 r-00 ° -18 -19 -20 14 The profile 3 is only a delta function at 10 km. Due to the way that the program calculate the integral, it results that the has a factor of 10~ 3 m ~ 2 / 3 for its delta value different from the original one. The box of integration for this value has 1 000 m length. The program converts data without taking this in account. Chapter 6. Data Analysis 48 Cj} profile 4 (with noise in the simulation) (DS) -12 -13 cf -14 -15 -16 15 Chapter 6. Data Analysis 49 Cl profile 4 (Lagrangian parameter) -11 - ! " '1 1 I i l i i i i I i i i i I -Lagrangian Parameter ft - 0.005 ft - » " 0.05 a • - * # ! i 0.5 • O o • * I 1 - • _ * 2 - 3 — a * 3.5 - # o 5 - a 10 - - » 50 # ft o 500 0 - o A ft * — a ° 5 ft - A i 5 - i i 1 , , • i * , i i i i i i -0 1 2 3 4 log h (m) 16 Chapter 6. Data Analysis 50 eg profile 4 (with noise in the simulation) (SA) n — 1 1 1 1 r Amount of noise - 0.01 4 0.02 D 0.05 0.1 * 0.2 * 0.25 0.3 0.5 -12 -13 -14 -15 2 log h (m) 17 The profile 4 is a flat step with a cut off at 10 m and a delta function at 10 km. In the DS plot, the cut off is clear for low noise level (0.2-0.5 %). For higher noise, the graph is a slope. The delta function value is around 10~ 1 7 m ~ 2 / 3 as oppose to 10~ 1 2 m ~ 2 / 3 . The lagrangian parameter chosen is 3.5 for the SA simulations. The profile decreases quickly at 10 m for cases with low noise level. For the part before 10 m, the C% value is higher then the original value 10- rn -2/3 Chapter 6. Data Analysis 51 Cg profile 5 (with noise in the simulation) (DS) -12 -13 o -14 -15 -16 18 Chapter 6. Data Analysis 52 -15 -20 Cg profile 5 (Lagrangian parameter) " I ' I I I ' • . -i i i i | * . : » : * • Lagrangian Parameter - " 8 - 0.005 0.05 8 s 0.5 8 „ 1 A X it O 1.5 « * ° • I : • 2 ' 3 - 5 i 3.5 5 10 50 500 i i i I . i i - I 1 " t I I i i 1 0 1 2 3 4 log h (m) 19 Chapter 6. Data Analysis 53 eg profile 5 (with noise in the simulation) (SA) -12 I " I , , , , . | , , i i " " « ! ' * , | , , , , | Amount of noise - 0.01 " 0.02 -14 -8 ° 0.05 0.1 -• t o • o • 0.2 * 0.25 ' * 0.3 -16 -0.5 • -18 - --20 - --22 I • i i 4 0 1 2 3 4 log h (m) 20 Chapter 6. Data Analysis 54 cn2 profile 5 - 10% noise (5000 MC simulations) c o log(Height) 21 The profile 5 is an inverse exponential in function of the height. For this profile, the value should be around 1 0 - 1 2 m ~ 2 / 3 at 1 m, 1 0 - 1 2 , 4 m - 2 / 3 at 10 m and 10~ 1 6 , 3 m ~ 2 / 3 at 100 m. For the DS simulation, for noise level bigger than 1 %, the value is higher than the expected value for the lower height. The curve that we should obtain is smoothed to a slope. The lagrangian parameter chosen for the SA simulation is 1.5. At low noise (0.01-0.1), the profile follows the expected curve. For higher noise, only the tail is over-evaluated by a factor 10~ 2. The Monte Carlo simulation for 5 000 iterations obtained has a C ^ value a little too high by a factor 10° 5 for low height and a factor 10 for high elevation. Around 10 m, we obtain the expected value. Chapter 6. Data Analysis 55 Cg profile 6 (with noise in the simulation) (DS) -12 -16 I I I , I I I I I * s I ' ' ' , | , I I I | 6 8 -— a - - • • • s > Amount of noise j 5 -- 0.01 * " 0.02 * • i ° 0.05 — o i t ° 0.1 o X « -0 o 0.2 D a a " I I ! I I I I I I I I 1 I 1 * 0.25 0.3 o 0.5 t , I , i i i I 0 1 2 3 4 log h (m) 22 Chapter 6. Data Analysis 56 Cg profile 6 (Lagrangian parameter) -14 -15 - I 1 1 1 " 1 i | i i 1 1 1 -1 - * i -1 -ft 1 I Lagrangian Parameter - 0.005 l * 0.05 -• O 0.5 1 2-. I o * . 3 -• o • 0 3.5 5 i a » 10 ft o 50 - i _ i o * • * a i 500 1 ! -0 1 2 log h (m) 3 4 23 Chapter 6. Data Analysis 57 eg profile 6 (with noise in the simulation) (SA) -12 -13 -14 -15 24 The profile 6 is similar to the profile 5 except that a delta function is added at 10 km. The results for the DS simulations are similar to those of profile 5. The delta function value has a factor 10~3 for the same reason described for the profile 3. As the noise increases, the evaluation of the delta function is lower than expected. The lagrangian chosen is 3.5. For any amount of noise, the profiles have similar shape. The profile is higher at low altitude as well as heights over 100 m. Chapter 6. Data Analysis 58 Cg profile 7 (with noise in the simulation) (DS) 1 k i 1 1 1 1 1 1 , 1 1 1 1 1 1 • * _ * o fi J * $ • * fl a * • * * 6 i o x fl * o - • • B * ft x * * fl o - o x Amount of noise fl • X • * 0.01 - 0.02 - * 0.05 • * o 0.1 * 0.2 fl * * 0.25 • - 0.3 0.5 - * * ' i l l , 1 1 1 • I , ' , , 1 I I I I I 0 1 2 3 4 log h (m) 25 Chapter 6. Data Analysis 59 Cg profile 7 (Lagrangian parameter) -12 - I 1 1 1 1 1 t , , , i | _ -A I - 8 # ft * 8 • -* -* o • a Lagrangian Parameter * •b * - 0.005 • 0.05 * • 0.5 * o 1 * * O • - 2 o • * 3 0 o 0 5 I • • 1 , 1 I 10 50 500 • I 0 1 2 3 4 log h (m) 26 Chapter 6. Data Analysis 60 eg profile 7 (with noise in the simulation) (SA) 1 , i . i | . , • i | , - * * i i i | , , , , | - * * « * * ° 2 * * . * * * * ~ 0 ° '8 * * * * ; 5 » 6 . b s 6 ; - - » ; s S - ' . » Amount of noise - - 0.01 4 0.02 D 0.05 . ° 0.1 0.2 * 0.25 0.3 -0.5 , , , I , , , , I 0 1 2 3 4 log h (m) 27 Figure 6.4: Simulations done with downill simplex (DS), simulated annealing (SA) and Monte Carlo methods The profile 7 is more complex with a variation proportional to the height and inverse exponential to the height. The DS plot is unable to get a profile close to the original one. The lagrangian parameter chosen is 5 so the value at low altitude is closer to the initial profile. The noise affects the low altitude C ^ value by a factor 10 for the most noisy signal for the SA simulation. Chapter 6. Data Analysis 61 The graphs in figure 6.4 show in the first case the simplex downhill method, followed by the two graphs of simulated annealing method. The first of the two showed different results without noise for different Lagrangian multiplier while the second one is for the best Lagrangian parameter at different noise levels. From these results, we notice that reducing the temperature of the system with the simulated annealing algorithm helps to find a better minimum and recover with a better accuracy the initial profiles. For the downhill simplex method, the addition of noise to the simulation reduces the recovery process. For simulated annealing, the simulations with noise tend to recover the profile with a better accuracy. The right choice of the Lagrange multiplier helps to direct the simulation to the correct range of values for Cjy. In the case of the simulated annealing method, the step function in the model tends to be simulated as a slope which increases with the amount of noise in the model. In the Monte Carlo method, results tend to be similar to the simulated annealing method. However, it was not always possible to obtain a profile when there was too much noise in the simulations. Nevertheless, results shown the possibility to implement Monte Carlo in order to recover profiles. In the case of the profile 3 representing only one delta function at 10 km, we observe a linear tendency for the points below 200 m. This is due to the program to attempt to model Cff(h)dh. 6.7 Turbulence profiles from data From the covariances measured at the Liquid Mirror Observatory, we have esti-mated the C% profile. Figure 6.5 and 6.6 show the profiles calculated from these data using different techniques and compared with Tokivinin et al. models. In the litterature, severals measurements has been done for different sites. Tokovinin and Travouillon presents results at Cerro Pachon observatory in Chile [16]. They found that the C ^ profile for the ground-layer can be rep-resented by a decaying exponent over 20 to 40 m. For a bad seeing, this decay-ing region can be extended to 100 m. They suggested a model for the average optical turbulence profile as a combination of two exponent C2N(h) = Aexp(-h/h0) +Bexp{-h/h1) (6.9) with A and B are in units of 1 0 - 1 6 m - 2 / 3 , ho and hi in m. A typical profile has the main characteristics of having strong turbulence near the ground and producing a seeing of 0.5 arcsec. At altitudes above 200 m, the C ^ value is about 10~ 1 6 and continues to decay slowly. In the results from our instrument, the turbulence is higher at lower altitude but similar around 100 m. In the second case (figures 6.7 and 6.8, Tokivinin et al. equation is used to compare with our results setting A = 3000, B = 30, ho = 30 and hi = 500. This second attempt gives a model that follows more our results. Although, the slope is bigger at lower altitude and decreases around 100 m. Here, only the smoothest curves are compared with the Tokivinin models. Chapter 6. Data Analysis 62 Figure 6.5: August 21s* C% profiles recovered with downhill simplex compared with Tokivinin et al. bad seeing profile Both plots are similar. On the other hand, no analysis could be done with Monte Carlo. This might be due to the non smoothness of covariances calcu-lated. As seen on the profiles recovered from data, the C^ value oscillates be-tween 10~ 5 and 1 0 - 2 0 m - 2 / 3 for height between the ground and approximately 200 m. We would expect a smoother curve for this range of measurements. For the smoother tests, a A C gain of 2 and a D C gain of 4 has been applied during the acquisitions (tests 5,6 and 7). The noisier ones have a A C or D C gain lower. During observations, we had to track manually. Also, the lower photodiodes may have received some interference from trees in the area. It is possible those intereferences cause the noise in the results. From the Cjj profile, the Fried parameter is also calculated for these measurements. The Fried parameter has been calculated for the smoothest result. For the site of the observatory tested, a seeing of 2 arcsec is expected. As seen in the figure 6.9, the Fried parameter is lower and smoother for the smoothest profiles (test 5, 7 and 9). A good site would have a Fried coherence length of 10 cm, while a bad site is characterized with a 5 cm value. The lowest Chapter 6. Data Analysis 63 CjJ donnees (simulated anneling) 08/21/05 log h (m) Figure 6.6: August 21s* Cfj profiles recovered with downhill simplex and simu-lated annealing compared with Tokivinin et al. bad seeing profile value of Fried parameter are rejected by their profiles which is oscillating a lot according to our simulations. Those profiles also presents a high value for the Allan variance and x 2-Chapter 6. Data Analysis 64 Chapter 6. Data Analysis 65 Chapter 6. Data Analysis 66 Fried parameter SA August 2005 T 0.2 -- test 2 " test 3 D test 4 - ° test 5 ? - * test 6 o 0.15 - * test 7 parameter -* test 8 * test 9 parameter -X X X X X A. 'test 10 •a 0.1 o V-• • • • • t i , ^ . . . i i , ' ^ ' 4 , i I I I I 0 1 2 log(Height) Figure 6.9: August 21 s' ro Fried parameter 67 Chapter 7 Conclusion We have built and tested a lunar scintillometer array. With it, we obtained a set of data at an astronomical site. This instrument has shown its capability to obtain a turbulence profile of the atmosphere below 200 m. Even if some improvement must be done before the next observations, the instrument is able to measure scintillation from the moon. We obviously need more data on a longer period to have a good analysis of a site and better accuracy for the profiles calculated. Some work could be done on the recovery technique by testing Monte Carlo parallel tempering for example. Instead of decreasing the temperature, one could have different sets of Monte Carlo at different temperature and switch between systems after a certain number of iterations. In the future, the instrument could be set with another instrument as a complement and to provide independent verification. It would then also be possible to estimate the impact of the noise associated with the photodiodes and the electronics and to have a better understanding of our measurements. This instrument presents a good solution for the measurements of turbulence profiles within the lower kilometre of the atmosphere. 68 Bibliography [1] D. Blanco. Effects of air pressure, temperature, and relative humidity on seeing. Proceedings of SPIE, 4004:552-558, 2000. [2] R. Barletti, G . Ceppatelli, L . Paterno, A. Righini, and N. Speroni. Mean Vertical Profile of Atmospheric Turbulence Relevant for Astronomical See-ing. Optical Society of America Journal A, 66:1380-1383, 1976. [3] F . Roddier. The effects of atmospheric turbulence in optical astronomy. Progress in Optics, 19:281-376, 1981. [4] V.I. Tatarski. Wave propagation in a Turbulent Medium. McGraw-Hill , 1961. [5] P. Hickson and K . Lanzetta. Measuring atmospheric turbulence with a lunar scintillometer array. The Publications of the Astronomical Society of the Pacific, 116(826):1143-1152, 2004. [6] M . Sarazin and F . Roddier. The E S O differential image motion monitor. AAP, 227:294-300, January 1990. [7] V . Kornilov, A . A. Tokovinin, O. Vozyakova, A. Zaitsev, N . Shatsky, S. F . Potanin, and M . S. Sarazin. MASS: a monitor of the vertical turbulence distribution. Adaptive Optical System Technologies II. Proceedings of the SPIE, 4839:837-845, February 2003. [8] R. Avila, J . Vernin, and L . J . Sanchez. Atmospheric turbulence and wind profiles monitoring with generalized scidar. AAP, 369:364-372, April 2001. [9] R. W . Wilson. S L O D A R : measuring optical turbulence altitude with a Shack-Hartmann wavefront sensor. MNRAS, 337:103-108, November 2002. [10] M . Gardner. Mathematical games, sa, 218:121-124, 1968. [11] William H . Press, Brian P. Flannery, Saul A. Teulovsky, and William T . Vetterling. Numerical Recipes in C: The Art of Scientific Computing. Cambridge Press, 1992. [12] M . J . D. Powell. A n efficient method for finding the minimum of a function of several variables without calculating derivatives. Computer Journal, 7:155-162, 1964. Chapter 7. Conclusion 69 [13] S. Kirkpatrick, C . D. Gelatt, and M . P. Vecchi. Optimization by simulated annealing. Science, Number 4598, 13 May 1983, 220, 4598:671-680, 1983. [14] V . Cerny. A thermodynamical approach to the traveling salesman problem: A n efficient simulation algorithm. J. Optim. Theory Appl., 45:41-5 I, 1985. [15] N. Metropolis, A. W . Rosenbluth, M . N. Rosenbluth, A. H . Teller, and E . Teller. Equations of state calculations by fast computing machine. J. Chem. Phys., 21:1087-1091, 1953. [16] A. Tokivinin and T . Travouillon. Model of optical turbulence profile at crreo pachon. Mon. Not. R. Astron. Soc, 000:1-8, 2005.
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Lunar scintillometer array for atmospheric turbulence...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Lunar scintillometer array for atmospheric turbulence measurement Dumais, Marie-Andree 2006
pdf
Page Metadata
Item Metadata
Title | Lunar scintillometer array for atmospheric turbulence measurement |
Creator |
Dumais, Marie-Andree |
Date Issued | 2006 |
Description | The performance of large telescopes is directly related to the quality of the site. Of particular importance is the turbulence profile within the lower few hundred meters of the atmosphere. We have developed a new technique for measuring this profile, and the Fried parameter r₀, by measuring spatial correlations in the scintillation of moonlight. Our lunar scintillometer array uses 12 photodiodes to measure the intensity fluctuation covariance function, and from this we obtain the index of refraction profile C²[sub N] to a height of about 200 m with a vertical resolution of order 2%. This paper describes the instrument, the technique, and results of initial tests. |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2010-01-05 |
Provider | Vancouver : University of British Columbia Library |
Rights | For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use. |
IsShownAt | 10.14288/1.0085734 |
URI | http://hdl.handle.net/2429/17523 |
Degree |
Master of Science - MSc |
Program |
Astronomy |
Affiliation |
Science, Faculty of Physics and Astronomy, Department of |
Degree Grantor | University of British Columbia |
GraduationDate | 2006-05 |
Campus |
UBCV |
Scholarly Level | Graduate |
AggregatedSourceRepository | DSpace |
Download
- Media
- 831-ubc_2006-0032.pdf [ 3.49MB ]
- Metadata
- JSON: 831-1.0085734.json
- JSON-LD: 831-1.0085734-ld.json
- RDF/XML (Pretty): 831-1.0085734-rdf.xml
- RDF/JSON: 831-1.0085734-rdf.json
- Turtle: 831-1.0085734-turtle.txt
- N-Triples: 831-1.0085734-rdf-ntriples.txt
- Original Record: 831-1.0085734-source.json
- Full Text
- 831-1.0085734-fulltext.txt
- Citation
- 831-1.0085734.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0085734/manifest