Earthshine blog

For a given stack of 100 images, how large is the uncertainty on derived albedo?

We investigate this using one of the ‘good’ stacks of 100 images for the V band. We call the stack ‘good’ because it seems that repeated fits to it with slightly different settings of the fitting algorithm give similar results.

We split the stack into 33 3-image averages, 16 6-image averages, 8 12-image averages, 4 25-image averages, and 2 50-image averages. We also fit the 100-image average several times. We fit each of these multi-image averages and plot the results:
– so the conclusion is that flux constancy varies.
Top left: log-log plot of how the RMSE depends on number of frames in each average. the RMSE is calculated in the areas used for the fitting of the model to the observed image. The fitted line has slope minus a third. Top right: log-log plot of how the estimated uncertainty of the Albedo-fit (given my routine LMFIT) depends on number of frames averaged. The fitted line has slope minus a third. Middle left: log-log plot of how RMSE and Albedo uncertainty depend on each other. The line has slope near 0.9. Middle right: Fitted albedo against number of frames. The error bars are the formal estimates of the albedo fit uncertainty (given by LMFIT). Bottom left: How the standard deviation of the Albedo determinations (in middle right plot) depends on number of frames averaged. To avoid ‘overplotting’ we have offset N slightly by adding small random numbers but still maintain a clusteringa round the actual N value.

We see some interesting things in this evaluation of ‘intra-stack variability’: The uncertainty (both RMSE and Albedo uncertainty) ought to drop as one over the square root of number of frames averaged, if the noise causing the scatter is normally distributed and independent and there were no biases. As it is, the errors drop with a minus 1/3 slope indicating that something is holding back error-averaging. We also see some hint in this in the middle right panel: albedo seems to depend on how many frames we average – the frames were also averaged sequentially (3-averages formed from frames 1,2,3 then 4,5,6 and so on). We can test if the order of averaging is important by doing some bootstrapping and randomizing (TBD) but there is also the possibility that increased averaging allows some subtle biasing factor to emerge. Finally, we see that intrinsic error (i.e. that not due to choice of image stack but due to factors like noise within a given stack) can be brought below 0.2% when averaging just 20 frames. We also see that the formal estimates of albedo uncertainty are not that far from the scatter of the albedo determinations themselves (but might be off by a factor of 2 or so).

So, in summary, we have found that scatter can be brought down to the few-tenths of a percent level by averaging and using our fitting method and observed images. [Wahey! That is actually a major point for the Big Paper! Rejoice!] We find that there is some sign of bias causing a monotonic evolution of the fitted albedo with increased averaging, as well as a decrease in scatter that does not drop as expected for averaged independent noise.

We can test if the order of averaging is important by bootstrapping, and this will be done. It could help rule out effects due to some time-dependence of the properties of the images in the stack.

CONTINUED LATER:

We have now tested some more ideas: We test effect of alignment as well as effect of the order in which the images are selected for averaging inside the 100-image stack.

Upper left: red line and black points show the old result – sequentially selected frames are averaged, without alignment. Green points and line – the images are selected non-sequentially as in ‘bootstrap with replacement’ and are also not aligned. We see that the slopes of the lines are the same – about -1/3 but that the bootstrapped results have about 25% higher error than the sequentially selected.

This gave the idea that ‘something changes’ along the stack which results in increased error. It could be that the Moon drifts slowly across the sky so that mis-alignment is more evident in sums of frames selected far apart as opposed to in sequence. There could also be ‘other time dependent’ things going on – such as various CCD-related effects.

Upper right: to test the drift hypothesis above we next aligned the selected frames when bootstrapping. As we see there is no effect of doing this! So drift is not the problem – or is it? Alignment consists of sub-pixel shifts that may induce variance-changes due to non-conservative interpolation in the image frame. We have considered this in this blog (see discussion and further blog references HERE). Since the green data look the same in upper left and upper right we conclude that alignment (i.e. possible flux non-conservation) cannot be the issue.

Bottom right: This shows that the scatter on bootstrapped albedo determinations are now larger than for sequential access, and that levels between 0.2 and 0.3% can be reached, at best.

So why is it ‘bad’ to pick the images in random order, rather than sequentially?

Speculation department: If it turns out that alignment simply isn’t the issue then what? We have the shutter and we have the bias. In the above we assume that shutter works fine and that the bias is the same in all frames. We test here the shutter, by looking at the total flux of each image in the 100-image stack:
Top panel shows percentage deviation of the total flux in each bias-subtracted frame relative to the average of the whole stack. The standard deviation is 0.46%. Flux is not constant from frame to frame. Why? Second frame shows the maximum flux in each frame after smoothing. The two graphs are very similar. The last frame shows the mean of a dark corner in each image. The graph does not look like the other two – so it is unlikely that strong sky variations or changes in bias caused the changes in total flux.

We guess that the changes in flux is due to either changes in atmospheric transmission or shutter variability.

But how could flux in-constancy cause the problems we see with increased errors when bootstrapping? Apart from a few big dips in the total flux, most points seem to be randomly scattered wrt the neighbours. In that case sequential access as opposed to random access of image sequence would not matter, would it?

Hmmm.

There is the ‘wobble’ of the image as seen in a film we have on this blog – waves of ‘something’ wash across the lunar disc causing small transmission changes. Sequential images may then be more alike than images selected far apart, and alignment, being whole-image shift based, would not repair the problem due to ‘wobble’. A test suggests itself – look at fits to single images in a stack and the scatter these have. Then consider whether that is more or less than the scatter seen between n-frame averages of the two types.

CONTINUED EVEN LATER:

It was tested whether albedo fit scatter increased or decreased when integer pixel shifts and sub-pixel shifts were used: It turns out that sub_pixel shifts is better than integer pixel shifts in most cases, but that a straight sum (i.e. un-aligned) is best! NB: This is the case for the present choice of 100-image stacks. YMMD.

Below, we have investiagted the phenomenon of ‘Mickey Mouse ears’ – i.e. the tendency for theoretical lunar images to have too bright cusps so that the residual image ends up with a bipolar structure ressembling Mickey Mouse’s ears. We have investigated what happens to the ears if a semi-empirical BS model image is used, based on scaling and clip-selecting the BS from the observed image. We showed that the ears went away.

Now we investigate what happens to the quality of the fit – for our goal is to derive terrestrial albedos with minimal error (both scatter and bias; bias we can live with as long as it is constant; scatter is a pain).

We selected 10 images that we fit with both semi-empirical BSs and theoretical BSs and tabulated the values for fit-parameters as well as some of their formal fit errors (output from LMFIT).

In this plot we show the changes, in percent, of some key parametrs and their errors, in the sense that the purely theoretical result is subtracted from the semi-empirical result:

Change, in %, of some key fitting parameters. (semi-empirical – theoretical)/(½*(semi+theor))*100.

A: albedo, is larger by several percent in the semi-empirical (S) fit, than in the purely theoretical (T) fit.

dA: error on albedo reported by the fitting routine, is mainly smaller in S compared to T.

alfa: the power to which the PSF is raised is less by several % in S compared to T. This implies a broader PSF in the final fit in S.

ped: the pedestal added to the whole image has both positive and negative changes.

dX: the shift of the model image in X direction needed to align model and observed image. Mainly smaller in S than in T.

cf: core factor on the PSF – has larger values in S than in T. So the PSF peak is ‘pulled up’ in S relative to the wings.

contrast: lunar albedo map contrast, has smaller values in S than in T. This is consistent with the ‘pulled up’ core which ‘copies’ detail better.

RMSE: the root mean square error in the area of the image used for minimization, is mainly smaller in S than in T.

So, what have we learned? The two formal estimates of errors – dA and RMSE are both smaller in S than in T, so S would appear to be a better way to go.

Taken together – the smaller errors (on A, and RMSE itself) in S compared to T is encouraging, but we really do not know how much closer to the ‘true value’ for A we are, just that we are probably closer – at least in terms of scatter, not necessarily in terms of bias.

Using the semi-empirical method compared to the fully synthetic method adds very little in terms of computational overhead, so we will use S henceforth.

We note that the width of the distribution of albedos in units of the mean is 8 percent smaller in S than in T. The S method gives tighter distributions of albedo. Of course, the intrinsic (geophysical) scatter in albedo is present in this number too.

It seems our synethtic model images of the Moon have a distribution across the lunar disc that does not match observations. When we fit such an image to the ‘band’ near the photo-equator and subtract the fitted image from the observed image we get a residuals image with ‘Mickey Mouse ears’ – i.e. the residuals in intensity are bi-polar.

We are testing a semi-empirical method where the BS of the synthetic image is replaced by the observed BS, scaled to the right value. When such an image is fitted to the observed image we get no Mickey Mouse ears problem:

Left image: Residuals after a well-fitting synthetic image has been subtracted from the observed. Note the polar nature of the residuals – the cusps are clearly too bright while the area between the cusps is too dark. Right image: residuals when a well-fitting image, based on the observed BS, is subtracted. There is very little large-scale structure now.

The two fits deliver terrestrial albedo as one of the fitting parameters.
For the left image (i.e. the entirely synethtic model image) we get A=0.37298, for the semi-empirical image we get A=0.36311, that is a difference of 2.7 %, which is considerable, given our science goals of 0.1%.

Whether the fit has been improved where it counts is a different matter. The RMSE for the two fits are identical (RMSE=0.026).

The large change in fitted albedo is a (regrettable) feature of our fitting method – change something a little bit and the fits change by several percent.

We suggest that a meaningful test of fit quality will be to fit the ensemble of images with the two methods above and see if the scatter is less in one than in the other. TBD.

In our fits of model-images to observed images we fit a model with several parameters. One of them is a ‘sky pedestal’ which is a constant that is added to the model image before fitting in order to accomodate sky brightness that is not due to the scattered-light halo.

We include this term because we figure that the sky might be uniformly bright for several reasons:

1) moonlight could be reflected from the ground up to the sky and back into the telescope
2) airglow could brighten the sky
3) unresolved stars and zodiacal lights shine from the sky behind the Moon
4) light could be scattered inside the telescope and contribute a more or less uniform ‘haze’

We show here the magnitude plot of that pedestal term against lunar phase for our 5 filters:

Error: axis labels have been switched!

In order to interpret this plot you should know that the pedestal term is ADDED to the model image, so a positive value of the pedestal implies that the model is representing a sky brightness, and that a negative value of the pedestal implies that the model needs to subtract something before it can fit the observation well.

Notice that the pedestal term is negative for all but the VE2 filter. That is, only in VE2 is a well fitting model one that has a sky brightness contribution. For the remaning filters the model needs a subtracted term for the fit to be good.

Notice also that the VE2 pedestal is about 4-6 times larger than the other filter’s, in absolute terms.

Notice also that the term goes to zero for large phases (i.e. near New Moon).

So – turning to speculation now – what is going on?

Only VE2 behaves as expected: As the Moon gets brighter more light has to be removed. This could be Moonlight scattered from the ground, onto the sky, and back to the telescope [See the paper by Bernstein, R. a. R. ∼A. “The Optical Extragalactic Background Light: Revisions and Further Comments.” The Astrophysical Journal 666, 663–673 (2007). They model the moonlight reflected from the ground near a telescope.]. Or light scattered inside the telescope – but not airglow (why would it depend on lunar phase?). Nor can it be zodiacal light because although it depends on lunar phase when you observe the Moon in particular it has the wrong phase-dependence – towards Full Moon the ZL near the Moon is low; towards New Moon the ZL near the Moon is stronger (because the Moon is nearer to the Sun and the ZL is strongest near the Moon).

So – the pedestal for VE2 behaves as we expect widely scattered light from the Moon should behave.

By the way, let us call this light ‘ADL’ from now on – ambient diffuse light, to tell it apart from the light scattered in the halo, which has a strong profile with distance from the Moon.

But for the other filters the pedestal behaviour is opposite!

I think this could be a sign of the fitting process having a problem – the method has to specify negative pedestals for something else to work. That could perhaps also be the case for VE2 but to a smaller extent than the conventionally expected scattering: if it is bigger for VE2 a contribution could still be taken out leaving a positive pedestal for that filter.

Do we have any reasons to believe in stronger ADL for the VE2 wavelengths than the rest? The VE2 filter is of a different nature than the rest – perhaps it scatters more light? This would match the phase-behaviour.

What could the role of airglow be? While the zodiacal light and the galactic background and stars come from behind the Moon and are blocked by this, the airglow is all over the image. What happens when we fit such images with just one pedestal term? Can the procedure go wrong in the observed way?

What do we know about airglow? Does it depend on wavelength?

WORKING HYPOTHESIS:

Let us assume that: In the VE2 filter there is a lot of scattering, and the method of fitting deals with this by specifying a large positive pedestal. For reasons not yet understood, there is a tendency to make the pedestal a little too small – but this is not seen in VE2 because the scattering is so large. In the other filters the scattering is less and the error comes through as a small negative pedestal.

PLAN:

Work on understanding why the pedestal needs to be a small negative number for the other filters. Perhaps it is a consequence of some of the choices of frozen parameters in the model. We do freeze ‘core factor’ and ‘rlimit’, for instance. Optimally, the factors that need to be frozen should be done so at values that do not induce strange problems such as the present one.

Secondarily, find out if the small negative pedestal influences albedo determinations significantly. With the DS intensity being proportional to albedo we can see that the small negative pedestal values directly correspond to albedo biases in the 0.2% range.

If you fit each image in a 100-image stack with our model, and plot the histogram of albedos, you get this:

The dashed line shows the fitted V-band albedo of the 101st image, which is the average image made up of the aligned and averaged 100 images in the stack.

In fitting the above 100 images we found residuals that were ‘wavy’ – probably the waves of ‘something’ that is seen in the movie here. The residuals look like this:

If you fit the 100-image average image, you get this:

The residuals are much smaller – but also ‘wiggly’. The RMSE of the 100-image average fit is about 3.7 times smaller than the single-image fit’s. Can this, and the width of the histogram (1.8%), be used to form a scaling argument so that the uncertainty on the 100-image fit can be deduced?

I’d like to think that

error_fit_100_average = SD_histogram*(RMSE_100/RMSE_1)

which would give us an uncertainty on the fit of the 100-image average of 0.5%whichis quite a lot. In contrast, the standard deviation of the mean, SD_m=SD/sqrt(N-1), is 1/10th of the 1.8%, or 0.18%. This is more in keeping with the experience we have had with fitting ideal images with added noise.

The advantage of the fit of the 100 individual images is that you get the histogram and thus a measure of the actual uncertainty – when you fit the 100-average image you just get the fit – formal uncertainty estimates depend on knowing (independently) the noise on each pixel. One way around is to consider Monte Carlo Markov Chain fitting using the Metropolis-Hastings algorithm – this iterates the fit in a Bayesian framework and generates histograms of all fitted parameters, but is very very slow because MANY iterations are needed. Another way would be to bootstrap the single images, somehow.

—–

Here are the similar histograms for the same night from V, B and VE2. Again, the dashed line is the albedo found in the 100-image average image, while the histogram is made up of 100 single-image determinations:

We now have 3 stacks for which the mean of the histogram and the albedo of the mean image are very close – this encourages us to think that the uncertainty of the determination based on the mean image is similar to the standard deviation of the mean for the histogram. Let us test this by fitting even more stacks and their component images along with the mean image. It would be a way towards using the mean-image determinations as the ‘real thing’ and designate the variations we see from stack to stack as geophysical and not based in analysis method bias.

In this entry below we suggested 4 things to try to improve the accuracy and precision of our fitting method. A fifth idea is that byt fitting small, sharp-edged regions on the DS we are allowing small shifts in the image to ‘pull’ thge fit. A similar comment applies to the effects of atmsopheric turbulence. We have elsewhere shown that the lunar disc ‘wiggles’ when you look at an animation of 100 images. Some regions (5-10 arcminutes in diameter) shift sideways by what looks like a large fraction of a pixel, relative to other parts. If we use small patches to fit on quality will suffer if we fit across a wiggly patch’s edges.

Using larger patches might therefore be an idea, and we test that here, while keeping all other thinsg fixed. That is, we use the exact same pixels as before but instead of concatenating 11 row-averaged bands we average the 11 bands into one profile. Again, we use the set of 10 ideal images to which realistic noise had been added.

We get:

mean albedo: 0.28031400
std. dev. of 10 fits:0.00021
SD expressed as percentage of mean: 0.08 %

Compare this to what we found before, using concatenated strips:

mean albedo: 0.28051
std. dev. of 10 fits: 0.00030
SD expressed as percentage of mean: 0.11%

It would seem we have about 30% less spread in the values found. As before the known fixed albedo of this ideal image was near 0.2808, so we have moved slightly away from this, and the bias is now about 2 SDs.

The use of averaging the 11 bands into one (they follow the contour of the lunar disc) does seems to ensure that flux is not lost outside a narrow band. We still use radial bins (5 pixels wide) so flux can be moved from bin to bin that way. We hesitate to make larger bins, as the degrees of freedom (surely?) would suffer. We have about 28 radial bins.

So, we have less scatter but more bias! Not sure this is a way it’s worth to go. Alsao recall that thgese are ideal images – we are not sure what happens with real images. One day we might fit all images using the concatenated and averaged versions of the fitting method and see which has the smaller scatter.

We jave more suggestions from last blog entry to try:

1) Use more sky pixels
2) Use more DS pixels
3) Use BS edge
4) Average several images

4. We can always do and is independent of other method choices. 3 is scary, so I will try 1 and 2 next.

We conduct two experiments, to asses the abilities of the fit-a-model-image method.

Experiment 1: Same image, with differing noise.
=============================

We pick an ideal image, generated from JD2455945, and to this one image we apply a PSF convolution, and 10 times generate images with varying Poisson noise. We then fit these 10 images, and look at the distribution of results.

We find

mean albedo: 0.28051
std. dev. of 10 fits: 0.00030
SD expressed as percentage of mean: 0.11%

The actual value for albedo used to generate the ideal image was 0.28083594, which is 1 SD above the mean value found.

So, we can (again) confirm that in principle, we can achieve fits with small bias and small errors, as per our original science goal.

Experiment 2: Same image with noise, different fit starting guesses.
=========================================

We pick now one of the noisy artificial images and refit it 10 times – each tiome starting at a random guess near the supposed best guess. We pick our vatiation froma normal distribution with amplitude 0.03 and multiply the fixed starting guesses by 1+0.03*N(0,1) for the 8 starting parameters (there are 8 – 7 before – we added shift om image along y-axis also). These guesses were then allowed to converge and we collected the results. We found:

mean albedo=0.280183
S.D.=1.4181365e-05
SD in % of mean=0.0051

We see that the sactter is very small, i.e. that convergence is almost always to the same value. However, there is a bias in that the known albedo (see above) was not found.

The bias is (0.280183-0.28083594)=0.0006529 which is 0.23% of the correct value.

From this we lkearn that
a) convergence may not always be dependent on starting point – the same best solution will be arrived at, and that
b) the noise added to the image itself will have enough influence, even if small (i.e. realistic for averaging 100 images in a stack, and coadding 20 ‘lines’ in each of 11 ‘bands’) to cause a bias in this perfect-case example, where noise is only Poissonian.

The bias is a mere 0.2% of the correct value which we could lkive with given our science goal of 0.1, I guess, but reality adds more problems to our images than just Poisson noise. These effects appear to give us a typical spread of 1% – about 10 times what we saw in experiment 1, above.

As it is, we are using 1/10 of the pixels in the image for our fits so we have at most a factor of 10 more pixels to work on. This could in principle give sus a factor f 3 reduction in the uncertainty, if all pixels are equally effective at driving the solution. If all went well, our 1% errors coul dbecome 0.3% errors which woul dbe acceptable (just), in the absence of a bias.

What can we do along this line?

1) Many of the DS sky pixels are not being put to use – they coul dhelp set the sky-level of the fitted model better.
2) Not all DS disc pixels are being put to use, but there is a limit how far onto the disc we can go before BS halo starts being a too importnat factor, but perhaps we could double the number of pixels.
3) We only use the BS disc and BS sky to calculate the total flux which is used to ensure same flux in model and observation. We might be able to also fit some of the edge near the BS or just try to match the level of the BS sky near the disc – but the levels are much higher and our cirrent use of errors based on differences between absolute levels would become dominated by the BS pixels in driving the fit.
4) We can use more images from the same night – at least inside some period of time where the geometry of Earth and Moon does not change significantly – perhaps half an hour or so. In that time some 10-30 images could be gathered in reality. If several filters are wanted then perhaps 10 images in each filter (in practise it seems we could reach that on some nights, inside a half hour) and then repeat for other filters.

In summary, if erros were due to Poisson noise onkly, we could hope to double the number of DS disc pixels, use a lot more of the DS sky pixels, and perhaps work towards 10 images in each filter inside each half hour, giving us, potentially, the required improvements in uncertainty.

We can work on this!

There is a “however”, however: All the noise is not Poissonian – there are issues having to do with:

a) Bias subtraction – and airglow – and zodical lights: letting the ‘pedestal’ solve for these may be a bad idea – perhaps we should fix the pedestal at predicted levels in each image? This needs some wor, since we are sensitive to errors there too!
b) The atmospheric turbulence causes ‘wiggle’ in the images,. described on this blog, and this is a factor in driving the error in the fits.
c) The model we fit is not faultless either – the BDRF is a simple one and we see hints of some azimuthal-angle problem ain fitting along the edge of the disc. Model couild be improved.
d) Other parameters in the model might be fixed too – such as the ‘core factor’. Experiments are still needed to see if this is a good idea.

After exhaustive (-ing?) tests of fitting the same image over and over again under varying conditions, we now turn to fitting all images in two colours from one night. We take B and V images from JD2456045.

