ABSTRACT Title of Dissertation: SHORT-TERM VARIABILITY OF ATMOSPHERIC EXTINCTION DURING THE NIGHT, UNDER CLEAR-SKY CONDITIONS, INVESTIGATED BY BROADBAND STELLAR PHOTOMETRY Ileana Cristina Musat, Ph.D., 2004 Dissertation Directed By: Prof. Robert G. Ellingson, Department of Meteorology, Florida State University and University of Maryland at College Park This is a study about the possibility of determining the aerosol optical depth by using star broadband observations from a whole sky imager. The main difficulty in such measurements consists of accurately separating the star flux value from the non- stellar diffuse light, which is overwhelmingly present in the whole sky imagery. A correction method to solve this problem is found and the monochromatic extinction at the ground due to aerosols is extracted form heterochromatic measurements. A form of closure is achieved by comparison with simultaneous or temporally close measurements with other instruments. The accuracy and precision of the method are assessed: the total error is a combination of random error of measurements and systematic error of calibration and model and is between 2.6 and 3% rms. SHORT-TERM VARIABILITY OF ATMOSPHERIC EXTINCTION DURING THE NIGHT, UNDER CLEAR-SKY CONDITIONS, INVESTIGATED BY BROADBAND STELLAR PHOTOMETRY By Ileana Cristina Musat Dissertation submitted to the Faculty of the Graduate School of the University of Maryland, College Park, in partial fulfillment of the requirements for the degree of Doctor of Philosophy 2004 Advisory Committee: Professor Robert G. Ellingson, Chair Professor Theodore Rosenberg Professor Robert Hudson Professor Ferdinand Baer Professor Sumant Nigam ? Copyright by Ileana Cristina Musat 2004 Dedication To my parents, Elena and Mihail ii Acknowledgements My heartfelt thanks go to my adviser, Professor Robert G. Ellingson, for his scientific guidance during all my time here, and for allowing me to go somewhat away from meteorology, ? stargazing. His encouragements, sympathy, good faith were the critical mass that happily added to my commitment to this work. I also thank the members of my exam committee, and especially Professor Robert Hudson for an earlier discussion that inspired this work. Professor Ferdinand Baer is thanked for the suggestions that led to a better dissertation. I thank Professor Sumant Nigam for his scientific and human support. I also thank Mr. David Yanuk for his prompt assistance regarding the disk space and other things ?helper?. I would also like to salute the friendly backing from Dr. Gabriela and Eusebiu Catana, and to thank Dr. Corneliu Pop. My deepest gratitude goes to my parents, for their unconditional love and support. Data used in this dissertation were obtained from the Atmospheric Radiation Measurement (ARM) Program sponsored by the U.S. Department of Energy, Office of Science, Office of Biological and Environmental Research, Environmental Sciences Division. iii Table of Contents Dedication..................................................................................................................... ii Acknowledgements......................................................................................................iii Table of Contents......................................................................................................... iv List of Tables ............................................................................................................... vi List of Figures............................................................................................................viii Chapter 1: Introduction................................................................................................. 1 1.1 Aerosol variability and its importance for the climate forcing..................... 1 1.2 Aerosol optical depth .......................................................................................... 2 1.3 Aerosol optical depth variability during the night ? the main aspect addressed in this work ............................................................................................................... 5 1.4 Summary of results ............................................................................................. 8 1.5 Work outline ....................................................................................................... 9 Chapter 2: The nature of the aerosol optical depth determination for direct star irradiance measurement with the Whole Sky Imager. Differences and resemblances of approach as compared with astronomical techniques................................................... 1 2.1 Modeling the star irradiance. Aperture corrections for point-like sources ......... 1 2.2 Point source appearance on WSI images. Phenomenology of the profile widening.................................................................................................................... 4 2.3 Methods of image segmentation. Aperture photometry and point-spread- function approaches. The sky background................................................................ 7 2.4 Summary of Chapter 2...................................................................................... 14 Chapter 3: Direct star irradiance measurement with the Whole Sky Imager ............. 16 3.1 Error in star irradiance. Star signal-to-noise (S/N) and calibrated measurements (CM) accuracy ........................................................................................................ 17 3.1.1 Temperature of the CCD, dark current error, flatfield, bias subtraction error, read noise .................................................................................................. 20 3.1.2 Systematic errors due to absolute flux calibration and geometric calibration in laboratory........................................................................................................ 34 3.1.3 Star flux measurement error: Signal to noise ratio and calibrated measurement error .............................................................................................. 43 3.2 Sky background measurement error: Constant height layer?s contribution to the background sky and extended sources.................................................................... 45 3.2.1 Constant height layer range from zenith.................................................... 45 3.2.2 Background sky correction necessity......................................................... 47 3.3 Summary of Chapter 3...................................................................................... 48 Chapter 4: Star spectra, atmospheric absorption and scattering background, extinction linearization in wide band astronomy, the method of AOD determination from the WSI heterochromatic measurements. ......................................................................... 50 4.1 Wide spectral band measurement linearization in astronomical practice......... 50 4.2 The method of AOD determination for the WSI heterochromatic measurements. The Method of Solving the WSI Equation..................................... 57 4.2.1. Introduction............................................................................................... 57 4.2.2. The WSI equation ..................................................................................... 58 iv 4.2.3 The Method of solving the WSI equation.................................................. 61 4.3 Stellar spectra used in retrieval......................................................................... 68 4.4 Atmospheric scattered light from non-stellar sources ...................................... 73 4.5 Summary of Chapter 4...................................................................................... 78 Chapter 5: Results and uncertainties, conclusions and future work .......................... 80 5.1 Errors calculation ? Validation of the technique of aerosol optical depth calculation............................................................................................................... 80 5.1.1 Systematic error of the atmospheric kernel (SBDART sensitivity to non- solar spectral flux)............................................................................................... 81 5.1.2 Aerosol optical depth computation uncertainties....................................... 89 5.2 Typical Results.................................................................................................. 93 5.3 Result discussion, conclusions and future work ............................................... 99 Appendix A............................................................................................................... 103 Appendix B ............................................................................................................... 156 Bibliography ............................................................................................................. 191 v List of Tables Table 3.1 Parameters for the Gaussian fit to the dark signal for the generic darks provided in the calibrations frames for the years 1998-2003. ?????????25 Table 3.2 Empirical assessment of the dark current for the dates in the first column, from the condition that the electronic noise is 1% of the sky signal. This assessment shows the dark current range, and it is not used in the actual S/N computation?.....31 Table 3.3 The dark current profile as function of the CCD temperature and the frame bias for WSI camera known as Unit 10?????????????????...32 Table 3.4 The dark current profile as function of the CCD temperature and the frame bias for WSI camera known as Unit 4??????????????????.33 Table 3.5 Maximum value of the angular resolution per pixel [?sr/px], for each calibration and unit?????????????????????????..37 Table 3.6 The range (km), along two spherical layers of interest visible with the instrument?????????????????????????????46 Table 4.1 The absolute calibrated flux of stars of 0 mag. and 0.003 mag. at a monochromatic wavelength of 5556A and in the V filter, as measured of quoted in the literature. All fluxes are in units of 10-7 mW m-2 nm-1. References: (a) Oke and Schild, 1970; (b) Tug et al, 1977; c) Wamsteker, 1981; (d) Budding, 1992; (e) Cox, 2000??????????????????????????????..?69 Table 4.2 Example of a few of the stars for which the reddening was considered important to be taken into account???????????????????...72 Table 5.1 Principal molecular absorption bands which fall within the WSI passbands. The minimum, mean and maximum values of the absorber abundance and its vi corresponding optical depth near the center of each molecular band for ozone, water vapor and molecular oxygen. For comparison, the molecular scattering optical depth at the current wavelength is given in parenthesis, in the columns ?f? and ?h?, for the maximum abundance of O3 and H2O, and for the climatologic value of O2 concentration, respectively??????????????????????..87 Table 5.2. Total error for the AOD retrieval from random and systematic uncertainties, given either as standard deviation or, when not mentioned, as relative error. See text for the computation relations????????????????90 vii List of Figures Fig. 2.1 A WSI image towards the North Galactic Pole (near the constellation Coma Berenices, in the center of the image). Constant zenith circles shown: 50, 60, 70, 80 and 90 degrees zenith angle???????????.???????????..9 Fig. 2.2 A WSI image from which the stars have been subtracted using profile-fitting photometry. Constant zenith circles shown: 50, 60, 70, 80 and 90 degrees zenith angle???????????????????????????????.10 Fig. 2.3 A WSI image from which almost all the star residuals have been subtracted. The perfect circularity is due to the atmospheric scattering of the airglow. Constant zenith circles shown: 50, 60, 70, 80 and 90 degrees zenith angle??...?????11 Fig. 2.4 An average of star-subtracted images from 20000208, still exhibiting some trails. Constant zenith circles shown: 50, 60, 70, 80, and 90 degrees zenith angle?.12 Fig. 2.5 An elliptical aperture centered on a star position, with the big axis perpendicular to the direction towards the center of the image????????...13 Fig. 2.6 Transect through star subtracted images when no zodiacal light and the Milky Way are under the horizon. Observation represents the airglow van Rhijn enhancement function and the atmospheric scattering of the airglow. ?????...14 Fig. 3.1 CCD temperature variations during each run, during the nights of interest, in order, for the years 1998(upper left) to 2003(lower right); each graph is clearly labeled with the year for which the CCD operation temperature was represented. The dark current in each pixel is a strong function of CCD operation temperature??...23 Fig. 3.2: Histograms of the distribution of the dark signal, in ADU, for ?generic dark? median files provided with each calibration file, for the period beginning with the viii dates indicated. Each graph is clearly labeled with the date beginning with which the generic dark applies: Upper left: images beginning with 20011128; lower right: images beginning with 20010821 (to 20010827). See Table 3.1 for values of the Gaussian fit????????????????????????????..26 Fig. 3.3: Histograms of the distribution of the dark signal, in ADU, for ?generic dark? median files provided with each calibration file, for the period beginning with the dates indicated. Each graph is clearly labeled with the date beginning with which the generic dark applies. Upper left: images beginning with 20010703; lower right: images beginning with 19981020 (to 19990202). See Table 3.1 for values of the Gaussian fit????????????????????????????..28 Fig. 3.4: Histograms of the distribution of the dark signal, in ADU, for ?generic dark? median files provided with each calibration file, for the period beginning with the dates indicated. Each graph is clearly labeled with the date beginning with which the generic dark applies. Upper left: images beginning with 19981007; upper right: images beginning with 19971010; lower right: images beginning with 19970404. See Table 3.1 for values of the Gaussian fit?????????????????...29 Fig. 3.5 Angular resolution variation with the zenithal angle in the object space, for 4 of the calibrations performed by instrument maintenance during the interval studied. Each graph is clearly labeled with the beginning date for which the calibration applies. Each pixel receives radiation from a solid angle close to 34 ?sr????....38 Fig. 3.6 Angular resolution variation with the zenithal angle in the object space, for 5 of the calibrations performed by instrument maintenance during the interval studied. ix Each graph is clearly labeled with the beginning date for which the calibration applies. Each pixel receives radiation from a solid angle close to 34 ?sr??.??...39 Fig. 3.7 Angular resolution variation with the zenithal angle in the object space, for 5 of the calibrations performed by instrument maintenance during the interval studied. Each graph is clearly labeled with the beginning date for which the calibration applies. Each pixel receives radiation from a solid angle close to 34 ?sr?????40 Fig. 3.8 The geometry of a ground observation of a constant height (H) layer. The range is measured along the layer; the Earth?s mean radius is R. ???????...46 Fig 4.1 Johnson?s system of broadband astronomical filters (from left to right) UBVRI normalized transmission and the response curve of the Whole Sky Imager. The wavelength is given in ?.???????????.??????????51 Fig 4.2 Observed star irradiance variation with respect to the airmass during two nights in Oct. 1999 (Langley plots). Upper: Langley plots for 7 stars during the night of 19991005. Lower: Langley plots for 7 stars during the night of 19991010.??...57 Fig 4.3 Normalized spectra of ?main sequence? stars, each clearly labeled with the spectral and luminosity class in the upper left corner. From upper left to lower right, the star spectra of types B0V, A0V Vega (Alpha Lyrae), F0V, G2V (the Sun), K0V, M0V are represented with the highest continuous line. Also represented is the star spectrum convoluted with the WSI responsivity curve (lowest broader line). For comparison, the spectrum in the Johnson V filter, centered at 550 nm is also represented??????????????????????????.??..70 Fig 4.4 The spectral surface brightness near the Draco Cloud (17h 53m 50.7s; 66 deg 19 arcmin 7 arcsec). Integrated in the WSI filter: 38.32 10-7 [mWm-2sr-2nm-1]?75 x Fig 4.5 Airglow scattered by a Rayleigh atmosphere with scattering optical depth of 0.2, 0.1 and 0.05????????????????????????.??..77 Fig 4.6 Airglow scattered by an aerosol atmosphere with scattering optical depth of 0.2, 0.1 and 0.05????????????????????????.??..78 Fig. 5.1 Normalized spectra of two stars of contrasting types. Left: a main sequence type O9, e.g. Mintaka (Delta Orionis). Right: a supergiant of type M2, e.g. Antares (Alpha Scorpii). The star spectra are represented with the highest continuous line. Also represented is the star spectrum convoluted with the WSI responsivity curve (lowest broader line). For comparison, the spectrum in the Johnson V filter, centered at 550 nm is also represented???????????????????...??85 Fig. 5.2 The atmospheric beam transmittance as function of wavelength 400-1100 nm (continuous line) and the transmittance due to aerosol only (bold line crosses). The aerosol optical depth at 355 nm is, from upper-left to lower-right charts, increasing: 0.07 (19991005, 0530 UT), 0.15 (19991011, 0530 UT) and 0.18 (19991005, 1130 UT).???????????????????????????????..88 Fig. 5.3 Autocorrelation of the AOD extraction model output for two input series of parameters, for three wavelengths at the center of the WSI bandpass: 550, 600, and 650 nm. Stars of all spectral classes have the effective wavelengths in the whole sky imager filter within the interval [550, 650] nm??????????????....91 Fig. 5.4 Autocorrelation of the AOD extraction model output for two input series of parameters, for seven wavelengths not central of the WSI bandpass: 400, 450, 500, 700, 750, 800, 850 nm????????????????????????.92 xi Fig. 5.5 Aerosol optical depth obtained from the same set of measurements using two starting values has better autocorrelation at the center of the passband????..?92 Fig. 5.6 Autocorrelation of the Angstrom parameters alpha and beta and of the simulated AOD at 355nm (which is outside the WSI passband), for the night of 19980904, with two sets of first guess parameters as input?????????.?93 Fig. 5.7 Night of 19980926, total AOD. Left: WSI measured AOD at 600 nm (thick line) and the simulation of AOD at 355 nm (cross symbol) from WSI-inferred aerosol Angstrom exponent and turbidity coefficient. Right: Raman Lidar measurement AOD at 355 nm (cross symbol and line between points) from backscatter at the monochromatic laser wavelength and simulation with the WSI (diamond symbol)??????????????????????????????94 Fig. 5.8 The night of 20020307. The comparative presentation of the RL AOD at 355 nm (line) and the AOD at 355 nm simulated from WSI (diamond symbol) inferred aerosol parameters. The alpha exponent value is around the value of 2.3, which is characteristic for pollution aerosols with diameters between 0.1 and 1 ?m????95 Fig. 5.9 The night of 20011115: AOD at 355 nm from RL (line) and from WSI (diamond symbol)??????????????????????????95 Fig. 5.10 The night of 20011220: comparison between AOD at 355 nm from RL (line) and from WSI (diamond symbol)????????????????.?..95 Fig. 5.11 Alfa parameter during the following daytime, from the Sun photometer. It can be seen that the curvature of the spectral AOD curve of the aerosols is not constant: it depends on the pair of wavelengths between which is computed???.96 xii Fig. 5.12 Night of 20000110: RL AOD at 355 nm (line) and the AOD at 355 nm simulated from WSI inferred aerosol parameters (diamond symbol). The alpha exponent value is around the value of 2.5, which is not close to the AOS simultaneous value of around 0.5, which is characteristic for larger particles???????..?96 Fig. 5.13 Night of 19980519, a smoke advection episode. WSI simulated 355 nm AOD and RL AOD at 355 nm measurement. AOD is over 0.8????????...97 Fig. 5.14 The night of 20030329. RL 355 nm (line) AOD measurement and WSI (diamond symbol) simulated AOD at 355 nm (excluding two episodes with clouds around 07-08 Hours)??????????????????????????????..97 Fig. 5.15 Night of 20030329: Alpha aerosol exponent variation in time from WSI (columnar) and from AOS (ground)???????????????????98 Fig. 5.16 Aerosol Angstrom exponent alpha measured at the ground by the aerosol observing system (gray line) at the beginning of the night, as compared with the same parameter inferred from the WSI (black squares). The straight line in the AOS alpha value is a defect measurement. No other simultaneous alpha estimations exist??..98 Fig. 5.17 Night of 20011024: apparent 4 hours lag between the RL 355 nm AOD measurement (line) and simulated 355 nm AOD from WSI star measurements (diamond symbol)??????????????????????????99 xiii Chapter 1: Introduction One of the fundamental attributes of the Earth?s climate is the radiative energy balance between the incoming shortwave radiative energy from the Sun and the energy absorbed by the atmosphere and Earth and radiated back into space as infrared radiation. This balance may be perturbed by small changes that occur during geological times in the Sun?s incoming energy and the Earth?s orbital parameters (both, external forcings) and by small changes in the composition and physical properties of the atmosphere, ocean and land distribution (internal forcings) (Thomas and Stamnes, 1999). There are also changes whose time scales are much shorter, and they refer to variations in intensive/extensive physical parameters of some of the atmospheric constituents, such as changes in aerosol composition and loading, in trace gases composition and abundances, and cloud cover extent and albedo. Incorporating an accurate mechanism describing the internal forcing of the climate, which is due to the atmospheric aerosols is one of the main components of the current climate modeling effort. General circulation models, energy balance models and radiative-convective models have been used with success to compute the influence of aerosols on the spatial and temporal variance of the surface temperature and of the warming/cooling rates, provided that the temporal and spatial variability of the aerosols is accurately known (Lenoble, 1984). 1 1.1 Aerosol variability and its importance for the climate forcing Extensive sets of aerosol observations from instruments either ground-based, air- or satellite-borne are used to construct climatology sets that can characterize their spatial and temporal presence in relation with their type. For example, the models constructed from observations by Shettle (1979) for ?rural, urban and marine aerosol models? in the boundary layer and the ?tropospheric aerosol model? for the aerosols above the boundary layer take into account a certain composition mixture (for rural model, for example, this is a 70% water-soluble and a 30% dust-like composition) to establish the distribution of the particle density with respect to their radius (in this example, a sum of two log-normal distributions), and the weighted refractive index of the mixture. Thus, from the number density and refractive index, using the Mie theory of scattering from spherical particles, the optical properties of aerosol models relevant to the aerosol radiative effects are calculated: the extinction coefficient, the albedo of a single scattering, and the asymetry parameter, i.e. a weighted average of the aerosol scattering phase function (Shettle, 1979). The direct radiative effect of the aerosols on the climate refers to the changes introduced in the radiative field by the presence of the aerosols in a molecular atmosphere (?clear-sky? case); they can increase the albedo, scatter the solar energy back and thus cool the ground and the troposphere (Lenoble, 1984). The indirect aerosol effect refers to the modifications introduced by the presence of the aerosols within or near the clouds and consists in the increase of the number concentration of cloud droplets and in the change of the liquid water path and thus precipitation efficiency of clouds. 1 Different radiative effects due to the different composition of aerosols (mostly scatterers, like water vapor and dust, or absorbers in visible like soot etc.) have been reported from computations with atmospheric models, for stratospheric and tropospheric aerosols, as well as a change in the aerosol direct forcing when in presence of clouds. For example, a soot (heavily absorbing) aerosol above a cloud produces a direct positive radiative forcing of 1.29 W/m2, as compared with only 0.24 W/m2 in the no-cloud case, while a sulphate (mostly scattering) aerosol produces a direct radiative forcing of -1.92 W/m2 in clear-sky conditions, as compared with a -0.50 W/ m2 in cloudy sky conditions, thus has a less pronounced negative direct forcing effect in the presence of clouds (Highwood, 2000). 1.2 Aerosol optical depth A key physical parameter used in all radiative models is the optical depth of the aerosols. The sensitivity of the radiative shortwave models to the aerosol optical depth variability, as compared to all the other atmospheric parameters in a molecular (clear-sky) atmosphere, is the highest (Halthorne et al, 1997). The aerosol optical depth (or thickness), AOD, is an extensive aerosol parameter and represents the extinction (absorption plus scattering) coefficient integrated over a certain path in the aerosol layer and, when the integration is done on the vertical path comprising the entire atmosphere, it is a measure of the aerosol loading present in the observed column. Aerosol optical depth is directly measurable by sun photometers from the ground from the beam transmittance (Harrison and Michalsky, 1994; Bergin et al, 2000), by air-borne sun photometers aboard an aircraft (Hartley et al, 2000; 2 Redemann et al, 2000), or by satellites equipped with a high resolution radiometer from the backscattered radiation above an optically known (mainly, ocean) surface (Tegen et al, 1997). Aerosol optical depth may be obtained indirectly from in-situ nephelometer measurements of absorption and scattering coefficients (Hartley et al, 2000; Bahreini et al, 2003). Spectral aerosol optical depth can serve further to obtain intensive aerosol optical parameters, such as the vertical distribution of the refraction index (Redemann et al, 2000). Another indirect method for determining the aerosol optical depth is to calculate it as the product between the aerosol mass load, i.e. the column mass per surface area computed from transport models applied to known sources, and the specific extinction cross section; this method is suitable for comparisons with satellite AOD retrievals, because is gives rather a global AOD assessment (Tegen et al, 1997). In addition to the optical depth, the aerosols are characterized by intensive quantities quantifying their scattering properties, such as the single scattering albedo (fraction of scattered light in total aerosol extinction), the scattering spectral dependence (the Angstrom exponent, ? ), the scattering angular dependence (the phase function and the asymmetry factor), and the hygroscopic growth factor (a relative humidity multiplier relating size of aerosol with relative humidity of the air) (Ogren, 1995). The principal AOD observational technique used in the in-situ experiments mentioned above and better known as the ACE, TARFOX, etc. experiments, was the inference of the AOD from the solar beam transmission through the atmosphere from above the aircraft to the top of the atmosphere (TOA). The ground-based 3 monochromatic observations with the photometers also use the solar beam transmission to infer the total columnar AOD, according to the Beer?s law of extinction, from what is called ?Langley plots? (Harrison and Michalsky, 1994). For example, the sun CIMEL photometer data, obtained at the Atmospheric Radiation Measurement (ARM) program (Stokes and Schwartz, 1994) at locations known as the Cloud and Atmosphere Radiation Testbed (CART), are used to obtain total atmospheric optical depth at the wavelengths of 340, 380, 440, 500, 670, 870, and 940 nm from the slope of a straight line fit of the sun spectral irradiance measurements versus, roughly, the secant of the solar zenith angle. As the sun apparently moves on the sky during daytime, the CIMEL set of observations provide instantaneous values of the optical path which are then reduced to the (zenith direction) optical depth, and the monochromatic AOD is obtained from the total optical depth, by subtracting the contribution of the molecular scattering and absorption. Apart from satellite measurements, which comparatively have a larger spatial scale (and a lower spatial resolution), a high spatial resolution of AOD observation can be achieved only by in-situ and/or ground-local daytime observations. Moreover, the ground-based aerosol data, such as, for example, that produced by the direct and diffuse atmospheric radiation measured by a network of sun-sky scanning spectral radiometers are useful in validating the aerosol products obtained by satellite measurements and can be inverted to obtain the aerosol size distribution and refractive index (Dubovik et al, 2000). 4 1.3 Aerosol optical depth variability during the night ? the main aspect addressed in this work Due to the aerosol high spatial and temporal varibilities in both composition and optical properties, a comprehensive small-scale modeling effort requires continuous observation, even when the sun-photometers and radiometers cannot be used, i.e. during the night. Sometimes there are no clear daytime periods, but there may be clear periods at night, where observations might be made with other instruments to augment the daytime ones. Diurnal variability studies are not a novelty; for example, the outgoing longwave atmospheric radiation was studied using satellite measurements (Ellingson et al, 2002). The variability of the aerosol during the night can currently be known operationally only from lidar measurements. Raman Lidar measurements of the aerosol backscatter vertical profile are currently used at the ARM CART site for estimating the AOD during the night (Whiteman et al, 1992). This powerful non- passive technique is able to obtain high spatial and temporal resolution data for water vapor and aerosols. The aerosol optical depth is computed at the laser emiting wavelength of 354.7 nm. At night, the atmosphere is illuminated by thousands of stars that traverse overhead. If one knew the star irradiance and could measure the time dependence of the monochromatic illumination at the surface, it would be possible to perform instantaneous Langley estimates of AOD ? perhaps even better than in the daytime, because of the large number of sources and the wide range of stellar zenith angles. Such time-dependent stellar measurements are not routinely available. However, the 5 Atmospheric Radiation Measurements program obtains broadband measurements (0.4-0.85 ?m) of stellar irradiance at the ground using calibrated Whole Sky Imagers (WSIs) at their various source sites. ARM uses these data to estimate nighttime cloud cover. The angular resolution of these data is nevertheless good enough to identify many individual stars, even without a system to compensate for the Earth?s rotation. However, these instruments have not been used heretofore to estimate night-time aerosols, primarily because of the difficulties associated with tracking the stars and removing the stray light from non-stellar sources, and because the inference of AOD from broadband observations is non-trivial. Of course, a broadband passive retrieval cannot rival in accuracy an active method as that of the nighttime lidar. Even at high altitude astronomical observatories, where boundary layer aerosol and a good part of water vapor columnar content are not present, stellar magnitudes in the medium band filters of an usual photometric system (e.g. Johnson, Cousins, Geneva, Stromgren, Washington etc.) would have to be determined to the 3 rd decimal (1 millimagnitude) to be suitable to rival the lidar extinction measurements (Young et al, 1991). Nevertheless, even if the accuracy could not be a match with that of the lidar, the nighttime AOD from WSI data might be probably less costly to obtain when deploying the instrument in places where a lidar system does not operate. The main objective of this work is to provide a frame for studies of temporal variability during the night of aerosols, by developing a realistic way to extract aerosol optical depth variability using starlight. The instrumental data we use come as radiance fields from the Whole Sky Imager of the ARM, at the Southern Great Plains (SGP) site in Oklahoma. The instrument is used mainly for observing the cloud cover, 6 such that, by design, determining the aerosol optical depth during the night is not its primary task. The temporal scale of the observed radiance field is 1 minute, with a 6- minute interval between successive observations. The instrument has a broadband spectral responsivity, and can measure only heterochromatic star irradiances. Solar photometers and radiometers cannot function during the night, and the Raman lidar (RL) deployed at the site is the only other instrument, beside WSI, which can offer an assessment of the aerosol optical depth. The computational strategy for obtaining the AOD under starlight conditions and from a field of view comprising practically the entire celestial vault has much to do with usual star photometry accomplished at the ground-based astronomical telescopes, but it is also very different from that in the matter of separating the point source light contribution from the contribution of extended light sources of the night, like the airglow (AG), zodiacal light (ZL) and integrated starlight (ISL). The method is suitable for further use for the WSI deployed in locations where no other aerosol measuring instruments are available. Because a point inside an aerosol layer situated at, say, 5 km altitude is seen by the WSI at a range (i.e., distance from the zenith, measured on the spherical layer) of 253 km if the point is at the horizon, and at a range of 8.4 km if the point is at 50 degrees zenithal distance, the spatial scale of the WSI observations can be considered mesoscale- specific and thus within the limits of spatial coherence of AOD (Anderson et al, 2003). Within the current work, a closure experiment was also attempted by using simultaneous measurements with the RL of the monochromatic aerosol optical depth and the alpha Angstrom spectral exponent (defined above) for the aerosols at the 7 ground level - which is measured by the Aerosol Observing System (AOS), to simulate the star radiation at the ground level. It has been found that the alpha parameter measured at the ground does not always provide the best fit for the WSI heterochromatic data, an explanation being that, within the column, the aerosol spectral properties cannot be accurately characterized by the ground-captured aerosols. This result will be discussed in more detail in Chapter 5. 1.4 Summary of results Simple Langley plots, as those derived from monochromatic observations with photometers, do not yield an accurate aerosol optical depth. We develop a method for obtaining aerosol optical depth from heterochromatic data. Boundary layer aerosols and troposheric aerosols can be accurately obtained from the WSI heterochromatic measurements because the observational spatial and temporal scales with this instrument are within the respective coherence scales for aerosol optical depth. The degree of heterochromaticity introduced by aerosol contribution is small and its weight is more accentuated towards each spectrum?s center of mass. The AOD dependence on the wavelength obtained implies interpolation towards the centre of the bandwidth (500; 600 nm) and extrapolation towards the margins. The parameters describing the aerosol optical depth variation with wavelength, (AOD=? ? -? ), alpha and beta, are thus very sensitive to the extrapolation and are not as reliable as the AOD at the centre of bandwidth. We consider the AOD at around 550-600nm (or at what the centre value is) as the most accurate figure in our non-monochromatic result. 8 In astronomical parlance, the core of this work is obtaining calibration equations for standard stars and color equations for nonstandard stars in a specific instrumentation system. In radiometric parlance, the method in this work is that of small perturbations employed by Rayleigh, i.e. adding together several thin slabs to obtain a fit for a function whose spectral shape is known. The Rayleigh method pertains to a class of techniques known as ?renormalization? techniques ? the summing expansions in small parameters to produce a uniformly valid expansion. Mathematically, the method in this work relates to the perturbation theory for integro-differential equations, and consists of obtaining from a nonuniformly valid expansion, an uniform one, i.e. expressing the optical depth as an asymptotic series expansion which maintains accuracy properties when the series is truncated to a reasonable number of terms. Why wasn?t this done before at astronomical observatories: lack of a deeper atmospheric insight, which only a photometric/radiometric test-site can offer; observations made usually in astronomic narrower-than-WSI filters, in standardized historical filter systems; little background atmospheric scattering in astronomical narrow field?of-view images; no interest in uncompensated star tracking observations. 1.5 Work outline The present work is organized as follows: Chapter 2 includes a review of the star photometry techniques and their applicability for the star fluxes measured with the WSI. The model for the star 9 irradiance is developed by taking into account the point-spread-function (PSF) variability across the image, the undersampling introduced by the detector, the particularity of the wide field of view observation with regard to the background gradients of the sky radiance. Atmospheric scattering of non-stellar sources is assessed as the main difficulty in star absolute flux computation. Chapter 3 presents the instrumental uncertainties arising from the charge coupled device (CCD) performance with regard to the physical conditions of the experimental set-up, the knowledge of the optical system transmission, the scale of the sky-projected image of the detection unit (pixel) and the systematic errors in the instrument geometric and absolute flux calibration. Chapter 4 describes the additional data regarding the stellar instrinsic spectra at the top of the atmosphere, the atmospheric absorption and scattering, and the model for computing the aerosol optical depth from multiple observations of spectrally different stars with a heterochromatic detector. The night sky radiance originating from other sources acting as background for starlight is treated extensively. Results from the previous chapters are used in the Chapter 5 to obtain the aerosol optical depth for a little under 300 nights in which a sky free of clouds could be observed, spanning intervals from under 1 hour up to 11 consecutive hours each night. A comparison is made with the AOD determined with the Raman Lidar (which is outside WSI spectral range), and with the nearest daytime AOD recorded monochromatically by sun photometers. The reasons for differences between WSI and the other AOD measurements are explained through arguments regarding the 10 different approaches in investigating the nocturnal boundary layer and the spatial coherence scale for aerosols. 11 Chapter 2: The nature of the aerosol optical depth determination for direct star irradiance measurement with the Whole Sky Imager. Differences and resemblances of approach as compared with astronomical techniques The subject of star imaging through an astronomical optical system is extremely vast, and covers domains such as the point spread function (PSF) of an optical system, the pixel sampling of an image, the historical astronomical set of filters, the methods for image segmentation, techniques for star tracking, the automatic analysis of multiple star frames. In this chapter we will assume that such subjects are thoroughly known, and, in this context, we will only evaluate the Whole Sky Imager as an instrument for performing astronomical photometry. The differences and resemblances with more usual astronomical instruments will be evaluated. 2.1 Modeling the star irradiance. Aperture corrections for point-like sources The star irradiance observed at the ground, F(z 0 ), or the ?direct radiance?, as it is called in atmospheric radiation studies, is the part of the star flux which is transmitted through the atmosphere: F(z 0 )=F(TOA) exp(-? /?), (2.1) 1 where F(TOA) is the star flux at the top of the atmosphere, ? is the optical depth of the atmosphere and 1/? is the geometrical path length (Thomas and Stamnes, 1999). Because of atmospheric scattering, the star might appear surrounded by a halo containing the star light scattered within a small angular distance from the original star image. In conformity with the causes that produced this scattering, the approaches to the star measurement are different: ? in astronomy, where the star is observed from above the boundary layer aerosols (2-5 km), the observer will try to include the halo, because it is caused mostly by the imperfections in the telescope and the dust between the optical components; a synthetic aperture with variable radius is thus employed to obtain ?growth curves? for the star flux, to yield all the halo starlight within the measured flux (Stetson, 1990b; Howell, 1989); the growth curves show that the flux above the sky level increases with increasing the aperture around the star, until it reaches a saturation level, at about 5 degrees angular distance from the star center. In view of the small field of view and fine angular resolution of a telescope image, this may mean that the total flux is completely extracted from a 40 pixels radius around the star. ? in Sun photometry, the halo observed around the sun will be excluded from any computation of direct normal star irradiance, because it is produced by the atmospheric aerosol scattering the direct sunlight in a close-to-forward direction, - and this flux loss has been already accounted as extinction through the scattering included in the total optical depth ? . 2 ? in WSI photometry, one would want to exclude the halo around the star, and also extract the right level of the star flux above the background sky; that is, the star flux should be extracted with a radius kept at the exact ?dimension? of the star, and the sky surrounding the star should not contribute with spurious radiation originating from sources which are not this star currently under observation and scattered within the star subtended solid angle. As it can be seen, the nature of a WSI star image extraction is completely different from both the one in astronomy and in Sun photometry. The particularities of the whole sky imager detector and optical system are such that the data suffer from lack of the spatial resolution required in a common astronomical observation. For example, if the exposure time were a little longer, the stars that blend in the background would become visible, their profile would emerge from the sky background to crowd a field already full of crowded stars. The close pair in Sagitta constellation, used to check the eye separation powers, shows that the instrument separation powers are slightly less than that of the human eye. But the WSI dynamic range, expressing the instrument?s sensitivity to small changes in the irradiance of the source is instead extremely high, due to a silicon diode detector type (CCD of silicon, frontal illuminated, quantum efficiency around 40%). The description of the instrument, its electronic and optical components is available from the ARM website (Whole Sky Imager (WSI), 2001), and has been detailed elsewhere (Musat, 2000). 3 2.2 Point source appearance on WSI images. Phenomenology of the profile widening In this subsection we will evaluate the causes for which the spread of the star image, which is smaller when the star is at the zenith from the situation when the star is at lower elevation, is caused by turbulence, by the star movement during the exposure time or by the off-axis viewing through the optical system? Stellar irradiance (flux) is proportional to the solid angle spanned by the source as seen from the Earth. From interferometry studies it is known, for example, that the true diameter of Vega (i.e. above the atmosphere, and neglecting the interstellar absorption) is 3x10-3 arcsec; this is the diameter which is used to verify whether a model of the star interior and photosphere is correct. But in a telescope with an entrance aperture of 16 arcsec, the stellar diameter of the ?seeing? disk is around 2 arcsec, for exposures up to 10 minutes. The shape of a star (King, 1983) on a WSI image is the consequence of the nature of the optical system and the length of observation (60s), which smoothes the profile of ?seeing? due to the atmospheric turbulence. (The ?seeing? and the scintillation are the effects of atmospheric turbulence on the refraction index and density of the air, respectively.) The diameter in a 60 s exposure WSI image is around 3 pixels, i.e. 0.9 arcdegree, if observed in the middle of the field of view. An image is usually widen by the following effects: ? atmospheric turbulence (coherence time less than 40 msec, in visible). The coherence time is defined as the time scale for which the turbulence is stationary, such that the speckle pattern of a star is not blurred during the exposure. 4 ? differential refraction along the atmospheric path of the polychromatic sources (at most, tens of arcseconds widening) (Stone, 1989), ? the non-compensated apparent diurnal rotation of the star within the field of view of the instrument (at most, a displacement of the center of the star kept within a pixel), ? the partial pressure of the columnar water vapor in the line of sight (Fluks and The, 1992), the scattering of the source from molecules and aerosols within the atmosphere, creating the aureole, mostly within the ring around the 4 th pixel distance from the center of a bright star. The diameter dependence on the partial pressure of the water vapor in the atmosphere was actually tested for the WSI images, and no correlation was found; Fluks and The (1992) measured this dependence at an altitude of 2300m, where the water vapor may already have been condensed into an aerosol haze which widened the image. The light from the source might also suffer scattering inside the instrument optical components: the transparent dome might scatter an amount of the radiation, as well as the glass of the fish-eye lens. No losses of incoming flux through scatter are registered inside the fiber optic placed between the lens and CCD. The instrument response curve, computed at the calibration, accounts for all transmission losses from the top of the dome to the CCD plane. The important angular diameter of the source is thus not due to actually resolving the star disk through the instrument, which would mean knowing the angular distribution of the intensity in the source, modulated by the atmospheric and 5 instrument optical transfer functions. It is instead the result of superposition of several dispersion producing phenomena. The fish-eye lens is just the prime surface of the optical system and has the property that it gathers light from as far as 90-91 degrees off the optical axis and refracts it such that the refraction angle is much lower; the fish-eye lens is followed by a system of other lenses that introduce corrections, such as ? for example - an achromatic doublet, which brings the red and blue light into a common focus. The f- number is the ratio of the (effective) focal distance to the aperture diameter. The WSI fish-eye lens is a Nikkor-type, f/2.8 lens, with the focal distance of 8 mm. The fish-eye lens, for which the ray intercept curves and ray aberration (Aldering, 1991) charts are given by Miyamoto (1964) and by Smith (1990), has a under-corrected spherical aberration for the yellow (called ?secondary spectrum?), while the red and blue have a common focus; this leads to a difference between the height in the image plane of diverse monochromatic rays coming from an off-axis point. Miyamoto patented a fish-eye lens in 1964, which made possible the correction for the large lateral color aberration, by adding a cemented doublet between the meniscus lens and the pupil; this made it possible to obtain color photographs through this lens. Kumler (2000) measured a Nikkor f/2.8, 8mm focal length fish-eye lens camera with a lateral color at 90 degrees of 16 microns, which, for a lens radius of 6.15 mm, means only 0.26%. Because the WSI has a fiber optics taper which reduces the lens-created image dimension by a factor of 2.47, the 16 ?m of maximum lateral color aberration are in the final image reduced to 6.5 ?m, which is only 1/3 of the lateral dimension of one CCD pixel (19 ?m); the lateral color decreases from this 6 value (1/3 pixel at the horizon, when imaging the sky vault), towards the center of the image (zenith direction). The most important optical aberration is the radial distortion, an increasing of the image height with the radius (3.8% at 90 degrees field of view, Kumler et al 2000). The radial distortion in the image plane is a cubic function of radial distance in the image, and is equivalent to a change of magnification in the image plane (Mamoulis and Macdonald, 1997). This is a very important fact when projecting the constant-height airglow layer. The layer is projected with a pincushion hyperbolic distortion. The aberration- due variations of the PSF across the image are still small, and are due mainly to the sphero-chromatism in the system. The observed star PSF is the result of multiple convolutions with ?degrading? factors: diffraction, ray and wave aberrations, ?seeing?, tracking error, pixel sampling, pixel non-uniformity due to differential absorption in the dielectric and in the metallic conductor layers built on its surface etc., and in the infinite limit, tends to a Gaussian function, conform the central limit theorem (Harvey, 2000). The observed PSF encloses the remainder of the original source energy. 2.3 Methods of image segmentation. Aperture photometry and point- spread-function approaches. The sky background Both aperture photometry and profile fitting of a PSF were applied to the WSI frames. 7 Detailed accounts of these comprehensive methods are given in : Stetson, 1987; Howell et al, 1989; Sterken and Manfroid, 1992a; Gullixon, 1992; Mighell, 1989; Allard, 1998. A point spread function fitting implies the detection of the star?s centroid in the image, the placing of a Gaussian-shaped surface with the top exactly over that center and the fitting of its two dimensional width until the volume of the Gaussian is equal (or within a preset error) to the star flux value. The resulting Gaussian widths are usually not constant across the WSI image. Because of the radial symmetry of the spherical aberrations, the two-dimensional width changes with the radial distance in the image (i.e., with the zenith angle), and is given by a skewed Gaussian. The sharp definition of the star image would favor aperture photometry, with the exception of a few instances when the system was defocused and the separation between the stars was extremely low. The sigma-width (i.e., the full width at half maximum/2.3548) of the Gaussian profile is small for the WSI images: only 0.45 pixels, and sometimes, when the seeing is better, even 0.39 pixels. The limit for undersampling is 0.85 pixels. The stars as seen by WSI are thus undersampled (Stetson et al, 1990), and a good PSF fit would imply ?masking? (Buonnano et al, 1989; Howell et al, 1996); this is the procedure of finding the star center inside the pixel which shows the maximum radiance value, without actually being possible to resolve underpixel details observationally. Also, another solution, suggested by Mighell (1989), was that a way to recover the Nyquist sampling frequency from under-sampled data is to renounce altogether the Gaussian function and use directly a ?digitized? (matrix) PSF, with twice the resolution of the 8 image. This means renouncing the advantages brought by having an analytical, easily integrable and interpolable function as the PSF. Examples of PSF fitting and aperture photometry applied to the WSI images are given in the Fig. 2.1-2.6. In Fig. 2.1 the only extended source of the background is the airglow, because star zenith in the image points towards the North Galactic Pole ( NGP ). Fig. 2.1 A WSI image towards the North Galactic Pole (near the constellation Coma Berenices, in the center of the image). Constant zenith circles shown: 50, 60, 70, 80 and 90 degrees zenith angle. The Milky Way is slightly under the horizon, and the solar zenith angle is much over 90 degrees (?dip?). Consequently, both the integrated starlight (i.e., non- resolved stars) and the zodiacal light are not present, and the subtraction of the stars let only the atmospheric scattering of the airglow to be visible (Fig. 2.2 and Fig. 2.3). 9 Fig. 2.2 A WSI image from which the stars have been subtracted using profile fitting photometry. Constant zenith circles shown: 50, 60, 70, 80 and 90 degrees zenith angle. 10 Fig. 2.3 A WSI image in which almost all the star residuals have been subtracted. The perfect circularity is due to the atmospheric scattering of the airglow. Constant zenith circles shown: 50, 60, 70, 80 and 90 degrees zenith angle. When the extended sources (but airglow) are missing from the vault, an average can be made from a set of star-extracted images (Fig. 2.4). Sky radiance can be then fitted with a zenith (i.e., radial) function, to obtain the airglow plus the atmospheric scattering of the airglow. The atmospheric scattering of the airglow is a sum of Rayleigh scattering, including the multiple scattering that finally reaches the observer, and the Mie scattering, which is more ?localized?, in the sense that peaks strongly in a small solid angle near the source?s direction. 11 Fig. 2.4 An average of star-subtracted images from 20000208, still exhibiting some trails. Constant zenith circles shown: 50, 60, 70, 80 and 90 degrees zenith angle. An example of an elliptical aperture is given in the Fig. 2.5. This sort of aperture for finding the sky median value around a star position can minimize the weight of radial gradients of the atmospheric scattered airglow. 12 Fig. 2.5 An elliptical aperture centered on a star position, with the big axis perpendicular to the direction towards the center of the image. Fig. 2.6 presents a transect through star-subtracted images, when no zodiacal light is present and the Milky Way is under the horizon. The observation represents the airglow van Rhijn enhancement function and the atmospheric scattering of the airglow. The conclusion is that the WSI offers an image of the sky that has a background non-uniformity due to the airglow, zodiacal light and the unresolved stars of he Milky Way (?the integrated starlight?), plus their scattering and absorption by the atmosphere. 13 Fig. 2.6 Transect through star subtracted images when no zodiacal light and the Milky Way are under the horizon. Observation represents the airglow van Rhijn enhancement function and the atmospheric scattering of the airglow. 2.4 Summary of Chapter 2 In the analysis presented in this chapter, the nature of the star image, as seen by the WSI instrument, was found to be a combination of dispersive phenomena in the optical system and atmosphere, and a limited resolution. Their summary is as follows: ? optical PSF is due mainly to distortion and spherochromatism; ? atmosphere is further attenuating the PSF encircled energy and is broadening the PSF through extinction and refraction, respectively, which are in turn affected by the atmospheric turbulence, producing scintillation and ?seeing? phenomena; ? optical PSF is also broadened by the heterochromatic nature of the observation; 14 ? total instrumental PSF suffers from the poor resolution (i.e., spatially under- sampling of the star image, and low separation between stars), due mainly to the coarse pixelization. In essence, it has been shown that this is an instrument with low resolution, high background light and coarse pixelization. Both aperture photometry and PSF fitting photometry have been used to obtain star fluxes from WSI nighttime images. The advantage of using elliptical aperture for median sky computation consists in both the maximization of the number of pixels like the ones under the star, and minimization of the sky radial gradients. A median sky was also constructed for several nights, by replacing the stars in the images with the median or mode sky values nearby. Thus an image of the airglow and atmospheric scattered airglow was determined, from which a dependence on the zenith angle has been extracted. For this analysis, nights in which the zodiacal light and Milky Way light are not present on the vault have been selected. The frames taken during such nights have been stacked and a mean sky value was inferred for each pixel. Some of the images, taken when the aerosol load was very low, are in fact registering the airglow source zenithal dependence, and, when compared with other nights, can give an indirect hint about their atmospheric scattered airglow and aerosol loading. It should be noted that in the following chapters, the high background signal from non-stellar sky sources, mainly the airglow, is found to be the principal obstacle in the accurate measurement of the atmospheric attenuated star flux. 15 Chapter 3: Direct star irradiance measurement with the Whole Sky Imager Sky monitoring with the WSI, presented in the previous chapter, allows the measurement of the star irradiance as attenuated by atmosphere absorption and scattering, on a background of atmospheric scattered light from extended sources of the night sky. This chapter presents the errors associated with star flux retrieval due to the physical principles on which the measurement is made and to the instrumental calibration. Geometrical calibration and absolute flux calibration of the WSI are addressed. In the last section, the WSI appearance of a constant height of airglow is assessed. Here are some of the questions to which we will provide an answer in this chapter: ? Is the variation of radiance, seen from image to image in the sky background, a real phenomenon, i.e. is it a recording of the airglow variation between the takes, or it is just the view of the dark current due to the thermal surface currents generated in the silicone by the incomplete cooling of the instrument? ? Is the star irradiance acquired on the CCD obtained as if projected from the sky dome at the corresponding zenith angle, or is the attack angle on the lens aperture significantly smaller? ? Is it possible to determine the star magnitude to the decimagnitude accuracy, to the centimagnitude or even to the millimagnitude accuracy? ? Was the flatfielding performed in the laboratory, or was it done using a median of frames viewing the sky? 16 ? Is there a shift in the stars color introduced by the processing method used on the raw frames? For example, if the flatfielding of the CCD was done using sky frames, then this operation introduced a shift in the star color. 3.1 Error in star irradiance. Star signal-to-noise (S/N) and calibrated measurements (CM) accuracy To finally compute the errors in aerosol optical depth retrieval, one must know the quality of the star flux measurement with the WSI. The errors in the star flux retrieval are due to the errors in the frame acquisition and manipulation, and are caused by the dark current subtraction, readnoise, digitization noise, flatfield division and bias subtraction for the CCD frame. These errors will be assessed next. There are instrumental uncertainties introduced in the final result by each of the processing operations performed on the image and by the conditions of instrument operation. There are also calibrations errors of the instrument. We will compute the measurement error in WSI star fluxes by fully propagating the diverse errors in the WSI formula. Because the instrument was designed to measure primarily the cloud fraction, the absolute calibration stability is ?less critical? (Marine Physical Laboratory website), while when making absolute measurements of stars, the absolute calibration stability for at least during 10 hours of night is fundamental. We will discuss each term in the rms error of star flux, which is the inverse of the signal to noise (S/N) (Newberry, 1991; DaCosta, 1992; Massey, 1992): 2/1 2 22 exp )])(1([ PgntnN n n nN N N S freaddarksky sky star starstar star ++++++ = ? , (3.1) 17 where ? N star , in electrons, is the total number of photoelectrons produced in the CCD by the source extended over n star pixels. ? N sky is the number of photoelectrons produced in semiconductor by the background (or ?sky?) diffuse sources, and calculated as the mode of the distribution of n sky pixels surrounding the object. ? n dark (in electrons/pixels/second) is the dark current, composed out of thermally generated electrons within the semi-conductor. ? t exp is the exposure time. ? n read (in rms electrons/pixel/read) is the readout noise, i.e. the number of electrons retrieved in the electric signal, and which are due to the random fluctuations introduced by the on-chip amplifier, the A/D circuit and from the output electronics. ? g (in electrons/ADU) is the gain of the analog-to-digital converter, which transforms voltages (real numbers) into digits (natural numbers). ? (? f ) 2 = 0.289 and is a factor which gives the error introduced by the A/D converter. ? P is the error introduced by the flatfield division and bias subtraction, which are considered very small. The signal measured in radiance units [mW/m 2 /sr/nm] may be converted back to photoelectrons, N star, through the equation: 18 N star = I star t exp / C, (3.2) where C is the calibration constant specific for the combination of neutral and density filters used in measuring I star . The denominator (the noise) is a combination of Poisson and shot noises. It is known that the photons arrival is a process that obeys the Poisson distribution law, such that the standard deviation is given by the square root of the number of events: ? (N star ) = (N star ) 1/2 (3.3) ? (N sky ) = (N sky ) 1/2 (3.4) ? (N Dark t exp ) = (N Dark t exp ) 1/2 (3.5) The last two terms in the denominator (N Read and the last term) are shot noises, or white noises, which are pure noises, and will not appear as square roots in the equation. The standard deviation of the measurement of the star: ? = 1 / (S/ N) star , (3.6) is further used in the error analysis as the root mean square (rms) error of one radiance measurement. When the errors due to the CCD detection process are much smaller than the rms error of the source?s photons fluctuation, then equation (3.1) yields the maximum signal to noise ratio: (S / N) max = (N star ) 1/2 , (3.7) which is also true for a very bright source. 19 This reevaluation was prompted by the fact that the n dark (the dark counts per pixel) appears to be alarmingly high for the data analyzed lately (the year 1998), and thus the value of n dark = 0.35e-/px, indicated to me by the instrument constructor, is believed to be too low to be true. No value for the n dark actually used in the reduction is archived. 3.1.1 Temperature of the CCD, dark current error, flatfield, bias subtraction error, read noise It is important to understand that too high a dark current might render the signal to noise ratio for a star insufficient to perform photometry. We will answer here the following question: Is the flickering, seen from image to image in the sky background a real phenomenon, i.e. is it a recording of the airglow variation between the takes, or it is just the view of the dark current due to the thermal surface currents generated in the silicone by the incomplete cooling of the instrument? The dark current noise is a part of the error introduced in the calibrated image by the electronic noise. This noise expresses the uncertainty of photon arrival, which is a process whose parent population is Poisson distributed, such that the rms is given by the square root of the expected value of the sample. The electronic noise is a sum in quadrature of the rms readnoise (n read ) and dark current noise whose variance is given by the dark count (n dark [ADU/s], the average dark count rate times the exposure time t exp ) : 2/1 exp 2 )( tnnNoise darkreadelectronic += , (3.8) The dark current for this CCD (TH7895B, grade 1) is known from the literature to double for each 6 0 - 9 0 C increase in the CCD temperature. The dark 20 current occurs mainly in the surface, at the interface of the Si-SiO 2 layers. Authors differ in assessing the temperature dependency of the dark current, but in general the number of thermally generated electrons or, equivalently, the electrical current generated in the unit silicon area, (n dark ), is given by the following formula, where T [K] is the absolute temperature, k B is the Boltzmann constant, and A, B, ? are parameters: ))/(exp( TkBATn Bdark ?= ? [e- /px/s or pA/cm 2 ], (3.9) The parameter A is proportional to the dark current generated at room temperature, B is proportional to the band gap between valence and conduction level for Si lattice and ? is a constant, usually between 1.5 and 3. In order to avoid the necessity of cooling the CCD camera with liquid N 2 , technology permits operation in multipinned phase (MPP) mode. This makes it possible for a ?120 0 C normal operating temperature requirement to be increased to, say, ?40 0 C with respect to the ambient temperature, while keeping the dark current and the noise associated with it at very low levels. As quoted in the technical file of the CCD TH7895M-H, there is a 30 e-/px/s dark current generated in the standard mode operating at ?40 0 C, while only a 0.5 e-/px/s dark current in the MPP mode at the same temperature. We do not know the parameters present in the exact formula for the dark current in the WSI CCD, and we have only the following information: the actual temperatures at which the CCD was operated in the 1998-2003 interval is recorded within the image file. We also have the information that the MPP mode was not used for the WSI CCD, (Shields, Dec.2000, personal communication). 21 We will assess then the whole range of dark current dependence of temperature in an empirical way, from the few recorded temperature points embedded in actual images files and the quoted date for the nights when the dark frames acquisition took place. The CCD temperatures during the nights of interest for this study have been represented in the Fig. 3.1, for each year between 1998 and 2003. It can be seen that in 1998-2000 the CCD was operated above ?28 0 C, which signifies that the dark current and thus the electronic noise present in the star measurements might be important. The measurements in 2001-2003 are comparatively better, with an operating temperature mostly under ?35 0 C. This is the reason why we regard data acquired during the many clear sky nights of 1998-2000 as noisy because of the dark current. We will try to assess how well the dark current subtraction was done and if the images fall into the category of ?bad-unusable data? or just ?ugly-real data? (Massey and Jacoby, 1992). On 1998 calibrated images, for example, the dark current degradation of the image can be seen as high rate, random bright spots on the part of the frame that is not acquiring a direct sky image, and thus are unrelated to either any sky optical source or cosmic rays (cosmic rays would not have such high rate). Lower values of dark current and consequently of dark noise are obtained beginning with the year 2000. Many a clear night were observed with a properly cooled CCD device in that period. The CCD 7895, in MPP mode, is quoted as having Fig. 3.1 (next page) CCD temperature variations during each run, during the nights of interest, in order, for the years 1998(upper left) to 2003(lower right); each graph is clearly labeled with the year for which the CCD operation temperature was represented. The dark current in each pixel is a strong function of CCD operation temperature. 22 23 a dark current of 0.5e-/px/s at 40 0 C, ?one of the highest performance CCDs as of 1994? (Martinez and Klotz, 1998). No MPP is used for WSI images. We know that a good dark frame is obtained as follows: ? a stack of dark exposures taken at the same exposure time as the raw image (60s) and with the shutter closed; ? in each pixel, a median value per stack is obtained. No real time dark image (which is taken one every 6 minutes), is provided to the user. This image is not archived. Instead, a ?generic dark? image per calibration file is offered for the qualitative assessment. We have represented the 15 generic dark current data for this CCD usage, one per calibration file, during 1998-2003, in Figs. 3.2 - 3.4. The generic dark is given as a combination linear in the exposure time (60s). Across the CCD, in each pixel there are a random number of thermally generated charges. Then we will expect that the ?generic dark? frame, made of dark medians in each pixel, will be well represented by a Gaussian distribution of the dark counts per pixel. We performed the Gaussian fit and the parameters for the generic dark current are given in the Table 3.1. For the units conversion we have used a gain or 4.37e-/ADU. It is most certain that the bias (997 ADU) was not subtracted prior from the dark data, such that we would better name these data the ?dark signal?. 24 Generic Dark period (yyyymmdd) Peak[number events] Mean [ADU/px], Mean[e-/px] Spatial ? [ADU/px], ? [e-/px] Interval limits for the histogram [ADU/px] After 20011128 522.0 1015.18; 4436.34 1.56; 6.82 994; 1034 20011119-20011128 428.90 1022.32 4467.54 2.11; 9.22 1002; 1042 20011115-20011119 340.43 1025.05 4479.47 2.40; 10.49 1004; 1044 20010907-20011115 and 20010828-20010907 and 20010821-20010828 411.74 1027.40 4489.74 2.52; 11.01 1007.5; 1047.5 20010703-20010821 221.11 1057.43 4620.97 4.69; 20.50 1037.5; 1077.5 20000506-20010703 486.97 1019.59 4455.61 2.13; 9.31 999.6; 1039.6 19990302-20000506 342.06 1411.85 6169.78 30.31; 132.45 1262.4; 1562.3 19990217-19990302 465.33 1290.62 5640.01 22.29; 97.41 1141.0; 1440.9 19990203-19990217 401.54 1342.54 5866.90 25.81; 112.79 1193; 1493 19981020-19990203 390.84 1360.16 5943.90 26.53; 115.94 1210.6; 1510.5 19981007-19981020 392.82 1355.27 5922.53 26.40; 115.37 1206; 1506 19971010-19981007 410.64 1339.10 5851.87 25.26; 110.39 1189; 1489 19970404-19971010 370.98 1379.03 6026.36 27.98; 122.27 1229; 1529 Table 3.1 Parameters for the Gaussian fit to the dark signal for the generic darks provided in the calibrations frames for the years 1998-2003. 25 Fig. 3.2 Histograms of the distribution of the dark signal, in ADU, for ?generic dark? median files provided with each calibration file, for the period beginning with the dates indicated. Each graph is clearly labeled with the date beginning with which the generic dark applies: Upper left: images beginning with 20011128; lower right: 26 images beginning with 20010821 (to 20010827). See Table 3.1 for values of the Gaussian fit. Fig. 3.3 (next page) Histograms of the distribution of the dark signal, in ADU, for ?generic dark? median files provided with each calibration file, for the period beginning with the dates indicated. Each graph is clearly labeled with the date beginning with which the generic dark applies. Upper left: images beginning with 20010703; lower right: images beginning with 19981020 (to 19990202). See Table 3.1 for values of the Gaussian fit. 27 28 Fig. 3.4: Histograms of the distribution of the dark signal, in ADU, for ?generic dark? median files provided with each calibration file, for the period beginning with the dates indicated. Each graph is clearly labeled with the date beginning with which the generic dark applies. Upper left: images beginning with 19981007; upper right: images beginning with 19971010; lower right: images beginning with 19970404. See Table 3.1 for values of the Gaussian fit. It is important to note the change of the generic dark spatial variance of the signal in time: from values of 22-30 ADU before May 2000, to smaller values of 2.5 ADU after Aug. 2001 and to 1.56 ADU after Nov. 2001. This seems to be correlated with the better cooling of the device when each of the later series of closed-shutter 29 images was taken. It is also possible that a change of CCD occurred. The dark current is more uniform spatially across the CCD when the chip is properly cooled. We do not know the variance of each pixel in ?generic dark? with respect to the same pixel?s median dark current in the stack of dark frames. We know only the inter-pixel spatial variance across the CCD. So we can say nothing about the variance and error per pixel in the individual dark frames used initially. It can be assumed that in the period 1998-2000 a correct dark subtraction was made from the raw images, with darks taken at same high temperatures as the raw frames; nevertheless, no information was given about the dark noise introduced by the dark subtraction. We correlate then what other indications we have about the dark noise. The only indication we have about the value of n dark is that the Noise instrumental =1% Background Sky Signal, during the night (Shields et al, 1997). As stated above, this is the sum, in quadrature, of the rms readout = 6.9 e-/read and (n dark t exp ) [e-/px per 60s]. Consequently, we proceeded to reconstruct the dark current from the sky signal, to attach a realistic range to the unknown dark current in the critical year 1998. As can be seen from this reconstruction made in the Table 3.2, in the last 2 columns of the table the dark current for 60s exposure, expressed in [ADU / ps / 60s] should be in the same range, if the 1% sky condition and the generic dark condition used, respectively, in the two last columns were to be close. It is obvious that the dark current is much higher at ?27C than it was admitted by the 1% sky condition. 30 Image from T CCD (C) Sky [ADU/60s] and Sky [e-/px/60s] Maximum (n dark ) [ADU/px/s] and [e-/ps/s] Maximum (n dark t exp ) [e/ps/60s] from Sky Maximum (n dark t exp ) [ADU/px/60s] from Sky Median (n dark t exp ) [ADU/px/60s] from ?generic dark frame? histogram 2001 1209 -41 460; 2007 1.4; 5.9 355 81 78.7 2001 0718 -33 420; 1835 1.1; 4.8 289 66 264.1 1999 1210 -27 640; 2797 2.8; 12.3 735 168 1809.2 Table 3.2 Empirical assessment of the dark current for the dates in the first column, from the condition that the electronic noise is 1% of the sky signal. This assessment shows the dark current range, and it is not used in the actual S/N computation. To solve this miscorrelation, without a proper knowledge of the actual dark frames used, we would proceed in the following way, using all the experimental data provided to an external user of the data, like this author: 1. Acknowledge the date at which the generic dark file was averaged from multiple dark files and find the corresponding CCD temperature during that day (days). 2. Construct the dark current dependence on the temperature for each or the two cameras used during 1998-2003 period. 3. Interpolate the dark current values at the other temperatures, which are registered with each night of interest. 4. Account for a camera change sometimes between 20000415 and 20000505. Proceeding as stated above, it is found that the part of the dark current, which is dependent on the exposure time, changes in median value as well as in the spatial variance across the chip. We found similarly that the change in the part of the dark current independent of exposure time (which is called ?bias?) is much smaller, 31 varying little around about 1000 ADU. The results of fitting a Gaussian distribution to the darks are presented in the Tables 3.3 and 3.4. Consequently, the corresponding temperature dependence of dark current can be assessed for any other operating temperature and a current value was included in the S/N computation. Calibration File (yyyymmdd) Generic Dark recording date, Unit 10 T CCD ( C ) Bias; ? Bias [ADU] n dark t exp ;? dark [ADU] After 20011128 74 darks from 20011203 and 20011204 -39?-40 999.56; 0.12 15.62; 1.60 20011119-20011128 59 darks from 20011124 and 20011125 -36?-38 1000.23; 0.13 22.09; 2.16 20011115-20011119 41 darks from 20011117 and 20011118 -36?-37 998.64; 0.17 26.41; 2.46 20010907-20011115 and 20010828-20010907 and 20010821-20010828 55 darks from 20010911 and 20010912 -37 1000.62; 0.15 26.77; 2.56 20010703-20010821 62 darks from 20010710 and 20010712 -30?-31 1003.51; 0.18 53.90; 4.76 20000506-20010703 37 darks from 20000603 -43(?) 997.72; 0.17 21.86; 2.18 Table 3.3 The dark current profile as function of the CCD temperature and the frame bias for WSI camera known as Unit 10. Calibration File (yyyymmdd) Generic Dark recording date, Unit 4 T CCD ( C ) Bias; ? Bias [ADU] n dark t exp ;? dark [ADU] 19990302-20000506 34 darks from 19990703 -27 1017.43; 1.02 394.34; 30.57 19990217-19990302 34 darks from 19990225 -28(?) 1005.10; 0.88 285.46; 22.54 19990203-19990217 32 darks from 19990215 -28 1007.24; 0.87 335.24; 26.08 19981020-19990203 18 darks from 19981030 -28 1012.76; 1.01 347.34; 26.86 19981007-19981020 56 darks from 19981010 and 19981012 -28 1006.77; 0.74 348.36; 26.64 19971010-19981007 37 darks from 19980701 -28 1007.14; 0.64 331.78; 25.48 19970404-19971010 47 darks from 19970915 and 19970916 -28 1009.95; 0.70 368.94; 28.21 Table 3.4 The dark current profile as function of the CCD temperature and the frame bias for WSI camera known as Unit 4. 32 The bias, or the signal obtained at zero exposure time, is a pedestal in the number of free conductors in the CCD, and is computed by averaging at least 100 bias frames, (which are taken with a zero exposure time). The bias has usually a highly stable value, ?at least during one night? (Gullixon, 1992). For this CCD, the bias indicated by the constructor is 997e-. As can be seen from the Tables 3.3 and 3.4, the value obtained from our computation ranges, for the two different cameras, from 997e- to 1017e-. The flatfield is the pixel-to-pixel quantum efficiency variations, which occur from non-uniformities in the CCD production and in the design of the frontal surface (for example, the metal gates on each pixel?s surface), which is illuminated when this CCD is in use. The error introduced by the flat fielding is considered very small. The error given by the pixel-to-pixel variation of the quantum efficiency is reflected by a variation in the gain across the CCD. A histogram of elements of the flatfield matrix (not shown) provided to the users in the calibration file and which is already normalized to unity, shows almost no dispersion, such that the gain has little uncertainty, and the CCD is considered good for absolute photometry (Gilliland, 1992). The read noise, which is a shot noise that occurs when the electrons acquired in each quantum well are shifted along the pixel columns towards the CCD margin to be read (i.e., counted), has been taken as 6.05 rms e- for unit 4 and 6.9 rms e- for unit 10, while the gain is, respectively, 4.9e-/ADU and 4.4e-/ADU. 33 3.1.2 Systematic errors due to absolute flux calibration and geometric calibration in laboratory 3.1.2.1 Absolute flux calibration It is essential to understand that the calibration of the frames was not made, as in astronomy, by using twilight sky images or images of the inside of the telescope dome. The absolute flux calibration was made in the laboratory, by imaging a uniform radiance source through the WSI camera, obtained by illuminating a Lambertian surface with a lamp of standard spectral irradiance. The irradiance of the image is known to be proportional to the radiance of the object (Smith, 1990). The absolute calibration of a frame acquired with a fish-eye lens has special features (Voss et al, 1989). When the instrument is uniformly illuminated by a constant density source, the lens produces a decrease of image irradiance with the zenith angle, which must be corrected with a roll-off function (?RollOff?), with the shape ? for the fisheye lens - of an inverted cosine function. Consequently, in WSI images, for example, the increase in atmospheric scattered airglow with zenith angle is already corrected for any diminution it might have had from the lens aperture- decreasing surface. What we see in the calibrated images is indeed due to a phenomenon in the object space and is not an artifact of the lens. The absolute flux calibration operation produces a calibration constant, which is specific to the filters (chromatic and/or neutral density), which are currently used. During the calibration, the measurement, M [ADU], of the radiance reflected by the plaque, is performed broadband in the central pixels of the instrument. The roll-off function correction, RollOff, is then applied to obtain the digital counts in each pixel 34 of the image, besides the central ones. The calibration constant is then computed as a relation that transforms the analog to digital units in the image into physical power units: A[Wms/m2/sr/nm/ADU]=? /? ?E(? )d? [W/m2/sr]t exp [ms]/{M[ADU] RollOff} / ?R(? )d? , (3.10) where ? E is the irradiance of the standard spectral lamp used to illuminate the Lambertian plaque, ? ? is the plaque reflectance, ? R is the spectral response function of the combined optics spectral transmission and CCD?s quantum efficiency, ? t exp is the exposure time for the calibration process, ? M [ADU] is the radiance falling on the camera in the central 3x3 pixels, expressed in digital counts, ? RollOff is the selective attenuation with the zenith angle in the object space, explained above. After the calibration constant is known thus from laboratory data, the newly acquired frames, expressed in photons per pixel or in ADU (digital counts), are flatfielded, bias subtracted and dark current subtracted, as explained in the previous section, and then can be expressed into radiative (power) units by a simple multiplication with the absolute flux calibration constant. A calibrated measurement (CM) in a pixel is then related to the radiance per pixel of the sources in the object space: 35 C.M.[W/m2/sr/nm/pixel]=A* Image Counts[ADU/pixel] / t [msec], (3.11) where t is the exposure time for the image acquisition, which registers ADU counts per pixel, and A is the calibration constant. In essence, the calibration transforms the image irradiance into the object radiance at the top of the camera. No astronomical source is involved in this process, contrary to the practice in astronomy. The flatfielding in practical astronomy is described by Chromey and Hasselbacher(1996). When using twilight sky images, the effect of small sky gradients on star flux extraction is important. 3.1.2.2 Geometric calibration The error in the projected area of a pixel on the sky is a matter of the accuracy with which the geometric calibration has been made. The geometric functions, relating a pixel position on the CCD with the area on the sky corresponding to it, are served to the user in the calibration file, and quoted with an accuracy of 0.25 pixels or 0.1 deg. This accuracy has been attained not by imaging markers in the laboratory (Kumler, 2000), but by a procedure involving star field imaging (MPL website). The geometric functions are a fit between the Cartesian position of a pixel in the frame and the zenith and azimuth of the corresponding direction in the object space. The set of 4 geometric functions, determined by star position technique for several consecutive days when the sky is very clear, are used for many subsequent nights, until the small vibrations in the instrument moving parts (the trolley for Sun and Moon occultation, for example) renders the calibration outside the limits of geometrical precision and a new calibration is required. There are 15 such calibrations performed during the period 1997-2003. 36 The solid angle in object space subtended by a pixel, ? , can be inferred to be given by: ? = d? [sr] / dA[px], where dA could also be expressed in physical surface, (19?m) 2 . Simple spherical projection geometry yields: ? = (sin ? /r) (d? /dr), where r [pixels] is the radial distance on the CCD from the optical axis and ? the zenith angle (or object angle with respect to the optical axis). Each pixel receives radiation from a solid angle close to 34 micro-steradians. The solid angle viewed by one pixel is maxim at 30 degrees zenithal distance, and drops by 3% at 85 degrees and with 8% at 90 degrees. The actual values of the maximum are given in the Table 3.5. The angular resolution variation with the zenithal angle in the object space, for all of the calibrations performed by instrument maintenance, is presented in Figs. 3.5 - 3.7. Unit 10: max ? Unit 4: max ? After 20011128 34.46 19990302-20000506 34.22 20011119-20011128 34.44 19990217-19990302 34.21 20011115-20011119 34.65 19990203-19990217 34.25 20010907-20011115 and 20010828-20010907 and 20010821-20010828 35.53 19981020-19990203 34.18 20010703-20010821 36.76 19981007-19981020 34.30 20000506-20010703 36.10 19971010-19981007 34.25 = = 19970404-19971010 34.25 Table 3.5 Maximum value of the angular resolution per pixel [?sr/px], for each calibration and unit. The angular resolution is found to vary with the changing of the lens performed in the spring of 2000. The angular resolution function is an inverted parabolic function. The precision of ? measurement is thus regulated by the precision of the angular resolution, d? /dr, and of the angle ? measurement. As theta here is taken from one of the geometric functions in the calibration file, computed from a polynomial fit, its precision reflects the quality of the accurate determination of the star center position 37 from WSI images and ultimately the low separation power (angular resolution) and distortions present in the optical system (particularly the lens). Fig. 3.5 Angular resolution variation with the zenithal angle in the object space, for 4 of the calibrations performed by instrument maintenance during the interval studied. Each graph is clearly labeled with the beginning date for which the calibration applies. Each pixel receives radiation from a solid angle close to 34 ?sr. Fig. 3.6 (next page) Angular resolution variation with the zenithal angle in the object space, for 5 of the calibrations performed by instrument maintenance during the interval studied. Each graph is clearly labeled with the beginning date for which the calibration applies. Each pixel receives radiation from a solid angle close to 34 ?sr. 38 39 Fig. 3.7 Angular resolution variation with the zenithal angle in the object space, for 6 of the calibrations performed by instrument maintenance during the interval studied. Each graph is clearly labeled with the beginning date for which the calibration applies. Each pixel receives radiation from a solid angle close to 34 ?sr. 40 Calibration functions for the night of 20011120 have been recalculated, as those provided by the calibration file were wrong. An algorithm was designed to compute the pixel corresponding to the zenith direction, and then the rotation angle between the true North direction and the margins of the CCD plate. Usually the margins are aligned to directions N-S and W-E. 2792 stars positions (of not all different stars) were detected and centroided from the images taken during the 20 Nov. 2001, and their azimuth and zenith angles were calculated from the right ascension and declination for the Epoch 2000.0. Assuming that the plate was kept horizontal, it has been found that the zenith vertical corresponded to a center at (x R , y R ) = (253.9, 256.6) pixels, and that the angle of true North was 2.454+/-0.001 degrees rotated to the West with respect to the CCD margins: [x R y R ] T = R (-? ) [x y] T , (3.12) where R the rotation matrix of angle (?? ), and (x,y) the pixel coordinates in a N-S aligned CCD (512x512 matrix). For the large number of degrees of freedom, the probability of the fitting functions is as follows: 90% for a chi-square per degree of freedom of 0.9657, 75% for 0.9818 and 50% for 0.998. It is believed that the CCD was rotated in the enclosure for the purpose to be used by the instrument mentors to study the pixel-to-pixel responsivity variations. These were the fits performed to the data regarding the radial projection of the zenith: ? a sinusoide (Kumler et al, 2000): r(pixels)=a sin [b ? (deg.)], a = 363.0804+/-2. 10-5, b=0.008+/- 1.09, ? 2 /(N-2)=0.809, when the error in r determination is taken as 0.75 pixels. The probability of the fit is 99.5%. The high probability indicates rather that the error in r has been over-estimated; we use the values 41 of chi-squared per degree of freedom to recalculate the uncertainty in r as [0.75 2 x 0.809] 1/2 , i.e. 0.67 pixels. ? a cubic polynomial: ? = a0 + a1 r +a2 r 2 +a3 r 3 , a0=-0.095+/-0.05, a1=0.341+/-0.001, a2=(-5.2+/-1.3)10-5, a3=(6.81+/-0.33)10-7, ? 2 /(N- 4)=0.9840, with error in zenith angle of 0.3105 degrees. The probability of the fit is around 70%. Statistically, this is interpreted as follows: in repeating the experiment, one should obtain in 70% of cases a fit that is no better than this. ? a function suggested for another lens, at the MPL website: r(px)=b0+ R 0 (b1 ? + b2 ? 2 ), b0=-0.0193+/-0.0377, b1=2.981+/-0.003, b2=(-5.78+/-0.76) 10-4, ? 2 /(N-3)=9.9. Probability of the fit is 0.01% (poor), when uncertainty in r is taken 0.25 pixels, and probability is 50% when the uncertainty is taken 0.82 pixels. Consequently, the geometrical functions were then assessed using the cubic polynomial fit above, with rms in zenith angle of 0.308 degrees, which are necessary for a fit with a chi-square of 1. The next night, the lens was restored to the N-S aligned position, and the sample could not be increased with other stars. The main error of such geometrical calibration functions is due to the error in the star centroiding for star centre position determination, and it is ultimately traceable to the low resolution of a fast lens (i.e., with small focal length). The relative error in the surface projected on the sky by a pixel is a sum in quadrature of the relative errors in radial resolution, zenith angle and dimension of a pixel. The highest value is that of the radial resolution error, given by the goodness of 42 the fit of the calibration functions, 0.70%; the errors in zenith determination and in radius are 0.1% each. The total relative error in pixel-to-sky projection is the 0.72%. 3.1.3 Star flux measurement error: Signal to noise ratio and calibrated measurement error Star flux is retrieved with an rms error equal to the inverse of the signal-to- noise ratio. For photometric purposes a 1% rms error in flux is necessary (Young et al, 1991). We have retained the star flux measurements that have a S/N above 100. We would like to determine the variance of the calibrated measurement (CM), defined in the first section of this chapter. We have already computed the root mean square error of the star irradiance measurement (F [ADU]). To recall the main results, the error in star irradiance was found to have three main components, which are, in order of importance: Poisson noise from sky and star, and camera error (camera noise and digitization error). The camera noise for unit 4 and unit 10 used in the period 1998-2003 was also computed in subsection 3.1.1. Depending on CCD temperature, an appropriate dark current value was used consistently. The image of sky background is the most error prone component. To limit this photon noise, we consistently used an elliptical aperture for median sky computations. Each ellipse has the axis ratio 5, a thickness between 8 and 12 pixels and the major axis perpendicular to the lens diameter passing through the star?s center. This should compensate for the radial gradients in the sky background, increasing from zenith to the horizon. The final error in star irradiance, as the inverse of the S/N ratio, is, for our measurements, under 0.01, but we have also used sometimes measurements with 43 0.02-0.03. These are stars that do not rise enough in the sky at the latitude of the observer, consequently are seen at low elevations, where the sky background is bright from van Rhijn enhancement effect (the airglow emitting layer is thicker towards low elevations), and from atmospheric scattered light of the extended sources (airglow, integrated starlight and zodiacal light). The airglow (AG) has a variable component (Leinert et al, 1995), which has a short period detectable in the WSI sky background. When a star with the visual magnitude 2, with typical value of irradiance (1000*10-7 mW m-2 sr-1 nm-1), or 31*10 3 photons (or, equivalently, electrons), in 60 seconds, and a S/N=100 is considered, then the formula for the propagation of errors in the calibrated measurement is, ? (C.M.)/(C.M.) = [? 2 (I (ADU))/I 2 (ADU) + ? 2 (A)/ A 2 ] 1/2 . Recall that I(ADU) is the star irradiance converted back to counts and A is the calibration constant. The rms of star irradiance, in ADU, is less than 0.03 for a S/N of 33, such that the relative error in star irradiance, ? (I (ADU))/I (ADU), is less than 10 -6 . For the gain, we know the value for both unit 4 and unit 10, and we inspected the variance of two flat images and two bias images provided in the calibration files. The flat matrix ? already normalized to unity, does not exhibit a Gaussian distribution, that is, the variance is extremely low. We concluded that the variance in the gain constant may be taken zero. From the above equation it results that the ?calibrated measurement? precision is limited not by electronics, or by the successive transformation applied to the raw image, but by the calibration constant, whose relative error is ?better than 0.5%? (MPL and WSI websites): ? (CM) / CM < 0.5%. 44 3.2 Sky background measurement error: Constant height layer?s contribution to the background sky and extended sources 3.2.1 Constant height layer range from zenith The polar coordinates along a spherical atmospheric layer, situated at a constant height, and having a negligible thickness as compared with the radius of Earth can be constructed with respect to the zenith point of an observer on the ground surface. The polar radius is measured from the point in the layer where the zenith direction for the ground observer intersects the layer surface. This radius measured along the spherical layer is called ?range?. The range variation on the layer depends only on the zenith angle of the current observed point, on the radius of the Earth (R) and on the height (H) of the layer, in the geometry presented in the Fig. 3.8. 45 Fig. 3.8 The geometry of a ground observation of a constant height (H) layer. The range is measured along the layer; the Earth?s mean radius is R. To obtain the range dependence on the zenith angle, range (? ), we have inverted the relation of ? (range, H, R) given in Taylor and Garcia (1997), to obtain: range =(R+H) arccos { R/(R+H) sin 2 ? + (1- R 2 /(R+H) 2 sin 2 ? ) 1/2 cos ? }. (3.13) Applying this formula to diverse layers, we obtain the following useful assessment of range visible with the WSI, which is given in the Table 3.6. Height of layer [km] Observation zenith angle ? [deg.] Range (? ) [km] Height of layer [km] Observation zenith angle ? [deg.] Range (? ) [km] 5 (tropospheric aerosol limit) 50 8.4 96 (O I airglow emission) 50 113.2 87.1 88.0 87.1 835. 90 252.5 90 1115. Table 3.6 The range (km), along two spherical layers of interest visible with the instrument. 46 As can be seen, when looking at the horizon, one is observing a range of 253 km on a 5 km height layer, but a range of over 1100 km on a 96 km height layer. 3.2.2 Background sky correction necessity For point sources there is a fraction of the light which can remain outside the signal-maximizing aperture; this fraction ? if one is interested in the star absolute flux measurement - is retrievable by aperture growth techniques (Stetson 1990, Howell 1990), and represents the part of the direct signal which is scattered outside of the solid angle constructed by the synthetic aperture around the line-of-sight. Usually, for non-absolute (i.e., relative) photometry of the stars, it is permitted to measure the partial flux present in an aperture or in a synthetic point-spread-function, provided that the fraction that is left unmeasured is constant for all the stars across the field of view of the instrument. The absolute photometric measurement of a point source is thus made of a direct part (flux) and a halo/aureole part, scattered by imperfections of the optics, dust in the interior/ exterior of the optical system and by aerosol scatterers. Molecular scattering is isotropic (phase function roughly does not vary with the scattering angle; i.e., there is some dependence in Rayleigh scattering, but it is not large by comparison to aerosols), while the aerosols scatter anisotropically and the phase function, on average, is 30 times higher for forward scattering than that for 40 degrees scattering angle, and 200 times higher than that for scattering in right angle or scattering backwards (Thomas and Stamnes, 1999). Because the halo appears in the immediate neighborhood of the star, and it is natural to consider it the result of the anisotropic scattering by the aerosols. Do we wish to include this part in the direct radiation received for a star? No, because it is attenuated (missing) when the 47 scattering due to the aerosols enters in the measurement equation as aerosol extinction optical depth. When one wishes to measure extended sources (like the airglow), one should correct for the signal scattered out of any observational aperture, that is still counted, because it is replaced by flux originating in other parts of the extended source, which is scattered into the solid angle of the observation. We should recall these considerations in the Chapter 4, where the accuracy in the determination of the sky background around the star in WSI observations is assessed. 3.3 Summary of Chapter 3 In this chapter we have described the flux calibration for point source measurements with the WSI, the geometrical calibration of the WSI and the background sky noise contribution in the point source measurement. The operations performed on the raw images were briefly explained and their respective contributions were computed pointedly for this instrument?s operation history, during the time frame of interest for this work: ? dark current (or thermal noise) subtraction; ? CCD reading of the CCD and the shot noise introduced by it; ? flatfield division to compensate for the pixel-to-pixel quantum efficiency variations; ? bias zero exposure subtraction of the electrons which are not generated by an image. Taking into account such particularities of the data, it was found that the absolute flux calibration of a star measured with the WSI is, in general, less governed by the random and shot noise of the instrument electronics, - although this produces a 48 contribution that was computed for each operating temperature-, but by the absolute flux calibration constant, which is determined from laboratory calibrations which are done periodically. Also, a geometrical calibration was performed for a case when the original calibration provided by the instrument?s mentors was accidentally wrong. The geometrical calibration was performed using the radiance images themselves, through a star field technique; the results were fitted with functions that are recommended in the literature for this type of lens. The solid angle seen by a pixel was also found to slightly vary with the zenith angle, or, equivalently, with the distance for the center of the chip. Finally, the contribution of the background sky radiance from the extended sources to the WSI image was examined in terms of non-uniformity across the observable zenith angle range. Depending on the height of the layer producing the atmospheric emission or scattering, the layer?s specific appearance on the images was examined in terms of the positional relevance. The background pixel-to-pixel non- uniformity introduced by this lens-specific footprint variation with the zenith angle is the main cause for the necessity of a background sky gradient correction, which will be further explained in Chapter 4. 49 Chapter 4: Star spectra, atmospheric absorption and scattering background, extinction linearization in wide band astronomy, the method of AOD determination from the WSI heterochromatic measurements. The extinction computation methods are as old and reliable as the ?art of photometry? (Howell, 1992) itself. In this chapter, such methods for solving the extinctions equations for stellar sources will be reviewed. The important set of practices from astronomical photometry is examined for the suitability to be used with the whole sky imager stellar extinction data, and found lacking. The principal constraints are to be found in the wide passband of the unfiltered silicone diode, and the non-photometric behavior of the extinction at a low altitude site. A different method of solving for extinction is proposed, and the assumptions under which it can be used are emphasized. At the end of the chapter, additional quantities useful in the problem, like the stellar spectra library for AOD retrievals and sky background scattered light correction are discussed. 4.1 Wide spectral band measurement linearization in astronomical practice Astronomical photometrists have perfected an observational system that affords the comparison of any two star flux measurements, whatever the ?weather conditions? under they were made. (Through the ?weather condition? at an astronomical site is it understood at most the water vapor haze, the occasional stratospheric dust and very rarely a pollutant aerosol.) The systems of filters 50 (Johnson, Str?mgren, Geneva, Washington) through which stars are observed are ?wide band?(width < 70 nm) or ?narrow band? filters (width < 30 nm) that are sampling the ultraviolet-visible-infrared spectrum. Catalogs of ?standard? (non- variable) stars for each such system are used as calibration comparisons to calibrate other stars through differential photometry within the same field-of-view. Each system of filters should have suitable catalogs, with stars of various colors and spectral type, and distributed in space in a way that at least a few standards will be found near any direction in which the telescope points. The current practice is that, after the current (?program?) star and the standard stars have been measured under similar extinction conditions (due to small angular distance between them on the sky), the ?filter system color transformation? equation is applied to the program star, to obtain its flux (or, in logarithmic scale, its magnitude). Fig 4.1 Johnson?s system of broadband astronomical filters (from left to right) UBVRI normalized transmission and the response curve of the Whole Sky Imager. The wavelength is given in ?. 51 An example of the color transformation equation is the transformation of visual magnitude obtained from a standard star measurement, which serves to calibrate the instrumental system: V std = v inst + a 0 + a 1 (B-V) + k ext X+ a 2 time+ a 3 (time) 2 +?, (4.1) where the V?s are the instrumental (just observed) and standard (from catalog) magnitudes, (B-V) is the color excess of the star in blue filter as compared to the visual filter, k ext is a (monochromatic) coefficient of extinction for the V filter, X is the airmass, and the ?a? s? are constants to be determined from the fit. The precision of the measurement is given by the residual between the left-hand side and right-hand side. The equation, with all constants determined through fitting (including the extinction coefficient), is to be applied then to finding the visual magnitude (or visual flux) of another star, as already mentioned. From this brief example it is obvious that the 70 nm wideband filters have been treated as liniarized around their mean or effective wavelength, and the magnitudes are treated as quasi-monochromatic. The historical systems of filters have been studied for orthogonality and found they do not span a Hilbert space (Young, 1991; Young, 1992a), i.e., they not form a system of orthogonal, linear independent vectors. Their cvasi-monochromatic treatment of magnitude and extinction is based on the Taylor expansion of the heterochromatic equation (Golay, 1974; Young, 1992b): ?? + ? ? +?= ? ?? ??? ????? ? ? Rd XF dXFd XFFRd TOA TOA TOA ]... |)))(exp(( |/)))(exp(( 5.01)[)(exp()( 0 0 22 2 00 , (4.2) where 52 ? the left-hand side is the heterochromatic measurement at the ground; ? F is the spectral star flux at the ground, F TOA is the flux at the top of the atmosphere; ? R is the filter spectral response function; ? ? is the second order moment of the response curve of the filter; ? ? 0 is the mean wavelength of the filter. ? An essential part in handling the extinction through linearization is accomplished by the hypothesis that the second term in the square parenthesis is much smaller than 1, and then taking ln (1+x) ~ x. This is equivalent to assuming that the extinction, spectral response curves and the spectrum are known to the second order derivatives and beyond, or the higher orders are dropped and a linearization is obtained. Young (1992b) has demonstrated that, even for heterochromatic filters with a passband of only 70 nm, there is a necessity to take into account higher order terms to achieve a high precision in photometry. Young also takes into account a supplementary filter, the atmosphere, and shows that the convolution with it brings new terms into the Taylor expansion, which cannot be neglected. The conclusion is that the linearization around the mean wavelength, for a heterochromatic observation, does not yield an equation of the type (4.1), from which both the extinction and the ?a? constants of the instrumental system to be precisely determined. Because the current systems of filters are a way to universally compare the star observations, they are not likely to be abandoned, in spite of the deficiencies. 53 After rearranging the terms in (4.2), and applying the assumptions mentioned before, the equation obtained is: m 0 (z 0 ) = m 0 (TOA) ? A X + B X 2 , (4..3) where m 0 (z 0 ) is the heterochromatic magnitude at the ground at the mean wavelength of the filter, m 0 (TOA) is the same quantity at the top of the atmosphere, X is the relative airmass, and A and B are two factors given by: A= ? {1+? 2 (? ?/? ) (F?/ F) +0.5 ? 2 (? ??/? ) }, (4.4) B= 0.5 ? 2 ? ?, (4.5) where the optical depth ? and its first and second derivative (? ? and ? ??, respectively), as well as first derivative of the star spectrum, (F?), are computed at the mean wavelength of the filter. This set of 3 equations is used by Rufener (1986) to obtain the ground extinction value from observing a pair of stars through the same airmass, one ascending and one descending (?M and D? method). It is a remarkable way in which this method handles heterochromaticity by linearization, by assuming that the optical depth is mainly due to molecular scattering, - which is a power law of wavelength (? - 4 ), with well-behaved derivatives. This is true enough at the observatory altitude of 2500m. The strategy for extinction measurements at the astronomical observatories is to observe at some intervals during the entire length of the night, beside the program stars, a set of chosen standard stars within the field of view of the telescope; compiled sets of standard stars with magnitudes fainter than 15-20 are available in advance, for example, the Landolt set; the set should be well represented in each field of view and 54 it must be evenly sampled in the filter system in use with respect to the color and luminosity. As the WSI does not perform deep field observations, the set of standard stars it needs is rather within the bright stars set (brighter than 6). A list of such standards, of stars for which the current knowledge is that they have a constant energetic output in the band of interest, has been compiled from the literature. Excluded from such a list are variable stars, either intrinsic variables or eclipsing variables, like, for example, Algol (Beta Persei), which has a period around 12 hours. (Its variability is nevertheless very well recorded by the WSI). As stated before, the extinction is spectrally, temporally and spatially dependent. ? By narrowing the field of view, the spatial variability of the extinction can be excluded. ? By observing in a multicolor system of filters, the spectral dependence can be accounted monochromatically and Langley plots can be drawn. ? By making only differential measurements of program stars relative to standard stars situated at close angular distance, the temporal variability of the extinction can be surpassed. In contrast with such techniques available for astronomical use, the WSI star irradiance measurement cannot be linearized on the premise that the extinction does not suffer a variation in time, or that it has only a linear drift in time. As seen in the previous section, the hypothesis that the extinction is the same at the beginning, middle and the end of the night would allow one to extract the extinction coefficient 55 as the slope of the star irradiance with respect to the airmass. An improvement is the hypothesis that the extinction ?drifts in time? (basically drops during the night) would allow one to add an increment proportional to the time interval, and with negative sign (Reimann and Ossenkopf, 1992a; Reimann et al, 1992b). This may be well enough for the stratospheric aerosols, which provide most of the extinction at astronomical observatories, and whose variability often can be ascertained only from several years (or at least seasons) of data (Sterken and Manfroid, 1992b). Rare injection of volcanic-originating aerosol in the stratosphere produces persistent and slowly decreasing stratospheric aerosol loading that has been studied over decades at the high altitude observatories (Burki et al, 1995). At a lower observational site like the ARM SGP CART, the assumption of a ?slow drift? of extinction during the night can be sustained only for a few nights a year, as the Raman Lidar measurements prove. For such rare nights with ?frozen extinction?, Langley plots technique can be employed over an ensemble of stars to obtain a mean value of the optical depth per night. The relative airmass value is fit for the mean 1972 ISO Standard Atmosphere density, and with a mean air refraction index at 700 nm wavelength (Kasten and Young, 1992). An example is given in the Fig 4.2. The next step is to allow a more complicated variation of extinction in time, and with the wavelength. 56 Fig 4.2 Observed star irradiance variation with respect to the airmass during two nights in Oct. 1999 (Langley plots). Upper: Langley plots for 7 stars during the night of 19991005. Lower: Langley plots for 7 stars during the night of 19991010. 4.2 The method of AOD determination for the WSI heterochromatic measurements. The Method of Solving the WSI Equation 4.2.1. Introduction In an exchange of letters between Maxwell and Rayleigh, the creator of the electromagnetic theory pointed out that, if the properties of the atmosphere were 57 better known, then the light of the stars could be used to determine the atmosphere density. Indeed, in the Rayleigh formula that links the scattering cross section per molecule, as inferred from Maxwell?s electromagnetic theory, with the squared index of refraction of the medium, - in the denominator appears the number of chaotic distributed scattering molecules (Jackson, 1999). Is there any value in revisiting such an old idea? The atmospheric properties are much better known today than in the lifetime of the two great physicists. In this work we revisit such formulas and others, regarding the atmospheric opacity to the extraterrestrial shortwave radiation, with the aim to characterize the opacity variations with a temporal length scale of minutes. Experimentally, when observing the stars, there is a precise time span and spatial ? angular extent for which an instrument can obtain meaningful readings. In our experimental setup, the lower temporal threshold is given by the instrument sensitivity to low intensity for a short exposure times and the upper temporal threshold by the confusion of the star position in the case of long exposures in a non- tracking instrument. Similarly, spatial detectability thresholds are given as following: the lower by the ?seeing?, i.e. the time averaged values of incoherent superposition of images due to the turbulence; the upper threshold by the apparent blending of stars trajectories, already mentioned. 4.2.2. The WSI equation Recall that the heterochromatic star irradiance measurement calibrated measurement (CM) is written as: 58 ? ? +? = ??? ??? ????????? dRpxsr dXttRclassf FnmsrWmI rel scat v abs v TOA )()/( )],()),(),((exp[)(),( ][ 2 5556 112 , (4.6) where ? I [W/m2/ sr /nm] is the calibrated measurement of the star (CM), ? f TOA (? ) is the normalized spectrum for the spectral Morgan-Keenan and luminosity classes classification of the star, ? R (? ) is the spectral responsivity curve of the instrument, ? the monochromatically transmission through the atmosphere along a slant optical path is expressed as the a multiple X of the (vertical) optical depth, ? ; ? ? [sr /px] is the solid angle subtended by the star pixels on the vault; ? F 5556 [W/m2/nm] is the spectral monochromatic flux of the current star at 555.6 nm. The equation describing the measurements is a Fredholm integral of the first kind (FK1). The FK1 equation above would be a linear FK1 equation when the unknown were to be considered the spectral flux of the star, or the product between the star flux and the instrument response function, which are both only wavelength dependent. In this case, the kernel is the atmospheric transmission, which is both time and wavelength dependent. It is indeed possible to assemble the measurements on a straight line, as the variations in atmospheric components are minimal, but would this be a Langley plot, which is by definition about source irradiance at a single wavelength? 59 In the WSI case, both the stellar spectrum and the responsivity function of the instrument are considered known and can be introduced in the kernel, along with the atmospheric transmission due to molecular absorption and scattering by the uniformly mixed components and the transmission known of more variable components, such as ozone and water vapor. The residual part of the atmospheric transmission, due mainly to the aerosols and probably, in small proportion to the absorbing gases whose current value differ from the climatologic mean, is the unknown functional of both the time and wavelength. Inside this unknown residual part might also appear the opacity due to the thin cirrus clouds interfering with the star observation. Rearranging the terms in the Equation (4.6), we will write the calibrated measurement of the star within the following product: 5556 / FRdIy m ? = ?? , (4.7) where y m stands for normalized ?measured? quantity, and is expressed in wavelength units. We included the absolute monochromatic star flux at 5556A in this term, with the purpose of having a compact formula when propagating the errors. The simulated (model) term is then the remaining integral of spectrum, response curve and transmissions: ? +?= ????????? dtXttRclassfy rel scat v abs v TOA )]),(()),(),((exp[)(),( . (4.8) The explicit notation for the non-linear FK1 equation is: ? ?= ??????? dtXtimeclassspectralKy rel )]),(()(exp[),,_,( , (4.9) where K(? ,t,? (t),class) =f TOA (? ) R (? ) T basic (? ,t,? (t)), (4.10) is the kernel, and 60 T var (? ,t,? (t))=exp[-?? (? ,t,? ) X(? ,? )], (4.11) is the transmission due to variable components of extinction during the night. The dependence on the zenith angle, although not an independent variable, was written above explicitly to emphasize that the unknown transmission is not a linear function of time. The problem we must solve is now of the form: m yy = , (4.12) where the equality plays in the least mean square sense of minimizing the difference (error) between both sides of the equation. 4.2.3 The Method of solving the WSI equation It is known that solving a FK1 equation, which is an ill-posed problem, in which a small perturbation in the input data, the right-hand side of (4.12), can lead to large perturbations of the solution (Kythe and Puri, 2002). For example, the FK1-type equation is usually obtained in atmospheric spectroscopic studies, when the monochromatic limb scanning from satellite observes the atmospheric radiation coming from a diverse range of depths and pressures. Knowing that the pressure broadening of diverse emission/absorption lines takes place for a certain pressure range, there is a one-to-one correspondence between monochromatic radiation received and the pressure level at which the active molecule emitted or absorbed it. The vertical profile of diverse active molecules can thus be computed. The FK1 problem was solved through methods like singular value expansion (Miller, 1979), truncated singular value decomposition, generalized singular value 61 decomposition, iterative methods (Siske, 1973; Strand, 1974), Phillips-Twomey regularization (Twomey, 1963) etc. These methods have as their basis a minimization of errors simultaneously with the smoothing (regularization or filtering) of solution oscillations, the so-called ?L-curve? (Hansen, 1992; Schimpf and Schreier, 1997). First, it is noteworthy that, were it not for the heterochromatic WSI measurement, the problem would be simpler. If one were to acquire a set of monochromatic measurements {y m (? )} at different times during the night, with, say, a spectrometer of negligible bandwidth, then the problem would reduce to a set of equations of the form: y m (? ) = exp (-? ? v (? ) X rel (? )), (4.13) where the sum is done after each component (gas, aerosol) involved in the light attenuation. Even assuming the atmosphere as stationary, the set of equation is non- linear in the parameters ? v , the optical depth (OD) of each atmospheric component, such that one of the non-linear minimization methods should to be used, for example, a Levenberg-Marquardt algorithm (Oikarinen et al., 1996). After OD computation, by knowing the extinction cross sections of the each of the separate channels of interaction, one can find the columnar amount for each component, i.e. the path integral of component?s density, provided that the measurements are made at wavelengths where the spectrometer is sensitive enough for the respective component. The present WSI measurements, which are not monochromatic, cannot be treated in the same way. In fact, a Langley line, while still possible to be drawn through heterochromatic measurements, would not have the usual significance. 62 Indeed, the logarithm ln (y m ), as a logarithm from an integral, would not be linear in the relative airmass X (? ), thus not easily separable into the classical form, ln (y 0 )-? X, - where we have dropped both subscripts for vertical optical path and relative airmass, respectively. Heterochromatic measurements that strictly obey a Langley plot are thus either prone to high errors, or refer to an average over the passband of a constant-in- time transmission through the atmosphere. Put in other words, the slope of a heterochromatic star irradiance vs. airmass, expresses in fact a double average, in time and over combined filter-source spectral curve, of the total extinction cross section; this averaged optical depth was determined in a previous work, as a prelude to this dissertation. The method that is used in this work to solve for the small variations in optical depth due to the aerosols and ozone, which are thought to be the cause of the short- time transparency variations during the night, is an iterative method first proposed by Hannah and Brown (1991). It implies successive approximations beginning from a first guess solution, and is using only quadrature. The pair of independent variables suitable for the current problem is the pair time-wavelength. To understand the importance of this pairing, we should refer next to the pressure-wavelength pair of variables in satellite spectroscopic observations. As stated before, the monochromatic limb scanning from satellite is based on observing the atmospheric outgoing radiation, for which there is a one-to-one correspondence emission wavelength- pressure, such that the profile of number concentration can be computed. 63 The attractive part of the iterative Hannah and Brown method for the current problem is the fact that, for an integral over wavelength, the method actually has a physical significance, amounting to continually re-calculating the effective wavelength of a slowly changing filter, or, equivalently, the center of mass of a slowly shape-changing integral. This continually adjusts the current step transmission with a new factor, which will be used to calculate the next step?s effective wavelength of the filter, and a new multiplicative correction for the transmission. When the ?saturation? in the correction is attained, then the sought transmission value, suitable for the current filter, has been constructed. An error control parameter, constructed simply from the successive corrections efficiency, can keep an account of the errors. Here are the details of the scheme: For a previous step (n-1), the transmission correction is C n-1 (t k ), interpolated to C n-1 (? j ). Then the n-th step of the new corrected filter?s effective wavelength is ? eff (t k ), while the new correction is C n (t k ), computed as: ? ? ? ? = ? ?? ? dFC dFC t n n keff 1 1 )( , (4.14) and ? ? = ?dFC ty tC n km kn 1 )( )( , (4.15) where F=(f TOA R T basic ) C 1 ?.C n-2 is the step (n-2) kernel, and FC n-1 is the current kernel obtained in the n-1 step, and FC n-1 C n is the kernel obtained in the n-th step, and will be used for the next one. 64 Again, for the next step, the term C n (t k ), which can be called ?the iterative ratio?, will have to be interpolated from the time domain to the wavelength domain C n (? i ). Besides a first step being the quadrature recommended for computing the integrals, Hannah and Brown devised their method around a second important step, - the interpolation/extrapolation of the correction coefficient (in our case, the transmission correction), explained as follows. They recommend the use of ?any approximation method, but least squares, rational or spline approximations might be natural candidates?. For our problem, this implies passing from: ? the set {C(t k )}, in bijective relation with respect to the set {? eff (t k )}, to ? the set C(? i ) of correction values, interpolated/extrapolated into the kernel?s initial wavelengths grid {? j }=[400,850] To obtain the new, interpolated, correction values for transmission, {C (? j )}, which is an exponential functional in ? (? ) and X (t), and a non-classical unknown for a FK1 problem, the following scheme was thought best. All measurements {y m (t k )}, taken at the same time t k , even for different filters (different spectral classes stars), were to be reduced with a unique optical depth, ? (? eff , t k ), and a corresponding known airmass function X (? (t k )). This interpolation condition has the significance that, in the average, it is possible to interpolate and extrapolate from an effective wavelength in the stars ?colored filters? to the entire range of the wavelengths of the total filter, [400, 850] nm. 65 Knowing the optical depth functional dependence on the wavelength will allow the delicate step of interpolating or extrapolating the transmission correction to follow naturally within the constraints of a physically significant law. Iterative unfolding or ?ratio method? is then a method that relies on the suitability of smoothing the data and interpolating the data for the next step, without losing too much sensible information. The constraint is in fact provided by the requirement that the wavelength dependence of the AOD is functional with a power form on wavelength. This is the same as assuming that the aerosol sizes obey a Junge distribution when integrated over the whole vertical column and for 3 orders in particle radius (from 0.01 to 10 microns). Kythe and Puri (2000), and Hanna and Brown (1991) recommend a one point Gaussian quadrature as the first guess. That is the same as taking the first guess of lambda effective of the compounded filter as being the lambda effective of the same filter with a kernel from which the aerosol extinction part is missing. There are a few other possibilities for choosing the initial value of the effective wavelength and the ratio ?C 0 ?: to construct these two values from a priori climatologic data of Angstrom parameters alpha and beta. The alpha exponent is given as a monthly mean for the site, in order, from January to December: 1.8, 1.9, 2, 1.95, 1.8, 1.78, 1.7, 2.1, 2.05, 2, 2, 2.1, 1.9, or as hourly 4-year mean for the nighttime: increasing in 0?12 hours local time from 1.8 to 1.9 (Sheridan and Ogren, 1999; Sheridan et al, 2001). After the algorithm has converged, the solution is found as the product of the successive iteration ratios ?C i ?, noted as C. The relative error of the algorithm is given 66 by (C-1)/C. When at the step n, the ratio C n ~1, the solution has converged. The AOD can now be computed as the optical depth producing the product of successive approximations of the transmission in many different colored filters (stars). From the solution at each computation wavelength, the parameter alpha is calculated, as the same between any two wavelengths. But it is thought as representative only for the real alpha coefficient for AODs between 500 and 650 nm. On this point we will comment in Chapter 5. The alpha Angstrom parameter is related to the particle size and represents here a ?mean size? that emerges from a smoothing (alpha can be allowed only within believable values [0,3] or [0,5], and then taken as a mean) and an interpolation (within Junge distribution of aerosol sizes). The turbidity parameter beta is also smoothed with the median. The smoothing is a filtering that allows to be taken into account only the alphas that are within a realistic range of values for aerosol transmission. It ought to be clear that the quantity determined by this scheme applied to real data is a residual optical depth, that is, it might include not only the aerosols, but also very thin clouds. This is a most common occurrence in the WSI data; there is almost never a perfect clear-sky day, there is at least one direction and path ? usually at low elevations - for which a star is observed through a thin Cirrus, with low opacity. We did not eliminate the cases in which extremely thin clouds appear to be present at the horizon. It is known that the form of the spectral AOD dependence, the power law, has been proved often suitable for other absorbers, for example it has been fitted to the water vapor (Michalsky et al, 1995), in the form 67 T= exp [-const (X abs ) b PWV b ], (4.16) where T is the transmission of the water vapor layer having a column abundance of ?PWV? (in gr/cm2) and seen through a length of path of X abs . 4.3 Stellar spectra used in retrieval. The stellar library of observed spectra used in this work is the ?UBVLIB? library compiled by Pickles (1989); each spectral interval is a weighted mean of multiple spectral observations of same type stars, published by different authors. The important feature of this library is that the multiple sources of the measurements, which are often at different resolutions and spectral samplings, are combined, filtered, checked for conformity with the spectral energy distribution from the Gunn and Stryker (1983) observational atlas. All spectra are normalized to unity at 5556 A (a wavelength near the mean wavelength of the visual filter, 0.55 ?m). The library contains all spectral Morgan-Keenan (?MK?) types (O, B, A, F, G, K, M) and all luminosity (I supegiants, II bright giants, III giants, IV subgiants and V dwarfs ?main sequence) types, for normal stars (solar composition) and for metal-weak and metal- rich stars. For the optical range, there are 6 spectral atlases combined (when overlapping), for better precision. Details of dereddening of the observed spectra are given in Pickles (1989). The 131 flux-calibrated spectra encompass the interval 115-2500 nm, in steps of 5 A. The library is used here to obtain synthetic normalized calibrated flux in the particular filter of the WSI instrument, and in the atmospheric filter. For use in the 68 atmospheric model, the spectra have been degraded to 5 nm steps, which is sufficient for the LOWTRAN model in the visible and near infrared. The spectra are also useful, in connection with the absolute calibrated flux for the standard star Vega (Oke and Schild, 1970; Tug et al, 1977; Wamsteker, 1981; Budding, 1992; Cox, 2000), to obtain monochromatic (555.6 nm) fluxes for the other stars observed. Table 4.1 gives the absolute flux values for a standard A0V star with a magnitude of either 0 or 0.03 (Vega) and the error 1-2% associated with the measurement by different authors. The spectra in Fig 4.3 are given as example of the spectra of dwarf stars (luminosity class V), for MK spectral types O to M, and are clearly labeled with the spectral and luminosity class. Source Flux of 0mag. star at 5556A Mean flux of 0mag. star in V filter (5504A) Flux of 0.03mag. star at 5556A Mean flux of 0.03mag. star in V filter (a) 0.348 0.364 0.336+/-0.002 0.352 (b) - - 0.347+/-0.001 - (c) - 0.364+/-0.002 - - (d) - 0.354 - - (e) - 0.375 0.344+/-0.002 - Table 4.1 The absolute calibrated flux of stars of 0 mag. and 0.003 mag. at a monochromatic wavelength of 5556A and in the V filter, as measured of quoted in the literature. All fluxes are in units of 10-7 mW m-2 nm-1. References: (a) Oke and Schild, 1970; (b) Tug et al, 1977; c) Wamsteker, 1981; (d) Budding, 1992; (e) Cox, 2000. Fig 4.3 (next page) Normalized spectra of ?main sequence? stars, each clearly labeled with the spectral and luminosity class in the upper left corner. From upper left to lower right, the star spectra of types B0V, A0V Vega (Alpha Lyrae), F0V, G2V (the Sun), K0V, M0V are represented with the highest continuous line. Also represented is the star spectrum convoluted with the WSI responsivity curve (lowest broader line). For comparison, the spectrum in the Johnson V filter, centered at 550 nm is also represented. 69 70 Monochromatic fluxes at 556 nm were constructed from the catalog visual magnitude, the spectra and the absolute calibration of A0V star flux. The error introduced by this step is discussed below. A supplementary check of the result was made with other spectroscopic measurements in the range 370-950 nm (Petford et al, 1985). The magnitude catalogs give a visual magnitude with an error of +/- 0.01mag (Jaschek and Jaschek, 1987), which is approximately a 1% precision in flux. For all computations we adopt the monochromatic value of flux at 555.6 nm for Vega as (0.344+/-0.002) 10-7 mW m-2 nm-1. Visual magnitude is defined for each star on a decimal logarithmic scale of flux ratio in a well-known standardized ?V? filter. We have computed the flux in the ?V? filter by integrating star normalized fluxes from the same catalog. Inferring a monochromatic magnitude from the visual magnitude and the absolute flux of Vega is a step usually avoided by photometrists, who work with the extinctions of colors. Here this is an unavoidable step. The reader not familiar with catalog magnitudes should know that they are given with a precision of 0.01mag (Jascheck and Jascheck, 1987). With this information, the computation yields the relative error in the flux value at 555.6 nm for any star of 0.92%. Interstellar reddening was allowed for, following the total extinction curve described in Fitzpatrick and Massa, 1999. The ratio of selective to total extinction in the V band, R (V), was taken as 3.1: R (V)=A (V)/ E (B-V), (4.17) where A(V) is the interstellar extinction in the visual band and E(B-V)=E(B)-E(V) is the color excess for the visual and blue filters of the Johnson system. Recall that the 71 color excess E (V) is defined as the difference between the star intrinsic magnitude V 0 , related to the energy emitted in the respective passband at the surface of the star, and the star observed magnitude V, which has been reddened by the interstellar dust scattering: E (V)=V-V 0 . (4.18) The interstellar transmission was written as a factor to the intrinsic flux: F observed (? )= F 0 (? ) dex (-0.4 A V f ism (? ) / 3.1), (4.19) where dex is the decimal exponent, A V is the interstellar extinction in the visual, f ism is the spectral mean extinction curve from a spectral energy distribution model of the galaxy. Outside the galaxy plane the ISM extinction decreases as the cosecant of the galactic latitude (Neckel and Klare, 1980). The stars were individually studied for distance (from astronomical paralax) and color excess, and the catalog of interstellar extinction for known lines of sight towards the Milky Way ?bulge? was consulted from Neckel and Klare, 1980. The interstellar extinction was considered for stars that are reddened (having E (B- V)>0.1), for whom the ISM extinction in the visual filter (A V ) is known, and which are at a distance of over 100 parsecs. Table 4.2 presents some of such stars (the list is longer). Star A V Deneb 0.36 Mirfak 0.34 CmaDelta-25 0.29 Table 4.2 Example of a few of the stars for which the reddening was considered important to be taken into account. 72 4.4 Atmospheric scattered light from non-stellar sources The equation of the night sky brightness I sky is given by: scaEBLDGLZLISLAGsky IeIIIIII diff +++++= ?? )( , (4.20) where the sources airglow (AG), integrated starlight (ISL), zodiacal light (ZL), diffuse galactic light (DGL) and extragalactic background light (EBL) are attenuated in the atmosphere by a diffusion opacity ? diff , and scattered by the molecule and aerosols, giving the contribution I sca (Leinert et al, 1998). The airglow constitute more than 50% of the total sky brightness, it is correlated with the solar radio flux at 10.7 cm wavelength, and has a variation between nights and within a single night which is not linear in time; besides an approximately 11-years cycle of variation, which has been measured in all AG components, a seasonal component for the AG variation in time has been established (Leinert et al, 1995). In the absence of the atmosphere, the airglow profile with the zenith would be given by a simple van Rhijn function, i.e. a geometrical enhancement of the emission layer thickness. That is, the AG increases with the zenith angle towards the horizon. In the presence of the atmosphere, the airglow is single and multiple scattered by air molecules, giving a supplementary component, I sca R (AG), which is directly proportional with the molecular optical depth (practically, with the air number density), and increases towards the horizon (Staude, 1975). The airglow is also scattered by aerosols (principally single scattered, due to the enhanced asymmetry of the aerosol scattering phase function), producing a supplementary component, I sca M (AG). 73 The other two extended sources present in the WSI images, the zodiacal light and the integrated starlight, have a very spatially localized presence. The ZL is the sunlight scattered by the interplanetary dust, and has a Fraunhofer spectrum (practically, the spectrum of the Sun), and is present, in the solar system, only in the inner part of the asteroid belt orbit. As seen by the WSI, the ZL appears as a cone with the apex towards the zenith direction and the base towards the sun dip. The atmospheric scattered ZL is then spatially limited, and the aerosol-scattered part of the ZL is further spatially localized by the fact that the aerosol scattering phase function peaks in the forward direction. The integrated starlight is also spatially localized towards the plane of the Milky Way, where most of the non-resolved stars can be found. ISL can be inferred for every galactic latitude and longitude from star counting algorithms, by integrating over star atlases and positional catalogs the fluxes of all the known stars fainter than the last magnitude stars resolved by a given instrument; the integration can also be made over a reliable spectral energy density model of the Galaxy. ISL has also been measured in the red by Pioneer 10 mission, which found that the integrated starlight towards the galactic centre is double than that towards the galactic poles. For the WSI images, the ISL consists of the stars fainter than roughly visual magnitude 5. In the Fig 4.4 (after Turnrose, 1974) the spectral surface brightness of the sky at Mt. Palomar was measured in a dark part of the sky, near the North Ecliptic Pole (i.e., on the direction perpendicular to the solar system?s plane). Assuming the same spectral distribution to be accurate for the SGP Oklahoma, -which may indeed be true, if one assumes no scattered light towards the zenith direction from sodium-based 74 night lamps and city light pollution, and integrating over the Mt. Palomar spectrum, we obtained a value of 38.32 10-7 mW/ m 2 / sr /nm as the mean spectral radiance in the WSI filter. This is very close to the zenith value recording during particular dark nights in the WSI images (33 10-7 mW/ m 2 / sr / nm), when the aerosol loading, and thus the Mie scattered airglow, were low. Other nighttime spectra of the dark sky brightness are available, with better resolution, from Mt. Hopkins and Kitt Peak in Arizona (Massey et al, 1990; Massey et al, 2000), but were not used here. Fig 4.4 The spectral surface brightness near the Draco Cloud (17h 53m 50.7s; 66 deg 19 arcmin 7 arcsec). Integrated in the WSI filter: 38.32 10-7 [mW m-2 sr-2 nm-1]. The theory of the correction for the high background at large zenith value is as follows. The radiance in each pixel is known as characterizing the incoming brightness in the solid angle subtended by the pixel. The main assumption of the sky background computation for a star object is that the ?On/Off? technique can render the star flux as the (atmospheric attenuated) value above a sky background. There are though two instances in which the pixel-to-pixel variation of the sky brightness makes 75 it impossible to calculate the value of the sky brightness in one pixel from that in the neighboring pixels. When the (projected) pixel surface on the emitting layer strongly increases with radial distance (effects of pincushion distortion), the AG emission within a pixel has a radial increase. A central pixel sees the layer perpendicularly and the layer has, say, 1 km 2 surface, but an equal dimension pixel from the CCD margin is seeing a much larger surface on the layer, which has the same specific intensity (surface brightness). In the same pixels, the atmospheric scattering of the airglow obeys a similar pattern: it is radially amplified at a rate that makes it difficult to reconstitute the ?sky? value in one pixel from the value in the neighboring pixels. The sky gradient within one pixel is unresolved, because the coarseness of the pixel resolution. Consequently, the sky value for the absolute star flux computation in the pixels adjacent to a star has too different (usually higher) a value that the sought pixel under the star, and the star flux ? extracted as an ?On/Off? - will be too small. This ?engulfing? phenomenon is therefore a sum of two effects: ? the amplification of the sources? emission and scattering with the increasing depth of the emitting and scattering layer. This is the enhancement with the zenith angle of the emitting AG per pixel (van Rhijn) and of the atmospheric scattered AG (Figs. 4.5-4.6). ? the radial increase in the surface brightness per pixel, which is directly proportional to the pixel footprint (i.e., surface seen on the layer) increase. The part of the sky that is, in any pixel, above its neighbor?s value for background sky, can be computed from the sky value of the airglow at the zenith of 76 the place, enhanced for the AG and scattered AG at the current zenith angle, and further amplified by the differential increase in the pixel footprint. It is essential to understand that this is not an instrument calibration flaw, but simply an inter-play of the emission/scattering phenomena, viewing geometry, and coarse pixelization. In the Chapter 3 an equation was established for the transformation of any pixel position into polar coordinates measured on the layer: (curved) range and azimuth. The surface of a pixel is modified by distortion into an area between four ? roughly- hyperbolas, which can be computed. Fig 4.5 presents the normalized intensity of a molecular scattered light from an extended, unit, uniform source as a function of the zenith, for values of the optical depth from 0.05, 0.10 and 0.20 (after Leinert, 1998). Fig 4.5 Airglow scattered by a Rayleigh atmosphere with scattering optical depth of 0.2, 0.1 and 0.05. The aerosol atmosphere scattering of a unit extended source is presented in Fig 4.6 (after Leinert, 1998). It can be seen that the sky extended sources produce an 77 atmospheric scattering, which has an increasing radial gradient. The WSI increasing footprint effect will further amplify the radial gradients of the sky in the image. Fig 4.6 Airglow scattered by an aerosol atmosphere with scattering optical depth of 0.2, 0.1 and 0.05. 4.5 Summary of Chapter 4 In this chapter we described the method for obtaining the aerosol transmission at the effective wavelength of various observed stars at each observational time, from the heterochromatic measurements. For all simultaneous pairs of observations, the aerosol parameters alpha and beta are determined. Cases when both stars have the same effective wavelength are excluded (i.e., they have a common spectral type). For each time, the method involves performing repeated iterations until the correction in optical depth becomes negligible (under a preset value of the convergence error). The choice of the spectral star catalog and the subsequent computation of the star irradiance in the instrument?s filter were also described. 78 The last part of the chapter explained how to compute the background sky level for zenith angles close to horizon (greater than 70 degrees). This is a correction applied to the median sky level computed from a pixel elliptic aperture around a star (explained in Chapter 2), and is specific for the WSI image at high zenith angles, in which the bright night-sky extended sources engulf the starlight; it is not a correction required in narrow-field astronomical observations, sheltered from high sky gradients. 79 Chapter 5: Results and uncertainties, conclusions and future work This chapter is dedicated to the assessment of instrumental and systematic uncertainties in the retrieval method of the aerosol optical depth from the whole sky imager star measurements. A second part presents actual results of the analysis of approximately 300 nights in which a sky free of clouds could be observed, at least for several hours, comparisons with Raman Lidar retrievals at a wavelength outside the WSI range, and comparisons with sun photometer successive results. 5.1 Errors calculation ? Validation of the technique of aerosol optical depth calculation In this subchapter we will evaluate the uncertainties associated with the aerosol optical depth calculation from the broadband measurements with the Whole Sky Imager of stars pertaining to various known spectral classes and luminosities. As the method to infer the aerosol optical depth combines theoretical assumptions and actual measurements, the random error and systematic error are closely intertwined. For example, in Chapter 3 it was shown that the calibration constant introduces a systematic uncertainty in the absolute flux measurement, while the electronics introduce a shot noise in the star flux measurement, and that the night sky background introduces a Poisson noise in the star signal. Also, Chapter 3 presented the errors introduced in the experiment?s ?pipeline?, beginning with the original CCD frames acquisition. In Chapter 4 it was shown that the stellar magnitude and interstellar extinction data have a limited precision. The use of an atmospheric model based on 80 LOWTRAN model can also limit the accuracy (Subsection 1). Subsection 2 will discuss the performance of the iterative unfolding routine for aerosol depth computation and the inherent loss of accuracy because the instrument is broadband. 5.1.1 Systematic error of the atmospheric kernel (SBDART sensitivity to non- solar spectral flux) The SBDART model, which is used in the atmospheric treatment with the exception of the aerosols, is quoted as within 1% of the direct observation of shortwave direct-normal irradiance observed with a normal-incidence pyrheliometer (Ricchiazzi et al, 1998). Based on a LOWTRAN model, which uses a k-distribution band model (Thomas and Stamnes, 1999), SBDART has a 20 cm-1 spectral resolution, which at 400 nm is equivalent to a 0.3 nm resolution, at 500 nm, to a 0.5 nm resolution, and at 850 nm, to a 1.5 nm resolution. The molecular absorption part of the model includes absorption cross section data for ozone, precipitable water vapor and 10 other molecular species which are given in average (climatological) concentrations, and are, in decreasing order of concentration: CO2, CH4, N2O, CO, NH3, SO2, NO, HNO3, NO2, NO3. The well-mixed gases - with constant or almost constant density vertical profile ? are: O2, N2, CO2, CH4, CO, and N2O (McClatchey et al, 1972). As emphasized by Thomas and Stamnes (1999), the current transmission models use the HITRAN data base of lines obtained from laboratory measurements of molecular transmission for given temperature and pressure conditions. A band model like MODTRAN is less precise than a line-by-line code in treating the molecular 81 absorption, but in the shortwave domain, which is of interest for WSI, the error is very small, 0.5%, for the integrated flux over 200-870 nm (Halthore et al, 1997). The most variable components of the clear-sky atmosphere are the ozone and water vapor abundances and the aerosol loading. ? The total column for ozone is obtained from the daily measurements of the TOMS instrument aboard Nimbus satellite, around 17-18 UT (11-12 Local Time). The ozone abundance is estimated from 3 pairs of wavelengths in the UV region, where the ozone is strongly absorbing. The ozone total column abundance, for the days considered, has values between the extremes of 246 DU and 428 DU (DU-Dobson Unit; for the days of 19981115 and 19990416, respectively). Using the climatological value instead of the current value, the uncertainty is 20% (Halthore et al, 1997). The TOMS ozone value used was obtained by looking under a certain viewing angle at the backscattered solar radiation, for solar zenith angles that are consistently below 64 degrees at the point of observation. Dahlback et al (1994) have shown that stratospheric aerosols (at 15-20 km height) can lead to an overestimate of 10-100 DU for the ozone determined from TOMS UV measurements made for solar zenith angles over 70 degrees, while the errors are under 10 DU for smaller zenith angles. In this thesis, a conservative value for the error in the column abundance retrieval by TOMS will be taken as 25 DU. ? The total precipitable water vapor (PWV) in the column can be provided by either the balloon sounding or by a microwave radiometer (MWR), which has a one-minute resolution. A mean value of the total precipitable water from the 82 MWR instrument around the hours when the balloon sounding is performed is used as scaling in the total PWV in the sounding (usually, at the hours 0230, 0530, 0830 and 1130 UT). The PWV for the nights of interest cover a range from 0.27 g/cm2 to 4.86 g/cm2 (for 20020125 and 19980801, respectively). The standard deviation of the MWR PWV measurement is usually very small, around 2% (+/- 0.006 g/cm2) and more rarely reaches 9% (for example: 0.055 g/cm2). In turn, the MWR-scaled balloon sounding PWV profiles agree with simultaneous PWV measurements by lidar to within 1%-2% (bias) (Ferrare et al, 2000). ? The pressure, temperature and relative humidity vertical profiles retrieved by the balloon sounding, completed, when necessary, by the AFGL standard profiles for midlatitude summer or winter (McClatchey et al, 1972), provide the data for the computation of the air density with a temporal resolution of 3 hours (when all 4 soundings exist). When a smaller number of nocturnal soundings exist, or when any of all soundings are missing, the vertical stratification of the air is assumed to be that given by the nearest night sounding. In this latter case, the total columnar water vapor continues to be taken as the instantaneous value measured by the MWR. The resolution of LOWTRAN and MODTRAN in the visible is very close, because the HITRAN database of lines ends practically at 435 nm (Thomas and Stamnes, 1999). The MODTRAN sensitivity to ozone and PWV has been assessed in the literature for the flux of a G2V star source (the Sun) and found to be 0.41% and 83 0.51% flux relative uncertainty, respectively. The combined uncertainty of the solar total flux due to these two uncorrelated absorbers is then 0.65% (Halthore et al., 1997). The variations assumed as characteristic rms uncertainties of ozone and water vapor were 28% and 10%, respectively. The stars in the lower metallicity range - O, B, A, - which have a spectrum with more weight towards lower wavelengths, are comparatively more sensitive to ozone absorption, while the G, K, M-type stars are more sensitive to water vapor absorption. Nevertheless, the stellar flux computed with LOWTRAN has the sensitivity to ozone and water vapor within a 1% relative uncertainty. An example of the actual shape of the spectra for O and M stars at the top of the atmosphere is given in Fig 5.2. The spectra (Pickles, 1999) are de-reddened observations of many same-type stars, which are normalized to 555.6 nm monochromatic fluxes. The spectra represented are not yet convolved with the atmospheric filter. It can be seen that their effective wavelength is different enough to be a good discriminating factor in the aerosol computation, even for a detector that observes only the integral flux. To assess the relative diverse spectral contribution by the diverse attenuators in the atmosphere to the extinction of the ground received flux, an example of the comparative strength of absorption opacity versus scattering opacity in the passband of the WSI follows. 84 Fig. 5.1 Normalized spectra of two stars of contrasting types. Left: a main sequence type O9, e.g. Mintaka (Delta Orionis). Right: a supergiant of type M2, e.g. Antares (Alpha Scorpii). The star spectra are represented with the highest continuous line. Also represented is the star spectrum convoluted with the WSI responsivity curve (lowest broader line). For comparison, the spectrum in the Johnson V filter, centered at 550 nm is also represented. In the visible and ultraviolet part of the spectrum, the energy of the photon is high enough to permit electronic bound-bound or bound-free transitions. The vibrational-rotational spectrum of dipolar molecules displaces the electronic transition frequency, which will occur at other frequency than the original electronic transition frequency (Kyle, 1991). For example, molecular oxygen electronic transition band at 13120 cm-1 (762 nm) is also displaced to the frequencies, in the visible, of 11564, 12969, 14525 and 15902 cm-1, thus has bands at the wavelength of 865, 771, 688 and 629 nm. The absorbing bands in visible for water vapor are roughly given by the intervals of [633; 667] nm, [700; 740] nm (?alpha band?) and [794; 870] nm (?0.8 microns band?). The ozone absorption in the visible is a quasicontinuum, which obeys the Beer?s law in the [440; 750] nm ?Chappuis? band, with a center in the 580- 605nm region (Anderson and Mauersberger, 1992). While the ozone opacity has a 85 lower maximum than that of the water and oxygen molecules, its influence encompasses a broader band of wavelengths. The changes in atmospheric optical depth near the wavelength center of the absorption bands mentioned above, and which fall inside the WSI passband, are given in the Table 5.1. It is useful also to compare the molecular absorption depth near the center of the respective bands with the molecular scattering depth in the same region, when the ozone abundance or the PWV are minimum, mean and maximum for the data analyzed with the WSI. The actual dates when these events occurred have been given above. The Rayleigh (molecular scattering) optical depth at the current wavelength is also given in parenthesis, in the columns ?f? and ?h? of the Table. The scattering optical depth is higher than the ozone highest optical depth around 605 nm, and higher than the oxygen optical depth in the bands centered near 629 nm and 771 nm, but is less important in the water vapor absorption bands range, and in the other molecular oxygen bands. 86 Band (nm) Wave- length near the band center (approx.) Lowest abun- dance observed of current absorber Optical depth (at the center of band) Abundance near the mean observed value Optical depth at the center of band Highest abun- dance observed of current absorber Optical depth (at the center of band) A b c d e f g h [440;750] 580?605 nm O3 246DU 0.032 318DU 0.043 428DU 0.055 (0.078) [633;667] 655nm H2O 0.216cm 0.007 2.99cm 0.038 4.77cm 0.054 (0.049) [700;740] 718nm H2O 0.216cm 0.071 2.99cm 0.342 4.77cm 0.473 (0.033) [794;869] 821nm H2O 0.216cm 0.036 2.99cm 0.165 4.77cm 0.226 (0.018) [628;637] 629nm O2 - 0.236 g/cm/mbar 0.023 (0.056) - [687;699] 688nm O2 - 0.236 g/cm/mbar 0.104 (0.039) - [750;781] 771nm O2 - 0.236 g/cm/mbar 0.012 (0.025) - Table 5.1 Principal molecular absorption bands which fall within the WSI passbands. The minimum, mean and maximum values of the absorber abundance and its corresponding optical depth near the center of each molecular band for ozone, water vapor and molecular oxygen. For comparison, the molecular scattering optical depth at the current wavelength is given in parenthesis, in the columns ?f? and ?h?, for the maximum abundance of O3 and H2O, and for the climatologic value of O2 concentration, respectively. An example of transmittance values for two days with clear-sky is given in Fig. 5.2. The ozone and PWV values for the day of 19991005, at 0530 GMT and 1130 GMT, are 283.2 DU and 1.1 g/cm2 and 1.42 g/cm2, respectively, and for day 19991011, at 0530 GMT, the values are 277.8 DU and 1.85 g/cm2. 87 Fig. 5.2 The atmospheric beam transmittance as function of wavelength 400-1100 nm (continuous line) and the transmittance due to aerosol only (bold line crosses). The aerosol optical depth at 355 nm is, from upper-left to lower-right charts, increasing: 0.07 (19991005, 0530 UT), 0.15 (19991011, 0530 UT) and 0.18 (19991005, 1130 UT). 88 5.1.2 Aerosol optical depth computation uncertainties The errors in the AOD retrieval are given by the measurement and methodological uncertainties and have been detailed in the Chapter 3 and 4. Regarding the methodological errors, recall that the iterative unfolding algorithm, presented in the Chapter 4, finds a solution for a Fredholm equation of the first kind with a kernel whose variable part is the atmospheric transmission. The error in the a priori known part of the kernel has been computed in the previous section as the LOWTRAN sensitivity to the input of ozone and water vapor. The relative error in the solution of the algorithm that has converged is given by the relative deviation from unity of the last multiplying correction to the solution. The convergence is equivalent to adding to the AOD solution a last correction layer that is so thin that it has the opacity lower than the preset accuracy required for convergence. Instrumental uncertainties have been computed in Chapter 3 from equations for error propagation for the measured star irradiance and the calibration constant. In the same chapter, the errors due to the electronics have been computed for the operational temperature of the instrument. The uncertainty in calibrated measurement is governed by the uncertainty in the calibration constant, the relative error of which is 0.5%. Table 5.2 gives the estimation of the AOD error at the WSI filter band center (500-650 nm) from instrumental and methodological errors referred to in the previous chapters. 89 Random uncertainty Systematic uncertainty Absolute Flux calibration constant 0.5% (Ch.3) Star flux accuracy (includes error from dark current subtraction, readnoise, digitization noise, flatfield division and bias subtraction for the CCD frame) < 0.03 rms (Ch.3) Solid angle per pixel 0.71% (Ch.3) Pixel surface projection on constant height layer 0.22% (Ch.3) Background sky correction (for large zenith angle only) 0.26% Vega monochromatic flux 0.002(10-7mW/m2/nm) rms (Ch.4) Catalog Magnitude 0.01mag (Ch.4) Numerical integration 10-12 Monochromatic absolute flux (other stars) 1% (Ch.4) Star calibrated measurement (CM) 1.11%(Vega) 1.37% (other) (Ch.3) Relative airmass (<84 deg.zenith) < 0.065% LOWTRAN model sensitivity (without aerosols) 0.72% (Ch.5) Kernel 0.72% (Ch.5) Per step iteration (from CM & Kernel) (1.32-1.55)% (Ch.5) Convergence error 10-5 AOD at centre of bandpass (2.6-3.1)% (Ch.5) Table 5.2. Total error for the AOD retrieval from random and systematic uncertainties, given either as standard deviation or, when not mentioned, as relative error. See text for the computation relations. As stated in Chapter 4, initial AOD values serve as input for the method of unfolding. These inputs represent a first guess of the spectral curvature and asymptotic value. Two input spectral dependencies for the first-guess AOD were submitted for the 66 frames measured on 20011024, with the (alpha, beta) aerosol parameters being (2, 0.002) and (2.5, 0.020). The output values for the AOD are correlated with a linear correlation coefficient that is above 0.99 for, in decreasing order, the aerosol optical 90 depth at 550, 600, 500, 650 nm, and with a correlation coefficient above 0.955 for the rest of the depths. The graphs for this example are presented in Figs. 5.3 and 5.4. Fig 5.5 shows the AOD autocorrelation as a function of wavelength. Fig 5.6 presents the autocorrelation of the Angstrom parameters alpha and beta, and of the simulated AOD at 355nm (which is outside the WSI passband), for the night of 19980904, with two sets of first guess parameters as input. In practice, it has been observed that submitting low absolute values for the first guess of AOD, and with a small spectral curvature produces alpha and beta values in the second iteration step that are within 5% of the convergence value. This is to say that the second iteration of unfolding (expansion) is already close to the result of unfolding. Fig. 5.3 Autocorrelation of the AOD extraction model output for two input series of parameters, for three wavelengths at the center of the WSI bandpass: 550, 600, and 650 nm. Stars of all spectral classes have the effective wavelengths in the whole sky imager filter within the interval [550, 650] nm. 91 Fig. 5.4 Autocorrelation of the AOD extraction model output for two input series of parameters, for seven wavelengths not central of the WSI bandpass: 400, 450, 500, 700, 750, 800, 850 nm. Fig. 5.5 Aerosol optical depth obtained from the same set of measurements using two starting values has better autocorrelation at the center of the passband. 92 Fig. 5.6 Autocorrelation of the Angstrom parameters alpha and beta and of the simulated AOD at 355nm (which is outside the WSI passband), for the night of 19980904, with two sets of first guess parameters as input. 5.2 Typical Results Annex A shows the aerosol optical depths at 500 nm, computed, from the techniques discussed in the previous chapters, from star irradiance measured with the Whole Sky Imager. In these graphs, the WSI nighttime measurements are compared with the Sun Photometer daytime monochromatic (500 nm) measurements, which precede or immediately follow. The two measurements seem to agree in the values obtained near Sunrise and Sunset. It must be emphasized here also that the parameters alpha and beta of the aerosol spectral AOD from WSI are derived from the expression ?dln? / dln? which is considered as an accurate representation of the 550/650 nm pair of AODs, because the effective wavelength of stars is usually between 500 nm and 650 nm. The turbidity parameter beta is obtained from this alpha and the AOD at 600nm. As these 93 three AOD values have the smallest retrieval error, the errors in alpha and beta are kept to a minimum. Because the AOD in the WSI wavelength range are not measured directly during the night by other instruments, we cannot perform a study of the model versus observation cases for the AOD. A secondary comparison is presented here, with the simultaneous Raman Lidar AOD measurements at the monochromatic wavelength of 354.7 nm. The Raman Lidar has a field of view towards the zenith direction only and it does not measure the same scattering air volume as the Whole Sky Imager. For the WSI, few stars cross the zenith direction. In the Figures 5.7, ?5.17, several nights of data are presented for the 355 nm AOD. The WSI-simulated and the actual RL measurements are very close. The agreement may be an expression of the spatial correlation scale of the aerosol optical depth and composition. Fig. 5.7 Night of 19980926, total AOD. Left: WSI measured AOD at 600 nm (thick line) and the simulation of AOD at 355 nm (cross symbol) from WSI-inferred aerosol Angstrom exponent and turbidity coefficient. Right: Raman Lidar measurement AOD at 355 nm (cross symbol and line between points) from backscatter at the monochromatic laser wavelength and simulation with the WSI (diamond symbol). 94 Fig. 5.8 The night of 20020307. The comparative presentation of the RL AOD at 355 nm (line) and the AOD at 355 nm simulated from WSI (diamond symbol) inferred aerosol parameters. The alpha exponent value is around the value of 2.3, which is characteristic for pollution aerosols with diameters between 0.1 and 1 micron. Fig. 5.9 The night of 20011115: AOD at 355 nm from RL (line) and from WSI (diamond symbol). Fig. 5.10 The night of 20011220: comparison between AOD at 355 nm from RL (line) and from WSI (diamond symbol). 95 Fig. 5.11 Alfa parameter during the following daytime, from the Sun photometer. It can be seen that the curvature of the spectral AOD curve of the aerosols is not constant: it depends on the pair of wavelengths between which is computed. Fig. 5.12 Night of 20000110: RL AOD at 355 nm (line) and the AOD at 355 nm simulated from WSI inferred aerosol parameters (diamond symbol). The alpha exponent value is around the value of 2.5, which is not close to the AOS simultaneous value of around 0.5, which is characteristic for larger particles. 96 Fig. 5.13 Night of 19980519, a smoke advection episode. WSI simulated 355 nm AOD and RL AOD at 355 nm measurement. AOD is over 0.8. Fig. 5.14 The night of 20030329. RL 355 nm (line) AOD measurement and WSI (diamond symbol) simulated AOD at 355 nm (excluding two episodes with clouds around 07-08 Hours). The alpha Angstom exponent is also measured during the night at the lower end of the column by the Aerosol Observing System (AOS); the sampled aerosol is dried to less than 40% relative humidity before entering the scattering chamber of the nephelometer (Peppler et al, 2000; Sheridan et al, 2001). An example of the (lack of) correlation between the AOS alpha and the WSI alpha is presented in Fig 5.15 and 5.16. 97 Fig. 5.15 Night of 20030329: Alpha aerosol exponent variation in time from WSI (columnar) and from AOS (ground). Fig. 5.16 Aerosol Angstrom exponent alpha measured at the ground by the aerosol observing system (gray line) at the beginning of the night, as compared with the same parameter inferred from the WSI (black squares). The straight line in the AOS alpha value is a defect measurement. No other simultaneous alpha estimations exist. It is difficult to make quantitative comparisons between the WSI and RL aerosols because the equation in alpha and beta Angstrom parameters used for the extrapolation is primarily valid between 0.5 and 0.6 microns where the stars have their effective wavelength within the WSI and atmospheric filter. Other instrument- 98 based techniques, such as the CIMEL sun photometer, use different alphas for different pairs of wavelengths (i.e., there is generally not a unique alpha for any given day), as shown in the Fig 5.11. Furthermore, on days when there is a substantial horizontal aerosol gradient, we expect to see time lags between the WSI and RL, due to the instruments viewing different portions of the atmosphere. An example of such lag is given in the Fig 5.17. In view of these difficulties, the agreement between the WSI and RL using WSI data extrapolated from the visible to the UV appears to be generally very good for most cases. We leave a more detailed investigation of the differences for future research. Fig. 5.17 Night of 20011024: apparent 4 hours lag between the RL 355 nm AOD measurement (line) and simulated 355 nm AOD from WSI star measurements (diamond symbol). Besides the few examples given here, the whole set of days, with both WSI and Raman Lidar measurements, is represented in the Annex B. 5.3 Result discussion, conclusions and future work It has been demonstrated in this work that the WSI instrument can be successfully employed to obtain an independent evaluation of aerosol optical depth during the night. The WSI-inferred aerosol optical depth?s precision can be kept very 99 good, for a photometric measurement, but we do not know its absolute accuracy. The only simultaneous measurements, from Raman Lidar, are actively scanning a different atmospheric volume than the WSI. Nevertheless, we have found instances when the two AOD?s measurements are correlating very well, as in the examples given earlier in this chapter. The WSI instrument has obtained a wide value range for the aerosol optical depth variations in nighttime; the simplistic view that the extinction is constant during the night or at most drops linearly from dusk till dawn does not offer the complete picture of the actual variation. The WSI AOD determined for a certain moment are accurate because the observed aerosol airmass is within the spatial coherence distance of < 200 km computed by Anderson et al (2003). The 5 km height layer seen at a zenith angle of 87 degree is placed at 88 km range from the zenith. The maximum distance measured at constant height between stars observations through this layer, which occur for stars on the opposite ends of a meridian, is then a little under 180 km, thus within the limit indicated. There are other experiments that endeavor to combine insight about aerosols from diverse instruments, with different results. Bergin et al. (2000) attempted to simulate the columnar aerosol properties from the ground aerosol properties when the boundary layer is well mixed (local noon). The mixing height was determined form the change in the sign of the vertical temperature profile. The aerosol scattering and absorbing cross sections, determined from AOS at an ambient humidity of 20%, were extrapolated using a hygroscopic growth factor for aerosols corresponding to the mean relative humidity in the aerosol 0-2 km layer. The simulation gave a synthetic 100 AOD, which was 40% less than the columnar AOD measured with the sun photometer, and the discrepancy was attributed to the aerosol loading above 2 km. Neither succeeded to simulate the total AOD based on the extinction coefficient measured at the surface; the coefficient is only poorly correlated with the AOD, and the synthetic aerosol optical depth is again underestimated with 30%. Taking advantage of the discriminating powers of the lidar zenith-direction observations, the significant vertical scale of the variability can be separated from the horizontal scale. It has been observed that the extinction-to-backscattering ratio (S a ) is inversely correlated (0.72 linear correlation coefficient) with the alpha exponent, that is, S a is directly correlated with the size of the aerosol particles (Ferrare et al, 2000). Extinction-to-backscattering ratio is known as a function of altitude and time, and the value in the vertical is not constant either in time, or in space. In Eulerian representation, with the aerosol layers drifting above the lidar location, the variation in S a was interpreted as a particle size spatial variability over both horizontal and vertical directions. Our analysis has emphasized the temporal dimension of the WSI measurements. The assumptions we have made were largely based on an a priori common vertical optical thickness for the volumes of air seen in each star?s direction. The next step would be to allow for a directional variability, in the sense that we should not assume that two stars having the same zenith angle (and seen through a same formal relative airmass), but at different azimuths, are attenuated similarly. One of the difficulties of such approach would be the fact that the star sources for which extinction is computed are not evenly distributed, or even conveniently distributed 101 over the sky, to endeavor correlations of extinction values over smaller azimuth ranges. There are of course the directions towards the Milky Way disk that offer the bulk of the star presence, and the directions towards the galaxy poles where the stars are scarce. The solution to this unevenness is to simulate the scenery of the extended sources in each pixel, and use the entire array of pixels to study the beam transmittance. But this would assume that one knows precisely enough the extended sources at the top of the atmosphere (ISL, ZL, etc.) or inside the uppermost emitting layers of the atmosphere (AG). In spite of difficulties, this can be a way in the future to determine the directional aerosol variability. 102 Appendix A The Sun photometer CIMEL monochromatic measurement of aerosol optical depth at 500 nm (diamond symbol), and the Whole Sky Imager aerosol optical depth at the same wavelength (square symbol), for a selection of dates from the calculations. Each graph is clearly labeled with the date of observation. 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 Appendix B Aerosol optical depth measured by the Raman Lidar at 355 nm (line), and the Whole Sky Imager aerosol optical depth simulated for the same wavelength (diamond symbol), for a selection of dates from the calculations. Each graph is clearly labeled with the date of observation. It should be noted that this wavelength is outside the WSI spectral sensitivity. 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 Bibliography Alard, C., R.H.Lupton, A method for optimal image subtraction, Astrophys.J., 503, 325-331, 1998. Aldering, G.S., G.D.Bothun, A fast focal-reducing camera for CCD imaging, Publ. Astron.Soc.Pacific, 103, 1296-1313, 1991. Anderson, S.M., K.Mauersberger, Laser measurements of ozone absorption cross sections in the Chappuis band, Geophys.Res.Letters, 19 (9), 933-936, 1992. Anderson, T.L., D.S.Covert, S.F.Marshall, M.L.Laucks, R.J.Charlson, A.P.Waggoner, J.A.Ogren, R.Caldow, R.L.Holm, F.R.Quant, G.J.Sem, A.Wiedensohler, N.A.Ahlquist, T.S.Bates, Performance Characteristics of a high sensitivity, three wavelength, total scatter/backscatter nephelometer, J.Atmos.Ocean.Technol., 13, 967-986, 1996. Anderson, T.L., D.S.Covert, J.D.Wheeler, J.M.Harris, K.D.Perry, B.E.Trost, D.J.Jaffe, J.A.Ogren, Aerosol back scatter fraction and single scattering albedo : Measured values and uncertainties at a coastal station in the Pacific Northwest, J.Geophys.Res., 104, 26793-26807, 1999. Anderson, T. L., R.Charlson, D.M.Winker, J.A.Ogren, K.Holmen, Mesoscale variations of tropospheric aerosols, J.Atmos.Sci., 60, 119-136, 2003. Andrews E., P.J.Sheridan, J.A.Ogren, In-situ Aerosol Profiles over the Southern Great Plains CART Site, in 11th ARM Science Team Proceedings, Atlanta, March 19-23, 2001. Ansmann, A., F.Wagner, D.Althausen, D.Muller, A.Herber, U.Wandinger, European pollution outbreaks during ACE 2: Lofted aerosol plumes observed with Raman lidar at the Portuguese coast, J.Geophys.Res., 106 (D18), 20725- 20733, 2001. AOS Aerosol Observation System, 2001: http://www.arm.gov/instruments/static/aos.stm 191 Bahreini, R., J.L.Jimenez, J.Wang, R.C.Flagan, J.H.Seinfeld, J.T.Jayne, D.R.Worsnop, Aircraft-based aerosol size and composition measurements during ACE-Asia using an Aerodyne aerosol mass spectrometer, J.Geophys.Res, 108 (D23), 8645, 2003. Bergin M.H., J.A.Ogren, R.N.Halthore, S.Nemesure, S.E.Schwartz, Aerosol optical depth estimates based on nephelometer measurements at the atmospheric radiation measurement program Southern Great Plains Site, on-line at www.arm.gov/publications, 1997. Bergin M.H., S.E.Schwartz, R.N.Halthore, J.A.Ogren, D.L.Hlavka, Comparison of aerosol optical depth inferred from surface measurements with that determined by Sun photometry for cloud-free conditions at a continental U.S.site, J.Geophys.Res., 105 (D5), 6807-6816, 2000. Bucholtz, A., Rayleigh scattering calculations for the terrestrial atmosphere, Appl.Opt., 34(15), 2765-2773, 1995. Budding, Edwin, Introduction to astronomical photometry, Cambridge Univ. Press, 272p., 1993. Buonanno, Roberto, G.Iannicola, Stellar photometry with big pixels, Publ.Astron.Soc.Pacific, 101, 294-301, 1989. Burkholder, J.B., R.K.Talukdar, Temperature dependence of the ozone absorption spectrum over the wavelength range 410 to 760 nm, Geophys.Res.Letters, 21(7), 581-584, 1994. Burki, G., F.Rufener, M.Burnet, C.Richard, A.Blecha, P.Bratschi, The atmospheric extinction at the E.S.O.La Silla observatory, Astron.Astrophys.Suppl.Ser., 112, 383-394, 1995. Chromey, F., D.A.Hasselbacher, The flat sky: calibration and background uniformity in wide-field astronomical images, Publ.Astron.Soc.Pacific, 944-949, 1996. Cimel Sun / Sky radiometer, 2001: http://www.arm.gov/instruments/static/cspot.stm 192 Cox, A.N.(ed.), Allen?s astrophysical quantities, 4 th edition, Springer, 719 p., 2000. DaCosta, G.S., Basic Photometry Techniques, in Astronomical CCD observing and reduction techniques, Publ.Astron.Soc.Pacific, 23, S.B Howell (ed.), 90-104, 1992. Dahlback, A., P.Rairoux, B.Stein, M.Del Guasta, E.Kyro, L.Stefanutti, N.Larsen, G.Braathen, Effects of stratospheric aerosols from the Mt.Pinatubo eruption on ozone measurements at Sodankyla, Finland in 1991/92, Geophys.Res.Letters, 21 (13), 1399-1402, 1994. Dick, J., C.Jenkins, J.Ziabicki, Design fundamentals of algorithms for photon- counting systems, Publ.Astron.Soc.Pacific, 101, 684-689, 1989. Dubovik, O., A.Smirnov, B.N.Holben, M.D.King, Y.J.Kaufman, T.F.Eck, I.Slutsker, Accuracy assessment of aerosol optical properties retrieved from Aerosol Robotic Network (AERONET) Sun and sky radiance measurements, J.Geophys.Res, 105(D8), 9791-9806, 2000. Elliot, J.L., E.W.Dunham, R.L.Baron, A.W.watts, S.E.Kruse, W.C.Rose, C.M.Gillespie, Image quality on the Kuiper airborne observatory. I. Results of the first flight series, Publ.Astron.Soc.Pacific, 101, 737-764, 1989. Ellingson, R.G., The ARM Instantaneous Radiative Flux (IRF) Working Group in 2000: A Summary of Accomplishments, Strengths, Weaknesses and Ideas for Future Activities, 15 p., ARM web site, 2001. Ellingson, R.G., M.Ba, 2002, A study of diurnal variation of OLR from the GOES Sounder, J.Atmos.Ocean.Tech., 20, 90-98. Ferrare, R.A., D.D.Turner, T.P.Tooman, L.A.Heilman, O.Dubovik, W.F.Feltz, R.N.Halthore, Characterization of the atmospheric state above the SGP using Raman Lidar and AERI/GOES measurements, Tenth ARM Science meeting proceedings, San Antonio, Texas, March 13-17, 2000. Fitzpatrick, E.L., D.Massa, Determining the physical properties of the B stars. I. Methodology and first results, Astrophys.J., 525, 1011-1023, 1999. 193 Fluks, M.A., P.S.The, On the correction of stellar spectra for the loss of radiation during its passage through the Earth?s atmosphere and through the spectrograph-slit, Astron.Astrophys., 255, 477-489, 1992. Garcia, F.J., M.J.Taylor, M.C.Kelley, Two-dimensional spectral analysis of mesospheric airglow image data, Appl.Optics, 36 (29), 7374-7385, 1997. Gilliland, R.L. and T.M.Brown, Time-resolved CCD photometry of an ensemble of stars, Publ.Astron.Soc.Pacific, 100, 754-765, 1988. Golay, M., Introduction to astronomical photometry, Reidel Publishing Comp., 364 p., 1974. Goldsmith, J.E.M., F.H. Blair, S.E. Bisson, and D.D. Turner, Turn-key Raman lidar for profiling atmospheric water vapor, clouds, and aerosols, Appl.Optics, 37, 4979 ? 4990, 1998. Gullixon, C.A., Two Dimensional Imagery, in Astronomical CCD observing and reduction techniques, Publ.Astron.Soc.Pacific, 23, S.B Howell (Ed.), 130-159, 1992. Gunn, James E., L.L.Stryker, Stellar spectrophotometric atlas, 3130 < ? < 10800 A, Astron.Astrophys.Suppl.Ser., 52, 121-153, 1983. Halthore, R.N., T.F.Eck, B.N.Holben, B.L.Markham, Sunphotometric measurements of atmospheric water vapor column abundance in the 940-nm band, J. Geophys.Res., 102 (D25), 4343 ? 4352, 1997. Halthore, R.N., S.E.Schwartz, J.J.Michalsky, G.P.Anderson, R.A.Ferrare, B.N.Holben, H.M.Ten Brink, Comparison of model estimated and measured direct-normal solar irradiance, J.Geophys.Res., 102 (D25), 29991-30002, 1997. Hanna, O.T., L.F.Brown, A new method for the numerical solution of Fredholm integral equations of the first kind, Chem.Eng.Sci., 46 (10), 2749-2753, 1991. 194 Hansen, P.C., Numerical tools for analysis and solution of Fredholm integral equations of the first kind, Inv. Probl., 8, 849-872, 1992. Harrison, L., J.Michalsky, Objective algorithms for the retrieval of optical depths from ground-based measurements, Appl.Optics, 33, 5126-5132, 1994. Hartley, W.S., P.V.Hobbs, J.L.Ross, P.B.Russell, J.M.Livingston, Properties of aerosols aloft relevant to direct radiative forcing off the mid-Atlantic coast of the U.S., J.Geophys. Res., 105 (D8), 9859-9885, 2000. Harvey, J.E., Krywonas, A., A system engineering analysis of image quality, in Current developments in lens design and optical system engineering, R.E.Fisher, R.B.Johnsonm W.J.Smith, W.H.Swantner (Eds.), Proc.SPIE, 4093, 379-388, 2000. Heitzenberg, J., R.J.Charlson, Design and applications of the integrating nephelometer: a review, J.Atmos.Ocean. Technol., 13, 987-1000, 1996. Highwood, E.J., Effect of cloud inhomogeneity on direct radiative forcing due to aerosols, J.Geophys.Res., 105 (D14), 17843-17852, 2000. Howell, S., Two-dimensional aperture photometry: signal-to-noise ratio of point source observations and optimal data extraction techniques, Publ.Astron.Soc.Pacific, 101, 616-622, 1989. Howell, S., Introduction to differential time-series astronomical photometry using CCDs, in Astronomical CCD observing and reduction techniques, Publ.Astron.Soc.Pacific, 23, S.B Howell (Ed.), 105-129, 1992. Howell, S., B.Koehn, E.Bowell, M.Hoffman, Detection and measurements of poorly sampled point sources imaged with 2-D arrays, Astron.J., 112 (3), 1302-1311, 1996. Jackson, J.D., Classical electrodynamics, 3 rd ed., 808 p., 1999. 195 Jaschek, C., M.Jaschek, The classification of stars, Cambridge Univ.Press, 413 p., 1987. Johnston, P.V., R.L.McKenzie, J.G.Keys, W.A.Matthews, Observations of depleted stratospheric NO2 following the Pinatubo volcanic eruption, Geophys.Res.Lett., 19 (2), 211-213, 1992. Kasten, F., A.T.Young, Revised optical air mass tables and approximation formula, Appl.Opt., 28 (22), 4735-4738, 1989. King, Ivan, Accuracy of measurement of star images on a pixel array, Publ.Astron.Soc.Pacific, 95, 163-168, 1983. Kumler, J.J., M.Bauer, Fisheye lens designs and their relative performance, Proc.SPIE, 4093, 360-369, 2000. Kythe, P.K., P.Puri, Computational methods for linear integral equations, Birkhauser, 508pp., 2002. Leinert, Ch., P.Vaisanen, K.Mattila, K.Lehtinen, Measurements of sky brightness at the Calar Alto Observatory, Astron.Astrophys.Suppl.Ser., 112, 99-121, 1995. Leinert, Ch., S.Bowyer, L.K.Haikala, M.S.Hanner, M.G.Hauser, A.-Ch.Levasseur- Regourd, I.Mann, K.Mattila, W.T.Reach, W.Schlosser, H.J.Staude, G.N.Toller, J.L.Weiland, J.L.Weinberg, and A.N.Witt, The 1997 reference of diffuse night sky brightness, Astron.Astrophys.Suppl.Ser., 127, 1-99, 1998. Lenoble, J., A general survey of the problem of aerosol climatic impact, Aerosols and their climate effects, H.E.Gerber and A.Deepak (Eds.), 279-294, 1984. Mamoulis, P., J.Macdonald, Geometrical Optics and Optical Design, Oxford University Press, 354 p., 1997. Martinez, P., A.Klotz, A practical guide to CCD astronomy, Cambridge University Press, 243p., 1998. 196 Massey, P., C.Gronwall, and C.A.Pilachowski, The spectrum of Kitt Peak night sky, Publ.Astron.Soc.Pacific, 102, 1046-1051, 1990. Massey, P., G.H.Jacoby, CCD data: the good, the bad and the ugly, Astronomical CCD observing and reduction techniques, Publ.Astron.Soc.Pacific, 23, S.B Howell (Ed.), 258-284, 1992. Massey, P. and C.B.Foltz, The spectrum of the night sky over Mount Hopkins and Kitt Peak: Changes after a decade, Publ.Astron.Soc.Pacific, 112, 566-573, 2000. Mattila, K., Has the optical extragalactic background light been detected?, Astrophys.J., 591, 119-124, 2003. McClatchey, R.A., R.E.Fenn, J.E.Selby, F.E.Volz, J.S.Garing, Optical properties of the atmosphere, Section 14 of the Handbook of optics, M.Bass (Ed.), McGraw-Hill, 1972. McKenzie, R.L., P.V.Johnston, C.T.McElroy, J.B.Kerr, S.Solomon, Altitude distributions of stratospheric constituents from ground-based measurements at twilight, J.Geophys.Res., 96 (D8), 15499-15511, 1991. Melfi, S.H., D.Whiteman, Observation of low-atmospheric moisture structure and its evolution using a Raman Lidar, Bull.Am.Meteorol.Soc., 66, 1288-1292, 1985. Michalsky, J., J. Liljegren, and L. Harrison, A comparison of sun photometer derivations of total column water vapor and ozone to standard measures of same at the Southern Great Plains Atmospheric Radiation Measurement site, J. Geophys. Res., 100 (D12), 25995 ? 26003, 1995. Mighell, Kenneth J., Accurate stellar photometry in crowded fields, Mon.Not.R.Astr.Soc., 238, 807-833, 1989. Miller, G.F.: Fredholm equations of the first kind, Numerical solutions of integral equations, Delves, Walsh (Eds.), 1979. 197 Miyamoto, K., Fisheye lens, J.Opt.Soc.Am., 54, 1060-1061, 1964. MPL Marine Physical Laboratory, 2004: http://www-mpl.ucsd.edu/ Musat, I.C., The Use of the ARM WSI to Estimate the Atmospheric Optical Depth at Night, Scholarly paper., 33 p., 2000. Neckel and Klare, Interstellar extinction catalog, Astron.Astrophys., 42, 251-283, 1980. Newberry, M.V., Signal-to-noise considerations for sky-subtracted CCD data, Publ. Astron.Soc.Pacific, 103, 122-130, 1991. Ogren, J.A., A systematic approach to in situ observations of aerosol properties, in Aerosol Forcing of Climate, R.J. Charlson and J.Heintzenberg (Eds.), pp.215- 226, 1995. Oikarinen, L., H.Saari, K.Rainio, J.Graeffe, H.Astola, Star-pointing spectrometer for measurements of atmospheric ozone, Proceedings SPIE, 2830, 224-235, 1996. Oke, J.B., R.E.Schild, The absolute spectral energy distribution of Alpha Lyrae, Astrophys.J., 161, 1015-1023, 1970. Peppler, R.A., C.P.Bahrmann, J.C.Barnard, J.R.Campbell, M.D.Cheng, R.A.Ferrare, R.N.Halthore, L.A.Heilman, D.L.Hlavka, N.S.Laulainen, C.J.Lin, J.A.Ogren, M.R.Poellot, L.A.Remer, K.Sassen, J.D.Spinhirne, M.E.Splitt, D.D.Turner, ARM Southern Great Plains Site Observations of the Smoke Pall Associated with the 1998 Central American Fires, Bull.Amer.Met.Soc., 81(11), 2563- 2591, 2000. Petford, A.D., S.K.Leggett, D.E.Blackwell, A.J.Booth, C.M.Mountain and M.J.Selby, Measurement of stellar integrated flux in the wavelength range 370nm- 950nm, Astron.Astrophys., 146, 195-198, 1985. Pickles, A.J., A Stellar Spectral Flux Library: 1150-25000 ?, Publ.Astron.Soc.Pacific, 110 (749), 863-878, 1998. 198 Redemann, J., R.P.Turco, K.N.Liou, P.B.Russell, R.W.Bergstrom, B.Schmid, J.M.Livingston, P.V.Hobbs, W.S.Hartley, S.Ismail, R.A.Ferrare, E.V.Browell, Retrieving the vertical structure of the effective aerosol complex index of refraction from a combination of aerosol in situ and remote sensing measurements during TARFOX, J.Geophys.Res., 105 (D8), 9949-9970, 2000. Reimann, H.-G., M.Boehm, W.Pfau, Numerical simulation of color equation for the uvby system, Astron.Nachr., 310 (1), 41-47, 1989. Reimann, H.-G., V.Ossenkopf, Correlation between atmospheric extinction and meteorological conditions, in Stellar Photometry, IAU Colloquium, 136, C.J.Butler, I Elliott (Eds.), 229-233, 1992a. Reimann, H.-G., V.Ossenkopf and S.Beyersdorfer, Atmospheric extinction and meteorological conditions: a long time photometric study, Astron.Astrophys., 265, 360-369, 1992b. Ricchiazzi, P., S.Yang, C.Gautier, D.Sowle, SBDART, a research and teaching software tool for plane-parallel radiative transfer in the Earth?s atmosphere, Bul.Am.Met.Soc., 79(10), 2101-2114, 1998. RL Raman LIDAR, 2001:http://www.arm.gov/instruments/static/rl.html Roberts, W.J., E.K.Grebel, Heterochromatic extinction. II. Dependence of atmospheric extinction on tellar temperature, surface gravity and metallicity, Astron.Astroph.Suppl.Ser., 109, 313-328, 1995. Rufener, F., The evolution of atmospheric extinction at La Silla, Astron.Astrophys, 165, 275-286, 1986. Schimpf B., F.Schreier, Robust and efficient inversion of vertical sounding atmospheric high-resolution spectra by means of regularization, J.Geophys.Res., 102, 16037-16055, 1997. Schroeder, D., Astronomical Optics, 2 nd ed., Academic Press, 478p., 2000. 199 Sheridan, P.J., J.A.Ogren, Observations of the vertical and regional variability of aerosol optical properties over central and eastern North America, J.Geophys.Res., 104, 16793-16805, 1999. Sheridan, P.J., D.J.Delene, J.A.Ogren, Four years of continuous surface aerosol measurements from the department of energy?s atmospheric radiation measurement program Southern Great Plains Cloud and Radiation Testbed site, J.Geophys.Res., 106 (D18), 20735-20747, 2001. Shettle, Eric P., Models for the aerosols of the lower atmosphere and the effects of humidity variations on their optical properties, AFGL-TR-79-0214, 94 p., 1979. Shields, J.E., M.E.Karr, T.P.Tooman, D.H.Sowle, S.T.Moore, The whole sky imager? a year of progress, Session Paper, on-line at http://www- mpl.ucsd.edu/people/jshields/, 1997. Siske, P.E., Iterative unfolding of intensity data, with application to molecular beam scattering, J.Chem.Phys., 59 (11), 6052-6060, 1973. Slusser, J., K.Hammond, A.Kylling, K.Stamnes, L.Perliski, A.Dahlbach, D.Anderson, R.DeMajistre, Comparison of air mass computations, J.Geophys.Res., 101 (D5), 9315-9321, 1996. Smith, Warren J., Modern Optical Engineering: The design of optical systems, 2 nd ed., McGraw Hill, 524 p., 1990. Staude, H.J., Scattering in the Earth?s atmosphere: Calculations for Milky Way and Zodiacal Light as extended sources, Astron.Astrophys., 39, 325-333, 1975. Sterken, C., J.Manfroid, Astronomical Photometry, Kluwer Academic Publishers, 272 p., 1992a. Sterken, C., J.Manfroid, The evolution of atmospheric extinction at La Silla derived from measurements in the Str?mgren photometric system, Astron.Astrophys., 266, 619-627, 1992b. 200 Stetson, P.B., L.E.Davis, D.R.Crabtree, Future development of the DAOPHOT crowded-field photometry package, Proceedings of the Conference, Tucson, AZ, Sept. 6-8, 1989. San Francisco, CA, Astronomical Society of the Pacific, 289-304, 1990a. Stetson, P., On the growth-curve method for calibrating stellar photometry with CCDs, Publ.Astron.Soc.Pacific, 102, 932-948, 1990b. Stokes, G.M. and S. E. Schwartz, The Atmospheric Radiation Measurement (ARM) program: Programmatic background and design of the cloud and radiation test bed, Bull.Amer.Meteor.Soc.,75, 1201-1221, 1994. Stone, R.C., An accurate method for computing atmospheric refraction, Publ.Astron.Soc. Pacific, 108, 1051-1058, 1996. Strand, O.N., Theory and methods related to the singular-function expansion and Landweber?s iteration for integral equations of the first kind, SIAM J. Numer.Anal., 11(4), 798-825, 1974. Tegen, I., P.Hollrig, M.Chin, I.Fung, D.Jacob, J.Penner, Contribution of different aerosol species to the global aerosol extinction optical thickness: Estimates from model results, J.Geopgys.Res., 102, 23895-23915, 1997. Thomas, G.E., K.Stamnes, Radiative transfer in the atmosphere and ocean, Cambridge University Press, 517p., 1999. Thorne, D.J., P.R.Jorden, N.R.Waltham, I.G.van Breda, Laboratory and astronomical comparisons of RCA, GEC and Thomson CCDs, in Instrumentation in astronomy VI, Proceedings of the Meeting, Tucson, AZ, Mar. 4-8, 1986. Part 2. Bellingham, WA, Society of Photo-Optical Instrumentation Engineers, p.530-542, 1986. Tug, H., N.M.White, G.W.Lockwood, Absolute energy distributions of Alpha Lyrae and 109 Virginis from 3295 to 9040 A, Astron.Astrophys., 61, 679-684, 1977. 201 Turner, D.D., R.A.Ferrare, L.A.Brasseur, Average aerosol extinction and water vapor profiles over the Southern Great Plains, Geophys.Res.Letters, 28(23), 4441- 4444, 2001. Turnrose, B.E., Absolute spectral energy distribution of the night sky at Palomar and Mount Wilson Observatories, Publ.Astron.Soc. Pacific, 86, 545-551, 1974. Twomey, S., On the numerical solution of Fredholm integral equations of the first kind by the inversion of the linear system produced by quadrature, J.Assoc.Comput.Math., 10, 97-101, 1963. Voss, K.J. and G.Zibordi, Radiometric and geometric calibration of a visible spectral electro-optic fisheye camera radiance distribution system, J.Atm.and Oceanic Technol., 6, 652-662, 1989. Walker, G., Astronomical observations, an optical perspective, Cambridge University Press, 347 p., 1987. Wamsteker, W.S, Standard stars and calibration for JHKLM photometry, Astron.Astrophys., 97, 329-333, 1981. Whiteman, D.N., S.H.Melfi, R.A.Ferrare, Raman Lidar system for the measurement of water vapor and aerosols in the Earth?s atmosphere, Appl.Optics, 31, 3068- 3082, 1992. WSI Whole Sky Imager, 2001: http://www.arm.gov/instruments/static/wsi.html Young, A.T., R.M.Genet, L.J.Boyd, W.J.Borucki, G.W.Lockwood, G.W.Henry, D.S.Hall, D.P.Smith, S.L.Baliunas, R.Donahue, D.H.Epand, Precise automatic differential stellar photometry, Publ.Astron.Soc.Pacific, 103, 221-242, 1991. Young, A.T., High precision photometry, in Stellar photometry, IAU Colloquium 136, C.J.Butler and I.Elliot (eds.), 81-92, 1992a. Young, A.T., Improvements to photometry. V.High-order moments in transformation theory, Astron.Astrophys., 257, 366-388, 1992b. 202