A New Method for Aerosol Measurement using Wide-field Photometry
Jan Ebr, Sergey Karpov, Ji?í Eliášek, Ji?í Blažek, Ronan Cunniffe, Ivana Ebrová, Petr Jane?ek, Martin Jelínek, Jakub Juryšek, Dušan Mandat, Martin Mašek, Miroslav Pech, Michael Prouza, Petr Travní?ek
DDraft version January 11, 2021
Typeset using L A TEX default style in AASTeX63
A New Method for Aerosol Measurement using Wide-field Photometry
Jan Ebr, Sergey Karpov, Jiˇr´ı Eli´aˇsek,
1, 2
Jiˇr´ı Blaˇzek, Ronan Cunniffe, Ivana Ebrov´a, Petr Janeˇcek, Martin Jel´ınek, Jakub Juryˇsek, Duˇsan Mand´at, Martin Maˇsek, Miroslav Pech, Michael Prouza, andPetr Tr´avn´ıˇcek FZU – Institute of Physics of the Czech Academy of Sciences, Na Slovance 1999/2, Prague 182 21, Czech Republic Institute of Theoretical Physics, Faculty of Mathematics and Physics, Charles University,V Holeˇsoviˇck´ach 2, 180 00 Prague, Czech Republic Nicolaus Copernicus Astronomical Center, Polish Academy of Sciences, Bartycka 18, 00-716 Warsaw, Poland Astronomical Institute of the Czech Academy of Sciences, Friˇcova 298, Ondˇrejov 251 65, Czech Republic (Received September 2, 2020; Revised January 8, 2021; Accepted)
Submitted to AJABSTRACTWe present a new method to measure the vertical aerosol optical depth (VAOD) during clear nightsusing a wide-field imager – a CCD camera with a photographic lens on an equatorial mount. A series of30-second exposures taken at different altitudes above the horizon can be used to measure the VAODwith a precision better than 0.008 optical depths within a few minutes. Such a measurement does notproduce any light and is thus suitable for use at sites where other astronomical instruments are located.The precision of the VAOD measurement depends on laboratory calibration of spectral properties ofthe system and of the response of the camera electronics to varying illumination levels, as well ascareful considerations of details of stellar photometry and modelling of the dependence of measuredstellar fluxes star color and position within the field of view. The results obtained with robotic setupsat the future sites of the Cherenkov Telescope Array show good internal consistency and agreementwith the simultaneous measurements from a Sun/Moon Photometer located at the same site.
Keywords: photometry — atmosphere INTRODUCTIONThe observation of natural light sources outside of the atmosphere allow for the measurement of the integral opticaldepth in the direction of the source without production of any artificial light (as in the case of laser-based methods).Subtracting the contribution of molecular scattering and absorption, which can be derived from meteorological data,the aerosol optical depth can be obtained. The Sun, as a bright and stable source, has been commonly used for thispurpose and photometers have become a standard tool in aerosol measurements during daytime (Karsten 1997); morerecently, the Moon has been used in the same way to extend the measurements into at least some nights, with somelimitations on precision due to limited photometer sensitivity and difficulties in precisely describing the variations inlunar brightness (Barreto et al. 2019). Using stars as reference sources extends the measurement capabilities to nightswhen the Moon is not visible as well as to arbitrary positions on the sky, but brings its own set of challenges. The stars,being much fainter than the Sun or the Moon, require a more sensitive optical detector, which must be kept properlycalibrated (P´erez-Ram´ırez et al. 2008; Baibakov et al. 2015). Moreover, if the ability to measure the optical depthin an arbitrary direction is required, the problem of insufficient accuracy of available stellar photometric cataloguesbecomes apparent.
Corresponding author: Jan [email protected] a r X i v : . [ a s t r o - ph . I M ] J a n Ebr et al.
The need for an absolute calibrated instrument can be avoided using the well-known Langley approach (Shaw1983) which relies on the known dependence of optical depth on the zenith angle for the case of horizontally uniform(stratified) atmosphere and allows the simultaneous determination of the calibration constant of the instrument andvertical optical depth (VOD) τ by observing a star at different zenith angles, that is, at different values of airmass A (which is defined here as relative to the airmass in the vertical direction). The main limitation of this method whenusing a single celestial target, such as the Sun, Moon, single reference star or field of stars, is that VOD variationsoccur much more quickly (minutes) than the Earth’s rotation causes the chosen reference source to traverse a largespan in zenith angles (hours). Thus, we propose a method where we measure a large number of stars at a large spanof airmasses in a short amount of time (a few minutes) using a series of images (which we call a “altitude scan”)taken with a wide-field imager. The VOD, and, after the subtraction of the molecular components, the vertical aerosoloptical depth (VAOD), is then extracted from a fit as a function of the airmass giving the difference between themeasured and expected brightness of the stars. The expected brightness is based on a homogeneous all-sky catalogueand a calibrated model of the optics and the detector (which can also be derived from the same data). Note that thisvalue is sometimes referred to in literature simply as the aerosol optical depth (AOD) and the dependence on airmassis implied.The method to determine the VAOD from the airmass-dependence of extinction is in principle straightforward as thelatter is expressed by a known analytical formula, but a significant instrument calibration effort is necessary to reach aprecision better than 0.01 optical depths (hereafter we will refer to optical depth values as dimensionless). In particularit is necessary to carefully consider any systematic measurement error correlated with the brightness of the stars, suchas any non-linearity of the detector or of the sky image processing methods, because those can systematically bias thefit.This paper is organized as follows. In Sec. 2 we describe the FRAM (F/Photometric Robotic Atmospheric Monitor)telescopes and laboratory calibration measurements necessary to properly process the acquired data. In Sec. 3 wedescribe the procedure to obtain a star extinction model as a function of the airmass including the molecular andaerosol components and the previous laboratory measurements. In Sec. 4 we discuss the details of the photometry(i.e. the conversion of the raw detected star signal to star brightness) and its uncertainties. In Sec. 5 we quantifyvarious systematic uncertainties. Finally we briefly discuss the comparison of the results with those from a Sun/MoonPhotometer and possible future prospects. FRAM TELESCOPES2.1.
Purpose, design, and operation
The FRAM design stems from the need for atmosphere characterization in astroparticle-physics experiments, whereultra-high-energy cosmic rays or very-high-energy gamma photons are detected indirectly by collecting the faint (fluo-rescence or Cherenkov) light created during the passage through the atmosphere of the extensive air showers (cascades)started by those particles. The light from these showers can travel from a few kilometers up to tens of kilometersthrough the atmosphere before reaching a detector. The amount of light detected can thus be significantly impacted bythe transparency of the intervening atmosphere, hence the need for accurate characterization of atmospheric extinction(Gaug 2017; Keilhauer 2019).There are many options in choosing the physical setup for a wide-field aerosol measurement, depending on furtherapplications of the instrument and the availability of complementary atmospheric data. In this paper we present resultsbased on the data from the three prototypes of the FRAM telescopes proposed as a part of the atmospheric monitoringsystem (Janeˇcek et al. 2017) for the future Cherenkov Telescope Array (CTA) gamma-ray observatory (Acharya et al.2013). Two FRAM prototypes (Prouza et al. 2019) have been installed at the future CTA South (CTA-S) site in theAtacama desert in Chile near the European Southern Observatory site Cerro Paranal and one prototype at the futureCTA North (CTA-N) site at Observatorio Roque de los Muchachos in La Palma, Canary Islands, Spain. The dataused in this paper has been taken by the three prototypes during 2017–2019 as a part of the preliminary CTA sitecharacterization campaign (Ebr et al. 2019a).The main requirement for the CTA FRAM (Janeˇcek et al. 2017) is to detect changes in the transparency across thewhole field of view of the Imaging Atmospheric Cherenkov Telescopes (IACTs) installed at CTA site in a single image,thus necessitating a 15 ×
15 degree field of view to include a significant margin (Ebr et al. 2017a). The CTA FRAMdesign is based on modified version the of the previous FRAM telescope (Prouza et al. 2010; Pierre Auger Collaboration2020) installed at the Pierre Auger Observatory (Pierre Auger Collaboration 2015) for several applications (Ebr et al. erosol Measurement using Wide-field Photometry × ×
36 mm 16-megapixel chip to obtain the required 15 ×
15 degree fieldof view. The cameras are equipped with photometric Johnson-Cousins BVRI filters, but only the B-filter data areused in this work. The B filter has its maximum transmission around 420 nm and is located at the short-wavelengthend of the spectral transmissivity of the lens used (see Sec. 2.2 and Fig. 1), and avoids almost all molecular absorptionbands. Therefore the molecular contribution to the VOD using this filter is almost completely described by Rayleighscattering and thus depends only on the integral air column above the telescope, which is easily available from globalweather models and has little variation with time. It would be possible to use a different selection of narrowband filtersso to avoid the atmospheric absorption bands altogether – or, for any other reason simplifying the analysis. Howeverthe use of standard Johhson-Cousins filters facilitates the comparison with various astronomical catalogues and alsoallows further use of the data for astronomical purposes – even though the catalogue we use, Tycho2, incidentally usesa slightly different bandpass and thus some conversion is still necessary.The CTA FRAM prototypes, deployed for the purposes of atmosphere characterisation on the site, take altitudescans from horizon to zenith in 7 images with 30 seconds of exposure per image most of the nights; some of the CTA-Ndata has been taken with a higher density of images close to the horizon and a gap in small zenith angles for sometime due to a software problem, but this does not seem to induce any systematic effects in the results.For all images the exposure start time is synchronized to UTC time scale using the NTP protocol to an accuracybetter than 1 second and all images in a scan are taken with precisely the same exposure length. Due to a softwareissue, for some time we had an occasional problem of the telescope mount still slewing during the start of an exposure.This may cause visible trails from bright stars in the image, but those are sometimes hard to detect visually, andvirtually impossible to detect algorithmically due to the variations in sky background and noise. The main effect ofthis problem was that the actual star images then corresponded to less than 30 seconds of exposure (as some of thelight was lost in the faint trail during the movement); if this happened to the image taken at the lowest altitude, thefitted VAOD would be biased by about 0.1, causing visible outliers in the time series of measured VAOD. Affectedimages can be identified from FITS metadata stored in the image files, specifically from the difference of requestedand reported telescope coordinates at the start of the exposure and all affected scans have been filtered out fromthe data. This experience shows the importance of registering as much technical information as possible during themeasurements. 2.2.
Spectral response
To model the atmospheric extinction as a function of airmass for a star observed with a telescope, it is necessary toknow the spectral response of the system R ( λ ), which is the normalized probability density of detecting a photon ofa given wavelength (as the camera signal is proportional to the number of photons, not, for example, energy density).The properties of each element of FRAM’s optics (lens, optical BVRI filters, a protective covering UV filter and theCCD camera) were measured in a laboratory. The spectral transmittance of the BVRI and protective filters wasmeasured using a Perkin Elmer Lambda 850 spectrophotometer for different incidence angles. The results show thatincidence angle has negligible impact on spectral transmittance. As the photon incidence angles are relatively smallin our telescope FoV we could neglect the angular dependence of the Fresnel losses. The spectral transmittance of thetelephoto lens varies with the distance from the optical axis due to different thickness of the lens elements, but weverified that its spectral dependence is negligible. The measurement setup uses a wide wavelength range light source Ebr et al. (LOT ARC light source with 350 W Xe Hg bulb) to cover the whole transmittance window of the lens. The light sourceis connected via optical fiber to a collimator illuminating a small region of the lens aperture. After passing through thelens, the light was collected at the focal point by the integrating sphere Labsphere 3P-GPS-053-SL and delivered viaoptical fiber to the Avantes AvaSpec Dual-channel Fiber Optic Spectrometer. All work was done in a dark laboratoryto exclude stray light. A reference spectrum from the lamp was first obtained without the lens, then the lens wasmeasured both centered (aligned with the optical axis) and at several offsets perpendicular to the optical axis to scanthe transmittance for different radial distances. The lens was tested with the iris fully open at f/2 . -0.2 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 200 300 400 500 600 700 800 900 1000 1100 t r an s m i ss i on / Q E wavelength [nm]lenslens mod.B filterB filter mod.CCDCCD specscover filtercover filter mod. Figure 1.
Measured spectral transmittance for the optical elements and quantum efficiency of the CCD sensor used in theCTA-N FRAM setup (points and dashed line). The modified and manufacturer supplied forms used for the calculation of R ( λ )are shown as solid lines. Uncertainties of the individual measurements are not shown. The relative quantum efficiency of the camera CCD chip was measured using the same broad band light sourceconnected using an optical fiber to a monochromator that allows to select the wavelength with specified full widthat half maximum (FWHM). The output light of the monochromator is connected to a variable attenuator and thendivided equally into two (the splitting ratio is wavelength independent). One of the fiber ends is coupled with theThorlab PM100USB calibrated photodiode (the accuracy specified by the producer is better than 1%). The second erosol Measurement using Wide-field Photometry T . The precision of the intensity of thelight on the output of the fiber is given by the accuracy of the calibrated photodiode, the precision of attenuator andfluctuation of the light source (monitored via the PM100USB control unit). The total intensity of the light illuminatingthe monitoring photodiode is (1000 ±
50) pW.Acquisitions were performed on the whole spectral interval from 300–900 nm. However, at both ends of the wave-length interval the reliability of the measurements decreases, we thus use a smooth approximation where high signalfluctuations are observed, see Fig. 1. The relative spectral quantum efficiency of the cameras was measured only ina limited range and with a limited sampling in wavelength; we use the manufacturer specifications instead, with thelaboratory QE measurement serving as a validation (at least in the relevant bandpass of the B filter); in any case, theQE of the camera varies slowly enough within the B filter transmission region so that its effect on the spectral responseis overshadowed by that of the other optical components. The overall R ( λ ) is obtained simply by the multiplicationof the contributions of individual components. The effects of the uncertainties in the determination of R ( λ ) on theaerosol measurements are explored in detail in Sec. 5.2.3. Detector non-linearity
The dark- and bias-subtracted CCD signal is often assumed to be a linear function of the incoming light, at leastfor values well below the maximum dynamic range of the camera. However initial photometric analysis of individualframes acquired by several FRAM telescopes, as well as the spurious dependence of measured VAOD on sky background(Ebr et al. 2019b) have demonstrated the presence of an unaccounted effect in the data analysis, which looks like tobe due to significant detector non-linearity, especially in a low light level regime (see Fig. 2). Therefore we developeda dedicated measurement routine, applicable both in the laboratory (Karpov et al. 2018) and on the remote telescopelocation, in order to study and characterize the detector response linearity.In order to characterize the linearity properties of G4-16000 cameras under specifically prepared experimental proto-col we used a dedicated calibration stand in FZU – Institute of Physics in Prague. The stand utilized an uncalibrated405nm LED light source and a screen placed in front of the camera illuminated in a smoothly non-uniform way so thatevery acquired frame covers a range of intensities simultaneously. The measurements at the actual telescope locationwere organized during day time with a closed dome, using scattered light and a section of the inner surface of thedome as a non-uniformly illuminated screen. Such a setup was chosen as it was possible to implement the calibrationmeasurements without physical access to the remote telescope site.The acquisition protocol consisted of obtaining a number of images with different exposures and then studying thedependence of signal level on exposure time. To do so, we acquired a series of “light” measurements, with shutterfully open, each immediately followed by a “dark” measurement of the same length with closed shutter to subtract thebias level of the images. To exclude any possible bias related to changes in the light source intensity, environmentalparameters and camera electronics, we randomly sampled the exposure interval from log-uniform distribution between0.1 s and 300 s. Moreover, to monitor the stability of the light source, after every “light and dark” image pair weacquired a similar pair of “control” images with an exposure interval arbitrarily fixed to 10 s. In the laboratory werepeated the acquisition sequence with different amounts of light reaching the detector. In the telescope experimentsetup, the intensity of light was constantly changing due to day/night cycle. Every run contained several hundreds tothousands of images acquired over long time intervals, from several hours up to several days.To speed up the frame read-out, we acquired only a 1024 × , and computed the mean signal values over these regions to be used in the analysis.We define the detector non-linearity N as the ratio between the measured signal level I obs and the “expected” valuedirectly proportional to the exposure time T exp (corrected for shutter effects as described in (Karpov et al. 2018)), N ( I obs ) = I obs / ( F T exp ) , (1) This procedure is done once for every continuous sequence of frames using in a single medium-intensity image on which sufficiently large(having more than 10000 pixels) sub-regions are identified so that the intensities within them do not vary by more than 1%.
Ebr et al. -0.6-0.4-0.2 0 0.2 0.4 0.6 5 6 7 8 9 10 11 12 13 M ode l - i n s t r u m en t a l [ m ag ] -0.6-0.4-0.2 0 0.2 0.4 0.6 5 6 7 8 9 10 11 12 13 M ode l - i n s t r u m en t a l [ m ag ] B [mag]
Figure 2.
Photometric residuals after fitting the instrumental star magnitudes measured on a single CCD frame with low meanintensity level with a model including Tycho2 (Høg et al. 2000) catalogue magnitudes, color term corrections to account for adifference in photometric system, and low-order spatial polynomial to fix residual uncorrected vignetting. Upper panel – originaldata residuals, lower panel – the residuals after additional correction for detector response non-linearity. Outliers correspondmostly to blended stars where the aperture photometry we use performs badly, as well as to intrinsic photometric uncertaintyof the Tycho2 catalogue and a difference of photometric systems. where F is a scaling coefficient proportional to the incoming-light flux level defined it in such a way that the meannon-linearity is equal to unity over the interval of intensities between 300 and 3000 ADU (analog-to-digital units). Inall our setups the intensity of illumination varied slowly over time, either due to the variations of the light sourceintensity with temperature in the calibration-stand setup, or due to changes of the in-dome illumination related to themotion of the Sun and clouds in on-telescope setup. Therefore, we reconstructed the per-frame values of F using therecorded “control” images all having constant exposure of 10 seconds as F ∝ I obs , /N ( I obs , ) , (2)where I obs , may be interpolated at any given moment, assuming smooth variations of incoming light, and the N term accommodates for the fact that “calibration” frames are also affected by the detector response non-linearity.Finally, we may iteratively solve Eq. (1) together with Eq. (2) to reconstruct the detector non-linearity as a functionof measured signal level. Fig. 3 shows an example result of a single detector characterization, first done directly ona remote telescope site, and, later after the return of the camera back to Prague, in a dark room environment usingCCD overscan regions to better mitigate the variations of a bias level. It is clear there that both measurements providethe same non-linearity curve, despite the difference in experiment conditions and illumination intensities. The cameraresponse is significantly non-linear, with measured signal level directly proportional to the exposure time (manifestingas a horizontal segment on the plot) only in a small part of the dynamical range. Instead, the non-linearity curveshows a characteristic shape consisting of several (two or maybe even three) log-linear segments with different slopes,with non-linearity being the most significant in the low intensity region. erosol Measurement using Wide-field Photometry M ea s u r ed / e x pe c t ed M ea s u r ed / e x pe c t ed Measured signal (bias corrected) [ADU]
Figure 3.
Results of non-linearity measurement for the same CCD camera (serial number 6069), first in the remote telescopelocation (upper panel) and after on the dedicated calibration stand in laboratory (lower panel). The saturation effects areobvious above approximately half of dynamic range ( ∼ ∼ − ). Non-linearity is defined asthe ratio of actually measured signal to the expected one. Every point corresponds to the mean value of a signal over large(tens of thousands of pixels agreeing on the incoming light intensity to 1%) regions of a single frame. The mean values of thenon-linearity are over-plotted in logarithmically spaced bins as violet dots. The outlier points from the mean curve mostlycorrespond to time intervals of rapid changes of incoming light intensity, which could not be properly reconstructed from“calibration” frames, and are insignificant for the determination of the non-linearity. The detector response is non-linear overthe whole dynamic range. The shape of non-linearity curves is similar among all cameras we tested, with a progressive improvement in thosemanufactured more recently, i.e. having higher serial numbers (see Fig. 4). They all show a characteristic structurewith a slope change at around 1000 ADU, which may either be tabulated, or approximately fitted by separate log-linearsegments. Finally, the signal on the detector may be linearized by dividing the observed bias and dark subtractedvalue with N as: I linear = I obs /N ( I obs ) . (3)The lower panel of Fig. 2 shows the effect of a linearization of the data on photometric residuals.2.4. Data processing
The raw CCD frame is first corrected by subtracting dark counts (due to the dark current in CCD chip) and biaslevel (due to A/D electronics) signals. For most purposes, subtracting a dark frame taken at the CCD temperature(or rather, to avoid introducing noise, a median of many dark frames) can satisfactorily remove both dark countsand bias, leaving only a small residual offset due to the bias drift with electronics temperature and history. However,while the dark signal depends only on the (stabilized) temperature of the CCD chip, the bias signal depends on thetemperature of the camera electronics, which are outside of the cold chamber and therefore the bias signal variesboth with environmental conditions and the history of its load. The bias frame is determined either using short darkexposures (the bias does not vary with exposure time, while the dark current for a short exposure is effectively zero) or
Ebr et al. M ea s u r ed / e x pe c t ed Measured signal (bias corrected) [ADU]Dome 06029 - AugerDarkroom 06029 - AugerDome 06069 - CTA-S0Darkroom 06069 - CTA-S0Darkroom 06132 - AugerDarkroom 06149 - CTA-NDarkroom 06204 - AugerDarkroom 06205 - Auger
Figure 4.
Results of non-linearity measurement for a number of CCD cameras used by the FRAM telescopes installed atthe future CTA sites as well as at the Pierre Auger Observatory (Pierre Auger Collaboration 2020). The cameras are basedon the same CCD chip and were mostly measured on a dedicated calibration stand under dark room conditions. There is aprogressive decrease in non-linearity with increasing serial number, which may correspond to the gradual improvements of themanufacturing process and electronics used in the cameras. just by scaling a normal dark frame according to the current camera electronics temperature – as the bias componentis so much larger than the dark, these methods give equivalent results (see Sec. 5). The non-linearity correction (NLC)described above turns out to be extremely sensitive to any residual offset and so the precise bias for each image mustbe subtracted. Originally, this was done using the observed correlation between the bias level and the electronicstemperature, determined over a very large number of dark frames at constant CCD temperature acquired over manymonths of operation (the dark current at the nominal operating temperature − ◦ C is very small compared to the biasvariations for these cameras). Later, the firmware of all our cameras was modified to read out the so-called overscanregion: the area of the CCD chip outside of the nominal imaging area. This area is physically masked from lightillumination and thus provides information on the dark+bias signal simultaneously to any camera acquisition. Theaverage (median) of the overscan region is then subtracted to the relative CCD acquisition frame thus correcting forany drift of the bias with temperature.With the bias and dark components subtracted, each image is linearized using Eq. (3) and then divided by a flat-fielding image produced from the combination of many images taken during dusk/dawn in the direction of the skywith the lowest expected brightness gradient to obtain uniform illumination – this is a standard procedure to correctfor vignetting (which is significant with the lenses used) and variations in detection efficiency across the CCD chip.Astrometric calibration of the images is then performed using the local installation of Astrometry.net (Lang et al.2010) code, which provides a WCS solution typically accurate to better than 1 pixel over the most of the field ofview. The point sources are detected in the calibrated image using the SExtractor software package (Bertin & Arnouts1996), and then identified with stars from the Tycho2 catalogue (Høg et al. 2000) if their predicted position on theimage lies within a short distance (typically 8 pixels) from the position of the detected object centroid. Any pair of erosol Measurement using Wide-field Photometry DATA MODEL3.1.
Model fitting
In general, the measured (instrumental) star brightness m inst is related to the catalogue value m cat , B by the equation(Ebr et al. 2017b) m inst M = m cat , B + Z i + g ( A M , A A , A O , k i , B − V ) + f ( B − V, x, y ) , (4)where f ( ... ) describes corrections due to the different response of the system to stars of different color index (asexpressed by difference between the Tycho B and V magnitudes ) and in different parts of the field of view, and g ( ... )is a model of extinction as a function of the airmass, the color index of the star and the actual value of the VAOD(described in detail in the following section); Z i is the momentary absolute calibration constant of the system (so-calledzeropoint) and M is a correction constant for any possible non-linearity introduced in the data during the photometricprocess. Note that Z i also automatically absorbs any effect due to the uncertainty of the overall normalization of thecatalog fluxes as discussed for example in (Bohlin 2007).The atmosphere is composed of several components with different vertical profiles – for each of these components,the airmass A is a different function of the altitude of the star above the horizon. In our case, we consider g ( ... ) asa function of three different airmass values: A M describes the molecular atmosphere according to (Kasten & Young1989); A A describes the tropospheric aerosols using the formula for water vapor from (Karsten 1965), as the relevantvertical distributions are similar according to (Wehrli 2008); A O is the airmass for ozone accroding to (Komhyr 1980)used in a similar vein as an approximation for the vertical distribution of stratospheric aerosols. The formulae forwater vapor and ozone should in principle be also used for the respective components of the molecular atmosphere,but in the passband given by our B filter, the total contribution of water vapor and ozone to molecular extinctionis only 0.001 and thus the effect of their different airmass formulae is negligible. Each airmass value is calculatedfrom the average altitude of the star above the horizon (corrected for refraction) – a more precisely correct way wouldbe to take the weighted mean of airmasses at several points during the exposure with the weight being given by thetransmissivity of the atmosphere at that airmass, but the difference with respect to a simple average of altitude is lessthan 0.0003 in airmass.The correction for the B − V color term in f ( ... ) is taken as a third order polynomial, but only stars with B − V < . B Johson = B T − . B T − V T ) by up to 0.015 mag for B − V ∈ [ − .
1; 0 .
8] and deviates more foreven bluer stars, which are however rare in the sky.In principle, the spatial dependence of f can reflect residual imperfections of normal flat-fielding correction but anadditional spatial dependence is introduced by the photometric aperture method, in fact the fraction of the total starlight captured by the photometric aperture is strongly dependent on the star’s PSF which, for wide-field setups, is inturn strongly dependent on the star’s position on the chip. Ideally, the spatial dependence would be a simple function ofthe star distance from frame center, but despite our best mechanical efforts, each hardware setup shows a deformationof the pattern of the residuals in some direction, likely due to a combination of non-orthogonality of the lens mountand imperfect adjustment of the orthogonality of the CCD sensor to the optical axis. To find a “photometric flat-field”correction map we first take a set of scans and process it without any further position-dependent correction and thencalculate the residuals of the fit for all stars as a function of the position on the chip and smear this data with a In the following, B − V will always refer to the difference between the magnitudes of a star in B and V filters in the Tycho photometricsystem , i.e. B T and V T (van Leeuwen et al. 1997). Ebr et al. f during the fit). t au_ F R A M B-V 0.00095 0.001 0.00105 0.0011 0.00115 0.0012 0.00125 0.0013-0.4 -0.2 0 0.2 0.4 0.6 0.8 t au_2 B-V
Figure 5.
The coefficients τ FRAM and τ as obtained by fitting a star extinction for individual synthetic stellar spectra as afunction of the spectra B − V color index. Continuous lines are the linear functions used to approximate these two coefficientsas a function of the star color index. Molecular extinction
For light of a given wavelength, the attenuation in magnitudes as a function of airmass A is simply∆ m = 2 . I ground I space = 2 . (exp( τ ( λ )) A ) = 2 . τ ( λ ) A ln 10 ≈ . τ ( λ ) A = k ( λ ) A (5)where k = 1 . τ , is the extinction coefficient in magnitudes. For molecular extinction, we obtain τ ( λ ) using theMODTRAN software package Berk et al. (2014) based on a typical vertical profile of the density of the atmospherefor the FRAM locations from the ECMWF ERA5 data . We are in principle able to obtain the profile specificallyfor the moment of each observation, but the difference with respect to using the mean profile is less than 0.001 inVOD for CTA-N and less than 0.0003 for CTA-S, based on a test calculation using 60 different profiles across oneyear. For CTA-S, this variation may be underestimated in the global models due to the lack of regular vertical profilemeasurements in the area and is subject to current research (Mar´ın et al. 2017).In order to study the dependence of the molecular extinction as a function of the airmass and spectral color weutilize a library of theoretical stellar spectra (Coelho 2014), from which we extract the normalised spectral density erosol Measurement using Wide-field Photometry F ( λ ) for each stellar model, and using the FRAM spectral response as obtained in Sec. 2.2 for R ( λ ), we calculate: I ground I space = (cid:90) F ( λ ) R ( λ ) exp ( − τ ( λ ) A ) d λ ≈ exp( − τ FRAM A + τ A ) , (6)where the other terms in the Taylor expansion are negligible (amounting to difference less than 0.001 optical depthsfrom the exact value for most of the spectra), while the quadratic term is important as a linear approximation wouldcause a difference of up to 0.03 optical depths at the A = 4 for some spectra. Both τ FRAM and τ are to be understoodas functions of the spectral properties of the star. For each stellar model we also compute the color index B − V (usingthe conversions from (Worthey & Lee 2011)), which is a spectral parameter available from the Tycho2 catalogue forevery star. For the region B − V < . τ FRAM and τ can be well parametrized as linear functions of B − V (see Fig. 5), thus the model for molecular extinction in magnitudes as a function of airmass and color index is g M ( A M , B − V ) = k M A M + k C , M A M ( B − V ) − k , M A − k , M A ( B − V ) (7)(the k , M term being a small correction with an effect on VAOD estimation of the order of 0.001). We checkedwhether a more general function would improve over this parametrization, but did not find significant improvementas most of the remaining spread is due to other properties of the stars not reflected by their color index (see Fig. 5).This procedure assumes that all stars from the stellar spectra occur with the same frequency, which is likely not true;moreover the relative contribution of stars with different properties may vary across the sky – the impact of this onthe uncertainty of the VAOD measurement is explored in Sec. 5. For the analysis of data in the B filter, where thecontribution of absorption is small and the molecular extinction is dominated by the Rayleigh scattering which isdirectly proportional to the surface pressure, a linear correction accounting for the variations in the pressure can beadded where pressure data are available. 3.3. Aerosol extinction
The wavelength dependence of optical depth for general aerosols is in general assumed to be expressed by the form τ ( λ ) = τ ( λ/λ ) − α (8)with α ∈ [0; 2] being the ˚Angstr¨om exponent and λ is a reference wavelength. In our case, we put λ = 440 nm whichis the mean wavelength of the light detected by FRAM for a star of B − V = 0 when α = 1 after subtracting molecularextinction – this choice minimises the uncertainty of τ for the given experimental setup. Repeating the integration ofEq. (6) using Eq. (8) for τ ∈ [0; 0 .
2] (the typical range of conditions we find at the CTA sites) and α ∈ [0; 2], we findthat the tropospheric aerosol extinction is well described with the following model (with a single free parameter k A ) g A ( A A , B − V, k A ) = (1 + 0 . α ) k A A A − . α k A A A ( B − V ) − . α k A . (9)This indicates that in the case of unknown ˚Angstr¨om exponent α ∈ [0; 2], assuming α = 1 leads to a relativeuncertainty of the order of 1% in k A , which is small in absolute value (we test this estimate against data in Sec. 5).The quadratic term in this case is negligible for all but the largest aerosol loads, but must be kept to measure thosevalues precisely.Additionally to the tropospheric aerosols, there is, at any moment, also some amount of aerosols present in thestratosphere. For those stratospheric aerosols, we refer to the data from (Chouza et al. 2020) which indicate thatduring the period covered by our measurements, their amount was relatively constant, amounting to τ = 0 . λ = 420 nm using the measured α = 1 .
6. For such a small value all other corrections are negligible and the model issimply g O ( A O ) = 0 . A O . (10)The aerosol extinction model is a sum of the two individual contributions and the total VAOD is thus obtained as τ = 0 . k A +0 . g O is valid only in periods of low global volcanic activity– for other periods, the split between tropospheric and stratospheric aerosols may be more complicated. However,since we are measuring the combined value, the uncertainty of g O only affects the quality of the description of theangular dependence of the aerosol airmass, not the absolute value of the VAOD. When tested on the data used in thenext sections, we find that using A = A M for all components of the extinction leads to a shift of 0.002 in the calculated2 Ebr et al.
VAOD, which is small, but not negligible relative to other uncertainties and thus the distinction between extinctioncomponents should be kept. The uncertainty due to the varying split between tropospheric and stratospheric aerosolswill be even smaller, barring a major global volcanic event.The complete model of atmospheric extinction, obtained as the sum of equations (7), (9) and (10), is thus fairlycomplex. Previously, we have followed a simpler approach, where we not only used a single value of A , but alsosimplified the model to g ( A, B − V, k A ) = ( k A + k G ) A − k C A ( B − V ) − k , M A . (11)where k C (which cannot be just computed for general aerosols) is fitted as a global parameter for one set of scans andthe quadratic term for aerosols is neglected – such a simplification however differs from the more detailed model by upto 0.004 in VAOD when tested on the same data, which is already a significant contribution as can be seen in section5. Figure 6.
Examples of plots used to select clear scans for aerosol analysis from data taken at La Palma. Top left: a very clearscan (VAOD compactible with zero). Top right: a borderline acceptable scan with the best fit giving VAOD around 0.6. Bottomleft: several light clouds, unacceptable scan. Bottom right: large clouds, unacceptable scan. The colors of the points indicatethe mean RMS deviation of stars in an altitude bin from the preliminary fit, the green arrows show the altitude coverage ofimages (to show whether a gap is due to clouds or missing coverage) and stars in red were marked for rejection as outliers. Notethat these are the results of the prelimiary analysis used for scan selection and the spread of the values for individual stars islarger than in the final fit used for the determination of VAOD.
Selection of clear scans erosol Measurement using Wide-field Photometry f ( ... ), as well as M and all the parameters of g ( ... ) that canpotentially vary with time, such as the optical depths for the molecular or stratospheric components, as fixed, basedon a rough estimate or previously processed data (see below), and fitting a ( Z, k ) pair to the scan. These plots (seeFig. 6 for examples) are decorated with visual clues about the mean departures from the fitted value in individual binsas well as other information, such as the altitude coverage of the individual images and identification of stars rejectedby the preliminary fit as outliers, so that they can be inspected efficiently in rapid succession. Then the altitude scansdeemed “clear” are grouped into batches of 50–100 and fitted with Eq. (4) simultaneously, with free parameters beingpairs of (
Z, k ) for each scan while all other parameters are considered as global variables. After the fit converges, thestars for each scan are binned in altitude and outliers above 3 standard deviations in each bin are removed and thefitting is repeated. PHOTOMETRY4.1.
Error model
To compare different possible choices for the photometric extraction method, we need to associate errors withindividual star brightness measurements to be used in the fitting procedure. To this end, we have acquired a series ofseveral thousand images of the same region of the sky with one of the FRAM setups and extracted star fluxes from asub-series of 700 of these for which the atmospheric conditions were exceptionally stable. The image processing wasdone with a photometric extraction software (IRAF with aperture radius 4) with a fixed error of 0.1 mag assignedto every star measurement. Each image was processed using the full pipeline, including near-neighbor filtering andoutlier rejection. (
Z, k ) is fitted on each image – the values produced have large uncertainties, but the goal is testingthe photometry software itself, not the determination of the VOD from this data.For each star that is found and accepted by the pipeline in more than 10 images, the mean differences betweenthe measured and expected star model brightness and the relative RMS are calculated; the former is then comparedwith the error reported by Tycho2 catalogue for the B magnitude while the latter is compared with the mean errorthat the IRAF software associates to the measurements of the star brightness. A systematic difference between themean magnitude difference and the Tycho2 error is somewhat expected as there is additional uncertainty in thecolor correction and it can be modeled by adding a systematic error of 0 . B . That the RMS is larger than theerror reported by the IRAF software means that our measurements are affected by other sources of noise in additionto the uncertainty in the background estimation and the Poissionian fluctuations of the light flux that are alreadyconsidered by IRAF when calculating the error on its measurement. These may be related to insufficient flat-fieldingbut preliminary tests show that at least a part of the effect is due to the dependence of the flux on the sub-pixelposition of the star centroid (the stars’ PSF are too sharply peaked near the camera center, and the CCD we useemploys a microlens raster array in order to improve pixel fill factor). This additional error is modelled by adding anextra error term of 0.028 independent from the star brightness.Overall, we can associate to each star measurement the error as (Fig. 7) σ = σ , N + σ + 0 . + (0 . B ) (12)and similar results hold for SExtractor measured values. This formula can thus be used for estimating the errors foreach star brightness that are used as weights for fitting the extinction model on the altitude scans.4.2. Choice of photometric method
In order to select the best photometric method, we compare the results obtained from the analysis of 100 randomlyselected vertical scans from each of CTA-N (labeled N) and CTA-S sites (two FRAM setups labeled S-WF0 andS-WF1). The S-WF0 camera did not have the overscan readout activated, so its bias subtraction is based only ontemperature parametrization. However, the CTA-N and CTA-S sites are very distinct geographically, so the actualaerosol conditions may be different. During the optimization of the photometric process, some scans appeared to becontaminated by clouds and were removed. The results presented are thus obtained on sets of 91 (from CTA-N) and80 and 92 scans (from each of the CTA-S setups).Figure 8 shows the aperture photometry results for as function of the aperture radius in pixels for both SExtractorand the IRAF apphot package. There are various differences between the two codes, most importantly in the shape ofthe aperture (different approximations for an intersection of a circular aperture with CCD pixel grid, especially evident4
Ebr et al. -0.15-0.1-0.05 0 0.05 0.1 0.15 0.2 6.5 7 7.5 8 8.5 9 9.5 10 10.5 11 e rr o r [ m ag ] B magnitudeMean magnitude differenceTycho error_Bmeasured starsRMS in bin 0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 6.5 7 7.5 8 8.5 9 9.5 10 10.5 11 e rr o r [ m ag ] B magnitudeMagnitude difference RMSIRAF reported errormeasured RMS 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1 6.5 7 7.5 8 8.5 9 9.5 10 10.5 11 e rr o r [ m ag ] B magnitudeStar magnitude error model
Figure 7.
Results from the analysis of a long series of exposures of the same sky region. Top left: mean magnitude differencefrom model vs. the Tycho catalogue B magnitude with RMS (quadratic mean) for each bin shown. Top right: comparison ofthe actual RMS in our star magnitude measurements (green), compared to the one reported by the IRAF photometric software(violet), as a function of star B magnitude. Bottom: the resulting error model for individual star measurements in each skyimage. for smaller radii, see e.g. (Bajaj & Khandrika 2017)), and the method used to extract the background – IRAF employslocal background estimation using an annulus around the star, whereas SExtractor constructs a global backgroundmodel. For both codes we obtain the best agreement of the measured star magnitudes with our model for an aperturesize of 3–6 pixels. We find significant differences in the non-linearity (expressed by the coefficient M in Eq. (4)) whichvaries between FRAM sites and it is generally more pronounced in SExtractor photometry. Disregarding 1 and 2 pixelapertures which produce unstable results (especially for SExtractor), we find that the VAOD variation as a functionof aperture radius is within 0.0025, although there is a roughly 0.003 systematic difference in VAOD as measured byIRAF and SExtractor. We choose the photometry done with IRAF at aperture 4 as the reference method (all furtheruncertainties are evaluated with this selection).Figure 9 shows the dependence of the fitted VAOD on the upper selection cuts imposed on the star magnitude andairmass to be accepted for the fit procedure with respect to the default configuration (used also to produce Fig. 8),which is < . < VAOD SYSTEMATIC UNCERTAINTIES erosol Measurement using Wide-field Photometry -0.01-0.005 0 0.005 0.01 0.015 0.02 0 2 4 6 8 10 12 14 16 VA O D d i ff e r en c e w . r .t. I R A F ape r t u r e IRAF apphotmeanRMS -0.01-0.005 0 0.005 0.01 0.015 0.02 0 2 4 6 8 10 12 14 16SExtractormeanRMS 0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0 2 4 6 8 10 12 14 16 R M S o f r e s i dua l s [ m ag ] NS-WF0S-WF1 0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0 2 4 6 8 10 12 14 16NS-WF0S-WF1 0.98 0.985 0.99 0.995 1 1.005 1.01 1.015 1.02 0 2 4 6 8 10 12 14 16 li nea r i t y c o rr e c t i on aperture radius [pixels]NS-WF0S-WF1 0.98 0.985 0.99 0.995 1 1.005 1.01 1.015 1.02 0 2 4 6 8 10 12 14 16aperture radius [pixels]NS-WF0S-WF1 Figure 8.
Comparison between two photometric software methods – IRAF apphot vs. SExtractor as a function of the apertureradius. The plots in the first row show mean and RMS difference of the VAOD values determined from the overall set of scans,second row shows RMS of residuals for individual stars for all three FRAM setups and the last row shows the non-linearitycoefficient M . We have tested various possible sources of systematic uncertainty in the determination of the VAOD with our methodusing the same set of 263 scans from three different sites as in the previous section, or a subset. These include thefollowing effects: • background estimation, evaluated as the maximal difference for a set of changes in the parameters of the annulusfrom which the background is calculated (the thickness of the annulus set to 4, 12 and 16 pixels instead of thereference value of 8, and the inner radius of the annulus set to 6 and 12 pixels instead of the reference value of9) as well as the change in the “salgorithm” IRAF setting from “mode” to “centroid”, • CCD non-linearity correction, estimated by varying all the parameters of the fit applied to the laboratory databy 1 standard deviation,6
Ebr et al. m a x m ag -0.004-0.002 0 0.002 0.004 m ean s h i ft [ m ag ] m a x m ag -0.004-0.002 0 0.002 0.004 m ean s h i ft [ m ag ] m a x m ag -0.004-0.002 0 0.002 0.004 m ean s h i ft [ m ag ] m a x m ag -0.004-0.002 0 0.002 0.004 m ean s h i ft [ m ag ] m a x m ag -0.004-0.002 0 0.002 0.004 m ean s h i ft [ m ag ] m a x m ag max airmass 5 7 9 11max airmass 5 7 9 11max airmass 5 7 9 11max airmass -0.004-0.002 0 0.002 0.004 m ean s h i ft [ m ag ] aperture_radius=2 aperture_radius=4 aperture_radius=9 aperture_radius=15 I R A F N I R A F S - W F I R A F S - W F SE x t r a c t o r N SE x t r a c t o r S - W F SE x t r a c t o r S - W F Figure 9.
Mean difference in VAOD values obtained from fits with different selection cuts in star brightness and airmass (withrespect to VAOD results for cuts at 9.5 magnitudes in brightness and 8 in airmass) for different photometric aperture radii(columns) and different choices of photometric methods and setup (rows). erosol Measurement using Wide-field Photometry -0.4-0.3-0.2-0.1 0 0.1 0.2 0.3 0.4 360 380 400 420 440 r e l . d i ff e r en c e wavelength [nm] -0.2 0 0.2 0.4 0.6 0.8 1 1.2 300 350 400 450 500 550 t r an s m i ss i on [ % ] wavelength [nm]lenslens-10%lens+10%BB-10%B+10% Figure 10.
Left: comparison of various measurements of the transmissivity of the B filter for different integration times andlight fluxes. Right: the worst case model used to put an upper limit on the uncertainty due to the determination of R ( λ ). • bias subtraction, especially important when the overscan data are not available and thus tested on CTA-S WF0data by varying the calibration curve between the bias signal and the electronics temperature within the rangeallowed by the dark frame data, • “processing” uncertainty given by the VAOD difference obtained between two processing codes which differ inthe choice of reference dark frames as well as in whether they consider dark and bias signal separately, • the numerical “photometric flat-field”, estimated by comparing the results obtained by creating the flat-fieldingdata from two different halves of the test dataset, • the unknown ˚Angstr¨om coefficient, quantified simply by putting α equal to 0 and 2 in Eq. (9), • the uncertainty of the star magnitude error model, compared with the extreme case of same uncertainty for everystar measurement.These effects combined amount to a total uncertainty of less than 0.001 in the determination of the VAOD. Thus thedominant sources of uncertainty are the choice of photometry software and selection cuts on maximal airmass andmagnitude of stars in the fit and the determination of the spectral response R ( λ ) that affects the modeling of molecularextinction.The uncertainty related to the choice of photometry method may be estimated as the maximal difference for IRAFsettings between apertures of 3–6 pixels and SExtractor settings between 3–8 pixels, as those roughly show the minimalRMS of the model description. This uncertainty may be split into a correlated part (corresponding to the shift of themean value of the VAOD) of 0.005 and an uncorrelated part (calculated as the RMS of the differences of individualVAOD shifts with respect to the mean) of 0.003. Similarly for the dependence on the airmass/magnitude selectioncuts, the correlated part is 0.001 and the uncorrelated part is estimated as 0.003 (it is slightly larger for the smallestcuts in airmass and magnitude, but those significantly reduce the number of stars available for some scans and thusincrease the statistical uncertainty).8 Ebr et al.
Determining the uncertainty of the measurement of R ( λ ) is not straightforward. The key question for our purposesis the spectral dependence, not the absolute normalization of the transmittance, which has no effect on VAOD deter-mination from the extinction model. Since in our spectral calibration method the signal measured through the opticalelement is immediately compared to a reference measurement, most wavelength-dependent systematics should cancelout, with the exception of non-linearity of the light detector. We have thus carried out a series of measurements (usingthe B filter as an example) with different settings of integration time and light intensity and compared the relativechanges as a function of wavelength – see left panel of Fig. 10. This shows that the measurements agree within 10%down to 365 nm, where the differences get overwhelmed with statistical fluctuations as the transmissivity of the filteris negligible. The largest systematic effect on the calculation of the molecular extinction, due to its steep dependenceon wavelength, comes from the shape of the short-wavelength slope of R ( λ ). To estimate an upper bound on thisuncertainty, we thus assume a non-realistic worst case scenario, where the measurement shifts by 10% right at theposition of this slope (see right panel of Fig. 10) – even such a situation leads only to a shift of 0.001 in the themolecular extinction, showing that the precision of the methods used to determine R ( λ ) is sufficient for this purpose.Similarly, assuming that the lens transmissivity does not decrease to zero at 365 nm but follows exactly the (rathersmooth) measured data at least down to 350 nm changes the predicted molecular extinction only by 0.0007.The largest uncertainty in the extinction model is that its parameters are determined from a set of stellar spectrawhich may not correspond to the actual stellar population in the area of measurement. From Fig. 5 we estimate themaximal resulting uncertainty in the determination of τ FRAM to be roughly 0.002 in the worst case scenario wherethe the distribution of the measured stars is heavily biased towards a specific part of the spectral library used. Thisuncertainty directly translates to an overall offset in the value of VAOD.The uncorrelated systematic uncertainty should ultimately be summed in squares with the statistical uncertaintyof individual fits, which is 0.001–0.002 for all scans in the test sample, leading to the final estimate of the maximaluncertainty of 0.006 correlated and 0.005 uncorrelated.In Fig. 11 we present two plots that support such precision as realistic. The left panel shows the histogram of theVAOD values for the 263 test scans, starting sharply at the value of the stratospheric contribution (which is alwayspresent) – this feature is particularly sensitive to the correct modeling of the molecular extinction, which for mostvalues of VAOD measured constitutes the vast majority of atmospheric extinction at these wavelengths. In the rightpanel, we show the comparison of the FRAM VAOD values with those measured by a nearby Sun/Moon Photometerin lunar mode (Jurysek et al. 2017). The mean systematic difference between the FRAM and the photometer is 0.005with RMS of 0.007 – consistent with the quoted FRAM uncertainties despite the photometer measurements havinguncertainties of their own (albeit difficult to estimate), which indicates that, apart of possible effects that are correlatedfor both instruments, the FRAM uncertainty cannot be much larger. CONCLUSIONS AND OUTLOOKWe have shown that using wide-field photometry, we can measure the Vertical Atmospheric Optical Depth with aprecision (for a single measurement) better than 0.008 optical depths. This is a relatively cheap and technically simplemethod, especially in environments where the production of laser light at night is undesired, such as astronomicalobservatories. Because the determination of the VAOD with this method is independent on the absolute calibrationof the system, the method is suitable for long-term monitoring, as slow changes in detector response do not bias theresults.The method presented depends on the assumption of horizontal uniformity (stratification) of aerosols. On the selected“clear” scans, the fit gives values of χ per degree of freedom between 1.1–1.4, using the error model developed inSec. 4.1, indicating very good description of the altitude dependence of extinction. However, this might be due tothe fact that any scans performed during conditions of non-homogeneous aerosols is discarded as “cloudy”. This isclearly an area demanding further study, which could be also facilitated through coordination with other atmosphericmonitoring instruments that are or will soon be deployed near the future CTA sites. erosol Measurement using Wide-field Photometry f r equen cy VAOD -0.02 0 0.02 0.04 0.06 0.08 0.1 0.12 0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1 F R A M VA O D Photometer VAOD xx-0.005
Figure 11.
Left: histogram of VAOD values from the test sample of 263 scans from three different setups; values lower thanthe stratospheric contribution are highlighted in green (the binning has been chosen so that the value lies at a bin boundary).Right: comparison between the VAOD measured by the CTA-N FRAM and by a nearby Sun/Moon Photometer (in lunarmode), selected for observations less within 15 minutes apart; the y = x line is added as reference as well as the dashed lineshifted by the mean systematic difference between the values (0.005). ACKNOWLEDGMENTSThis work was supported by European Structural and Investment Fund and the Czech Ministry of Education, Youthand Sports (MEYS) within the projects CZ.02.1.01/0.0/0.0/16 013/0001403 and CZ.02.1.01/0.0/0.0/15 003/0000437)and by the MEYS projects LM2015046, LM2018105 and LTT17006. We appreciate access to computational resourcesprovided by Computing center of the Institute of Physics of the Czech Academy of Sciences (Adam et al. 2020) andsupported by project LG13007. This work was conducted in the context of the CTA Consortium.
Software:
IRAF (Tody 1986), SExtractor (Bertin & Arnouts 1996), MODTRAN (Berk et al. 2014), Astrometry.net(Lang et al. 2010) REFERENCES
Acharya, B. S., Actis, M., Aghajani, T., et al. 2013,Astroparticle Physics, 43, 3,doi: 10.1016/j.astropartphys.2013.01.007Adam, M., Adamov´a, D., Chudoba, J., et al. 2020, EPJWeb Conf., 245, 03034,doi: 10.1051/epjconf/202024503034Baibakov, K., O’Neill, N. T., Ivanescu, L., et al. 2015,Atmospheric Measurement Techniques, 8, 3789,doi: 10.5194/amt-8-3789-2015Bajaj, V., & Khandrika, H. 2017, Comparing AperturePhotometry Software Packages, Space Telescope WFCInstrument Science Report Barreto, A., Rom´an, R., Cuevas, E., et al. 2019,Atmospheric Environment, 202, 190 ,doi: https://doi.org/10.1016/j.atmosenv.2019.01.006Berk, A., Conforti, P., Kennett, R., et al. 2014, inProceedings Volume 9088, Algorithms and Technologiesfor Multispectral, Hyperspectral, and UltraspectralImagery XX; 9, 1–4,doi: 10.1109/WHISPERS.2014.8077573Bertin, E., & Arnouts, S. 1996, A&AS, 117, 393,doi: 10.1051/aas:1996164 Ebr et al. (cid:104) (cid:105) erosol Measurement using Wide-field Photometry21