We have before fitted ‘lines’ going from DS sky onto the DS disc itself, and then moved on to another position along the DS edge. Thsi helped us see that fits along the edge could vary systematically, and we have speculated that it was due to the BRDF-model used had some azimuthal angle-dependence that was different from reality.

Instead of getting several fits along the edge and combining these, we now choose to combine lines and make one fit, thus allowing the errors to balance out somewhat.

A typical fit, for 11 lines, looks like this:

The upper panel shows 11 concatenated ‘lines’ going from 50 pixels off the DS edge to 70 pixels onto the disc. The white curve is observations and the red is the model. The lower panel shows the difference between the two at ‘best fit’.

We have 8 parameters in the model that together describe such things as the terrestrial albedo (the main goal of the exercise) and parameters that adjust the PSF, as well as the relative alignment of model and observed image.

In the above, the residuals are small, buts how systematic behaviour near the ends – going along the sequence of the zig-zags we are essentially sweeping along the edge of the DS – so something is not well at extremes of the sequence, but near the middle the residuals are very small.

We smooth the concatenated lines (and corresponding model line) before fitting – this helps reduce a problem we had before, namely, that repeated fits of the same image starting at slightly diffeent (1% different) starting guesses could finish up at slightly different ‘best fits’.

For all images in B and V for one night we can compatre the various quantities determined in the fitting:

From top to bottom, left to right: Albedos in B and V images against time of day (not really Julian Day, but fraction of day since that JD started); alfa of PSF wings against time; PSF core radius (here held fixed at 3 pixels); pedestal against time; image shifts in X and Y; lunar albedo contrast factor; PSF core factor; and RMSE against time. Black crosses are B-band, while green crosses are for the V-band.

Here we see a clear difference in B and V albedo, with some scatter. The PSF is behaving diffrently in B and V (core factor, pedestal and wing alfa). Lunar albedo amp contrast factor is unaffected by B and V. RMSE is bigger for B than for V. Image shifts are near 1 pixel in each idrection.

The above choice of a fixed core radius was prompted by noting that leaving it free to vary inside the range 2-10 would cause the fit to pin the core radius to the lower limit in V and the uppe rlimit in B. We argued that some of the factors core radius/core factor/ andpedestal were degenerate and chose to fix the core radius at a compromise value. We note the extreme values of the core factor parameter for B, compared to V. In essence, a large core factor implies that most of the flux is in the core of the PSF and little is in the wings. It is hard to say if this is what we expect for B vs V – the colour-dependence of light-scattering in optics (and atmosphere) is not understood by us. Yet.

The observations in B and V were not obtained at the same times (because we use filter-wheels) so it is not easy to combine the B and V albedo determinations. Some of the scatter seems to be correlated (first points in left top panel above), but the rest do not. The means and standard deviations and standard deviations of the means, and the SD in percent of the mean, for B and V are:

B: 0.3608 0.005 0.0012 1.3%
V: 0.2902 0.002 0.0006 0.8%
mean SD SD_m SD/mean*100

So we have scatter on albedo determinations of near 1%, like BBSO. Formal errors on the determination of albedo for a single image depend on being able to estimate the uncertainty of the data. As we clearly have some systematics going on (plot at top of this entry) this is not possible.

For the night JD 2455945 (the subject of our paper) we have fewer good images, but get is:

B: 0.37007500 0.00020615528 0.00011902381 0.055706352
V: 0.28093333 0.00023094011 0.00016329932 0.082204595
mean SD SD_m SD/mean*100

The B and V here are ALBEDOS, not magnitudes. Converted to magnitudes the Balbed/Valbedo ratios on the two nights have Johnson B-V colours that differ from one night to the other by 0.06 mags.

Thus we are told that one night’s earthlight was 0.06 mags bluer than the other. This is several times our estimated observing precision.

Is that a lot or a little? Our EGU 2014 poster tells us which …

We have been looking at the ‘pedestal’ which is a part of the model we fit to our observations. With a black sky and perfect bias and dark-current removal, the pedestal should be 0, and the fitted model should have a halo that sloped off from the Moon and perfectly explained all brightness seen on the sky.

In reality, there is background galactic light, zodiacal light, airglow and any terrestrial light (Moonshine reflected from the ground!!), to allow for. To this end we have operated with this pedestal in the model, which was a vertical offset of the model. In looking at fits we notice that the pedestal is degenerate with the halo-parameters in such a way that they tend to lie along a line when a lot of fits are compared, each started at slightly different points – even for the same image repeatedly fit. The use of this pedestal therefore leeds to a bias in the determined albedo. Which is not good.

We therefore have tried to evaluate expressions for the additional sky brightness so that the pedestal could be fixed at an assigned, physics-based, value and not interfere with the fit itself.

In doing that we had reason to revisit and elaborate on Henriette’s error budget analysis. In particular, we performed the analysis allowing for error terms in DS and BS fluxes, dark currents, and exposure times. We introduced the uncertainty on the superbias factor and see now that it is the dominating factor when area-means are considered.

For area-means of 20×20 points, the error due to the superbias factor is as large as the error due to Poisson statistics on the DS. For smaller areas than 20×20 the DS Poisson error dominates.

Since levels of sky brightness are of the same order as the contribution due to factor uncertainty, or larger, we see that we have a problem.

We could try to estimate sky brightness better, but as we realize that fitted albedo varies, for a single image treated in different data-reduction ways (bias removal, and stacking), we need to also re-consider the use of ‘anchoring’ of the model halo in the DS sky part of the image. This is what BBSO achieve by fitting a line to the DS sky and subtracting it from the image.

Improvements in estimating sky brightness could also be pursued, but all constant contributions will disappear in ‘anchored’ procedures. The point is to not let the fit be biased by such terms.

We continue the investigation into how well albedo can be determined by looking at justa single image and applying various methods to it.

We compare the ‘old’ method to a new method. The old one consists of plotting row-averaged slices that start on the DS sky and proceed onto the DS disc. We place the slice at 9 different locations up and down the edge of the lunar disc. The ‘new’ method consists of averaging the whole image (and fitted model) into one row by averaging along columns. We do two runs of the old method – varying the distance we go onto the DS, and one run of the new.

Black crosses indicate the albedo found by fitting along 9 slices ono the DS disc going from 50 columns away from the edge to 100 columns beyond that. Red crosses indicate the same but going to 50 columns onto the disc. Green crosses indicate results found when column averaging the whole image but starting from 9 different starting guesses of the fitting parameters. The dashed and dotted lines give means and standard deviation of the mean for all trials.

We see that the three methods give a spread of values, with few repetitions. The standard deviation of the Mean is about 0.25% of the mean, but the standard deviation of the distribution is about 1% – single determinations will therefore vary a lot, while the mean over many trials can give a mean with small error.

We notice that the green crosses – which are results of identical trials in the ‘collapse the whole image into one row’ method but starting from slightly different starting guesses – all lie below the mean and have some scatter – so the new method is biased wrt the old one and the methdo itself is finding ‘best fits’ that depends on starting conditions. That the convergence is not to the same number each time show that the LMFIT method used reaches a stopping criterion. We have chosen max 500 iterations.

Neither of the above findings are reassuring – fitting method choices determine the outcome. Possibly, averaging over a large number ffo different methods would give us a result that is stable and near the truth.

In summary, we are looking at method-dependent biases of half and whole percentages, and errors of means near 0.2%.

It is not clear if the bias comes from lunar albedo-map and reflectance-model choices, or are stochastic in nature. THe effect of stochasticity is seen from the green crosses where all things are fixed, except the starting guess for the iterated model-fit.

We explore how well albedo can be determined from just a single image, varying the fit approaches.

A single 100-image stack of the Moon was selected, aligned (or not, see below), bias-reduced and averaged and then subjected to our model fit. The fit selects 9 21-row wide ‘strips’ cutting across the DS edge of the disc so that we have 50 columns of sky and then 100 columns of disc to fit on. The 21 rows are averaged in observation as well as model.

Applying this method we get 9 estimates of the terrestrial albedo. We selected three methods for centering stack-images:

1) Align the images in the stack once allowing for sub-pixel shifts, then average
2) Align iteratively, refining the alignment wrt an updated reference image, then average
3) Do not align at all, just average

In this figure we show the albedo determined from the three sets:
The determinations are shown in sequence, separated by vertical dot-dashed lines – method 1,3 and 2 (in that order). The overall mean and standard deviation is shown

We see systematic deviations: all Method 1 values are below the mean, all method 3 values are above the mean and method 2 is a bit above and a bit below.

The mean of these values is 0.3047 with a standard deviation of the mean of 0.0014 (0.45% of the mean), which is not too bad! However, we could clearly do better if the method-dependency was not present. How to remove this?

Methods that shift images by non-integer-pixel amounts risk being non-conservative – i.e. the total and area-fluxes are not conserved. We have elsewhere in this blog estimated how bad the problem can be [see http://iloapp.thejll.com/blog/earthshine?Home&post=369 also see http://iloapp.thejll.com/blog/earthshine?Home&post=299]. We found that it was slight, means can only change by much less than 0.1% during a sub-pixel shift of the image.

Not aligning images at all risks comparing ‘smeared’ observations with ‘knife-sharp model images. The smearing occurs due to image-wander during exposure and can be as much as a pixel or two.

Revision of our alignment methods seem in order, given the above. Possibly shifts by whole pixels only? This is at least conservative away from edges.

Note that method 3 (the middle set) has the least internal scatter. For this set the standard deviation of the mean is just 0.2% of the mean. That is pretty strange, actually, since this is the ‘do not shift images before coadding’ method. Hmm. Then again, perhaps it is not strange, and we are simply seeing that the likely most fuzzy image has the most stable fit – as the strips select different parts of the light-and-dark lunar surface to fit on, contrast causes the fit to wander if the image is sharp. But on the third hand, our method now allows for ‘contrast scaling’. Hmmm, mark 2.

Shifting and co-adding involves two issues – the interpolation required at a given non-integer pixel shift may not be flux-conservative (our tests seem to show this is only a small problem). The shift itself may be poorly estimated. Our shifting method is based on correlation – it does not look specifically at such things as ‘are the edges still sharp’ – i.e. Sally Langfords Laplacian edge-detection method. The use of edge-information only is not clearly the best method since it is based on very few image pixels and thus suffers from noise. Perhaps a hybrid method can be envisaged – “correlations and maintaining-sharp-edges”?

Since the start of this project we have worked on a method that fits suitably-convolved amd scaled theoretical images of the Moon to real observed images of same. Below we show that the method has improved, by fitting a profile from the same night with two methods from 2013 and now.

This is a pdf file showing the three plots:

As you can see, the present method leaves residuals that are structure-free – this is not the case with the older methods.

We are therefore approaching a situation where all data can be fit pretty much perfectly by a multi-parameter model.

The model uses these fitting parameters: a vertical offset representing sky brightness, a shift of the image along rows, a parameter for the width of the PSF core, a parameter that vertically scales the core inside that width, a parameter that sets the slope of the wings, a factor that scales the lunar alnbedo map contrast, and the terrestrial albedo.

Since the results still show signs of a dependency on just which part of the Moon we are fitting on, we conclude that the BRDF functions we use – Hapke 63 – are not capturing the right azimuthal-angle dependence of the reflectance. This area is now where improvements are to be found.

If we succeed we shall actually be contributing to Hapke-type work, by providing empirically validated BRDF functions.

At MLO the aerosol load of the atmosphere is measured and a database of hourly, daily and monthly average values for the absorption cofficient, in three optical bands, exists. Here is the link:

http://www.esrl.noaa.gov/gmd/dv/data/index.php?site=mlo&category=Aerosols&frequency=Hourly%2BAverages

The data cover the period in 2011+2012 where we observed..

Our fitting procedures have produced estimates of the power a PSF has to be raised to in order for the convolution of model images with the resulting PSF to produce a ‘good fit’ at the edge of the DS disc in lunar images.

We compare our ‘alfa’ values with AERONET absorption coefficients, and get this plot:

We used daily average AERONET values, which goes some way towards explaining why the scatter for our data (the ‘Power’ alfa) is large while AERONET seems less messy. The overplotted red and green lines are two different robust regression fits. The slope is negative implying that large absorbtion coefficients correspond to small ‘alfa’ values. This matches how we empirically undersatnd the behaviour – on bad foggy nights (presumably with lots of aerosols in the air) we got broad PSFs with small values of alfa.

Repeating the above with hourly-average data the situation does not improve,

So our first interpretation of the above is that ‘alfa’ is only a poor descriptor of the aerosol load – ‘alfa’ is really just a dummy fitting parameter for us – alfa is the value that causes a good fit at the edge of the DS disc. It describes the shape of the PSF. Apparently thisis only weakly coupled to the total amount of absorption.

We must understand the AERONET data better first.

In our data-reduction methods, image-alignment is often involved. It is based on interpolation. An important questions i whether the mean is conserved during such alignment?

To test this we iterate the shift of an observed image, measuring the mean of a patch given by selenographic coordinates. We iterate many times – i.e. reuse the previously shifted image for next round – in order to get a clearer answer. Visual inspection shows the iterated image getting ‘fuzzy’ as the shifts recur. That is, image standard deviation suffers. How bad is the situation for the mean of an area?

The answer is ‘not very much’. Iterating the shifts 40 times we find thatthe per-shift loss in mean value over a realistically sized patch on the DS is -0.005 %. The mean in the patch was about 10 counts which is ‘bright’ for the DS. At other lunar phases the counts may be just 2 or 3 – still just a factor of 5 from the above.

We would never shift the same image 40 times, of course, but the iterated shifts above allows us to build up a change in the mean that can be fitted with a regression. We have learned that one shift of an image may ‘cost’ 0.005% of the mean flux. This is a lot less than the goals for accuracy we have set, which are ‘0.1%’.

Added later: an identical test, but using ROT to rotate an image bya random amount, iterated, showed that errors per rotation were almost always less then 0.01%.

We are interested in just how good the method we use is when it comes to measuring B-V in small areas on the DS and the BS of observed images. In this entry we could find an error as low as 0.005 (once we stopped measuring incorrectly). Is this thelower limit?

We now use synthetic images with realistic levels of Poisson noise added to get the answer. We use ‘syntheticmoon_new.f’ from Chris Flynn and input the ideal image from JD2455814 – i.e. the almost Full Moon – and coadd 100 images, each with realistic noise added, and folded with a realistic PSF (alfa=1.7). We do this twice, with different noise seeds – convert to instrumental magnitudes, and then subtract the two images, pretending they are B and V images and call the result our ‘synthetic B-V image’.

The resulting image is then measured in the areas designated by the agreed selenographic coordinates – those used in our little paper in Table 1. The errors found are:

0.0004 in the synthetic image
0.001 in the case of the real observed images

So, something is generating extra noise for us in the observed images, increasing the noise by more than a factor of two. This could be due to many things – real images are bias subtracted and flat fielded – synthetic are not; there is no readout noise in the synthetic images; image registration does not have to be performed (but could be simulated) in the synthetic images and we have seen variance do strange things during the necessary interpolations that take place during image alignment. There may be more issues to think of here.

This was for a situation with full illumination on both patches. What happens when one is in Earthshine only?

We will simulate this, but, for now, assuming the DS is 1000 fainter than the BS we would expect the error level to be sqrt(1000) larger in a DS-BS image, or about 0.012 mags. This is close to what we see (0.015).

Added later: Ona synthetic image of our lucky night (JD2455945) we find that the lower error limit on B-V BS-DS is

0.0021
and, for the observed image (see this entry)
0.005

So, again, we have a factor of two or so more nois ein reality than the best-case synthetic world predicts.

There is potential for improving our technique! For now we should report the above in the little paper. Fine tuning bias subtraction, flat fielding and image alignment – and understanding the role of image ‘wigglyness’ in alignment issues is needed.

Our observed images are a convolution of a point-spread function (PSF) and the Object (the Moon). We know from analysis of images of various sources – the Moon itself, Jupiter, and stars, that the PSF has the character of a power law along the lines of a core with wings that follow 1/r^alfa with alfa near 2.9.

It is in principle also possible to determine the PSF by de-convolution. Clasically we would expect that

Image = Object * PSF

Where the ‘*’ implies spatial convolution. Taking FT of both sides we get

F(I) = F(O) x F(P)

where F is the forward Fourier transform, and the ‘x’ implies multiplication. Rearranging and taking inverse Fourier transforms (f) we get:

P = f(F(P)) = f(F(I)/F(O))

We have I (the image) and if we had O we could calculate the PSF P. Usually P is very noisy because of amplification of small noisy signals at high frequencies due to the division above. We can average over several estimates of P, however.

For us, this is easy because our data come in the form of 100-image stacks. Unfortunately we do not know O, the object. But we can make synthetic images of O that are highly realistic!

We have applied the above procedure on an image-stack of the Moon, for which we also have a synthetic image. We average the 100 PSFs that are generated. The radial profile from the peak looks like this (black curve):

The blue line is a power law with alfa=2.1 and the red line is a power law with slope 0.3. For this particular night we know that a single power law (plus a narrow table-lookup core) of alfa=2.6 has been estimated for this image using forward modelling techniques where we fit parts of the scattered halo at about 100-150 pixels distance from the photo-centre of the bright side of the Moon.

The two estimates of the PSF are thus quite different, and I wonder why.

Possibly it is a sign that the edge-fitting can be performed well with a large family of PSFs. We know that the fit is excellent at the edge of the disc.

I suggest that a next step could be to test both PSFs or types of PSFs and see how they perform in various situations. How close to the observed image is the ideal image once it is convolved with the above profile?

Following up on this entry we test the Force method on ideal images. This should enable an understanding of the methods actual strengths and weaknesses.

We take an ideal image and fold it with two PSF’s – one with alfa=1.74 and one with 1.76. The 1.76 image is the ‘sharper’ of the two, so we now subject this image to a test where it is folded with various PSFs so that it becomes ‘less sharp’, subtract it from the 1.74 image and look at the difference, expressing it as a percentage of the known intensities (known from the ideal image it is all based on).

Sharper pdf here:

In the upper panel we see the cross-section of the ideal image (black line) and the 1.74 image (red) and the 1.76 image (blue). In the lower panel we see the percentage differences. Black curve is for the difference between the red and blue curves (errors on the DS larger tan 1% everywhere, growing fast as we approach the BS), while the green and red curves show the percentage difference between the 1.74 image and a sequence of images formed by folding the sharper, 1.76, image with PSFs using alfa=1.98, 1.96, … , 1.88. (top green curve has alfa=1.98, next down is red and has alfa=1.96, etc) The difference can, evidently, be less than 1% (for alfa=1.90) but not all over the DS at the same time – there is a tendency for the error to dive faster near the BS and in fact become negative. Near the DS edge to the sky it is possible to have errors in the 0.1-1% interval.

This test is realistic in that the PSF used is developed from real images; in that the procedure of modifying one PSF to another through the use of raising a canonical PSF to a power (alfa) is theoretically underpinned; and in that the brightness levels used in the ideal image are realistic.

The test is unrealistic in assuming that an image with one PSF can be transformed into an image with another PSF by convolution; the procedure of generating PSFs from a canonical PSF through exponentiation is probably only good in the wings of the PSFs.

The test shows that the method is not foolproof. Possibly it is better to first deconvolve an image and then forward convolve it in order to match PSFs – but this procedure needs to be tested on its own.

The test seems to show that for areas near the DS/sky rim a halo can be removed so that errors are in the few tenths of a percent range.

‘

After some deliberations we now return to the ‘edge fitting method’ and its results. There are 19 nights on which more than 17 B and V pairs of closely spaced images exist. We consider those B and V image pairs with negative lunar phases, and look at how the B minus V (albedos, not colours!) values are distributed against lunar phase. We see this:

pdf here:

The error bars are estimated on the basis of error propagation in the B and V albedos from the edge-fitting uncertainty procedure. We have added small offsets in phase to get separation of points. The above actually contains data from 19 nights. To better see things we calculate the nightly median value Of B and V albedo, and now plot these:

pdf here:

Here, error bars are given by the standard deviation of the nightly values. Each point is labelled with the JD number. We see three night-after-night sequences: 2456073-76, 2456016-17, and 2456045-47. Given the error bars, these sequnces follow the same pattern – a rise in B minus V (albedo) as you go from new moon (left) to half moon (-90 degrees). Since the three sequences painmt the same pattern we are confindet in saying that ‘Something Is Going On Here!’.

At the moment we do not know if this is geophysics in the shape of Earth albedo changing with phase – or some aspect of the halo (which increases from left to right above) being harder to model as phase grows.

Next it may be appropriate to study what the Earth actually looked like in the above three sequences. The negative phases selected for the above plots all correspond to the Sun illuminating Earth from Western Pacific/Australia/East ASia and westwards. As time passes on a given day the Moon sinks further in the West (as seen from Hawaii) and more of Asia contributes to the earthshine. On days in a sequence, at the same local hour, less of Earth is illuminated as seen from the Moon so the contribution to earthshine drops but is more and more ‘sicle-like’. SUmmarizing:

For the negative phases selcted we expect:

a) during a single day – contribution to earthshine by continental areas increases, ocea contributions decrease,
b) same time of day but consequitive days – earthshine contribution comes from areas closer and closer to the edge of Earths disc.

I think this means that – if clouds are evenly distributed – we should see reddening of the albedos during a sequence taken on one day, and a blueing in a series of days.

We see b)! Have we seen a)? Below is a plot of the B albedo minus V albedo values for each of the days in question plotted against JD day fraction.

sharper pdf here:

We see little evidence of any up- or down-turns in this. Perhaps 2456104 shows an upturn – i.e. the opposite of what we thought we’d see. The remaining narrow sequences above seem to show a slight downward trend as the day passes – i.e. reddening. Potentially, clouds dominated these days so that little of the surface was visible?

In understanding this this post, we now compare the results to those from the edge fitting method. We remind thereader that the Force method gives us the difference in B-V [magnitudes] between the BS and the DS, whil ethe dge-fitting method gives us albedos fro Earth in e.g. B and V bands. We look at the same range of phases in the two methods and repeat the relevant ploits here. First the edge-fitting method results:

Upper panel shows us B-band albedo minus V-band albedo against lunar phase (upper panel) and against the average alfa (B-alfa and V-alfa averaged).

We repeat the plot from the Force method here:

Here are the pdfs:

In the upper panels we see a slight dependence on lunar phase – B albedo minus V albedo as well as B-V [mags] do tend to rise towards Full Moon. Note that this implies opposite colours! The rise in B albedo relative to V albedo means that the Earth appears bluer as we approach Full Moon in the edge-fitting method, while larger B-V [mags] in the Force method implies a redder Earth as we approach Full Moon.

In the lower panels we see that the alfas used in the edge fitting method and the derived albedo differences have little relationship with each other, while the ‘incremental alfa’ needed in the Forcing method bears a strong relationship with B-V, as discussed before.

This seems to rule out the worst of the possibilities given in the Force posting – namely that if B and V albedos depended as strongly on alfa as B-V does then we would have a serious problem. In the present situation we are not in that position – but where are we then?

For both methods we seem to have outliers for phases less than -120 degrees. Let us take a closer look at these points. They are from night JD 2456073. The upper plot shows that the B and V albedos we could derive from edge fitting’ differs by .1 between the upper group and the lower group. From the Force method we see that B-V [mags] differs by 0.5 between these two groups.

The images involved have these names :

2456073.7452091MOON_B_AIR_DCR.fits
2456073.7472223MOON_V_AIR_DCR.fits
2456073.7555269MOON_B_AIR_DCR.fits
2456073.7658635MOON_B_AIR_DCR.fits
2456073.7758928MOON_B_AIR_DCR.fits
2456073.7781942MOON_V_AIR_DCR.fits
2456073.7862409MOON_B_AIR_DCR.fits
2456073.7882301MOON_V_AIR_DCR.fits
2456073.7964398MOON_B_AIR_DCR.fits
2456073.7983881MOON_V_AIR_DCR.fits

There are more data points for the edge-fitting method, in the plots above, than there are for the Forcing method. This is because more B and V combinations were picked tested for the former. The images picked were observed withing half an hour of each other. These images should be visually inspected.

[later:] Tried that – it gets messy: differences in the level-differences of the DSs in B and V can be spotted but do seem to depend on knowing the exposure times, and this is one of the things we do not think worked – better to trust the edge-fitting method, since it is a ‘common mode rejecting method’.

In this blog-entry we discussed a method to force halos to cancel by convolving the sharper of two images in a B/V pair, until their halos matched and could be cancelled by subtraction.

From the many possible B and V image pairs that we have, a number has been selected, and the method applied. The criterion used was that the linear slope on the DS disc should be within 1 sigma of 0. The result, for negative lunar phases, is this:

Here is a sharper pdf:

In the upper frame above, we plot the B-V difference between the BS and the DS. On average, B-V seems to be 0.4 bluer on the DS than on the BS. Given that the BS has B-V near 0.92 [Allen, 3rd ed.] we can convert the above to B-V for the DS by subtracting it from 0.92. That would give us B-V values for the DS near 0.5. Fred Franklin found less blue values since his B-V for DS was 0.17 below the BS value. That is at the lower limit of what we find.

We plot +/-1 error bars – they seem to get bigger as phase approaches Full Moon. We note that the scatter is larger than the errors – so something other than the errors we can estimate are causing the scatter.

Is there a trend in the above? By eye, it seems like there is – as phase approaches FM we have more positive ordinates – implying that the DS gets bluer with lunar phase.

The points fall in ‘columns’ corresponding to single nights – so something is causing a colour-change as the night progresses. All data are of course extinction-corrected so we doubt that the nightly colour-change is due to extinction. Is it a true color-change in earthshine?

We can compare the photometric B-V derived in the above way with B-V differences in albedo found by fitting the DS disc edge.

Now for the lower panel:

In the lower frame we plot B-V against the PSF-parameter alfa needed to bring the sharper image into the same state of blurriness as the less sharp image – so that halo cancellation could occur. There is a very clear relationship between alfa and the B-V achieved. I guess this is bad news? Nature (i.e. the B-V colour of Earth) should not depend on a PSF parameter used in the above convoluted way. Oh dear.

I still think it is worth comparing the present results with the ones based on fitting the DS edge. If we get a contradiction it may mean anything – but a strong agreement would imply that the ‘edge-fitting method’ has some gremlins hiding in it – i.e. that results are very dependent on those alfas.

It is also possibly revealing that the nightly points fall in columns [probably sequences once telescoped out to single nights – will check]. It is telling us something about things that change on an hourly basis – this includes airmass (i.e. extinction) and also the Earth-Moon-Observer aspect – so, still interesting!

Measuring B-V on the DS in B and V images with different PSFs: The Case For Forcing them!

We have detected some image-pairs where the B and V halo appear to cancel. We have plenty of other image-pairs where they evidently do not cancel. This is seen by structure on the sky – DS sky in particular as well as ‘slopes’ across the DS in the direction towards the terminator near the BS. We here investigate whether it is possible to ‘force’ the image with a narrower PSF to acquire a halo that resembles the broader image’s halo so that cancellation can occur.

If the differences between the PSF width-parameters (alfa) is small we expect to convolve the sharper image with a quite narrow PSF that only broadens the halo a little. We will use PSFs that are based on our ‘canonical PSF’ and simply raise it to a large number – this will generate a more and more delta-function like peak.

We identify centered pairs of B and V images that are close in time; identify the image with the narrower PSF – by using previously fitted alfa values; we then broaden the sharper image, while conserving flux, and inspect the difference between the B and V images. The results of such a trial is shown in the Figure:

Here is a less fuzzy pdf file of the above plot:

We see 15 panels showing a ‘slice’ through the disc centre, of B-V (magnitudes). Each panel has an increasing power (alfa) used to raise the canonical PSF to – at first we raise it to 1.5, then to 1.6, and so on. A smaller power will broaden the PSF a lot and thus give us ‘fuzzier’ images. We expect that when the image is too fuzzy we will see slopes in the B-V slice – just as when two observed images in a pair are subtracted and one has a broader PSF than the other. As alfa increases in the figure we reach a point where it is so large that the sharper image no longer is made appreciable more fuzzy – in the limit of large alfa it becomes an almost identical copy of the original image – hence, towards the end of the sequence of plots we essentially see what the slice across the original B-V image looked like.

In each image we notice the DS to the left and the BS to the right – both situated between the two vertical dashed lines, denoting edge of disc. We see a ‘level difference’ between the DS and the BS. We notice the slope of the DS – upwards for large alfa and downwards for smaller alfa – e.g. at alfa=1.7. Somewhere near alfa=1.8 and before alfa=1.9 the slope is horizontal.

We expect the real DS to be ‘flat’ [we need to check this against realistic models] so when alfa is near 1.85 it seems we have induced enough additional fuzziness in the sharper of the two images that the halos cancel and the level difference between DS and BS is what it is in reality, [We need to test this on images with known properties: Student Project!].

We notice also that the level difference depends on alfa, even when alfa has gone past the ‘now DS is flat’ point. This is a warning that incorrect estimation of
the limiting alfa value may give us spurious results for delta(B-V)! We estimate the slope on the DS (between vertical dotted lines) and look for the values of alfa where the DS slope is not statistically different from zero at the 3 sigma limit. We also do this for the BS. We find midpoints of the fitted DS line and the fitted BS line (red lines in uppermost plot) and take the difference bétween these, forming a sort of ‘+/- 3 sigma interval of confidence’ for the B-V difference between BS and DS. This is plotted here:

Here is the pdf:

We see thaat the generous +/-3 sigma interval allows quite a range in B-V – and also that if we do not estimate the right alfa in the above procedure – but make it too large then we will underestimate B-V.

For the 3-sigma limits we squint at the above graph and read off upper and lower B-V limits for the BS/DS difference: 0.5 to 0.7.

All of the above has been done for an average of an 11-row slice across the middle of the disc. The method could in general be extended to treat broader strips, or the whole ds – planes could be fitted instead of lines.
First, though, it is necessary to understand if the above is even correct! We need to look at artificial images of the Moon. These should have their own “B-V colours” then we fold the B and V images with seperate PSFs and try to force one to be as fuzzy as the other, and so on. The actual slope of the DS, when fitted with lines or planes could then be estimated and in fact used as the goal slopes for the above procedure. Work work work!

In V band, extinction is about 0.1 mag/airmass.

The Sun’s apparent brightness is -26.74.

0.1 mag of extinction means that about 10% of the light is taken from the sun and spread over the sky, about 21,000 square degrees. Half of this is scattered out of the atmosphere, half down to make the blue sky.

The integrated magnitude of the sky would then be of order -26.74 + 2.5*log10(10*2) = -23.5.. i.e. 3.2 magnitudes dimmer than the Sun.

Spreading this over half the sky (21000 square degrees or 2.7E11 square arcseconds) gives a reduction in surface brightness by 2.5*log10(2.7E11) = 28.6 magnitudes.

This yields a daytime sky brightness in V band of -24.24+28.6 = 4 mag square arcsecond.

If we have any daytime exposures with the ES telescope — we could check this and get numbers in other bands as well…

Is it possible to ab initio calculate the expected colour of earthshine? Earthshine, in the visual band, is due to diffuse scattering from clouds and the surface of the Earth, and Rayleigh scattering from molecules in the atmosphere. Can we simply calculate the expected B-V colour of Earth, seen from space?

If we assume that the diffuse scattering on clouds preserves the colour of the light (I think this is true – clouds look ‘white’) and if we assume that clouds cover most of Earth (a fair assumption given that the oceans are dark and cover 70% of the Earth) then we can make a model of Earth’s spectrum consisting of a fraction of the solar spectrum plus the complementary fraction of light ‘blued’ by Rayleigh scattering. This model assumes no ‘true absorption’ in the visible bands (also a fair assumption; ozone may actually absorb at UV wavelengths but the solar flux at those wavelengths is small).

By considering the Rayleigh effect and imposing flux conservation we get the following B-V colours for the spectrum of Earth:

f B-V
——–
0.0 0.63
0.1 0.51
0.2 0.40
0.3 0.29
0.4 0.20

where ‘f’ is the fraction of the beam that undergoes Rayleigh scattering. We have elsewhere estimated the B-V colour of Earth (here we do not mean ‘earthshine as measured on the DS’) to be 0.49 +/- 0.02. That would seem to imply that we have a correspondence with (crude) theory if 10% of the beam undergoes Rayleigh scattering.

Since most photometrically measured extinction coefficients are in the 0.1 – 0.2 mags/airmass range it does not seem unlikely.

In meteorological observations pyranometers are sometimes used to monitor the clarity of the air during the day – they typically measure the global radiation (all wavelengths, all directions – i.e. from the Sun and from the sky) and the ‘diffuse’ radiation – i.e. only the part of the light not coming directly from the Sun. Perhaps it is possible to estimate the fraction of light removed from the incoming beam, by scattering, using such, essentially bolometric measurements?

Using a SURFRAD pyranometer data set from Boulder CO, USA, I determine the noontime ratio of Diffuse to Global radiation, and plot it for each day of a whole year:

There is a great deal of variability! This is due to clouds. On perfectly clear days I guess that the ratio of diffuse to global radiation is at a minimum since the contribution to the diffuse is only from the blue sky while the decrease of the direct is also only the Rayleigh. Clouds tend to increase the former and decrease the latter giving a larger ratio.

The incoming beam contributes to the diffuse radiation and to the global radiation. The same amount of light scatters down to Earth as scatters into space so:

Incoming = Global + Diffuse

Then D/I = D/(D+G)

Looking at the plot it is clear that near the middle of the year you
can see D/(D+G) near 0.1. This is not exactly the same as saying ‘the
bolometric extinction coefficient is 0.1’ but we are close.

Problematic is the large amount of clouds – I shall look for a station that has less clouds!

It seems that Desert Rock Arizona is less cloudy:

Notice that we are getting D/(D+G) as low as 0.07. The altitude of Desert Rock AZ is 1007 m, while the station near Boulder, above is at 1/4 the altitude (213 m). Does that explain the different minimum ratio?

At MLO they publish a transmission coefficient which is approximately the same we are trying to look at – it would seem that at MLO the transmission can be 92-94% when there is no volcanic dust: http://www.esrl.noaa.gov/gmd/grad/mloapt.html

In considering B-V images, as here we are have to know how far the BS halo reaches onto the DS. We have looked at that before, here. We now revisit this issue. Here is a contour plot of the DS and the halo, with colours so that we can see how far the halo is likely to reach – we see it reach into the sky on the BS – how far do we think it reaches onto the DS?

Consider the yellow end of the greens, for instance – on the sky that contour lies 40-50 pixels from the BS rim – on the DS it lies adjacent to the red area which is the BS. The yellow-green contour therefore does not interfere much with most of the DS. The blue contours on the DS, howvere, are represented on the sky far away from the BS – so the blue areas on the DS may be interfered with by that part of the halo.

The above is V-band image. In B-V images we rely on much of the halo being similar in B and V and thus cancelling. The above noticeable gradient in DS brightness is absent in the B-V image.

[Update: See note at end!]

In this post Chris showed that the night of 2455945.1xxx has halos, in B and V, that seem to be at a constant distance across the DS. This implies that the ‘gradient problem’ is absent on this night – perhaps because the night was very clear or the haloes in B and V almost identical. [Inspection of the ‘alfa’ for the profiles shows that V images had alfa in the range 1.547–1.549 while B images had alfa in the range 1.544–1.547. These values seem close but halo shape is very sensitive to alfa, so analyze first!]. It is possible to generate a B-V image for this night using the only good images for B and V (in the sense of ‘being on Chris list of good images’, as explained elsewhere). We used one B and one V image:

2455945.1776847MOON_V_AIR_DCR.fits
2455945.1760145MOON_B_AIR_DCR.fits

both bias subtracted but not flat-fielded. Exposure times were taken from the image headers.

This is the result:

The colours are of course chosen arbitrarily, to aid excitement and imagination! For this image we corrected for extinction using

Vinst(idx) = -2.5*alog10(Vim(idx)) – Vam*0.1 ; kV=0.1
Binst(idx) = -2.5*alog10(Bim(idx)) – Bam*0.15 ; kB=0.15

and calculated B and V images from the B and V instrumental images, in an iterative manner:

BminusV=Vinst*0.0+0.92 ; BminusV is an IMAGE
for iter=0,10,1 do begin
V = Vinst + 15.07 – 0.05*(BminusV)
B = Binst + 14.75 + 0.21*(BminusV)
BminusV(idx)=B(idx)-V(idx)
print,iter,mean(BminusV(idx))
endfor

Here, ‘idx’ is a pointer that selects for all pixels on the lunar disc – both DS and BS.
Vonvergence was swift and independent of initial value. Convergence was to B-V=0.925.

A plot of a ‘slice’a cross the disc above is here:

We see a BS B-V value near the dashed line at 0.92 and a DS value between 0.75 to 0.8.

The colour image reveals some structure on the DS. It seems we see a faint flatf ield pattern? This should be looked at. On the BS colour differences eem to be lunar in origin and coul dbe compared to the literature on this subject.

The deep ‘cut’ in the middle is due to some small-distance difference in the B and V halos, we think.

Until we are told otherwise we think this is the first map or image showing the colour of the Dark Side of the Moon, which is also an album by Pink Floyd.

Images in other colours may be fothcoming.

Another version of the image above, here with a legend is here:
On the night of JD2455945.1ish the Earth was iilluminated over South America and shone on the Moon. Seen from the Moon the Earth looked like this, from the Fourmilab Earth Viewer:

So, mainly Pacific Ocean but also also all of South America. Need some satellite images of this to see how the clouds were distributed!

Note added later: The (far) above is unfortunately quite sensitive to how well the images are aligned. This needs to be looked at before anything else.

Here is a nice link to images of the
colour of the Moon obtained in sunshine and with simple DSLR cameras –
i.e. jpeg files: http://www.datarescue.com/life/kepler/moontests/raw_vs_jpeg.html

We have previously considered B-V images of the Moon. This was done with ‘raw’ images – that is, images where the halo had not been removed. Since we have the BBSO linear method implemented and since it does clean up the DS we can also calculate B-V images for the Moon based on these.

We have 55 pairs of B and V images thata re close in time (about 1 houror less apart). Using the standard star calibration relationships that Chris worked out fropm NGC6633 standard stars, we can convert images to instrumental magnitud eimages and from there to calibrated B and V magnitud eimages. We also corrected for extinction since the images were not obtained at the same time.

SInce the calibration relationships depend on B-V we have to assume some B-V values and iterate (Chris solves algebraically). The iterations converge quickly. We use onlythe brightest pixels in each image – i.e. the pixels delineating the BS – for calculating the mean B and mean V values needed to update B-V in each iteration.

The values for BS B-V that we converge to have this distribution:

The mean B-V=0.989, and the S.D.=0.019. The accepted value – e.g. Allen (4th ed), Table 12.14, gives B-V=0.92 (van den Bergh observations?). We therefore have a significant discrepancy. It should probably be noted that our values come from phases near 90 degrees, while the Allen values may be from ‘Full Moon’ conditions.

If we accept the above B-V (BS) values at face value we can continue:

The Sun has B-V=0.642 (Holmberg et al, MNRAS 2006). One reflection off the Moon reddends this value by 0.989 – 0.642 = 0.347. This value will also apply to earthshine that is observed after one reflection off the Moon, even if it is the DS. [We ignore here any colour-dependencies in the lunar surface mare vs highlands!]

If we can estimate the B-V of earthshine as seen on the lunar surface, we can work backwards to what the B-V of that light was before it struck the Moon – it will be the observed value minus 0.347.

Before trying this we need to understand to which degree the use of BBSO linear images, as opposed to ‘raw’ images, has helped us observe the true colour of the DS – has an important amount of the BS halo been removed from the DS?

We generate centered B-V images and plot the average of 20 rows across the middle of the images:

We see two panels – each panel is the result of using a fixed B image and two different V images – all three taken a short time apart. The black curve is the run of B-V values in the ‘raw’ image – that is, the image where no effort has been made to remove the BS halo. The red curve is from images cleaned with the BBSO linear method. The deep jag in the middle is the BS/DS terminator. The DS is to the left of this and the BS to the right. Since the BS is not altered by the BBSO linear method the red curve covers the black curve on the BS.

On the DS we see that cleaning the image has resulted in a slight reddening of the DS – it was ‘too blue’ in the red images.

We also see that the ‘linear gradient’ in B-V across the DS is unaltered qualitatively by the cleaning of the image. Why?

If we push on, ignoring the not-yet-understood gardient, and assume that the part of the DS closest to the sky has an un-polluted B-V value then we can calculate the colour of earthshine before it strikes the Moon, as explained above. First we extracted DS B-V values for that part of the DS disk that is to the left of 90% of the vertical columns on the disk. These values were on average 0.29 +/- 0.05 below the BS value.

If the BS value is given the canonical B-V=0.92, then we have a B-V for the DS of 0.63.

Franklin (JGR 72, no 11, p.2963-, 1967) measured B and V repeatedly on the DS. The difference between his mean B and his mean V is 0.64. We are close, but we are worried about scattered light!

Subtracting the effect of reflection once on the Moon brings us to the value for B-V of earthshine, before it strikes the Moon, that is, as it would be seen in space:

B-V_ES = 0.28.

There is one published B-V value for earthshine, based on Mariner II data in the 1960s. The paper is http://adsabs.harvard.edu/abs/1964JGR….69.4661W by Wildey. Unfortunately I cannot make head or tail of that paper!

Playing a bit more with the above, we can consider the effect of Earth on light – the Sun has B-V=0.642 when it strikes Earth. If the earthshine has B-V=0.28, then the bluing effect of Earth is 0.28-0.642 = -0.36 in B-V.

The PSF of the VE2-band Arcturus data deviated markedly from the other bands (next to previous post). We check on this further by looking at the PSF profiles of the stars in NGC6633, which were used to calibrate the filters (in this post).

We plot the central flux normalised PSFs in B, VE1, VE2 and IRCUT versus V (shown in green in each panel). The PSF is formed from 95 different stars spread over the frame in each case: e.g. this shows the stars in the IRCUT image.

The PSFs look like this:

Same plot as above but in log-linear form:

None of the PSFs are as sharp as the V band one — they all have more
flux at greater distance from the center. Interestingly, the VE2 band
profile is not discrepant, so the conclusion is that the IRCUT data for
Arcturus may simply have been out of focus.

It would seem the lesson from this is that the PSF can vary quite a bit — but note well that this is in the core only — as this is not a test of the extended power-law wings at all, as we are only probing to a radius of about 4 pixel (~30 arcsec) with these data.

In the previous post, we noticed that the VE1 and IRCUT filters have very similar profiles to the V band — at least for the data we obtained pointing at Arcturus.

Here we show a colour map (V-VE1, on an arbitrary magnitude scale: properly calibrated colours will come later).

It’s clear that the deep artifacts around the edge of the moon are much smaller now, compared to the B-V colour maps we have been producing to date.

Data used:

2456015.7558321MOON_V_AIR_averaged.fits 2456015.8108611MOON_VE1_AIR_averaged.fits

This shows that the core of the PSF of the V and VE1 filters are quite similar — and the power law tails as well — not just for Arcturus but for lunar images too.

This is a followup to the previous post on the V and B band PSFs.

Here is the stellar profile of Arcturus, shown in B (blue) and V (green). Clearly the B band light has a different core to the PSF than V.

Arcturus data were taken on the night of 2011-03-22. The files used are

align_stacked_2455643.4800437Arcturus-B-FILTER.fits align_stacked_2455643.4891739Arcturus-V-FILTER.fits

Oddly, the V band data are affected by what we called “shutter bounce” — a non-axisymmetric feature to the right (along crows in the CCD) of the PSF, but the B band data show no sign of such a bounce (the PSF profile can be exceedingly cleanly centered, unlike V band). No explanation for this for the moment!

Here are the other bands, all compared to V (green symbols):

PDF version here:

Here are the same profiles, this time with the data averaged into radial bins. The solid green line is the V band profile in each case.

VE1 and V are very similar — but this is what we would expect. IRCUT is not too far off V either.. the B and in particular the VE2 filters differ significantly from V. The difference between VE2 and all the other filters is huge, so we need to look into that next!

It is clear from the colour maps we have been making of the moon, that the B band PSF must be quite a bit broader than the one in V-band. We have only measured the V-band profile with any precision to date: this is a first go at the B-band profile.

The images used are these:

2455945.1776847MOON_V_AIR_averaged.fits

and

2455945.1760145MOON_B_AIR_averaged.fits

both taken at an airmass of ~1.3 – so quite high in the sky (which is good).

A wedge shaped region is defined for each, as shown in this figure by the green vectors:

The (V band) image is on a logarithmic scale, with counts/pixel shown along bottom. The red circle marks the fitted (by eye) position of the moon.

Centers and radii are :
V band center (195, 287), radius 147 pix
B band center (197, 291), radius 147 pix

Plot shows log radius (from center of moon, in pixels) versus flux (counts/pixel) in V band (green curve) and B band (blue curve).

It’s clear that the dark side of the moon (DS) is bluer than the bright side (blue curve offset by ~0.05 in the log on the DS).

Plotting distance from the edge of the bright side, we get this:

Focusing on the falloff of light off the brightside edge (coincidentally the B and V band fluxes are almost identical in the two images, so they are effectively normalised at the peak flux) — we see in the figure above that the blue light has considerably broader scattered light than the V band. (Note that the radial scale is now linear, not log).

Next plot shows the falloff of light on a log-log scale, as a function of distance from the lunar edge (less 3 pixels, in order to catch maximum flux).

Clearly blue light is more scattered than V!

Up to now, we have been using a single V-band scattering profile for the light — which is based on this technique of measuring the scattered light off the lunar edge, fitting a power law to the falloff at large distance (the power law tail of the profile, which has a slope of order -2.7 to -2.9 on very clear nights, at least in the V band).

The power law tail in B appears to have the same slope as in V — the B band light beyond ~ 10^1.5 = 30 pixels from the lunar edge out to the limit of the wedge (250 pixels) follows V band closely.

However, there is a large excess of light from the lunar edge out to about 25 pixels. This is messing up our colour maps, producing the dark (very blue) artifacts around the moon.

More work is certainly needed! A deconvolution method to extract the true PSF from full moon images, in all our bands, might be the thing to try next, rather than the “poor man’s” wedge method! (Note that in the wedge method, we correct for the fact that the moon is not a point source, we do not use the profile of the light in the wedge directly!)

We have some images of the SKE inserted in the ray path of our telescope, during Moon observations. We also have some images of the same from the BBSO telescope, kindly lent to us by Phil Goode.

BBSO images om two separate periods about a week apart. Top row shows histogram equalized imges of the BS (left and the DS (right) with SKE superimposed over the BS. Lower row shows the same – BS short exposure on left and DS long exposure with SKE on right.

Histogram-equalized image of the Moon with SKE inserted from our own telescope. SKE runs from upper left to lower right and is opaque above that line.

Our image shows a halo around the Moon that stretches behind the SKE – hence we know that halo is formed to a large degree in the secondary optics! On the BBSO images we do not see this as clearly – either because the BBSO telescope is of a different design or because the exposures with SKE (right column of images) covers much more of the BS. Assuming the latter we understand why there is no halo behind the BBSO SKE – the BS was effectively blocked and no bright light entered the secondary optics and could not cause a halo. On our image the SKE is inserted experimentally only, leaving a large part of the BS to shine into the secondary optics.

So, we learn that the halo is not from atmosphere or primary optics alone – apparently a large fraction of it comes from the secondary optics!

We also see the importance of having an SKE. While the BBSO group uses the ‘linear BBSO method’ to remove scattered light over the DS they have a much smaller problem than us because the halo from the BS is nowhere near as strong as ours is!

We now see how terribly important the SKE is.

I have been looking at the use of GPUs versus CPUs for our scattered light analysis. We need to be able to convolve artificial Lunar images (outside the atmosphere) with the instrument PSF. GPUs offer a considerable speed advantage.

First look (CPU versus GPU):

Upper left panel: artificial lunar image outside the atmosphere.

This artificial image is then convolved with a 2-D Gaussian-like PSF
which has fat (powerlaw) tails, and which closely reproduces what we see
in real data.

Upper right panel: convolution using 2-D FFT code running on a CPU

Lower left panel: convolution using 2-D FFT code running on a GPU

Lower right panel: the ratio of the two methods, i.e. the ratio of the two previous panels

There is a lot of structure in there, mainly images of the lunar
crescent turning up in different places — at a level of about 0.1% of
the intensity.

IMPORTANT: the CPU code was written in double precision, whereas the GPU was in single precision.

(The above reproduces with more explanation an earlier post)

Notes: The CPU code calls the FFTW3 libraries from Fortran (Dec’s ifort compiler is used), just using the standard Fortran to C wrappers provided with FFTW3. The GPU code is in written in CUDA.

Second look (CPU only, single versus double precision):

The plot above shows the ratio of the single precision CPU versus
double precision CPU (i.e. no GPU results shown on this plot).

There is similar structure in the ratio — and at about the
same level as the GPU tests gave, i.e. discrepancies at the level of a
few x 0.1% of the intensity.

Third look (CPU in double precision, renormalisation)

In this plot
we compare CPU double precision, applied to the ideal Lunar image, and without “min/max
renormalisation”. (Min/max renormalisation means scaling the input image so that the smallest value in the frame is 0.0 and the largest value is 1.0).

The ratio panel of the two convolutions (bottom right) shows noise only, and at a
very low level — 1 part in 1E7. Highly acceptable!

Fourth look (CPU, single precision, renormalisation)

This plot shows the same as the previous one — but with single precision rather than double. The artefacts are back, at the same old level of a few x 0.1%!

Thoughts:

We might already be able to conclude from the above that double precision FFT/CPU is robust (negligible
artefacts), but that a single precision CPU, or a single precision GPU,
produces similar sized (few x 0.1%), and thus slightly worrying, artefacts.

But I need access to a double precision GPU to test this. Hope to do so next week!
The acid test will be the results of comparing double precision on a CPU to a GPU.

We do have many problems to contend with – but now and then we are confirmed in what we set out to do: Here is an example of the scattered-light reduction of an observed image (2456034.114etc) where the EFM Method seems to be doing well compared to BBSO linear. We compare to a synthetic model image generated for the moment of observation (i.e. the libration and geometric factors are representative of the observing situation):

Top panel: slice across the Moon at row 296 so that Grimaldi near the edge is transected. The black curve is the synthetic model unconvolved – it is in units of W/m². The red curve is the scaled and offset profile along the same row of the EFM-cleaned image; and the blue curve is the BBSO-linear cleaned image identically scaled and offset.
Middle panel: detail of the above.
Bottom panel: difference between EFM and ideal and BBSO-linear and ideal, along row 296, expressed in percent of the ideal value.

The observed data were scaled since the units are different; they were offset because the EFM model did not have 0 value on the DS sky. The BBSO-linear, being ‘anchored in the DS sky’ did have a 0 value in the sky. The same scaling factor and offset was used on EFM and BBSO, for comparison, however – hence the little offset on the DS sky at right.

What do we see? Well, the EFM-cleaned (red) line follows the synthetic model quite nicely between column numbers 325 and 370 while the blue line (BBSO-cleaned) diverges all the time. Near the DS edge (columns 370-385) the synthetic model is higher than EFM.

What does it mean? The BBSO-linear method has better removed the flux on the sky – it is designed to do that, while the EFM is designed to minimize the residuals squared in a mask on the sky around the Moon. This implies that the BBSO-linear method, in the present case, would be better than EFM on the lunar disc near the DS sky. As we move further onto the DS disc the BBSO-linear method will fail more and more since the halo is not linear with distance from the edge – the method underestimates the amount of scattered light on the disc between the edge and the BS. We do see how the blue line diverges more and more; we do see the red line cling closer to the ‘true value’ (assuming the synthetic model is ‘true’) across the disc, before it too fails nearer the halo and the BS. The behavior nearer the edge may be a consequence of how we model the synthetic images – we have to use a reflectance model to make the synthetic images – and the angles of emergence and incidence corresponding to ‘near the edge’ is not one for which the reflectance model is inherently very good. The models we use are based on the simple ‘Hapke 63’ model. We speculate that the disparity between observations and synthetic models in the columns 370-385 is due to model inadequacy.

Anyone using a reflectance model as simple as the Hapke 63 model will encounter the above problem if they try to use pixels near the lunar disc edge – the natural thing to do, if the halo is removed using the linear method is to use edge or near-edge pixels.

Hence the problem of reflectance modeling and the inadequate linear method become coupled! Our EFM method seems to be a way around this obstacle – allowing use of disc areas further from the edge where simple reflectance modelling is adequate – hence we should be getting more reliable terrestrial albedo data. One day.

In the EFM-method we determine the alfa values for every image. Is there a link between the alfa value for one filter and the rest in an interval of time? It is our understanding that alfa is determined by the amount of scattering in the optics plus the atmosphere. We therefore expect that on ‘bad’ nights the alfa values will tend to move in the same direction. We investigate this here.

We find all alfa values in all EFM-treated images. We sort them into half-hour bins. We calculate all the alfa values in each bin and plot the results. Below is a pdf showing all the plots between some filters, at different ‘zoom-levels’. The image shows the last zoom-level, highlighting the dense ‘clump’ of points:

There seems to be a general agreement that the alfa values are correlated – bad nights (i.e. ‘broad PSFs’) occur in all filters at the same time. Since VE1 is just about identical to IRCUT the scatter seen above means that the fitting routine is unable to make a perfect match – or that observing conditions, during the half-hour bins used, changed.

Using 15 minute bins does not improve matters:
I therefore suspect that the fitting method does not find the best fit each time.

We compare the B-V values on the DS of images from JD2456034: we look at images only exposed to bias-removal (‘DCR’ images) and images cleaned with EFM and images cleaned with the two variants of the BBSO method: linear and log. We show a ‘slice’ across the disc at row 256:

Top left: black is the B-V slice from the DCR images; red is the B-V values from the BBSO-linear cleaned images. The second graph in the upper left panel is the difference between the B-V values (BBSO-linear images minus cleaned image). The vertical dashed lines show the edge of the lunar disc and the start of the BS in column units.
Top right: same, but for BBSO-log method.
Bottom left: same but for EFM images, as shown elsewhere on this blog.

The results are quite different – the BBSO-linear method has given us fairly constant B-V values across the disc – they are about .05 mag below the DCR values (i.e BBSO-linear are bluer than DCR values).

The BBSO-log values seem completely unrealistic.

The EFM-cleaned values also look a bit unrealistic in that there is a spatial dependence on the magnitude of the B-V relative to DCR – in a way that looks like a remnant effect of the halos. The Delta(B-V) value changes sign across the disc.

The night JD2456034 is very close to New Moon and we know from other results that this is when the BBSO-linear method is likely to work best (the halo being small). Since there is a phase-dependency in the overall results for BBSO methods over and above what EFM shows we know we cannot universally use BBSO results. On the DS, towards the edge, the BBSO-linear method should be very good – it is ‘anchored in the sky’, unlike the present EFM method, and therefore should be unbiased near the edge. Our EFM method, at the moment, only minimizes the square of the residuals on the sky, inside a mask.

Near the sky, the DS B-V values in the BBSO-linear and EFM-cleaned images differ by about 0.04 mags. We should look at EFM methods that also ‘anchor in the sky’ and see what we get then.

The present largest worry is not the B-V offset but rather the dependence on position on the disc that the EFM method shows.

What does the literature tell us we should expect B-V to be under earthshine?

At the moment I only know of Franklin’s 1967 paper. It gives B and V values of the earthshine.

http://adsabs.harvard.edu/abs/1967JGR….72.2963F

The average B-V seems to be +0.64. Can we find any other information on B and V?

In this entry on the blog: http://iloapp.thejll.com/blog/earthshine?Home&post=272 we discussed the Laplacian method of estimating earthshine intensity by removing the influence of the scattered light halo. We have tested it on a good image from near New Moon (i.e. lots of bright DS, very little interfering BS). Consider this plot:

In the top two rows we see the observed image, the EFM-image and the Laplacians of these. In the second row second plot there are three additional curves, in black. The centre one shows the estimate of where the disc edge is, based on the estimates we have of disc centre and disc radius. The other two lines are parallel outliers at + and – 40 pixels from the edge. In the third row we show first the coaddition of all rows along the edge, using the estimates of the edge and the outliers as start and stop – this allows a coaddition that centres on the disc edge. Second in third row is the same thing based on the Laplacian of the observed image. We see the signature of the double derivative of an edge. In the last row we show the same as in the third row but for the EFM-cleaned image and its Laplacian.

The application of the EFM method removes most of the halo and this is consistent with the slight lowering of the slope after the ‘step’ in the first panel of fourth row, compared to first panel of third row. We extract the size of the ‘jump’ as well as the distance between the minimum and the maximum of the Laplacian edge signatures. The jump size is estimated as the difference between the mean of the first third of the curve and the last third. We collect the four estimates of the earthshine intensity – two jump heights and two estimates based on the Laplacian signature, and plot them below for 4 differently rebinned versions of the original images – unbinned; binned 2×2; binned 4×4 and binned 8×8. This is done to ‘sharpen’ the edge in case the edge is ‘fuzzy’ so that the 3-point numerical operator used in constructing the Laplacian will interact correctly with the edge.

We see color-coded graphs. The estimates of earthshine intensity based on the observed image and the EFM-cleaned image are black and blue, respectively. The estimates based on the ‘signature size’ in the Laplacian of the observed and EFM-cleaned images are red and orange, respectively.

The differences are quite large – several percent, and there is a dependence on the rebinning of the image – even in the case of the step method. The results that ‘ought’ to be the best is the jump-size one based on the EFM-cleaned one since the halo is absent and because the jump size method should have less statistical error-of-the-mean than the signature size method since more pixels are involved in estimating the former. The error-in-the-mean for the jump estimates for the unbinned image is of the order of tenths of percent – clearly dominated by systematic effects.

The above suggests that there are large method dependencies.

The ‘jump estimator’ may be victim of the slopes seen near the jump – especially in the EFM case – with such slopes the mean before and after the jump is dependent on just how you estimate the jump. Possible improvements is to extend a linear regression from the first half or third of the points, before the jump to the midpoint and do the same from the other side and calculate the ‘jump size’ as the difference between the two extrapolations. The slope is still a concern, however. It is due to the lunar surface albedo variations. The Jump size itself is a function of earthshine intensity as well as lunar surface albedo and therefore on libration, as well as the geometric distance factors. Being estimated right at the edge the reflectance is a constant, leaving lunar albedo as a variable factor (through libration). This could be taken out with knowledge of the parts of the surface of the Moon right at the edge – i.e. via modelling.

The ‘signature size’ method has different challenges – it is based on fewer pixels for one thing. It also depends on the lunar albedo – and not on reflectance, provided the edge is used (and not features on the DS, as Langford et al do). However, the method depends on knowing that the image is in good focus. Bad focus will blur the edge causing a lowering of the signature size – for area methods such as the EFM we do not depend on a sharp edge but can estimate levels from regions chosen away from the edge. Of course, we still get blurring from lunar features overlapping.

The above study did reveal dependence on binning – but if blurring was being addressed the results are not encouraged as there was no convergence.

So many things to compare! Super student project! A study of artificially blurred synthetic images could be done, and ‘step’ and ‘signature’ methods compared.

In an interesting article Sally Langford, et al. [ http://adsabs.harvard.edu/abs/2009AsBio…9..305L ] describe an analysis method for earthshine images. Essentially, the Laplacian’s effect on images is to detect edges and the size of the ‘jump’ in going from e.g. sky to lunar disc in the observed image translates into an amplitude of a characteristic signature of the derivatives of an edge. This enables the large-spatial-scale features of the scattered light halo to be removed from the short-spatial-scale features of edges and craters.

Our method for removing effects of halo is to construct a forward model of the halo, based on the BS of the Moon (caught in the same image as the DS – Sally’s FOV is smaller – about 20 arc minutes; our is 60 arc minutes), and subtracting it, leaving the DS almost halo-free.

How well does the Laplacian method work on our type of images? Sally’s images were well-exposed images of the DS alone – allowing large counts – 10’s of thousands; ours are co-addition of 100 images taken so that the BS is not overexposed, leaving the signal on the DS typically near 10 counts. One image with 10,000 counts should thus have the SNR of 100 of ours.

We took one of our good images – from JD 2456095 – near New Moon – and plot the profile across the image near disc centre in the image that we obtain using the forward modelling method (named ‘EFM’ in the plot and this blog) and the same profile from the Laplacian of the image. We average 20 rows.
Inspection of the columns (approx # 100-120) where the lunar disc edge is we see a clear ‘step’ up from the sky level in the EFM-cleaned image but there is no sign, above the noise, in the Laplacian of the same image. Inspecting the whole Laplacian of the EFM image we do see a faint signature:
The above image is ‘histogram equalized’ to show the feature. Our own method – also shown when histogram-equalized is here:

We see that the halo has been only partially removed on the right.

I think we would be hard pressed to extract a signal from the Laplacian of the image. It also bothers me that the whole earthshine signal is reduced to the value of the derivatives in just edge pixels. In our own image we have hundreds of pixels to measure on – in the Laplacian we get only a signal along the edge.

We should note that the remnants of halo on the right are much less evident in the Laplacian image, suggesting that there is a kernel of a good idea in the method. For reference, the Laplacian of the raw observed image is here (histogram equalized):

It does seem as if the Laplacian helps remove a lot of the halo, but also reduces the analyzable part of the image to the DS edge.

The Langford paper mentions both smoothing images and co-adding them – the former is done at the resolution of the worst image for a sequence. The Langford paper analyses features on the disk – not the edge.

Working out how the Laplacian is best used, on realistic images of the type we have, would be a good student project!

We have, below, studied the ‘structures’ or ‘bands’ that, in some images, stretch from side to side. One concern was that the structure was induced by some step in image reductions. We therefore compare raw images to EFM-cleaned images to see if the structure is present in both.

Visual inspection of histogram-equalized versions of raw and EFM versions of the images reveals the structure, although it is easier to see in the EFM image since the halo and most of the BS is gone:

Plotting the average of the first 20 columns of each image, and scaling them to each other shows that the structure better:

The black line is from the EFM-cleaned image, while the red curve is from the raw image.

This seems to rule out that the ‘structures’ are induced by image-reduction steps.

We remain suspicious of the possibility that the bright light from the BS itself somehow is causing the problem. It is not a reflection since the structure is a deficiency of counts. Only way to explain it is by some process subtracting counts from these bands – or lowering the sensitivity of the whole set of bands – perhaps some sort of non-linearity caused by strong light in one part of the CCD affecting the whole row?

We still need to map when the structures is present – we want to see if it appears in a timeline or is a function of the placement of the image in the frame. Inspecting the first 700 images of 2000 for the EFM frames shows no sudden dependency on time of this problem.

Idea: also extract the maximum counts of the image and use this as a parameter – we do have suspicions that non-linearity sets in far below 55.000 counts.

Question: Why is the structure along rows, while readout direction is along columns? What is special, in a CCD chip, about the rows?

More to follow.

We have extracted the mean value of alfa (the parameter that describes how broad the PSF for a given image is) for a given night, and the corresponding value of the extinction coefficient for that night. We have only few nights where the extinction could be determined.

The plot of one vs the other looks like this:

(download and look, etc)

We see that there is a tendency for alfa to be narrowly constrained, and that the extinction has a broader distribution. In general there is no strong relationship between the values, but if we ignore outliers and the effect they have on the regression (plotted as a red line) we see a general tendency for high extinctions and small alfa values to be related: For B it is quite clear. V would be clear but for the outlier, IRCUT also seems to be clear. VE1 and VE2 are all over the place. A broad PSF is given by small values of alfa. We expect broad halos (i.e. broad PSFs) on hazy or turbid nights – nights on which extinction also should be large.

Factors that determine scatter in alfa are things like image focus and how the nonlinear fitting routine determined it should stop. Physical factors include haze and thin cloudyness on the night in question.

Factors influencing the scatter in extinction include the actual regression: we used all nights with more than 3 observations used to determine the airmass vs extinction line, and for which the determinations from 3 regression methods agreed to a S.D. of less than 0.02. The three methods were – ordinary least squares, and two ‘robust’ methods (“ladfit” and “robust_poly_fit” in IDL). While doing the actual regression it was necessary to eliminate some outliers by hand. Physical factors include haze and cludyness and whether the halo around the Moon was well captured inside the image frame.

The relatively low value of alfa for the VE2 filter is still not understood.

At the moment we are displaying our results by providing plots of ratios of the DS/BS ratio in observations, to that in models. We do this chiefly to get rid of common factors – such as solar irradiance and distance-related geometry.

What remains in such a ‘ratio of ratios’ are the effects of:

1) The Earth’s actual albedo,
2) the model’s Earth albedo,
3) Earth’s reflectance (real vs modelled),
4) lunar surface albedo’s in the reference patches (real ratio vs model ratio)
5) effects due to the choice of lunar reflectance model.

We are really only interested in 1.

2 is an assumed value so that the results we get for terrestrial albedo are relative to that choice. We use a value of 0.31.

3 is a choice – we expect that any errors made in this choice will be seen as a phase-dependency in the results, and we can therefore control or at least understand it. Earth is more Lambertian than the Moon. The Earth has edge-darkening, the Moon has very little. We use a Lambertian model for an otherwise uniform Earth.

4 is observable, but only with difficulty – you need a good total lunar eclipse, then the albedos in the two patches on the Moon – or their ratio – can be measured. BBSO has done this. We have not (yet). In the model we make a choice, based on which lunar albedo map we use. As long as it is fixed the results will be relative to that choice. Perhaps we can use published images of total lunar eclipses to extract the ratio?

5 is a choice – we have models for Lambertian reflectance, as well as the Hapke 63 model and other, as yet untested, more advanced reflectance models. We expect that incompletenesses in these models will be seen as a phase-dependency in the results.

Of the above only some could possibly induce a phase dependency: 3 and 5.

We have reduced all data using both lunar models based on the Lamertian reflectance and the Hapke 63 model. We show them next – look here for a discussion of what we are actually showing: ‘ratio of ratios’ and all that:

(as before, plot needs to be downloaded as it does not show up well on this blog).

The plot is for 5 filters with the EFM method applied. The left column is the ratio of DS/BS in obs to the same in model, while the right column is DS/total in obs relative to models, where ‘total’ means the disk-integrated brightness of all source counts (plus a few stars that we can ignore!).

The first page is for the Lambertian lunar reflectance. The last page is for the Hapke 63 reflectance.

We notice that left and right columns are quite similar. We notice that the ‘jump up’ at angles corresponding to about 40 degrees from New Moon is much smaller in the Hapke 63 model results. Inspection of the files that correspond to the individual points in the ‘jump up’ and those next to the jump, reveals that the ‘jumped up’ points have a different processing history: they are the result of coadding and averaging single image sequences, while the rest are stacks of observations that were averaged. Why these should be different is unclear, as yet.

We conclude that there was an effect of lunar reflectance model on our results – and that Hapke 63 is better than Lambert. This is not surprising as the Moon is well-known not to be Lambertian in its reflectance.

So – we are beginning to see observations constrain theory!

There is still some scatter to account for, and we shall return to this mattrer, using the estimates of nightly seeing available from the measured alfa.

There is also some ‘slope’ in the result wrt phase – so the EFM method has a success rate that is phase dependent. WHile the second-best method (BBSO logarithmic; not shown here) has some ‘upturn’ towards Full Moon (center of plot), the EFM has an even slope down towards FM. Perhaps some empirical fine-tuning of the method will remove this problem too.

In a previous post (here) we have compared the Wildey and Clementine albedo maps. These maps are important for our synthetic modelling code since the albedo (along with reflectance assumptions and correct geometries) are the basis of constructing realistic model images, used in analysis. We can compare these two maps very directly by accessing both and plotting the mean albedo in boxes at common lon,lat positions:

Evidently, the Clementine albedos are lower than the Wildey ones by a factor of two for dark areas and by some tens of percent at brighter areas.

This may be an explanation for the discrepancy we have between observed morning/evening brightness ratios and modelled ones, which Chris Flynn put his finger on.
We seem to have [not shown, but material could be inserted] an observed difference (expressed in magnitudes) of 0.12-0.14 magnitudes at absolute phase 90 degrees between the morning and evening integrated brightness. In models, based on Clementine albedo and either Hapke 63 or Lambert reflectances the difference is more like 0.3 magnitudes. If the Wildey map is more correct than Clementine, in terms of the highlands/mare albedo, then using the ‘flatter’ Wildey albedos in the synthetic code would help on the morning/evening brightness issue by lowering that ratio in the models.

The issue with using Wildey instead of Clementine is that Clementine is a global map – Wildey only covers something like -89 to 89 degrees in longitude and something similar in latitude so there is no remedy for modelling under lunar libration. We could ‘scale’ the Clementine map, using the above relation and see what that gets us, though.

Before we do that, we should fully understand how our own synthetic code uses the Clementine map – I believe there is more to its use than merely being used as a lookup-table. Hans will be able to tell us about this.

A second-order robust polynomial fit to the above data (where points where Clem > Wild have been omitted) is:

W=0.067509410+C*0.50252315+C^2*1.3644194

Earlier we showed the performance of the linear and logarithmic BBSO method. Now we have added the EFM method.

As before, plots are not shown well on this blog, so please download this pdf file and look at it:

You see the same 5 panels on each page – one panel per filter. On the panels you see the behaviour of the ‘ratio of ratios’ against phase – that is, the DS/BS ratio in observations relative to a set of synthetic lunar model images.

First page: RAW data – i.e. no scattered light has been removed. Full Moon is at the centre of the plot so we see the effect of scattered light – the obs/model ratio is increasingly not 1 as we near FM.

Second page: the linear BBSO method has been appllied to images. We see a reduction in scattered light – points ‘move down’ towards the 1-line.

Third page: logarithmic version of BBSO method – a slight improvement is seen.

Fourth page – the new one: the EFM method – a quite large improvement is seen over the best of the others! Many of the selected points are now lying on a flat seqeunce, except points towards New Moon, and some outliers.

The EFM data has been ‘selected’ in the sense that the parameter alfa, determined from images, has been used to select ‘good cleanups’. A histogram of the detected alfa values was made and a notable peak found and only images with alfa in a narrow range were used for the above plot. Alfa between 1.67 and 1.73 were picked. The absence of VE2 points is due to this – the peak of the VE2 alfa distribution is between 1.60 and 1.66. We must investigate whether these are ‘good’ solutions. [added later: a brief visual inspection seems to imply that the 1.7 solutions are the good ones – not the majority, which for VE2, lies in the peak at 1.63. Hmm.]

Speculatively, we note the ‘turnup’ of points towards New Moon. We ahave earlier discussed that this may be due to some feature or failure of the synethtic lunar image model to correctly portray the Moon at large phase angles. On the other hand the ‘turnups’ now look less like a gradual sequence, and more like a ‘jump’ up. The ‘jumps’ occur near 40 degrees from New Moon – …. is this the ‘rainbow angle’? The jump seems large – the rainbow angle paper spoke not of factors of 2 and 3 and 4 but of percentages. On the other hand the phenomenon has never been observed, as far as we know. It would be nice to have points ‘outside’ the rainbow angle to see the jump go down again (if this is the rainbow, that is). We note here that as we get closer to New Moon the lunar sicle is narrow and it becomes increasingly important to place the ‘patch’ in which the BS brightness is measured – here is an opportunity to experiment with DS/total ratios, instead.

Work continues …

We evaluate the effect of two different scattered-light removing techniques on our data, comparing to no removal at all. We do this by considering the DS/BS ratio – and we look at the ratio of this ratio in observations to the same ratio in models – so it is a ratio of ratios, ok?

The graphics do not reproduce well on this blog so I put a link to a pdf file here: download it and read on.

There are 3 pages to look at, and on each page there are 5 panels – one panel per filter, one page per data set.

The first data set shows the ratio of ratios extracted from data where only a bias has been subtracted, plotted against lunar phase. Second page shows the ratio of ratios when the LINEAR BBSO method has been applied, and the last page shows the ratio of ratios when our logarithmic variant of the BBSO method has been applied.

Full moon is at phase 0.

We see that there is a scatter and there is a systematic dependence on lunar phase – the points ‘curve up’ towards full moon (the middle) and towards the edges (new moon).

Near Full Moon it becomes increasingly difficult to remove scattered light because the BS is closer and closer to the patch on the DS where photometry was extracted. The Eartshine is also weaker and weaker as you approach Full Moon because that corresponds to approaching New Earth.

Therefore a hypothesis was that the ‘curve up’ towards Full Moon was due to incompletely removed scattered light. Looking at the raw vs linear vs log methods it is evident that this ‘curve up’ is substantially reduced by application of either of the scattered-light removal methods. The linear method does well and the log method adds an increment of improvement.

I believe these plots are a demonstration that BBSO linear is not perfect – and that small improvements are possible. It is open to discussion what BBSO has done about this problem – light seems not to be completely removed in their method. How do they average their data to compensate? We should note that our image scales and observing methods (ND filters; coadd etc) are not identical.

Soon we will add a third method – the EFM!

By running the computers for one week we have been able to process all observed images (both singles and those that had been coadded from stacks) with the EFM method, and thus estimated the exponent, alfa, that the basic PSF must be raised to to model the halo around the Moon. Plotting the histogram of 9000 values for alfa we get:

We see a broad distribution from about 1.3 to 2.1. Inspection of the results reveals that only the images with alfa values very close to ~1.7 are any good. The lower values correspond to hazy or even partially cloudy nights. The values higher than ~1.7 have yet to be examined.

A typical (histogram equalized) image of the residuals for one of the images with alfa near 1.71 is here:

A ‘slice’ across the image shows this:

Upper panel is the slice – black being the image, red being the fitted halo; second panel is a detail view of panel one showing the DS, vertical dashed lines showing the limits to the sky on which the halo is fitted; third panel shows the difference between the black and red curves in panel two.

We see that the EFM method has been able to fit the sky so that it is essentially just noise – even on the BS (right) side of the disc, and that the DS has been revealed for a wide area onto the disc itself, only near the terminator is there a problem with the subtraction – the residual dips down.

The results for applying the linear BBSO scattered-light removal method to all data are further considered. Here we show data where the lunar disk is well centred and the radii extracted are between 131 and 149 (real range is 132-148) pixels. (Sorry about the image rotation!):

Rotate image and then use this caption: Results from comparing model images to observed images – models calculated for the observing moment but no halo generated. In the observations the halo has been attempted removed via the BBSO linear method. Column (1) DS counts to total counts for observation (black symbols) and same for Model in red symbols, against lunar phase; Column (2) observed DS to total ratio against airmass; (3) ratio of the DS/tot ratio for observation and models.Note that Full Moon is at 0 degrees phase.

We note that for lunar phases between New Moon and about 100 or 110 degrees the model and observed ratio of DS to total flux behave similarly. For phases closer to Full Moon the model DS/tot ratio is much smaller.

That is probably because near Full Moon the halo gets closer to the DS patch where the DS is measured: In the model no halo is present so the DS brightness is not polluted until the patch is actually part of the BS – in the observations the halo moves with the lunar terminator towards the BS and becomes harder and harder to remove so that it starts to pollute the DS – even if measured in images with some of the halo removed.

In the middle column we see a fair spread in observed DS/tot vs airmass. Most of the observations are for AM < 3, while a few have even higher airmasses. In the bulk of the observations we see no filter-to-filter consistent evolution of DS/tot ratio against airmass, and conclude that there appears to be no dependence on airmass; this is good as the whole point of simultaneous observing of DS and BS is that external factors, like extinction, affect both sides identically.

In the last column we see the ratio of the DS/tot number between observations and models. The behaviour noted from the first column, above, is evident here – the DS/tot ratio in models is much smaller than in observations for phases inside 110 degrees or so and the ratio of ratios goes up accordingly. For phases between 110 or 100 degrees and New Moon we see a more constant behaviour consistent with the observation made above that model and observations behave similarly. We see, however, a consistent pattern from filter to filter – there is some ‘curvature’ with phase. The inner part of this could be due to lingering halo effects, but the outer ‘curve up’ must be due to effects in the model, not in the observations. Near New Moon it seems that the DS/tot ratio for observations grows compared to the same ratio for models – the models’ DS is not as bright as in the observations near New Moon. We note that the DS model brightness is a function both of the reflectance description of the Moon and of the Earth.

Summary:

1) There is indication that halo-removal near Full Moon (phases between 0 and +/- 100 degrees) is increasingly incomplete.

2) There seems to be little or no effect of airmass on the DS/tot data.

3) There is some indication of a reflectance problem in the models for phases approaching New Moon – either in the Earth description or in the Moon.

Comments:

1) is good news – we thought we could only do removal up toi half Moon – the limit seem sto be 10 or 20 degrees beyond that, towards a fuller Moon. 2) is very good news – it opens up for our method to be applied at low altitudes, i.e. for small lunar phases near New Moon.
3 should be investigated by now redoing the whole reduction but with some simple and single change in the modelling – we cannot change the BRDF for Earth (it is always Lambertian), but we can change the lunar BRDF from Hapke 63 to Lambertian.

In these models the Earth was modelled as a ‘cloud free Lambert sphere’ with ocean/land contrast built into the single-scattering albedo. The Moon was modelled by the Clementine single-scattering albedo map with Hapke 1963 reflectance.

Residuals image of Moon minus model of BS and its halo, fitted to the sky (i.e. the EFM method). The green and blue colors are a few tenths on both sides of 0. The mask used covers the image except the disk of the Moon.

In this application of the EFM to an observed image of high quality we see that the residual sky is not even. Particularly on the BS (to the right) there are some streaks. A blob at the bottom also shows unevenness, as does the striations to the left of the DS.

This illustrates a shortcoming of the EFM. Since the PSF used is rotationally symmetric it is possible that the unevenness seen above are in the observed image – perhaps some sort of reflection effects in the optics. The EFM is a global method seeking to subtract all the scattered light with one model. The BBSO method is essentially local – it estimates the scattered light in a wedge off the disk. On the other hand the BBSO linear method almost certainly removes too little of the halo as you get nearer to the BS. The EFM method attempts to model all parts of the halo.

The linear BBSO (i.e. as done by BBSO) method and the modification that operates on logarithmed images can be compared:

A strip, 20 rows wide, across the image was made and the rows of that strip averaged, and the average plotted. Black is the original bias-subtracted image, red is theordinary linear BBSO method and blue is the log-method. We plot the absolute value of the fluxes to avoid problems with log.

Both the linear and the log method does well on the sky at left, but there are differences on the edge-near disk. Tests on synthetic images, where we know the flux to expect on the disk, has shown that the log-method is more accurate than the linear method.

The residual mean (about 2) is unsettling, but has to do with the use of absolute values and the averaging over several rows – inspection of the images reveal values distributed around 0.

In order to precisely remove scattered light from the observed images we need to know the centre of the lunar disc in image coordinates, as well as the radius. These numbers are used by the BBSO method, while the EFM method can work without them. Extraction of fluxes from designated areas of the Moon also requires knowledge of the disc coordinates.

We have, as described elsewhere, found a fairly good way to estimate disc centre and radius, and have more than 5000 images (singles or sums of stacks) with know coordinate estimates.

We can check on the quality of these data by inspecting the time evolution of the disc radius in terms of the lunar synodic period (27.322 days).

Plotting the detected radius against observing time modolu 27.322 we get, for the 5 filters:

There seem to be some outlier groups, as well as a general scatter. The scatter is on the order of 2-3 pixels while the outliers reach 5. These outliers can be identified and the relevant images inspected.

We fit a general sine curve to the data, and get:

Filter: B
Offset Amplitude Period
141.169 7.98152 27.6367
+/- 0.0108026 0.0166443 0.00208776

Filter: VE1
140.663 7.81091 27.5796
0.0122441 0.0199136 0.00187041

Filter: V
140.647 -7.47896 27.5736
0.0413606 0.0837630 0.0111588

Filter: VE2
139.845 8.51341 27.5771
0.0245651 0.0383531 0.00435187

Filter: IRCUT
140.744 7.98575 27.6255
0.0358718 0.0596396 0.00708377

The period is not close to the expected 27.322 days. We expect this is due to a poor fit (in turn due to the outliers). We identify the outliers. 108 images are found that have radius more than +/-3 from the fitted sine curve.

Upon inspection, it turns out that not many of the identified outliers are obvious ‘bad images’. The determination of radius and disc centre is therefore somewhat deficient.

Update:Thanks to help from Rong Fu it now seems we have a working Fortran 90 code for the EFM Method! The code is being tested against synthetic test cases and a performance report will be written. The code will be set up to run on the DMI Cray and with it we will be able to process much of the observations we have accumulated, and improve the paper significantly.The point of the new code is that it is no longer written in C (which is a bit of a mystery to us old-timers) and can now more easily be experimented with and optimized.

On the night 2456003 we happen to have a nice sequence of the Moon rising through a clear atmosphere and we can determine the instrumental alfa – i.e. the power of the PSF due to instrumental effects only.

Normally alfa depends on atmospheric and instrumental effects but on a steady night it is possible to gather data for several airmasses and then extrapolate to airmass zero.
The different colors are for different filters (more or less matching the wavelength – blue to red). The data were obtained from the fits of the EFM method.

We see a little bit of scatter in VE2 (red) but otherwise a remarkable agreement in intercept and slope. It would seem that the instrumental alfa is very near 1.724 for all filters. Since we have a ‘base PSF’ (with a power of 1.6) which we raise to alfa’th power to fit observed PSFs we see that the instrumental ‘real alfa’ is near 2.76. The diffraction limit is at 3. We are not saying we have a system that is diffraction limited – but any values of alfa larger than 3 would be impossible!

If we assume that the instrumental alfa is always fixed or decreases (large alfa means narrow PSF, small alfa means wide PSF) with time then measuring the alfa_instr on all good nights would allow us to perform some sort of quality check on the data for that night. Also, keeping an eye out for the largest alfa_instr ever will help us set better upper limits.

Such diagnostic testing as the above should be made a routine part of a future pipeline for reducing data.

Note: This is our second visit to night 2456003 – we have also looked at the same as the above but based on FFM, here. The value for alfa_instr is nearer 1.76 for the FFM, suggesting that there is method-dependence. However, FFM is currently our most troublesome method, of the FM methods.

A central tool in both our and BBSO’s analysis methods is the understanding of the lunar reflectance since it relates the observations to fluxes from the Earth and ultimately the terrestrial albedo. BBSO has measured and determined their own lunar reflectance function. Until now we have been using the Hapke 1963 reflectance model.

Kieffer&Stone 2005 have a reflectance model suitable for phases between Full Moon and Half Moon – these are based on terrestrial ROLO observations of the Moon.

From space SOLSTICE has imaged the Moon at lunar phases to 170 degrees and they show a preliminary reflectance function.

We compare now the irradiance expected from the Moon at a range of phases given our own eshine code (i.e. based on Hapke63), and the Kieffer&Stone model as well as the LLAMAS reflectance from SOLSTICE. The curves have been scaled vertically to coincide at intermediary phase angles.

The crosses are from the Kieffer&Stone irradiance model and should only be used from 0 to 90 degrees. The triangles are from our simulation code, and the red symbols are the LLAMAS data as hand-digitized from Figure 3 of the Holsclaw et al 2010 paper. The K&S model is wavelength dependent and has been interpolated to 282nm which is the wavelength of the LLAMAS data. Our own model has no wavelength dependence.

We see some disagreement on what the slope is as a function of phase for intermediary values, and large disagreement at large phases. The Hapke63 and LLAMAS values strongly disagree at phases from 120-170 degrees.

That is the phase-range in which we currently are limited to work (since the SKE does not work), so we have an issue related to data interpretation. I think what we see here explains why the Lambertian albedo we derive for Earth for lunar phases near New Moon can be above 1, as noted elsewhere in this blog. It is because the Hapke63 model gives a lunar irradiance that may be 2-3 times greater than the observed (LLAMAS) value for relevant phases. In the FFM method we scale the DS intensity to total flux in observations and synthetic (i.e. Hapke63-based) images and with synthetic irradiance exaggerated by so much the ratio to DS flux is low in synthetic images and we have to compensate by driving Lambertian albedo up.

So, we understand stuff!

For better use of the synthetic image code we need to implement a more realistic reflectance. Since LLAMAS data are only preliminary we have to wait for them to publish a full wavelength- and libration-dependent model. In the meantime we could perhaps roughly scale the Hapke63 reflectance – or something like that?

BBSO is somewhat ahead of us here, as they have derived their own reflectance model. We do not have resources to do that currently, although we have plenty of data to do it with. It is a suitable student project (PhD) or even work for a postdoc.

We have a result from night 2456003.

The 5 colored plot frames show the ratio of DS/total-flux for all filters using 4 different scattered-light reduction methods as a function of time. The two red curves are the BBSO method and our logarithmic variant of it. The two blue curves are the EFM and FFM methods. The solid curve in the frames with color is the value expected from a constant-albedo Lambert sphere. The last frame is showing the evolution of airmass with time.

The colored curves show the DS/total flux ratio divided by the mean of the BBSO method curve. The spread (+/- 10 %) is BIAS between the methods – they do not find the same answer.

Note that the slope of the observations match the expected theoretical slope.

The method with least scatter is the FFM, closely followed by the EFM. The two BBSO-like methods are not doing well, with our own logarithmic variation on the BBSO method doing least well.

Of interest is comparing VE1 to IRCUT since these filters are almost identical – the small wiggles in one set of curves do not re-appear in the other filters’ curves. This helps us understand what is geophysical signal and what is observational (and data-reduction method) noise.

VE2 is all over the place! This is probably related to the yet not understood intermittent focusing problem we have with that filter.

I find it very encouraging that the slope in the observations is so close to the theoretically predicted slope: The slope is there because during the observing period the phases of Moon and Earth slightly changed, which alters the illumination of the Moon by the Earth and Sun.

When we reduce more nights of data we will be able to begin to see if this slope is always there or changes sign with rising/setting Moon. If that is no problem we can begin to trust the signal we receive – even if the signal is biased by factors we do not yet understand.

This is a feather in our cap! We now know what nobody else knows: How much of the earthshine-based result is due to method and not geophysics. This will have consequences looking forward when enough data has accumulated and has to be interpreted in terms of climate processes.

I am sure that improvements in our analysis-methods will continue to aggregate and we may see the bias between the methods lessen.

We show results allowing a comparison of the EFM and FFM methods:

All observations for night 2456003 have been reduced to DS/(total flux) for both EFM and FFM method, and a direct albedo determination is also available from the FFM. The DS/total values are of course proprtional to terrestrial albedo and are plotted as crosses (FFM) and boxes (EFM) above. The line is the direct albedo estimate from the FFM method. Values have been normalized to a scale near 100.

The spread is near 5% but the standard deviation of the mean would be near 0.5-1%. If these numbers represented all data from one night we would be in line with BBSO claims, or slightly better.

Analysis of the BBSO method on these data awaits.

For the night 2456003 we now have both EFM and FFM results (and BBSO method results, to be analysed later).

We look at the mean difference and standard deviation of various quantities determined from the images by the two algorithms. The differences are calculated as percentages thus: 100*(FFM-EFM)/((EFM+FFM)*0.5). Data are for about 40 observations with a mix of filter repetitions – since these are differences the mixture of filters should not matter as same-filter images are compared.

alfa Mean: -0.3805 SD: 0.2027
a Mean: -66.64 SD: 0.6907
BS Mean: -0.0004525 SD: 0.0001109
TOTall Mean: -5.082 SD: 0.5449
TOTnoPed Mean: -0.01426 SD: 0.003551

DS23 Mean: 30.68 SD: 7.404
DS45 Mean: -7.415 SD: 9.557
x0 Mean: 0.9189 SD: 1.491
y0 Mean: 0.06020 SD: 0.3928
radius Mean: -11.19 SD: 1.802

The difference in offset a is imposed by the codes and can be ignored. The two totals are: TOTall=sum of pixel intensities in image, TOTnoPed=sum of image pixels minus the estimated bias and sky brightness pedestal (i.e. most close to ‘total original light in the image without sky’). The two DS values are calculated at 2/3 and 4/5ths of the radius on the DS of the disc. The difference is radius is due to user setting and is and error I can see now. Ups. More later!

We notice that disc center coordinates (x0,y0) have a bit of variance – the algorithm does not find the same lunar center and radius each time.

Fixing the codes so that radius is determined in the same way we rerun the analysis and get:

DS23 Mean: 21.70 SD: 5.806
DS45 Mean: 1.648 SD: 6.742
x0 Mean: 0.7158 SD: 1.765
y0 Mean: 0.1147 SD: 0.2926
radius Mean: 0.2838 SD: 1.524

The problem on radius is solved now. The other results are more or less the same – DS23 is significantly biased while DS45 is not (given the scatter).

What can we learn from this analysis? The error on DS and total flux is very small. DS45 is not biased, but the scatter is very large – that is, the EFM and FFM methods arrive at the same answer but both are subject to noise.

The ‘width’ of the halo we see is due to scattering in the atmosphere and the telescope. With data of the Moon rising it may be possible to extract information on how much of the halo is due to atmosphere and how much is due to the telescope. We first inspect data from night 2456003.

We do two things:

1) We run the FFM method on all aligned and coadded images thus getting values for alfa and airmass, and terrestrial albedo.

2) We calculate the fractional area in each image where the intensities fall in a specific fractile. We choose the fractiles between 0.1 and 0.3 % which turns out to be mainly in the BS halo. Doing this for all images gives us fractile area, or ‘halo size’, as well as airmass.

First we look at the contour plots from each image, in the B band. We have subtracted not only a bias frame but also a fitted linear surface – fitted to the corners of the bias reduced image.

We see 9 images of the Moon (each from the aligned and coadded images available at 9 different times in the B band on the given night).

We have extracted the area (counted in number of pixels) of the fractile between 0.1% and 0.3% (heavy dashed and next contour outwards) cumulated intensity. [Quite complex statistic that: 0.1% refers to the 99.9% brightest pixels of the image. 0.3% refers to ‘99.7% brightest pixels’.] We plot this, against other information extracted:

Upper Left: fractile area vs. time (Moon was rising). We see a smaller area as the Moon rises – less light is being scattered into that intensity fractile.
Upper right: Airmass and fractile area almost proportional.
Middle Left: FFM alfa vs time – the width of the PSF depends on where the Moon is.
Middle right: albedo vs airmass – the FFM-fitted albedo depends on the airmass: Not so good!
Lower left: albedo and alfa are dependent: Not so good
Lower right: fractile area and alfa are dependent: Not so good.

The last three results are not independent, of course.

We can learn a few things from this:

a) the FFM is not (yet) able to determine albedo etc independent of image quality – this is probably due to a degeneracy in the factors of the model we use – it seems that pedestal, alfa and albedo are coupled.

b) the nice relationship between fractile area and airmass, and alfa and fractile area promises an equally nice relationship between alfa and airmass. By linear regression we find:

alfa = 1.76 – 0.0036*airmass

which suggests that alfa (in the B band) is 1.76 at zero airmass – i.e. the contribution of the telescope in the B band to the PSF is a power law with power 1.76.

For all bands we get:
———————–
alfa_0 error
———————–
B 1.757 0.0010
V 1.763 0.0014
VE1 1.757 0.0006
VE2 1.745 0.0002 (one outlier discarded)
IRC 1.758 0.0004
———————–

Within +/- 0.001 it is possible to see no real color-dependence in the alfa_0 parameter. A consistent value for almost all bands is alfa_0=1.76. The VE2 PSF may be slightly broader with alfa_0 at 1.75. VE1 and IRCUT (which are almost identical filters) have very similar alfa_0 and errors.

We have on some nights by other means found alfa near 1.8.

The problem that albedo depends on the PSF power has to do with our model definition and fitting techniques and is open to improvement by adjusting fitting weights and so on. While it is too bad we are not ‘there’ yet with the FFM we see a tool for detecting when we succeed – namely, when the albedo no longer depends on alfa, and this information will be put to use.

From night 2456000 to 2456007 we have excellent data. Here they are (sorry about the sideways plot) for each band and each night. I show the ratio of the DS intensity in EFM-corrected images to the total image flux. Alfa was always near 1.7 so flux is not expected to be lost off the edge of the image. We have scatter at the 1% level and up.

We have a fairly close relation between observed and expected ratio if we chose the simple Schönberg phase law (blue curve), while we seem to be contradicting the more advanced Lommel-Seeliger phase law (red curve). Hmm? Moon and Earth modeled as Lambert spheres in this argument.

The larger scatter on the first and fifth nights in some filters indicate observing problems and not geophysics!
I think issues in the EFM method are: centering of the Moon in order to extract DS intensity; and choosing the pixels to use as source. On the first night problems are due to very low earthshine at this phase.

The phases or Sun-Earth-Moon angles covered by the 7 nights extend from 97 to 34 degrees. Unfortunately skipping the interesting 42-degree angle. Except for the brevity of the observing window on the last night I see no strong indication that we are getting particularly bad data at high airmass – some of these data points are at airmass 3-7!

Now then – the above to be repeated with the BBSO and FFM methods: all for the poster at the EGU.

When we use the BBSO or EFM methods for removing scattered light we end up with a number for the DS flux – usually in a box on the DS.

When we apply the FFM we end up with a number for the terrestrial albedo – based on all pixels in the image.

The precisions of the two approaches is not the same. In a test on synthetic images with added, realistic, noise, and realistic variability in the offset and PSF-parameter alfa we find that the

S.D. of the determined albedo is 0.12%, while the
mean error on DS intensity is 1.2 %.

This difference was seen in the Figures of the paper we submitted. The present results are a re-calculation of those simulations – it looked too good to be true in the paper, but I can conform that it is indeed the case. I speculate that it is due chiefly to more pixels being used to determine the albedo in FFM than the DS flux in the EFM.

Further calculation of the albedo from DS intensity would carry the 10-times larger error forward to the albedo!

“Eureqa” is software that searches for regression formulae that fit given data in the most parsimonious way – it not only tests which regressors are the best, but also searches a large space of possible functional forms. We have used Eureqa to hunt for solutions for SSA (i.e. albedo) in sets of data extracted from a large number of synthetic images.

A training set was generated from the averages in 64 equal-size boxes laid out over synthetic lunar images.
Each synthetic image was generated for the same JD, but each had its own realization of Poisson noise. One or 30 of these frames were averaged, and a set of 1000 such images were generated with random values of alfa, pedestal and albedo sampled from uniform distributions.

The distribution limits are:
alfa 1.5 to 1.85
ped 0.0 to 50.0
albedo 0.1 to 0.9

Eureqa was run until no further progress seemed possible.

Results:

formula from 30 frame averages:
“SSA = (199,9*V29 – 199,3*V48)/
(V53 + V47*V47 + 399,8*V29*V47*V47*V47 – V1*V29 – V1*V47 – 398,5*V47*V47*V47*V48)”

formula from single frames:
“SSA = 60,55*V28 + 171,8*V11*V28 – 0,006181 – 60,62*V32 – 159,8*V40*V40”

The VXX refer to the mean of the XXth box on the lunar image – counting from the corner and along rows, then columns.

Statistics 30 frames 1 frame
————————————————————–
“R^2 Goodness of Fit” 0.9999703 0.99935872
“Correlation Coefficient” 0.99998515 0.99967931
“Maximum Error” 0.0055585436 0.027577193
“Mean Squared Error” 1.6123362e-6 3.5070792e-5
“Mean Absolute Error” 0.00093508038 0.0043065976
“Albedo error” 0.32% 1.4%
————————————————————–

We see that averaging 30 frames gives better results than using single frames, as expected: We see the maximum error better by a factor of 5, the mean square error by 22, the mean absolute error by 5. Thus, the error in the fit scales about as you would expect – by the
SNR.

We see albedo errors of 0.3 – 1.4%. These are calculated using the Mean Absolute Error (i.e. best case choice – not conservative). These are all higher than our science goal of 0.1% so we must use more frames.

Discussion.

How would one use a system based on the above formulae?

Each formula is trained on centered images from a given point in time. It can be used only on images also centered and with the same image scale and from the same time. A training set can be generated for any point in time, and it is possible to ensure that the image scale is the same as for the observed image. Observed images may not be centered, however. It is probably more difficult to shift the observed image, because of edge
effects, than it is to determine the shift and then shift all synthetic images equivalently and train the formula on the shifted synthetic images.

While the use of a simple formula such as the one above would be extremely fast, it does not appear to offer an advantage in terms of error over our other methods.

The method is somewhat similar to Chris’ principal components method, but governs itself in finding the optimum boxes – i.e. areas on the lunar image to use.

The method quickly finds similar-quality solutions if started from a fresh set of images, but the formulae may differ in choice of detected variables. This implies that the possible-solution space has not been sufficiently sampled by the 1000 image training-set.

Generating larger training sets is possible – generating 1000 takes about 20 minutes, though – so it is not realistic to generate sets that are e.g. 100 times larger.

Training the formula takes 10-20 minutes using 1000 images and 2 CPUs on an ordinary PC workstation. Perhaps there is scope for much larger sets using the CRAY and multiple processors?

Our analysis of images is hampered by their numbers and size. Many of the pixels in each image represent sky or otherwise similar areas: In other words – the ‘degrees of freedom’ of the image is probably less (much less) than the number of pixels in the image. In order to explore the potential of methods that operate not in the image plane but in various transform planes I investigated just how many singular values are needed to satisfactorily reconstruct a typical lunar image from an SVD – a singular value decomposition. First I just tried to reconstruct an image completely, using all 512 singular values:

At upper left is the original (equalized) image. At upper right is the reconstructed image, and at lower left is the difference image. At 512 singular values the RMSE in % was very low – 0.00219 %. The RMSE fell below 0.1% at 420 singular values. The usual plot of the singular values sorted by size shows a ‘break’ near 100 singular values, suggesting that by using more than 100 SVs it is mainly noise that is being reconstructed. Inspection of the spatial reconstruction error as function of increasing number of SVs showed that at 80 SVs the added effort goes mainly into reconstructing the sky – i.e. noise.

I have re-fitted the synthetic images with the FFM method, but this time used weighted least squares to perform the fitting. The weight was 1/sigma in each pixel where sigma is the Poisson SD. The previous bias that was seen for single images appears to be gone now. The scatter for single images is also much smaller.

Here are the results of a test of the FFM method similar to the previous
one (below) but now based on 100-image stacks. The improvements are
appreciable: There is almost no bias (5e-5 or 0.02%) in the determination of albedo (top panel), the scatter is less (SD=0.00034 or 0.11%) and the range of useful lunar phases seems to extend to almost +/-120 degrees. The mean bias is 5 times less than the SD so in effect we have no bias now and are determining albedo to within our science goal. Dotted line shows +/-1% from true albedo value in top panel.

Note that the error bars are unrealistic – they are merely the formal values that the fitting method gives under the assumption that all data are statistically independent. We can get empirical error bars instead by doing some Monte Carlo or bootstrapping on a few of the results.

Test for alfa=1.6 remains to be done. More later!

For a full month of single images (alfa=1.6) we show how well the FFM does at recovering the terrestrial albedo from the noisy synthetic images.

The crosses show the solutions for the single scattering albedo used in the Earth-component of the model that generates the synthetic lunar images. The dashed line shows the true value (0.297). Bars on the points show the LM-estimated 1 SD. The two dotted lines show the 1% upper and lower levels from the true value. New Moon is at phase 0. The simulations were performed on single realistic image calculated using alfa=1.6.

The solutions are clearly biased from to the true value and the bias is significant at the several sigma level, but the bias is small compared to other methods’, such as our own implementation of the BBSO method. Solutions are inside the +/-1 % limits for phases roughly between -90 and +90 degrees fromNew Moon. BBSO quotes 1% levels of uncertainty on albedo for a single night. We get 1/3 % in single images.

If the realistic limits on lunar observability is +/- 40 degrees from New Moon we have significantly extended the range available with the FFM method, compared to the logBBSO method applied to stacks of 500 images which was seen to be more than 1% biased outside about +/- 45 degrees!

Here is a 3-panel figure showing additonally the detected values for offset and alfa:

We see some breakdown of the method near New Moon – but NM cannot be observed anyway so that is OK.
I think that the bias in alfa is causing the bias in albedo – so if we could fix that …

Similar plot for alfa=1.8 (incomplete – it crashed! Gahh!)

The Full Forward Method for determining not only earthshine intensity but also terrestrial albedo has been tested on, so far, a singe synthetic lunar image corresponding to a few days after New Moon.

By minimizing on the difference between the ‘observed image’ and trial models from Hans’ synthetic image generator over the whole image frame, for three values of PSF width, we can retrieve the terrestrial albedo to within 0.07% for alfa=1.6, 1.7 and 1.8.

This is below our science requirement, and gives confidence!

The model must be tested for images corresponding to a whole month now. The calculation is time-consuming since it is not programmed for the DMI Cray yet, so meanwhile the method runs in a mixture of Fortran and IDL on a workstation.

The fitting method can handle offsets in the image so to some extent an incorrect bias subtraction should be subsumed by this analysis method.

The method relies on Chris’ insight that a library of lunar images for a given phase but different terrestrial albedos can be generated from the linear sum of an ‘albedo is 0’ image and an ‘albedo is 1’ image. The weighting function between these extremes is simply the albedo and is one of the fitting parameters.

Testing remains also for how dependent we are on knowing the PSF – i.e. not just the exponent but also things like non-power-law shape.

This result brings to an end many weeks of worries I had that we would not in time be able to show that we have a method that is good, and now – after some more calculations overnight – we should be on line to write the results section of our paper!

Wahooo!

Some more analysis of scattered light.

I take three synthetic lunar images —

1) a thin crescent with no ES, i.e. BS only

2) ES only

3) BS + ES

Raw images shown on the left — with scattered light on the right. (All images shown at the same display levels, making the halo hard to see!)

The albedo adopted was 0.30 for the ES on (2) and (3).

All three are generated with the correct flux scaling for the components
i.e. image (3) is the sum of (1) and (2).

These three images are convolved with the PSF with alpha=1.8.
No noise is included in the final output images.

I then compute the amount of light which has scattered off the
original synthetic images and onto new pixels. (A list of all pixels which
had non-zero counts in them in the original images was kept and the
counts in these pixels in the convolved images compared to the total
counts in the convolved images. The normalisation of the convolution is set
to the total counts in the original images).

While doing this, I checked how much light scatters off the edge of the 512×512
frames and into the larger (3*512)x(3*512) frames used for the FFT. This turns out
to be ~0.04%, a good deal less than our target accuracy (0.1%)) for alpha=1.8, so this is not a major source of error.

The amount of off-moon scattered light is as follows

1) BS only : 4.2 percent
2) ES only : 0.58 percent
3) BS and ES : 3.9 percent

let’s call this the “halo light”.

these numbers highlight why getting the BS to ES normalisation correct
when generating forward scattering fits to data is important!

We want to measure the DS to BS ratio to an accuracy of 0.1 percent, so the approximately 4 percent of halo light has to be gotten right when figuring out the amount of BS light. The difference between BS only halo light and BS+ES halo light is of order 0.3 percent, 3 times more than our desired accuracy. This shows that the ES makes a non-negligible contribution to the halo light.

Those numbers were for alpha=1.8. For alpha=1.6, I get

1) BS only : 6.9 percent
2) ES only : 1.58 percent
3) BS and ES : 6.0 percent

Now the ES is making a very significant contribution to the halo light — 16 times more than the accuracy with which we would like to know the BS light.

The light lost off the 512×512 frame is small at 0.04% – half of our target error.
However, for alpha = 1.6, this lost light increases to 0.3 percent of the total.

Lost light off 512×512 array for centered thin crescent Moon

alpha lostlight
2.0 0.004 %
1.8 0.04 %
1.7 0.11 %
1.6 0.32 %

It looks like we need to take into account lost light off the 512×512 array
when we are doing the fitting for the BS, especially if the images aren’t well centered, as is presently an issue!

More added

In response to Peter’s comments:

We have both done tests on artificial data which suggest that we are able to measure alpha with an accuracy of +/- 0.01. I figured out the corresponding missing light off the edges of the image for this uncertainty in alpha:

alpha lostlight
1.80 +/- 0.01 0.038 -/+ 0.004 %
1.70 +/- 0.01 0.11 -/+ 0.01 %
1.60 +/- 0.01 0.32 -/+ 0.03 %

So if we really have alpha=1.8 and can measure it with an accuracy of 0.01, we are pretty safe — at least for this thin crescent which is well centered.

Peter asked: what if extinction and alpha are correlated? Great idea! We should certainly look into this!

I had another idea — what about if we made longer exposures of the moon which would saturate the moon but sample the halo light well all the way to the edge of the frame? Could we then measure alpha with appropriate accuracy where it counts — at the moment we are merely extrapolating a not very well sampled halo beyond the edges of the frame to estimate the missing light (in 30 ms exposures, the halo light is very close to the bias level at the edges of the frame, so is not well measured out there, and a bias error could lead to systematics in the estimate of alpha, perhaps larger than the error estimate adopted above of 0.01, which was obtained for rather idealised synthetic data.

The behaviour of the scattered light around the moon beyond the edges of the frame is still very much an open question — we didn’t acquire the data at the last full moon to settle this, but just got more questions!

More still: position of the Moon

Tested the effects of the positioning of the moon using the same setup as above

For 1 pixel shift of Moon in the x-axis, I got negligible changes in the amount of lost light.

For a 10 pixel (circa 70 arcseconds) shift (to the left along the X-axis for a crescent moon which is illuminated on the left side):

alpha lostlight
2.0 0.005 %
1.8 0.04 %
1.7 0.12 %
1.6 0.33 %

The changes to the fractional lost light are very small. So the good news here is that we are not very sensitive to the precise position of the moon if it’s near the center.

In order to better understand the various data-reduction methods, I test here a modified forward model. It consists of a test on synthetic images – this allows us to test what the error is relative to the ‘truth image’.

A synthetic ‘observed image’ is generated as the convolution of an ideal image with Poisson noise and an intensity offset as in a pedestal. This pedestal models the effect of small error sin bias subtraction.

Trial images are next generated and compared to the observed image – when a good match is found the winner model image is declared. We estimate errors in a box on the DS. We use a grid search in alfa for the PSF and in offsets. We CALCULATE a factor that scales the intensity so that once the offset and the PSF are chosen the resulting image is scaled so that the total flux is the same as in the observed image.

We used 0 bias or offset and we used 0.01 as the offset. We used alfa =1.6, 1.7 and 1.8. We allowed loops that bracketed these values.

The test was based on the RMSE between the observed image and the trial images. For the winning model image the error between the known ideal image intensity inside the box and the intensity inside the same box in the offset and scaled (but unconvoluted) image was calculated.

We found errors in the 0.0 to 0.75% range, seemingly independent of the value of alfa or offset.

The above was for one lunar phase near new Moon. We repeated the test on an image with a phase larger (like quarter moon). We found errors from 0 yo 1.85% with a tendency to be large for small alfa – i.e. broad PSFs.

The above shoul dbe repeated for all phases.

The results are similar to what we have found for the EFM method – but for that method results are available for all phases, so – … (more later)

On the above we note that the best solution found seemed to have biases in the alfas and the offset – nevertheless the best solution was fairly good, as shown.

So: the above test was witha method that unlike the EFM did not constrain the trial images to ‘fit through the sky’, but it did, like the EFM, constrain the total flux.

Here are some results from a sensitivity-tst of one of our scattered-light-removal techniques. The technique is a forward convolution method but is not the EFM method, but shares elements from that method.

The EFM method consists of these steps:

0) An observed image is at hand.

1) From the observed image an idealized ‘source’ is constructed from all the bright pixels of the image.

2) A model observed image is calculated on the basis of trial values of alfa (the power-law exponent in the PSF) and the source image in 1 by FFT convolution

3) The model image is now shifted and scaled in intensity so that the halo matches the sky part on the DS and so that the brightest part on the BS disc matches.

4) the RMSE between the model image and the observed image is calculated on an area of the image that includes large parts of the sky around the Moon.

5) The trial image that has the smallest RMSE ‘wins’

6) That winning image is subtracted from the observed image. Since the model image did not contain any DS pixels (see step 1) the difference between the observed image and the model image should be the DS.

That was the EFM method. The method I test now is based on the above but omits step 3. In the EFM method step 3 was included in order to ensure flux conservation while making the model ‘fit the sky’ (this is reminiscent of the BBSO method where the sky if fitted ensuring that the model is anchored in the actual sky level thus making the residuals there 0).

We show two plots. Each plot has 3 panels – they are:
Panel 1: The Moon image being analysed.
Panel 2: A slice through the lunar disc with three graphs – the red one shows the known solution (the tests here are performed on synthetic images where a halo is applied with a convolution, using aflfa=1.8, to simulate effects of atmosphere and optics); the black one is the ‘observed’ image in step 0 above; the blue curve is from the image in step 2.
Panel 3: The black curve is the ideal image showing the structure on the DS, the thick green curve is the difference between the blue and the black curve in Panel 2 – i.e. it should be the model estimate of the DS.
Panel 4: This is the percentage difference between the two curves in panel 3.

The difference between the panels on the left and the right is that on the left a small value has been subtracted from the image in step 0 to simulate an uncertainty in bias subtraction; on the right no such value has been subtracted. The value subtracted on the left is 0.01.

A typical uncertainty in bias level may be a number of that magnitude – the consequence is an error in the estimate of the DS flux of 5-10 %.

Henriette’s estimates of how well we can subtract bias levels can be used here to introduce more realistic values of the error (0.01 may be too large, we hope).

The above shows the sensitivity in a scattered-light-removal technique to small errors in bias subtraction. This was NOT the EFM method which anchors itself in the DS sky, but a similar analysis is required for that method before we can present it in a paper.

Note how the slopes of the DS sky in the above plots do not match, despite the same alfa being used in the generation of the test image and the model image. Remember that in a lin-log plot like this a power law is not linear. Inspection of the blue curve in panels 2 above shows that it does indeed curve – the different slopes on the DS sky may be due to the proximity of the DS flux itself (which is also convolved and contributes to the DS sky halo) in the black curve. In other words: the halo on the DS sky in a real image may not be entirely due to the scattered BS, which we have assumed in the above.

This suggests that we should test methods based on convolving entire synthetic images of the Moon for the modelling part – make up the model of more than just the BS as in step 1. If we make up the models of BS+DS images we will have to enter the intensity of the DS, or directly, the terrestrial albedo, as a parameter in the grid search.

This is actually a feature of ‘Chris Method’ that we have seen in this blog.

This is our next step, but first we want to know if the above simpler method (just loop over values of alfa, not alfa and albedo) can be shown to work better if we do as we actually do in the EFM method – namely, we ‘force a fit’ on the DS sky. Since the preliminary trials of the EFM has shown up some not-understood problems of their own we are still working on the matter.

In our tests of various scattered-light removal techniques we are using synthetic images to which various realistic things have been done – smearing by a PSF and ‘adding’ Poisson noise, etc.

Poisson noise should not be ‘added’ because this will increase the mean value of a region of an image. Instead, Poissonian noise should be drawn from the Poisson distribution with a mean at the value of the image pixel under consideration.

We have two codes to generate Poisson noise in images and we test these next. We take one synthetic (i.e. noise-free) image and define a region on it. In that region we calculaet the mean. We then generate Poisson-noisy images using our IDL code and a Fortran code. These are different in that the IDL code writes out floating-point images while the Fortran code writes out integer images.

In a loop over 1000 trials we generate noisy images from the original image and each tine calculate the mean in the predefined region. We set aside the percentage difference between the IDL code and the ideal image and the Fortran code and the ideal image.

Then we plot these percentage differences:

The two distributions are symmetric about 0 and practically identical. From this we conclude that both codes ‘add’ Poisson noise to the image correctly and that the writing-out of the integer resulting image in the case of the Fortran code does not cause any rundoff problems, even if the region inside of which we tested had a small mean value (about 3.9 pcounts). The regiion was 19×19 pixels large.

We do note that the spread, although similar, is quite large – SD is 4 and 6 counts respectively. This implies, unsurprisingly, that the number of boxes of size 19×19 that must be averaged over is in the hundreds if the standard error of the mean is to be at our science goal of 0.1%. We knew that – this is just a reminder …

Final note: the Fortran code uses a fixed seed to generate random numbers. A modification was needed to read in a random seed instead, from the command line. This is supplied when the code is called by the user.

Peter and Chris took images of the sky near the moon, going outwards from the Moon in 0.5 degree steps, in order to measure the powerlaw falloff of light using the MLO telescope.

Data taken on night of 2012-01-08.

We got good data for 4 positions, 1.0,1.5, 2.0 and 2.5 degrees from the moon,

WCS coordinates for the frames were found using astrometry.net. This is a great site!

Stars in the fields were found by searching Simbad for objects brighter than 10th magnitude. A small program was written to locate these stars and do the photometry based on the magnitude equations determined previously for the CCD, and their J2000 coordinates. This is quite nicely automated now. Bias frames were subtracted from each of the images and the magnitudes and sky flux of the known stars in the field computed. The magnitudes were compared to the known magnitudes and came out very nicely.

The images were taken at
POS3 : 06:27:24.313 +22:22:27.07
POS4 : 06:27:25.994 +22:51:53.62
POS5 : 06:27:23.910 +23:22:22.01
POS6 : 06:27:26.336 +23:51:34.44

Some example images with stars of known magnitudes:
POS6:

POS 3:

The Moon during all of this was at 06 27 36 +21 50 28.

Distances from the center of the moon were computed for each star, and the sky flux next to that star plotted as a function of Lunar distance on a log-log plot.

PDF version of this plot:

Upper plot shows the magnitudes of the stars found in all four fields versus their correct magnitude. Lines is the 1:1 relation.

Lower plot shows the brightess of the sky versus the distance from the Lunar center in log(arcsec).

Line is a powerlaw fit, and has a slope of -0.6.

This is very curious. We were expecting -1.5 to -2 or so. What is this extra light?

The surface brightness is about 16 mag/sqarcsec in the inner most position at 1 degree from the moon, dropping to about 17.5 mag/sqarcsec at 2.5 degrees from the moon. This is a very slow drop in the brightness. The Moon was in the Milky Way plane, but the luminosity of the Milky Way is about 20-21 mag/square arcsec. So it’s not the Milky Way we are seeing.

For the present this is a cosmic puzzle!

UPDATE

Improved coordinates of the Moon were computed from the times of the exposures by Peter and a better plot results. Here it is:

Slope of the halo is now -0.9. That’s a bit better. More work to be done on this at the next full moon.

On the night JD245524 conditions were good: Moon was setting in clear dry skies and CoAdd mode was possible. We obtained images in all bands, and present the reduction to DS/BS ratios using the EFM method here.

Scaled superbias was subtracted from all images. No flat-fields were applied. Images that were poorly centred in the frame were removed. Images for which the estimated PSF exponent was above 2 and all negative values of DS intensity were removed. On the remaining 171 images DS/BS ratios were extracted using two patches on the DS disc – near the rim and further from the rim: (2/3 and 4/5 times radius away from the centre of the disc).

We plot the PSF exponent determined and the B,V,VE1,VE2 and IRCUT bands and DS/BS ratios below:

We see that the exponents of the PSF (upper frame) are quite stable for each filter, but that the VE2 filter has a smaller best value (i.e. PSF is ‘thicker’ in VE2).

For the two patches on the DS disc (close to rim middle panel; further from rim bottom panel) we see quite stable DS/BS ratio values but that the value does depend on where the DS patch is located. The V filter is represented only by data near fractional JD = 0.77. We see a large scatter there: this was when the Moon passed directly behind the MLO antenna, which interfered with the images and extraction of the photometry. The considerable scatter near 0.77 should therefore be discregarded – the data should be ignored at that point. This also removes all V data from futher consideration.

For the remaining data we look at stability of mean values.
VE1 mean: 0.0083 +/- SD 0.0006 Zm: 77.8037
VE2 mean: 0.0073 +/- SD 0.0005 Zm: 67.6726
B mean: 0.0109 +/- SD 0.0006 Zm: 127.6579
V omitted……………………………………….
IRCUT mean: 0.0081 +/- SD 0.0005 Zm: 108.8674

Here, SD is the standard deviation of the data for a band while Zm is the ratio of the listed mean and the listed SD divided by the sqrt of the number of points minus 1 – i.e. the ratio of the mean and the SD of the mean. The number of points in the different bands are 33, 25, 60, 11 and 51 for VE1, VE2, B, V and IRCUT respectively.

Now, these SNRs need to be MUCH bigger to reach our Science Goals of SNR=1000. If the errors are randomly distributed and independent we need 60 times more data-points in the best case (B).

We need to discuss why we have so few data :

a) There are only about 150 good data points (i.e. images). We actually had 210 observations. This number was reduced by the need for images witha well-centred lunar image: This requirement comes from needing similar images to give the EFM method a level basis of data to work on – if the Moon is sometimes in one side and sometims in the other side of the image frame the weight given to thedata during the location of the best fitting PSF varies and leads to systematic effects – i.e. bias. So by centering the Moon better we can have more data to work with.

b) By far the largest factor limiting the number of images is, at present, the faulty shutter: We may have taken more than 1000 images on that night, but most of them failed because the shutter stayed shut. This is despite some improvement received when the Mylar blanket was installed.

c) The amount of data is also dependent on the time interval available: it necessarily starts at sundown for these types of CoAdd images and stretches until the Moon sets – perhaps 2-3 hours at best. The present interval is only 1.2 hours long. If time was not wasted by fiddling with the dome control time could be saved.

d) As we observe in a cycle with equal numbers of images assigned to each filter we are limited by the length of the longest exposures. Change of this strategy could give us a larger number in some filters. Often it is the B filter that takes the longest to get by a factor of almost 10, compared to the broader and redder filters.

The above are the main reasons for the limited number of images we get on even a good night. In order theya re:

1) Wasted time due to shutter failing.

2) Limited observability of Moon in CoAdd mode.

3) Spending too much time on long exposures in B rather than getting many short ones in VE1 etc.

4) Wasted time due to dome and precise telescope pointing issues.

Of these there may be little to do presently about #1, and nothing to do about #2, but #3 should be considered and experimented with and #4 may be possible to do something about soon. I estimate that if 1 could be adressed we could have 5-10 times more data, while #3 and 4 could give us 50-100% more data if adressed.

Thus it seems we will likely not get the 60 times more data we estimate may be needed to reach SNR of 1000 by averaging short-exposure CoAdd mode images alone. This means we really need the Lund mode operating, and the SKEs. Perhaps that is where the investment of time is most valuable?

………….added later:

Here is the above plot with the ‘antenna interference’ points removed:

Note how well the VE1 and IRCUT (same filter, almost) line up.

I looked at a bright white LED instead of a bright green laser.

Imaged it with the Canon 35 mm camera as before.

The result is exactly the same scattered light power law around the source.

Did this to check that nothing funky is going on with using a laser as a point source!

I added a scattering plane to the setup – a translucent tupperware lid was placed directly in front of the camera lens. Imaged the LED with and without the scattering plane. The lid scattered light very nicely (left) versus naked LED (right)!

The red curve shows the scattered light profile, the black curve the naked LED. Light is scattered on all scales in this experiment. Can one have a regime where scattering in the optics is dominant and elsewhere where atmospheric scattering is dominant? Not sure at all about this — my feeling is if the atmospheric scattering has a shallower power law index than the optics scattering, then it will always win. Could scattering be scale dependent?

The power law index for the naked LED is -2.6, for the scattered version it’s -1.5. This latter value is very similar to what we get for the moon in the 35 mm camera shooting through air! The tupperware lid simulates the Earth’s atmosphere, presumably by accident.

Scattering through camera optics, air and tupperware lids all produce power law profiles. Is that not interesting? Scattered light follows a power law under a wide range of conditions!

Observed the full moon with a Canon 35 mm camera, eclipsed by a “sharp edge” (the front balcony of our house) and then moved the camera so that the Moon appeared from behind the edge and imaged it again. The moon was about 45 degrees above the horizon on a very clear sky. Sector was chosen to avoid the gum trees seen lower left

Image above is the UNECLIPSED moon close to but detached from the edge of the house. The radial profile was measured in the inner region of the sector shown, with the outer region defining the background sky for subtraction.

Image above shows the light from the moon, now placed just behind the roof’s edge. The radial profile was measured in a similar manner as before.

The idea is to test if the scattered light far from the moon is the same in both cases, i.e. is mainly due to the atmosphere.

The result is that it does indeed seem to be due to the atmosphere — the light falls off as a log-log power law exactly the same and, for two images with the same exposure time, the falloff has the same surface brightness on the sky.

The falloff at large R (> log(R)=2.1, or about 125 pixels (1 pixel is about 1 arcminute) is exactly the same in the two cases, going out to more than 1000 arcmin (15 degrees). There is some funny behaviour near the knife edge, with some missing light in the eclipsed profile relative to the uneclipsed one, but it was very difficult to estimate where the center of the moon was behind the knife edge — since the moon is eclipsed — and the behaviour is quite sensitive to the moon’s position. Tried to figure out where the moon was using stars in the background, but the field distortions in the 35mm lens were too severe to model with a handful of stars.

If we try this again it would be better to leave the camera fixed and wait for the
Earth to turn and eclipse the moon, rather than moving the camera. A sequence of images as the moon approaches the knife edge might make modeling where the moon ends up behind it easier.

UPDATE: The measurements on this post have now been superseded but more reliable data from NGC6633. See blog entry for July 1st, 2012!

On 23/12/2011 Peter and Chris observed M41, an open cluster in the Milky Way plane at 06:46, -20:40. It was a very clear night and we got good images (although off center) in all bands. Image below is a centered B band frame.

Stars identified from the WEBDA database for this cluster are shown below:

Data for these stars are here:

http://www.univie.ac.at/webda/cgi-bin/frame_list.cgi?ngc2287

The transformation to standard V and B are rather good. There is a small
colour dependence in the V filter (-0.08 mag/mag) and quite a large one in the B filter (+0.28 mag/mag). The scatter is quite low in the V transformation — 0.02 mag, and a little larger in B — 0.05 mag.

Transformations

Instrumental magnitudes convert to true magnitudes as follows:

* compute instrumental magnitudes
V0 = -2.5*log10(ADU/sec)
B0 = -2.5*log10(ADU/sec)

* instrumental colour
(B-V)_0 = B0-V0

* transform instrumental colour to Johnson B-V
B-V = -1.056293 + 1.529783*(B-V)_0

* transformation to Johnson V (scatter is a very small 0.02 mag)
V = 15.15 + V0 – 0.08*(B-V)

* transformation to Johnson B (scatter is quite large at 0.05 mag)
B = 14.46 + B0 + 0.26*(B-V)

These are for a (rather small) 2.5 pixel aperture – designed to avoid crowding by nearby stars. The correction for the aperture size is ~0.1 mag brighter.

There are no standard magnitudes for VE1, VE2 or IRCUT. I have searched for
similarities with V, B or I, since these are the bands we have external magnitudes for.

The following definition of VE1 magnitude gives a very close 1:1 match between V and VE1:

VE1 = 16.34 – 2.5*log10(ADUs/sec) + 0.06*(B-V) (scatter 0.04 mag)

so VE1 is functionally quite similar to V.

On the other hand, VE2 is very similar to I. With this definition for the VE2 magnitude:

VE2 = 14.22 – 2.5*log10(ADUs/sec)

one gets very close to a 1:1 relation with I, with a scatter of only 0.03 mag. No colour term either — so VE2 this is an excellent match to I.

We didn’t get any good frames in the IRCUT filter. We’ll have to try that some other time!

Summary : Our V and B are much like their standard Johnson counterparts.
Our VE1 is functionally more or less the same as Johnson V, and VE2 is very close to Cousins I band. The latter seems a bit odd. I haven’t seen the filter curves, so they’ll be fascinating in this context. What is the IRCUT filter going to show? More to follow.

Applying also the logarithmic BBSO method to the rainbow images, we get the following results:

We see much more scatter in the early part of the sequence, compared to the ordinary BBSO method. This is both the part where we had the lowest fluxes (shortest exposure times) and (we presume) a proximity to the ‘rainbow phase’ – so are we sure if we see scatter due to rainbowing or just low exposures?

As earlier we expect that scatter rises when erroneous additive levels are removed – but do we have a way to decide if ‘too much’ has been removed? The log method is able to ‘remove too much’ while the linear method is more likely to underestimate how much is to be removed.

The larger total/DS ratios in the early part indicate lower Earthshine intensity. Remember that the rainbow is supposed to be brighter than ‘non-rainbow’ Eartshine.

It is probably realistic to assume that problems have arisen in the early part of the sequence due to low exposure times, not because of low actual fluxes.

Shall we consider co-adding frames into bins – every connected set, for instance?

The Moon was observed going through the ‘Rainbow angle’ on Dec 21 2011. Observation started at Moonrise and stopped at sunrise a few hours later. All filters were in use and good exposures were obtained with minimal dropouts due to shutter stick.

Bias was subtracted and the scattered light removed using the BBSO method. The total counts for a 23×23 pixel box on the DS was extracted in each image as was the image total flux. The ratio of the total flux to DS box is plotted in the figures above in the upper panels. DS box and total are plotted below. To the left are the results when scattered light has been removed; to the right when it hasn’t.

The step up in total flux near the start is due to adjustment of the exposure time used.

The increase in scatter that appears to occur when scattered light is removed may be quite realistic as the total/DS ratio will be more stable when the DS total contains an (erroneous) contribution from as-yet unremoved bias or sky levels. When these contributions are removed the true scatter of the DS data becomes evident. Note that, e.g. the cyan points lie on a straighter sequence in the corrected image than in the uncorrected implying variability in the scattered light, which was quite realistic for that night.

The observed downward trend in the total/DS ratio is consistent with a gradual increase in earthshine intensity.

The sequence of observations started as the Rainbow angle had already been met. The Moon moved from red to blue side of rainbow.

We next add the expected evolution of the Earthshine intensity for this night. We use the synthetic model and generate models for the ‘cloud-free’ case – i.e. The Earth has no clouds but does have surfaces with different albedo: land, ice, water – and the ‘uniform’ case whjere Earth is just a Lambert sphere with a phase function.

The solid line shows the evolution of Earthshine intensity for a Lambert sphere; the dashed line shows for a cloud-free Earth. The two vertical dotted lines delineate the observing window. The large dip correspnds to the (cloud-free) Pacific ocean.

We see that Earthshine is predicted to increase during the observing session if the phase-function rules – but decrease if the surface albedo rules.

At the observing time the Moon was at -19 degrees declination so the sunglint was also southernly (Sun was near -23 degrees in decl). The sunglint position does not provide all the earthshine but is at the centre of the area that contributes. For the observing time we see the sunglint moving off South America and into the South-Eastern Pacific:

We can use the modelling code to make a map of the areas that contributed to the earthshine at the moment of observation.

This map corresponds to UTC 1 hour on that day. As the sunglint moves further West the dark Pacific comes to dominate the average albedo.

It is quite possible thatthe Pacific was not cloud free at these times and therefore helped to increase the average albedo. From http://goes.gsfc.nasa.gov/goeswest/fulldisk/vis/ we can download full disc images from the GOES West satellite at 3 hour intervals.

USA is at upper right. South-Eastern apcific at lower right – Sun is just rising on Hawaii ner athe middle of the image. While therea re some large swirls of clouds in the South Pacific it does not seem to be enough to cause the growth in average albedo observed. A quantiative calculation would have to be performed.

Summary: We have observed a decrease in the total(DS ratio which implies a rise in the DS intensity – i.e. increasing Earthshine. We have shown that Earth’s phase function was growing at the time and thus contributing more light to the DS of the Moon, but we see that the rotation of the Earth is bringing darker areas into view thus counteracting the phase-function tendency.

The noise model is improving rapidly. Below is a simulation of a lunar image (right) from a synthetic moon provided by Peter; and the corresponding real image (left).

Images are equalised histograms with the same scaling of log(flux).

See www.astro.utu.fi/~cflynn/logbetter.jpg

Until now I had been using Gaussian sampling for the noise, but this fails badly when only a few photons are expected in a given pixel — while this is not an issue for most of the moon (ES>10 counts/pixel in typical exposures, and BS is of course very bright) — it matters a lot when modeling the halo far from the moon.

Replaced the Gaussian noise model with a combination of a Poisson model for low count rates and Gaussian for high count rates (>25 counts per pixel). The Poisson code is both very slow and inaccurate at high count rates, so a combined model was necessary.

The Poisson code comes from here:

http://www.netlib.org/random/zufall.f
Credit : W.P. Petersen, IPS, ETH Zuerich.
(I should check it’s public code).

More on the method to get a Poisson sampled count in regions where the expected number of photons is < 1 to follow later!

MORE stuff on this:

There is a timing hit doing Poisson noise, it’s quite
long winded.

For 30 summed images the total run time is about 40 seconds
on my laptop with Poisson noise. for 30 images with Gaussian
noise, the run time is about 2 seconds. The FFTs only take about
1 second of this. So there is presently a huge performance hit. I think this can be greatly improved.

Important note: 30 summed images using Gaussian noise look
very similar to 30 summed images using Poisson. This is
because 30 images of Gaussian sampling simulate Poisson
sampling rather well. It’s very different though when talking
about a SINGLE frame:

ABOVE : GAUSSIAN noise (left) versus POISSON (right): for a simulated lunar image in a single 30 ms frame. The very low light levels are much better modeled in the Poisson frame on the right, where the flux level drops to a few photons per pixel — and very much better when the level drops well below a photon per pixel.

Gaussian (right) versus Poisson noise (left) — but for 30 coadded 30 ms images. The differences cannot now be seen by eye, since mathematically the samplings are now more or less the same thing.

I have now written my own subroutine to do Poisson noise, and it seems to reproduce very well the results of the code acquired from http://www.netlib.org/random/. It’s slightly slower though; thought I could improve on the run speed easily, but I was wrong!

Need to still check for small systematic differences between expected mean number of photons and the numbers actually returned by the Poisson routine, since we need to worry about this at the 0.1 percent level. More on that later.

We remove scattered light from synthetic images and test the BBSO method in 4 ways.The synthetic images were generated with alfa=1.6. 1 synthetic image was used.

The BBSO method was first applied in the standard way – i.e. linear least squares fit to all pixels in 5 degree cones on the sky off the DS disc. Then we fitted to log-intensity, converted the correction into linear values and subtracted from the linear image. The two methods were applied to a patch on the DS disc 2/3rds and 4/5ths from disc center – i.e. close to the rim and not so close to the rim.The figure shows the percentage difference between the cleaned up image and the known value, as a function of the day of month. New Moon is near day 13. The absolute difference is shown but all corrected intensities were larger than the true intensities.Red symbols show the results for the standard method (i.e. linear regression), blue is for the log-method. Crosses show the results at 4/5ths from disc center and boxes show for 2/3rds from disc centre.We see that results with errors below 10% are available for only about 5 days. We see that the log-method is somewhat better than the linear method, and that best results are found close to the rim.We also see an asymmetry – as thelunar phase passed Half Moon the enalysis switched to the other side of the lunar disc. Apparently a 50% improvement in the derived intensity can be had by using the log-intensity method, but mainly on one side of the Moon. The reason is unclear. As images were generated ideally – i.e. without libration and without distance variations a fixed coordinate for the patches in which the corrected intensity was evaluated always accesses the same pixels in the synethtic image – except when we switch to the other side of the disc after New Moon.

My last post was on techniques to measure the ES using artificial images as input.
Three boxes are defined on the image which are mainly sensitive to alpha, ES and threshold. Box 1 (off the BS limb) turned out to be completely dependent on alpha alone — and one can get a very accurate alpha measurement from just measuring the photons in the box (see plot below). The scatter in alpha is only 0.01 in tests of the method based on semi-realistic synthetic images made by coadding 10 x 30ms images together and properly including Poisson noise. Error in the bias level or registration of the moon and overall normalisation of the image is NOT yet included. Having determined alpha, it takes a little care to then compute the amount of light due to BS scatter in the box on the ES limb, and subtract it. Doing this, one then gets the actual ES level in the box. This correlates very well with input albedo and can be empirically turned into an output albedo (see plot below), with a scatter of 1% — our target accuracy. So it looks to me as if 1% is achievable! Yay! However, there are quite a few systematic effects not yet included here, all of which will shift us hither and thither. So that’s yet to be done!

Basically I think this method is the same as Peter and Jacob are pursuing at the moment — simply trying to fit away the unwanted scattered light in the ES box — except that I do it “empirically” by building a calibration from a range of synthetic data with different alphas and albedos — whereas they are working from a single frame — which is what we will actually have to do in practice.

Anyway, all of the above gives considerable oomph to the idea that the glass is 90 percent full!

I’ve developed an kind of ’empirical’ method for measuring alpha and ES on a frame, I am working from 30 co-added 30 ms exposures of a fake moon, so the S/N is rather good.

Here’s how it works:

Three boxes are defined on the lunar image as shown below

Image shows an example synthetic lunar image (30 co-added 30ms exposures)
for some value of alpha and earth albedo (threshold is identically zero in
the tests done to date).

Box 1 : (far left) designed to be mainly sensitive to alpha. (11×21 pixels)

Box 2: (middle) mainly sensitive to ES level on DS. (26×21 pixels)

Box 3: mainly sensitive to threshold level — i.e. the bias (51×21 pixels) (not analysed here).

A set of synthetic lunar images is created for a range of values of albedo and alpha, chosen randomly in the range albedo = [0.0, 0.6] and alpha = [1.6,2.0].

For each alpha and albedo pair, the total photons are measured in boxes 1 and 2.
I then ’empirically’ calibrate relations between the photon counts in the boxes and alpha and albedo.

Above is the plot for the counts/pixel in Box 1 — it correlates extremely well with alpha and has no measurable dependence on albedo. The lower plot implies alpha can be measured with an accuracy of ~0.01 from the counts in in box 1. (Note well that alpha here is the power to which the PSF is raised — it is already a power law like r^-1.6, and it is then raised additionally to power alpha).

The counts/pixel in box 2, selected to be mainly albedo sensitive, are indeed very well correlated with albedo. There is a very small dependence on alpha (lower panel). Taking into account both dependencies results in recovering the albedo with an accuracy of about 0.03 — compared to a typical albedo of 0.3, for an accuracy of 10%. We require 1% accuracy! Not sure what’s driving the 10% figure at this stage. Is it merely Poisson?

We can improve things by using more of the frame — the boxes were selected pretty arbitrarily just to get the ball rolling. The modeling is also pretty ideal at this stage as the error in the threshold is not yet included, and one has to have a good synthetic moon image to get started — any error in the synthetic moon (i.e. before scattered light and Poisson sampling is included) will bias the method.

Update: method now much improved — fitting functions work much better using log of the counts in the boxes, rather than the counts in the boxes (as I should have realised — we are dealing with a power law in the scattering!). This leads to great improvements, as discussed in the blog entry above.

Using synthetic images where the DS intensity is known, we have generated ‘fake’ observed images by convolving with a PSF and adding noise (Chris’s syntheticmoon.f used; alpha=1.8). We then forward modelled these images using the EFM method (i.e. using the ‘observed’ image BS as the source and allowing that to scatter in model images until the model images matched the input image).

Subtracting the model images left us with the DS corrected (the BS goes away, per construction). We extracted DS intensities only (adding the BS is a separate chapter of complexities!). We extracted DS intensities from the corrected image, from the synthetic observed image, and from the known input ideal image. We plot these images as a fucntion of the day number during one month.

The black symbols show the DS intensity in a 15×15 patch at 2/3 of the radius for the synthetic observed images (i.e. the ones that suffer from realistic scattered light from the BS); the blue symbols show the known intensity from the ideal pre-scattering images; the red symbols are from the EFM-corrected images. The period of New Moon is excluded because our method fails there (there is not enough crescent to make a good model of the halo from). Near Full Moon is inaccessible since the DS is so small then and it becomes difficult to extract the DS intensity (nowhere to put the patch!).

Annotations in the upper panel show the choices of various settings of the analysis code.

Bottom panel shows the difference between true and corrected DS intensity expressed as percent of the true value. In the upper panel points inside the two sets of dotted lines have been selected for analysis in the lower panel – the mean absolute error is printed.

We should next apply the BBSO method to these images.

We see that this method also reaches the 1% error level but that the number of days during which this level is attainable is smaller than with the FME method.

For the night of 2455864 we have tested the effect of using a 5×5 patch on the sky vs an 11×11 patch. This is the patch in which the average is taken and through which average the model curve is forced to pass. We had noted that a 5×5 patch was small compared to some of the wiggles in the sky background so we tested the use of 11×11 patches.

We also test the effect of using 64-bit rather than 32-bit code – we do this with ‘no change of compilers or compiler flags’.

Finally we tested the effect of taking smaller steps in alfa. The figure summarizes these three experiments:
The difference in detected alfa was small when we went from 5×5 to 11×11 patches. The average change was near 0, the standard deviation apeared to be about 1 alfa step (i.e. 0.001), the max and min of the change were -0.6 and 0.4 but only 2% of the 384 images suffered a change in alfa larger than 0.01 in absolute value (i.e. outside the plotting range). We conclude that there is a tiny indication of a low-bias in determined alfa when using 5×5 patches, so we will use 11×11 from, now on.

For the 32-vs-64 bit analysis we found that for step size in alfa of 0.001 the change in going from 32 to 64 bit was most often 0. A small fraction of images had a +/- change of 0.001 in alfa – i.e. the step size.

The effect of smaller steps was negligible.

We next looked at changes in the determined DS intensity.

For each filter and for each of the three modes of the EFM code run on the CRAY we extracted the median DS intensity in a box 11×11 pixels large positioned 2/3 of the radius on the dark side. We also extracted the BS from an 11×11 box positioned at the center of gravity of the image in the observed image.

Run1 was the test for the effect of using a 5×5 vs 11×11 patch on the sky and the BS to force the model through.

Run2 was the test for effect of changing from 32 to 64 bit code.

Run3 was the test of the effect of using smaller step sizes in alfa.

We found that:

For all filters there is a 6 to 7% effect on the BS/DS ratio of going from a 5×5 to 11×11 patch on the sky when fitting the model to the observed image.

The effect of changing from 32 to 64 bit code was negligible, as above – the largest change in BS/DS ratio was 0.02% in the VE1 and VE2 filters; the others were at the 0.007-0.009% level.

The effect of smaller step sizes when searching for the optimal alfa was a little larger – in the range from 0.19 – 0.9 %; largest in the VE2 filter.

We conclude that careful choice of the way in which the model is fitted to the sky and the BS is warranted.

We have used the ’empirical forward model’ method to reduce 384 images from night 2455864.

Using the DMI CRAY we stepped through 1000 values of alfa from 1 to 2 for each image and required the best possible fit on the sky part of each image.

We plot the best value of alfa found for each filter againts the time of observation:

We see that the values for alfa are narrowly grouped for most filters, except for the VE1 filter which jumps between two values. Unexpectedly we see that the blue filters have high alfas (i.e. narrow halos) while the reddest filter (VE2) has the smallest value of alfa (i.e. the broadest halo). This is opposite of what we might expect from physics, and the code plotting these data must be reviewed.

The scatter of alfa for e.g. the IRCUT filter is 0.0036 which is 0.2 %. Since we had 98 IRCUT values the standard deviation of the mean alfa is about 0.00036 which is about 3 times larger than the required accuarcy estimated on this page below, to achieve a DS accuracy of 0.1%.

Visual inspection of the fits show that improvements are possible, so more experiments with parameter settings is called for. The above were results from fits of the absolute residuals. Next we try fitting the relative residuals, i.e.

(obs – model)/obs

instead of

(obs – model).

We have now tested that – and it was not a good idea: fits got worse.

Henriette and I ran the BBSO linear extrapolation method for removing scattered light with two settings of a parameter that is arbitrarily set. We got output for many files with these two settings and can therefore compare the two results for a single image and measure the difference.

On the image below we see the relative differences between two pairs of two such processed images.

On the left-hand side of the upper image the white cones correspond to percentage differences in the 10-30% range. They gray or black areas are much lower.

The second image is doing better with most differences in the single-% range.

The images processed here are single exposures in CoAdd mode – so there is lots of noise on the DS and the sky. Performance of the BBSO method will probably improve with higher SNR.

We changed the processing parameter (cone width) from 8 degrees to 6 degrees. BBSO uses 5 degrees (in the JGR Qiu et al paper). We can now look at all images where we have these two processing results and make some summary statistics.

Below are two histograms of just that for 383 pairs of BBSO-method corrected images- the first shows the histogram of mean difference values (they are in %). The other is the histogram of the median difference value. Both measures of performance are centred over 0% – the mean difference is a broader distribution than the median – by about a factor of 10 (FWHM mean is 2% while FWHM median is 0.2%). The differences were calculated for all pixels (DS sky and DS itself, but excluding BS and BS sky) corrected by the BBSO method evaluated with the two settings as described above.

The result suggests that correcting the DS for scattered light using the BBSO method and extracting a subregion mean will give unbiased results with a few % scatter, while extracting the median of some corrected subregion will give the answer to a few tenths of a percent.

With various systems to simulate realistic observations at hand – Chris’ noise simulator and Hans’ ideal image simulator – it becomes possible to pose some questions, such as:

1) Is it better to reduce every
blessed image and then average derived results, or should images be
coadded and then reduced?

2) Knowing that the bias has a 1-count
amplitude, 20-minute period is it then a sensible
strategy to scale a low-noise ‘superbias’ to the observed but noisy bias?

3) What is the relationship between precision in alfa and scatter in
DS/BS ratio, in the presence of realistic noise?

A partial, idealized, answer to 3) is hinted at by considering the change in DS and BS intensities as a function of change in alfa:

DS change: -15.9935 %
BS change: 0.0893807 %
alfa ch: 0.995023 % [Note that PSF was normalized]

I.e. If we change alfa by 1% we get a 16% change in a typical DS point and a 0.1% change in a typical BS point. This suggest that we need to know alfa to an accuracy of 0.6% (alfa was 1.7) or 0.006. This would be the typical step-size of a grid-search, for instance, and the tolerance on any downward-descent search. Perhaps it is easy to get such accuracy with methods that ‘fit the sky’ – such as both the BBSO and our own forward methods.

How does the above change in the presence of realistic noise?

Now have a non-interactive version of the BBSO method for removing scattered light, written in IDL.

It takes about 10 seconds to ‘clean up’ the DS of one image.

The code is in SCIENCEPROJECTS/EARTHSHINE/ and runs by using a script:

./go_linearclean.scr filename.fits

a file is produced called

LINCLEANED_filename.fits

The method is non-intercative because the user no longer has to specify disc rim and disc center with a cursor. The estimation is done by edge-detection methods and an algorithm to fit a circle to three points in the plane, not on a line.

I have made an attempt to figure out what the noise levels are at various fluxes in individual frames.

28 images of the nearly-full moon were used. They were very carefully aligned (using subpixel rebinning and the iraf task imshift). The shifts were determined with a fortran program which centered the moon using center of light — this was written because no iraf task I could find could easily align such images without reference stars.

Below I show in green circles the regions that were used to measure the standard deviation of pixels in a 10×10 box on the lunar face. These were chosen because they are relatively flat and have a range of luminosities. All the pixels in the scattered light off the lunar edge were also used, in a sector starting from the center between the position angles 40 and 50 degrees (i.e. off to the upper right at 45 degrees).

The standard deviation of the selected pixel positions over all the 28 exposures was computed, and a plot made of the sqrt of the mean counts in the pixel versus the standard deviation of the counts.

If the noise is Poissonian, then these are related 1:1. The plot below shows that this is only approximately the case (solid line). A closer model is that the noise is about half of what is expected (dotted line). Why the noise is sub-Poissonian is open to speculation, but must be due to some slight smoothing we don’t account for (possibly iraf introduces some smoothing when rebinning the data for the subpixel shifting).

The points at the lower left are those from the scattered light off the lunar edge. The other points are the pixels in the 5 small patches indicated by the green circles in the previous figure.

It’s difficult to fill in more points because it’s hard to find flat regions of the appropriate luminosity. Flat regions are needed because the technique is very sensitive to slight misalignments in the images (if one examines regions where the luminosity is changing rapidly — such as near the lunar edge, craters, strirations etc on the lunar face — the standard deviation goes haywire).

Assuming Poissoninan noise in test data generated from synthetic moons appears to be a conservative option, as it is pretty close to the correct value and possibly an overestimate, so currently I’d recommend that.

Update and error correction: The ADU conversion for this chip is approx. 3.8. The statistics above are were done in ADU, not photons. This of course accounts for the factor of two (=sqrt(3.8)) problem above. Thanks to Peter for spotting this error.

Using photons instead of electrons gives a pretty close fit t the 1:1 relationship. We could try this for a range of lunar exposure times, using the same flat regions on the moon and the halo off the edge, to fill in the gaps for a range of luminosities.

Peter sent me images of the ideal crescent moon for a range of Earth albedos — 0.0, 0.29, 0.30 and 0.31. I made a map of the scattered light to see how much brighter the ES (Earthshine) is than the scattered light from the crescent. Over most of the dark side, ES does indeed dominate over the scattered component — which augers well for measuring its ratio relative to the BS (bright side).

Top left: ideal moon outside the atmosphere with NO earthshine at all (Earth albedo = 0) i.e. we see the BS (brightside) only.
Top right: scattered light from the BS only.
Bottom left: scattered light of the same ideal moon, but with the ES for albedo=0.30 included.
Bottom right: the difference between the top right and bottom left images.

Interestingly, one can subtract the modeled BS to obtain the ES over the entire lunar disc. The isophotes of the remaining light are then nicely round around the full, ES-only illuminated disc, as seen in the lower right panel. This might be a way to check that the solutions are coming out the way they should — i.e. fit the BS and scattered light from it only, remove it — and the rest should be the ES. Errors in modeling the BS will show up as asymmetries in the ES that remains. Perhaps we could think about solving for the ES directly, by fitting a lunar model with a BS only, and recovering the ES on the lunar face directly?

After reduction of the observations the color of the earthshine can be considered. It will depend on the type of surface that predominantly reflected the sunshine but should be even across the DS at any one time – apart from lunar surface details, of course.

Any systematic gradient in the color of ES across the DS is a sign that the scattered light removal was imperfect – so a test can be set up.

I have spent some time examining to what extent GPUs could help us speed up scattered light fitting of ideal moon images to the observed images, while visiting Swinburne in Melbourne.

Ben Barsdell (Swinburne) had some c++ code to convolve two fits files using an NVIDIA compatible GPU (GTX480) running CUDA. The FFT library is FFTW, the same as I have been using for CPU-based modeling of the scattered light.

Upper left: ideal lunar image (intensity is log scale).
Upper right: convolved with our best PSF using CPU based FFT.
Lower left : convolved with PSF using GPU.
Lower right: ratio of the two methods (seen in more detail below).

On a desktop my fortran code calling FFTW does the three FFTs needed —
ideal moon image, psf image, their multiplication in the Fourier
domain, and the inverse FFT for the final result — in about 1100
milliseconds (ms). (This excludes the time needed to insert the 512×512
images into 1536×1536 images to sufficiently reduce edge wrapping
effects, which brings the total runtime to about 3000 ms). So the FFTs
are taking about 300 ms each.

Using the GPU, we attained FFT
speeds of about 40 ms each, a speed of a factor of 10 or so. This is
typical of what Ben expects for such applications (he is doing a careful
study of the types of astronomical problems GPUs can be profitably
applied to, and where the bottlenecks typically lie.)

Ostensibly
this means we can speed up our light modeling code by about 10 times — possibly
quite a bit more because the overheads per modeled image can be reduced
quite a lot by careful programming, since we want to explore a large
parameter space, but don’t need to do the same overheads each time we
run.

VERY IMPORTANT: the CPU code did the FFTs in double complex precision, whereas the GPU was doing single complex precision.

We
compared the output images of both methods — for the case of scattered
light from an ideal moon with a power law fall off PSF with a slope of
about r^-2.8.

VERY IMPORTANT: there is significant structure left in the methods if we divide the CPU output by the GPU output, as shown int he image above.

The scatter in the mean about 1.0000 is 2.028E-4 and there is a frame minimum of 0.9936 and frame maximum of 1.003, so the deviation from unity in the ratio of the two methods is not worse than ~0.7% anywhere on the frame. There is certainly structure, and it looks like it might be too much for our purposes! (We would like this to be better than 0.1% at the very worst). We are checking if this has to do with the single precision used.

Unless this problem can be solved, it is not clear that the speed gain with the GPUs
is worth having!

Tested two methods for scattering light from an ideal moon. The Ideal image is shown upper left (log scale) with a completely dark sky, narrow BS crescent and ES.

Bottom left is the moon after scattering the light using our current PSF and a brute force method — simply add all the pixels in the ideal image multiplied by the PSF. Bottom right is the same, but using an FFT method instead. The scattered light at the edges of the frame is too high to simply use the 512×512 frame in the FFT, instead the idea moon is placed in the center of a 1536×1536 frame, the FFT carried out, and then clipped back to 512×512. Upper right is the ratio of these two methods. There is a lot of structure at a very low level — the display limits on the ratio are set at 0.9993 to 1.0009, and the ratio is with 0.1% of unity everywhere on the frame. This should be sufficiently accurate for our purposes, but we should double check!

The histogram of pixel values in the ratio image is shown below:

For reducing data so that the terrestrial albedo can be determined we have forward and inverse methods.

Inverse methods:

1) The BBSO linear sky-fit method. Yields a cleaned-up DS so that DS/BS ratio can be found and used for the inversion. Probably works best near the edge of the DS as the scattered light is not a linear function of distance except far from the BS.

2) Andrew’s clean-up method based on fitting a model image convolved with a variable PSF until it can be subtracted perfectly from the sky-part of an observed image. As a model of the Moon an observation – minus the sky part and the DS – is used.

Forward method:

3) A synthetic model of the Moon is used within a convolution framework to produce an image that ‘looks just like the Moon and its halo’. Since a parameter in the method is the terrestrial SSA it is directly determined in this way, but does depend on the same assumptions about time, geometry, lunar reflectance and so on as does the inverting methods. Depending on how much of the lunar disc is included in the fit there may be a difference in how the knowledge of lunar albedo and lunar reflectance influences bias.