Blog Image

Earthshine blog

"Earthshine blog"

A blog about a telescopic system at the Mauna Loa Observatory on Hawaii to determine terrestrial albedo by earthshine observations. Feasible thanks to sheer determination.

Effect of flatfields

Data reduction issues Posted on Jun 05, 2015 12:10

We have generated ‘super flatfields’ by considering all the flat fields we have – we have many exposures through each filter taken on sky, the dome inside and the hohlraum lamp. While the patterns are the same for the three types (mainly a roof-tile like diagonal pattern’ we have dips in the corners of especially IRCUT and VE2 (and VE1). It is therefore not obvious which flat field is the ‘correct’ one. Assuming that the ones with the least vertical relief must be the ones closest to reality we have adopted one superflat for each filter.

We need to investigate if there are other ways of testing which FF is best – and we have the possibility of ‘scamp’ and ‘sextractor’ to generate background fields, which we assume are some sort of smoothed residual after sources are removed. Using frames from NGC6633 we generate these background fields from ‘raw’ frames and ‘bias-subtracted and flat-divided’ cluster images.

A sample of the background fields, shown as contour plots, for one science frame is shown here:

The upper field is for ‘raw’ and the lower for ‘bias-subtracted and flat-divided’. The contours seem to have a mean near 18 for both – the contours evidently are not the same. The differences between the two backgrounds amount to 0.1 – 0.2 counts.

We should center the backgrounds on sky-coordinates and average and see if the fields are consistent given method (we assume so). TBD.

If the patterns are stable we have learned that applying bias-subtraction and flat-division can alter the levels by a few tenths of a count. This may be important for photometry of the DS. Presumably it is better to correct for the flatfield than not to, but that depends on whether the flat-processed fields are consistent. Perhaps warp can be used to generate a stacked field, and its background?

Added later:
Doing just that we get the average background frame for ‘raw’ and ‘bias-subtracted, flat-divided’ processing:

The relief in that field (took the most-overlapping 512×512 region) is from +0.1 to -0.02 — so some averaging has occurred, but there is still a one-tenth of a count effect of applying the flat (in, B, and when 14 images could be stacked).

So, we cannot rule out that applying flats have a large effect compared to DS photometry – i.e. 0.1 out of 1 is 10% and at best 0.1 out of 10 is 1%.

The effect is regional inside the frame and so we can well have an effect on the limited-region fittings we perform on the DS.

Might want to consider what the effect on the BS-halo, evaluated at the DS, is if the BS also has flatting-effects brought to it.

Basically, applying flats are necessary as long as they do not bias results due to bad production of the flat-field itself.

Two-way smear from overexposure

Data reduction issues Posted on Apr 27, 2015 10:41

Here is an example of a two-way smear. The red area is the overexposed region, and it has bloomed along vertical columns in the image. To the right of this bloom is a smear extending hundreds of pixels. What is that?

Read-out direction is top-down.

What is causing this smear to the right? Is it related to the ‘depression of the sky in a broad band left-right’ that we see in some images?

The intensity of the smear is at the 500 count level – i.e. it acts just like part of the halo from the BS. This has, surely some consequences for data reductions? Do images that have not bloomed have this halo? Is this why the bright edge of the BS-halo is always hard to model?

Can we see if the fitting of the BS halo is easier in images that have been modestly exposed, compared to those exposed right up to 50.000 or 55.000 counts?

Extinction two ways

Data reduction issues Posted on Jan 13, 2015 15:34

Atmospheric extinction is caused by absorbtion and scattering of light.

Usually photometry is performed on images of sources by collecting the flux from an area near the source and relating it to the airmass of that particular observing moment. In a Langley plot the relationship between magnitudes and airmass will be an (almost) straight line, and the extinction coefficient is the slope of that line.

As light passes through more and more atmosphere there will be more absorbtion as well as more scattering. Both remove intensity from the direct beam. Absorbtion just takes the light away (at least at optical wavelengths) while scattering redistributes the light near the source. If the method of measuring light from the source therefore allows re-capture of the light scattered – for instance by allowing for large collection areas in the image – then the effect of scattering will be reduced and a smaller extinction coefficient should resultr, as compared to the extinction calculated when scattering is allowed to remnove light from the calculation.

Thus, when considering an extended source like the Moon we have an opportunity to test the magnitude of scattering vs absorbtion by explicitly collecting all light near the Moon (for instance, by summing the flux in the whole image frame) and comparing to what we find when light is collected only from a small patch on the lunar bright side.

We have done this for 5 wavelength bands and a score of observing nights where conditions were right. We measured flux in each image in two ways – either as the total flux in the image frame or as the flux inside a small patch on the lunar brightside defined by selenographic longitude and latitude. Langley plots were constructed and slopes measured using regression, along with error estimates of the slopes.

Results are plotted here:

Langley extinction for whole-image fluxes plotted against extinctions determined from a small bright-side patch in each image.

If the image is rotated and inspected (sorry about that) we see that for the five filters extinction determinations often agree, despite the expectation that extinction determined from the small patch (allowin scattered light to escape) would be larger than the extinction determined from whole image. Only a few points deviate significantly from the diagonal line.

This implies that scattering at MLO was very low, compared to absorbtion – as scattering includes Rayleigh scattering thge plot for B should be considered carefully.

In a few cases the plots allows us to identify ‘bad nights’ when the extinction determinations do not follow the diagonal.


Data reduction issues Posted on Aug 30, 2014 11:43

We have commented earlier that for an otherwise well-fitted model image the ‘horns’ or cusps of the fitted image appear relatively too bright compared to middle-of-the-disc pixels. We investigate this now.

For an observed image and the corresponding best-fitted model image we extract certain pixels defined by the incoming and outgoing angle wrt to local normal on the lunar disc and compare the intensities in the observations to those in the model image. Here is an example:
The upper panel shows the observed image with certain pixels picked out (the even dark band along the terminator). These pixels all see the Sun under angles between 55 and 65 degrees wrt the local normal. In the lower panel we show the intensities of those pixels in the observed image (black) and the model image (red) plotted against the angle between local normal and the direction towards the Earth.

We see that the model ‘dips’ for large angles (i.e. when the direction towards Earth is almost paralell with the local surface) while the model values (red) are unwaveringly high. We also notice that the model values start low at low angles and rise with a modest plateau starting at 75 degrees or so.

If we repeat the above for other sun-angle bands we see the same pattern – model is ‘flat’, observations ‘dip’ when the angle to Earth is large, leaving the cusps always brighter in the model than in the observations.

Bulk extinction coefficients

Data reduction issues Posted on May 23, 2014 13:12

In this entry we estimated extinction coefficients for all data from each filter – that is, we did not study each night individually, but took all B, V etc data seperately.

We found, at that time:

B 0.15 mag/airmass
V 0.10 mag/airmass
VE1 0.08 mag/airmass
VE2 0.06 mag/airmass
IRCUT 0.12 mag/airmass

Since then we have eliminated some observations that we now know had problems of one sort or the other, and have the opportunity to re-estimate extinction coefficients. We now find:

B 0.18 mag/airmass
V 0.11 mag/airmass
VE1 0.06 mag/airmass
VE2 0.09 mag/airmass
IRCUT 0.05 mag/airmass

kB is 0.03 higher; kV 0.01; kVE1 0.02 lower; kVE2 is 0.03 higher and kIRCUT is 0.07 lower. The changes of +/- 0.02 are as expected given Chris’ analysis of extinction from single nights, but the change in kIRCUT is large – however, it is now more in line with kVE1: the two filters are almost identical, so that is a step in the right direction.

We found these extinction coefficients by plotting extinction-corrected flux against lunar phase and fitting a third-order polynomial. For trial values of the extinction coefficient

Large halo

Data reduction issues Posted on May 17, 2014 12:23

In some of our images sky conditions were such that the halo obviously was ‘too large’ for good results to be expected. We seek to eliminate these images now.

One simple measure of a halo that is “too large” is by looking on the BS and finding those images where the halo drops slower than in other same-phase images. We calculate and index of halo-narrowness based on the ratio of brightness 20 pixels beyond the BS edge, and the maximum brightness in the image. We call this the Q20 index.

On nights with ‘large halos’ we also often see a bright sky quite far away from the lunar disc, on the DS part of the sky. We quantify this by the mean brightness of the 10×10 box in the lower DS corner of the bias-subtracted image.

We plot now the phase vs Q20 and the sky brightness vs Q20:

Each row is for one of the 5 filters, left column is Q20 against the illuminated fraction of the lunar disc (i.e. phase). Right column shows sky brightness on the DS vs Q20. Q20 is the ratio of the BS intensity 20 pixels from BS edge on the sky and maximum brightness of the image.

We see that Q20 is quite constant against lunar phase, implying that we are correctly normalizing the halo. We see only a small drop in Q20 for (perhaps) V and VE2 filters as phase becomes very small.

We see some outliers where Q20 is high compared to most of the data at same phase.

We see (right column) that sky brightness and Q20 are related – for large values of both, but otherwise Q20 is a constant.

In fitting images to models – discussed abundantly elsewhere in this blog – we realize the need to concentrate mainly on images with as little halo as possible, so we wish to eliminate the images with a large Q20 value. Where shall we set the cutoff?

On the plots the dashed line shows 1.3 times the median Q20. We (arbitrarily) choose to get rid of all images with Q20 above this line.

We also note that Q20 is largest for VE2 and smallest for B – light is scattered further the longer the wavelength is. Does this tell us anything about the scattering mechanism? I’d venture that if the scattering was basedon the Rayleigh mechanism we should see more scattering in th eblue – wonder how Mie scattering goes? At least some of the scattering is in the optics and is fixed and unrelated to the atmsopheric effects – perhaps that can be confirmed by considering some optics?

What of the sky brightness?

We see that VE2 and V have higher sky brightness values than the other filters – VE2 so more than V. Most B, VE1 and IRCUT images have sky brightness below 1 (count), while VE2 has mos between 1 and 10. V has most between 0.,1 and 1 butr a fair fraction above 1.

For VE2 we spculate that the sky brightness is realted to long-wavelength flux from the sky – this could be from water vapour – so we suggest that VE2 sky brightness is compared to MLO weather data for humidity. Since the H2O may be high up in the atmsophere and not a ground level where the hygrometer is, this may not be easy. Satellite imagery for the nights in question could be brought to play on this issue. [Student Project!].

We are uncertain why V should have sky brightness issue if nearby B and VE1/IRCUT do not.

We arbitrarily suggest to also eliminate all images with sky brightness above 2 counts (10 for VE2).

Removing such images from ‘Chris list of best images’ we are left with 494 in number, and the list itself is available on request.

How bad is the Drag?

Data reduction issues Posted on Sep 22, 2013 11:04

During data-reduction we found the many ‘dragged images’. These are images where the shutter probably did not close before readout started. They look like this one (histogram equalized):
The ‘drag’ is from the BS but is probably present everywhere, but not visible from the DS. The BS has a maximum near the right hand limb where the counts are near 10000. The ‘drag’ is near 240 at most. There is no gradient in the drag along the direction of the drag:

The absence of a gradient in the drag helps us understand its origins – it is the effect of constant-speed readout. The readout speed is supposedly 1e-6 s per pixel (in the image header: this may be read from the camera but could be a fixed number entered at setup ….). Reading 512×512 pixels should take 0.26 seconds. The exposure time for this image was set to 0.0274s. Being a ‘dragged image’ we have no idea if this is right.

I do not quite understand the ratio of maximum drag value to maximum BS value. If the drag/BS ratio is 240/10000 and the BS exposure time is 0.0274 s the drag must have lasted just 240/10000*0.0274 = 0.00066 s. The camera is a frame transfer system – is frame transfer that fast? Perhaps the 0.26 s readout time is a slower line-by-line operation? This page mentions ‘a few msec’ for the frame transfer: . We have 0.66 msec, so – close.

If the shutter got into the act in mid readout the amplitude of the drag is affected by this. I think many dragged images we have seen are like the above one – dragging is seldom much worse than this.

Let us check other dragged images with higher exposures and see if the drag/BS ratio is constant or not.

The above was prompted by pondering whether we could have a ‘shutter-less’ camera. It would appear not – at least not if the drag is due to open shutter during readout. The readout would have to be 24 times faster than it is and I am not sure CCD cameras come with such fast readout. 24 MHz readout?

Colour woes

Data reduction issues Posted on Aug 19, 2013 12:57

In this post we commented on the problems we have with the requested and actual exposure times. It turned out NOT to be a good idea to correct the requested exposure times – certain colour problems got worse!

What happened was that B-V for the BS plotted against time of day showed a decided reddening towards the end of the Julidan day – this corresponded to the setting Moon. Airmasses were growing towards the end of the day.

When time of day and B-V were plotted we could see slopes using UNcorrected exposure times of 0.4 magnitudes/day. If we DID correct the exposure times as explained in the above link we got almost a whole magnitude per day in slope!

This got me thinking about whether the extinction corrections we perform are correct – they depend on knowing what the airmass is and knowing what the extinction coefficient is. If either, or both, are wrong we can expect colour-terms to pop up.

The airmass we have been using assumed surface pressure and 0 degrees C and 0 % relative humidity and a wavelength of 0.56 microns. More appropriate is the pressure that occurs at the MLO (we guessed it is 0.6*760 mm Hg), and a temperature of at least a few degrees C (we never observed in frost), and a relative humidity of 20%. Some of these quantities could actually be replaced by the available MLO meteorological data available – more on that later. We should also use an appropriate wavelength for B and V of 0.445 and .551 microns, for B and V respectively. However, using these better values will at most give an offset in B-V and cannot explain an airmass-dependent colour term. [or can it? gotta think about that one …]

Next we must consider the formula used for the airmass. Looking in the IDL code ‘’ provided by Marc W. Buie, we see it is the “cosine based formula derived by David Tholen” with an update 2003. The formula appears to be:
cz = cos(zenith[z])
am[z] = sqrt(235225.0*cz*cz + 970.0 + 1.0) – 485*cz
When we have a reference for this I shall return to the question of whether the approximate formula is good enough.

Finally there are the extinction coefficients. Since we see an increase in B-V with airmass it is either the V magnitude that is being over-corrected or the B magnitude that is being under-corrected. Experimenting, I find that kV=0.08 and kB=0.17, instead of the 0.1 and 0.15 that Chris found from standard-star observations, givea a colour-free B-V progression with time of day (at Moon-set – Moon-rise is too underrepresented to be useful here). Is it within the error margins that kV and kB should have these values? Old man Chris says he estimated the values 0.1 and 0.15 ‘by eye’ here but I think he also listed some regression results elsewhere, with errors – but I can’t find ’em!

Anyway, using the values 0.08 and 0.17 gives a mean BS B-V of 0.898 +/- 0.054 (error of mean +/- 0.004). This value is relevant when keeping the discussion about the “van den Berg value for B-V=0.92” in mind, link here. Our own average of vdB62 data gives 0.88 +/- 0.02 – within the errors we have, easily.

Core vs. Wing

Data reduction issues Posted on Jul 24, 2013 23:56

In using the fitting of model to observations at the edge of the lunar DS we have the following situation:

It is seen that the ‘edge’ of the observed image is not ‘sharp’ – it slopes over perhaps 3-4 columns on either side of the midpoint of the jump from sky to disc. The models often have problems matching that low slope and tend to be ‘sharper’ at the edge, at ‘best fit’. We fit +/- 50 columns on either side of the edge so more pixels on the sky participate in driving the fit than does the ‘edge-pixels’.

The shape of the edge is driven by convolution with the PSF – and the PSF core dominates in setting the sharp details. The slope of the DS sky halo is given by the wings of the PSF.

There is an opposition between fitting edge and halo slope – lower alfa would broaden the PSF and ‘lift’ the halo wings, increasing the slope at the same time as ‘dulling’ the edge of the model edge. With more pixels driving the fit at the sky halo the edge is bound to be ignored – thus too large alfa’s probably are the result. The pedestal added to the model image to match sky brightness is competing with alfa’s effect on the sky halo – part of the level of the sky halo is provided by the pedestal value while part is given by the effect of alfa on the PSF itself.


1) Exponentiation of the canonical PSF to power alfa may not preserve the shape of the PSF near the core in the required way.

2) More weight could be given to the ‘edge pixels’.

3) Slope of the edge is not just due to rounding by the PSF – there is also some ‘jitter’ from co-addition of 100 frames so alfa is used to account for more than the effects of the PSF.

4) The fitted pedestal is sometimes negative – it should probably always be positive due to sky brightness after bias subtraction, but the bias varies by 1 count due to the temperature effect and is scaled, as described elsewhere. The pedestal is often near -0.1 counts – are we sure that the bias has been correctly scaled to the level of 0.1 counts?

Is RON stationary?

Data reduction issues Posted on Apr 17, 2013 16:52

We are using co-add mode – that is, we take stacks of 100 images and co-add them to increase the SNR. Since there is RON (read Out Noise) in each frame – about 2 ADU/pixel, tests have shown – it is important to know whether a sequence of added frames have a mean that converges.

We now test this directly on the stack of 100 images: “2456015.7742682MOON_V_AIR”. We extract a small square of size (2*w+1) around regularly spaced points covering the image plane and calculate a mean subframe (averaging along the stack direction), the average value of that mean subframe (a scalar for each of n coadded subframes), and plot these. These are shown in the PDF file:

The top panel in each page shows where the point was chosen. w was 4 so subframes are 9 pixels by 9 pixels, centred on the point. Second panel in each page shows how the mean value of the average subframes evolves, along with estimates of the +/-1 sigma error bounds (calculated based on the value the series converges to, or at least evolves to, in the last 10 steps of the series; the known RON (about 2.4 ADU/pixel in all observations, estimated from bias frames differences), the number of pixels in each subframe and the number of subframes co-added). Last panel shows the mean value of each sub-frame used (along the stack direction).

What do we see? In page 1 a point on the sky has been selected and we see a mean value that evolves inside the expected error bounds, and we see that the subframes have mean values with some spread, but no trend.

Page 11 shows a point on the DS. Nice evolution.

Page 24 is on the BS and is very hard to explain. This holds also for points ‘near the BS’.

It seems we have little to worry about on the sky and on the DS, but on the BS we see strange evolution of the running mean! We have before touched on such subjects in this blog, when we tested effects of alignment of images – perhaps this problem was hiding inside the other problem?

[added later:] Chris asked a good question. Here is a plot of the running mean of the whole frame, as frames are taken from the stack. The value is expressed as a percentage deviation from the middle value on the curve – about 1378 counts. A small drift is seen at the 0.05% level from one frame at left to all frames at right. I think this could be due to sky variations – or slight drift of the Moon inside the frame causing ‘light to fall off the edge’. If drift is the explanation we may have an answer for why the mean in a much smaller sub-frame, when near the BS, drifts so much more – image brightness gradients are being sampled in a small subframe. Experiment at top should be redone but with drifts taken out. No large drift was evident sp I guess we are learning that a small drift can be very important!

This has an implication for how we measure the albedo – if we use ‘DS patch divided by BS patch’ we run a risk . if we use ‘DS patch divided by total flux’ we are much better off.

Further adventures of VE2

Data reduction issues Posted on Mar 12, 2013 16:11

We have been following the mystery of the VE2 filter. Ana Ulla asked a good question: “Do you see the same problems in images not of the Moon?”, so I checked that, using our copious collection of Altair images. I fit a Moffat profile to each well-exposed image of Altair and extracted the sky level. I plot it here as a function of which filter (0=B, 1=V, 2=VE1, 3=VE2 and 4=IRCUT) and as a function of time, for one night (JD 2455845):

In the top panel we see that VE2 (the brownish points) have two levels, and in the lower plot we see that the VE2 level switches from high to low after some time. (Bottom panel has fractional day as x-axis). We also see that the black points (B) increase their level (slightly) at the same point (0.87) where VE2 drops to the low level. Blue points (VE1) have some history like that too. V (red) is quite stable. Green is IRCUT.

If the excess in VE2 (and others) was due to a source of light, then we would hardly see a down-shift in the level. Ana’s idea was that thin cirrus could enhance the nIR band of VE2 (quite like the Johnson I band) – but that does not explain why B rises, unless that is a particular behaviour of cirrus clouds, which I doubt. I do not think anything above is consistent with a light source. I think this is what we have had to conclude from the other results presented in this blog – e.g. the inverse exposure time dependence and the inverse airmass dependence.

I am guessing when I say that I think it is something to do with electronics. The camera knows which filter is being used! (Not).

The VE2 Mystery

Data reduction issues Posted on Mar 08, 2013 21:25

During reduction of our data it has become evident that the VE2 filter gives us certain problems. A new problem is that VE2 exposures tend to have a ‘sky component’ that is much larger than in the B,V,VE1 or IRCUT filters. Here is a plot of the relationship between exposure time and the ‘pedestal’ – i.e. the offset that has to be subtracted from the VE2 image to bring the sky level down near to 0 (in bias-subtracted images, of course).

It does seem like the pedestal height is a decreasing function of the exposure time. Why is it any function of exposure time, when the other filters do not show this behaviour?

If there was a source of light in he image frame, then its level should INcrease with exposure time, surely? Here is the plot for the VE1 pedestal:

The level is much smaller than in the VE2 case – and it does not appreciably depend on exposure time.

What is going on?

VOM plots and ideal models

Data reduction issues Posted on Feb 12, 2013 11:41

In several posts we have considered the behaviour of variance-over-mean images (VOM). In the ideal world these should be 1 because the noise (once RON is subtracted) is Poisson-distributed. We have seen how this is hardly the case in observed images.

We now consider a test using ideal models with and without imposed image drifts. We generate 100 synthetic images with Poisson noise, and RON as well as a bias.

Below, each set of three panels show a cut across one image in the synthetic stack, then a slice of the VOM image before image alignment and then the same slice after alignment. On the left are images that were not drifted while to the right images were allowed to ‘jiggle around’ by a few pixels.

We see that VOM is 1 on the BS in un-jiggled images, but that DS and sky values fall below 1.

We see a HUGE effect on images that were allowed to drift.
Some of it we understand. Before alignment, VOM rises on structured areas of the drifted images because surface albedo variations are being introduced for a given pixel along the stack direction. The effect on the DS and sky is much less – perhaps because the Poisson noise is so large comapred to the signal variations. After alignment, VOM falls to slightly below 1 on the BS, except near the edge. On the DS and the sky, though, a large lowering is seen. So far it is not understood how this comes about.

Any strange effects seen in observed images will be all the larger since the images do not just drift but also ‘shimmy’ because of atmospheric turbulence.

The effects of aligning images to sub-pixel shifts is part of the above.

Let us learn from this that a noise model probably is hopeless to build in the presence of image shifts – despite realignment – and that sub-pixel interpolation is not a welcome added bother. We could just omit single images with the most drift and use a straight average of the remaining stack. In real images we have the option to not use stacks that have a l ot of drift – but we do not know the extent of the ‘shimmy’ for the remainder.

These realizations above have the most impact on our ability to interpret results that discuss the effect of alignment – alignment reduces some problems but probably adds some others.

Noise model

Data reduction issues Posted on Feb 08, 2013 12:37

It would be a Good Thing if we could explain observed images in terms of a ‘noise model’ – that is, a model that explains why the signal and the noise bears a certain relationship.

We consider here a stack of 100 raw images taken in rapid succession. The difference between adjacent images should only contain noise in that any objects shown in both images should subtract. The difference between two series drawn from the same poisson distribution but independent has a variance that is twice the variance in either of the series subtracted, so we apply a factor of 1/2 to the difference image to get the variance image.

We subdivide such a variance image into smaller and smaller squares, and in each square we calculate the variance. We generate an image of the ratio between the variance image an dthe mean original image, also subdivided In the case where everything is Poissonian this image should be a surface with value 1. In reality there will be noise in this surface – and there will be strcutures seen wherever imperfect obejct subtraction took place. Here is the result:

In the three rows we results from different subdivisions – into 8×8, 4×4 and 2×2 squares. In the leftmost column is a histogram of the values of the image and to the right is plot of the profile across the middle of the image.

We see, to the right, that variance over mean (VOM) is not 1 everywhere. We are evidently picking up a good deal of variance near lunar disc edges [The image corresponding to the above is a quarter Moon with the BS to the right and DS to the left, situated in the image field center]. We see that the sky manages to have VOM near 1 and that parts of the DS does this, but that most of the BS appears to have VOM>1. Even ignoring the peaks that are due to edges we see a value for VOM near 2 or more [also seen in the rest of the 100-image sequence].

Since we used raw images we have to subtract a bias level. We subtracted 390 from all raw image pixels. We applied the ADU factor of 3.8 e/ADU to all bias-subtracted image values. We then calculated the variance and subtracted the estimated RON (estimated elsewhere at near sigma_RON=2.18 counts, or variance_RON=4.75 counts²; consistent with ANDOR camera manual “8.3 electrons”).

The subtracted bias is a little small – the observed mean bias value is nearer 397, but if we use that value we get strange effects in the images – only a relatively low value of the bias mean gives a ‘flat profile’ for the sky in a slice across the image. This is one poorly understood observation.

We also do not understand why the VOM is nicely near 1 on the DS while it fails to clearly be near 1 on the BS – both areas of the Moon have spatial structure which is bound to contribute to the variance in the difference image during slight misalignments.
That is the second poorly understood thing.

Progress towards a ‘noise model’ is therefore underway, but there is some distance to go still.

How would non-linearity in the CCD influence the expectations we have from Poisson statistics?

Effect of image alignment on result quality

Data reduction issues Posted on Feb 07, 2013 10:42

We are using stacks of 100 images. These are acquired by the CCD camera in rapid succession – 100 images can be obtained in less than a minute. During that time there may be small motions of the telescope, and the turbulence in the air, as well as the slight change in refraction due to airmass changes may cause image drifts in the stack.

In one approach we ignored the possible drift and just averaged the 100 images. In the other approach we use image alignment techniques to iteratively improve the alignment of the stack images: First we calculate the regular average, then we align all stack images against that average image, then we calculate a new average image on the basis of the alignments and re-align all images against this average and so on. This procedure turns out to converge, and after 3 iterations we can stop and save the last average image.

We would like to know the effect of doing this on result quality. We therefore generated two sets of averaged images – the first using the simple first kind above, and the second the iterative method.

We estimate ‘errors’ in the bootstrapping way. That is, we extract MHM averages (mean-half-median, as explained elsewhere) of the DS intensity on raw and cleaned-up image patches and also estimate the statistical error on these values by bootstrapping the pixels inside the patch, with replacement. This bootstrap procedure gives us a histogram og MHM values and the width of this distribution is a measure of the ‘error on the mean’. We express this error as a percentage of the mean itself, for RAW images and images cleaned with the BBSO-lin, BBSO-log and EFM methods.

we now compare the results to see the ‘effect’ of performing alignment of stack images:

PSF Alignment RAW EFM BBSOlin BBSOlog
1: one without 0.73 1.25 0.74 0.79
2: one with 0.57 1.03 0.62 0.68
3: two with 0.58 0.99 0.63 0.69

Table showing errors (in percent of the mean). Lines labelled 1: and 2: show results for ‘one-alfa PSFs’ with and without alignment. The third line shows the effect of using a ‘two-alfa PSF’ and alignment.

We see that there has been a large reduction in errors by using alignment. Raw images improved by 20%, EFM images by 17%, BBSO-lin by 16% and BBSO-log by 14%.

The effect of using a two-alfa PSF on aligned images is small – indeed, all images except EFM experience a small increase in the error (probably not significant).

The effect of alignment on single-alfa PSFs is not investigated.

We conclude that alignment is a beneficial operation. We note that EFM has the largest errors but other arguments imply EFM is the better method to use – this is related to the stronger phase-dependence seen in non-EFM images, and is discussed elsewhere.

Stability of imaging in rapid sequence

Data reduction issues Posted on Feb 06, 2013 23:10

Our standard practice to date has been to take 100 exposures of the moon, each with an exposure time of 10 to 40 milliseconds, depending on the brightness of the moon. The exposure time is set so that the bright side attains about 50,000 counts maximum, in order to get good signal but not risk overexposure (saturation of the CCD wells).

We have looked at the differences between successive images in these 100 image stacks. A good night with a fullish moon was chosen (JD2456106)


and the relative difference between the first 50 consecutive pairs of images in the stack computed, and made into a movie — here

(this can be viewed with mplayer in linux)

An example of nine consecutive frames is shown below :

During the sequence of 50 exposures in the movie, the moon only moves about 1 pixel down and 1 pixel to the left, so drift in the position of the moon from one frame to the next is practically negligible.

White and black on the grey scale run from a -10% to a 10% difference with the previous image. Note that the two compared frames are normalised to the same total flux (this is typically a correction of ~1%).

We see plenty of structure in the difference image, which we interpret as being due to air cells — turbulence — above the telescope. The telescope size (3cm) and short exposure times (20 ms) mean that scintillation is very significant (because of the turbulence). In tests on bright stars with our telescope and 10 ms exposures, we were getting scintillation fluctuations in the fluxes of order 50% — very large — as one would expect (Dravins et al 1997a) (paper attached below).

Some of the structures look greatly elongated, and remind us of the “flying shadows” discussed by Dravins et al 1998 (paper attached below). To quote them : “Atmospheric turbulence causes “flying shadows” on the ground, and intensity fluctuations occur both because this pattern is carried by winds and is intrinsically changing”.

Analyse sequences of images in 100 image stacks, and make an estimate of the amount of turbulence present following Dravins et al’s work — for correlation with our fitted parameter alpha of the PSF — to see if the amount of turbulence in the frames has an effect on the PSF’s long power-law tail.


Dravins et al 1998:


Dravins et al 1997a (the scintillation figure is on p 186):


Dravins et al 1997b:


PSF with two parameters: Better?

Data reduction issues Posted on Jan 26, 2013 17:29

In this post: I proposed to fit more elaborate PSFs than we have been using.

I present here the results of fitting a 2-parameter PSF. The previous one was characterized by an exponent alfa, and had the form

PSF ~1/R^alfa.

The new one allows two exponents – one regulating the slope of the PSF up to a cutioff point in terms of radius, and the other parameter taking over outside that point. Thus we have 2 new parameters to fit – exponent and cutoff point.

Using such a PSF and evaluating success on the parts of the image designated by the mask – and using a mask consisting of sky on both sides of the Moon, but no polar areas – i.e. a ‘belt across the image’ – we arrive at images with slight changes in the estimated halo on the DS of the Moon and larger changes in the area on teh BS-sky.

This plot summarizes what was found:

The top panel shows absolute value of residuals (obs-EFM-model) in a slice across an image after the halo has been removed (along with most of the BS, as is always the case in the EFM method). The PSF used was the new tow-parameter PSF.
The middle panel shows the result when the usual one-parameter PSF has been used.
The last panel shows the difference between the two panels, expressed as a percentage of the first one. The illegible y-axis label for the bottom panel shows 10^-2 at the bottom, next label is at 10^0

Note that absolute values of residuals are plotted above.

We see the DS to the left, the almost-erased BS to the right and the sky on either sides.
Clearly, the BS sky has been removed much better in the first image – i.e. with the two-parameter PSF. On the DS we see changes in the range from 1 to a few percent on the sky-ward part, rising to tens of percent near the BS.

This implies that a signifcantly different amount of scattered light has been removed using the two PSFs – but which is the better result? Judging by the RMSE per pixel on the mask in the two cases there has been a significant improvement in going from one to two parameter PSFs. The RMSE per pixel in the one-parameter PSF case is about 5 [counts/pixel], while the RMSE per pixel is 0.23 [counts/pixel] in the two-parameter case. Most of this change has evidently taken place in the BS part of the sky.

We need to know if there has been an improvement on the part of the image that matters – the DS!

The best fitted parameters were quite alike – near alfa=1.7. The best-fit cutoff point was near 30 pixels.

Best of the data – revisited

Data reduction issues Posted on Jan 23, 2013 15:24

Chris thoroughly investigated the quality of all data and arrived ta a list of the 534 best images:
I have collected all those images in ‘cubes’ consisting of raw observation, EFM- and BBSO-cleaned images, images for incidence and emergence angles and the ideal synthetic image and lunar longitude and latitude images – each ‘cube’ contains centred images for a given ‘stack’ of 100 observations. We thus have 534 ‘good coadded’ images with all secondary data that belongs to that moment. Each cube is about 16 Mbyte uncompress,a nd slims down to 5-7 Mbyte upon being bzip2-ed.

I inspected the correlation of these observed stacked images with their synthetic models and found that 8 of them are poorly correlated with their synthetic models. I inspected each case and found that the problem was either focus (the observed image is blurred) or a slight alignment problem.

Omitting those images we arrive at this list of 525 best images:

These are all available as ‘cubes’, upon request.

In the future we might wish to only work with these images, since alignment and image-reversal issues have been dealt with (the observed images swap East and West upon the meridian flip).

We can set up a Mark I dataset, to be updated later as we find issues or solve the remaning (small) alignment issues.

Lunar colour map

Data reduction issues Posted on Dec 12, 2012 06:00

We have used B and V images of the moon from JD2455859 to make a map of the B-V colour of the reflected light from a thin crescent moon.

We use the photometric transformations derived from our open cluster observations:

The B band image of the moon was shifted to the position of the V band image, and the colour computed from the fluxes in B and V directly from the pixel fluxes.

The map below shows the result:

The scale at bottom runs from B-V=0.1 to B-V=1.5.

The darkside (DS) has a B-V colour in the range 0.6 to 0.9 — whereas the crescent has B-V ~ 1.1. The measured colour of the BS is B-V ~ 0.9 (Van den Bergh, 1962, ApJ, 67, 147), so we may be a little too red. More work needed on this.

We are seeing roughly the colour difference expected between ES and the BS — ES should be a few tenths bluer in B-V.

There is a clear residual halo in the colour map, indicating that the scattered light around the crescent falls off differently in B and V filters (as we expect).

Very interestingly, the crater pattern becomes much harder to see in the colour map. The V and B images look like this (after registration) V on the left, B on the right:

The lowlands appear to be slightly redder than the uplands in the colour map, by a few tenths of a magnitude in B-V.

More work needed on this — this is just a progress report!

Technical info:

a) Images used were:
2455859.1313477MOON_V_AIR_DCR.fits and 2455859.1365256MOON_B_AIR_DCR.fits

b) The transformations used were (airmass is 2.6 for V and 2.7 for B, extinction coefficients for V and B are 0.10 and 0.15 respectively).

Vinst = -2.5*log10(sumv/exptimev) – 2.6*0.10 ! airmass = 2.6
Binst = -2.5*log10(sumb/exptimeb) – 2.7*0.15
V = Vinst + 15.07 – 0.05*bmv
B = Binst + 14.75 + 0.21*bmv

which implies that B-V = 1.35 * (Vinst-Binst) – 0.43

(the airmasses above were incorrect in the first version of this entry, they have been corrected now)

Scaled Bias

Data reduction issues Posted on Nov 23, 2012 09:06

One of the reduction steps performed has to do with the scaling of the bias images due to the thermostat-induced temperature variations in the CCD chip. This temperature variation causes a 20 minute period in the mean level of the bias with an amplitude of almost 1 count – thus of importance to our attempts to analyse extremely small signal levels.

We take Bias frames on both sides of all science exposures – one just before and one just after. If we were to just subtract the average of these frames from our science image we would be adding noise to the result – we therefore need to subtract a smoother bias frame. We have constructed a ‘super bias’ frame as the average of hundreds of bias frames – it is very smooth, but probably has a level that is unrelated to the actual level in each frame.

By scaling the super bias to the mean level of the average bias frames taken at the time of the science frame we get a scaled superbias that we can subtract – it has the right level and very little noise.

We need to understand how the scaling procedure performs, so we have extracted the scaling factor from the 5000+ exposures we have.

Top frame shows the factor on the superbias as a function of the sequence number, and the bottom panel as a a function of the Julian Day of the observation.

Most factors are near 1 but some stand out. 9 files have factors above 1.09 – their Julian days (integer part) are:
2455938 (7 images) 2455940 2456032 (one on each).
A list of the 209 images with factors over 1.08 is here:

The 4 unique JDs are: 2455814 2455938 2455940 2456032, with the majority of cases on 2455814.

These images should perhaps be inspected very carefully for problems with bias.

A close-up of the factors nearest 1 looks like this:

The best of data, the worst of data ..

Data reduction issues Posted on Nov 22, 2012 10:28

Using the methods described by Chris in the entry below – i.e. at

I have also selected for the likely good data by taking extinction into account and looking for linear sequences. This resulted in a set of images that I deem ‘good’. That list can be compared to Chris’ list.

Chris has 3162 ‘good’ images, while I have 2990. The cross between the lists finds 2273 instances on both lists. This list of ‘jointly agreed good images’ is here:

Methods for selecting good images differ slightly: Chris does not yet consider the alpha value found in the EFM method – I select for alphas in a narrow range near the mode of the distribution for each filter.

With converging selection criteria the list above would expand somewhat – perhaps best to keep the criteria a little different to avoid duplication of potential errors.

The joint list contains 49 unique nights. The distribution of images over nights is:

58 2455856
363 2455857
404 2455858
33 2455859
228 2455864
65 2455865
276 2455917
173 2455923
117 2455924
11 2455938
6 2455940
15 2455943
10 2455944
5 2455945
2 2456000
3 2456001
8 2456002
9 2456003
3 2456004
6 2456005
4 2456006
5 2456014
31 2456015
54 2456016
34 2456017
3 2456027
4 2456028
16 2456029
11 2456030
10 2456031
7 2456032
15 2456034
2 2456035
26 2456045
49 2456046
27 2456047
20 2456061
2 2456062
5 2456063
3 2456064
35 2456073
21 2456074
31 2456075
11 2456076
7 2456090
12 2456091
9 2456092
6 2456093
18 2456104

where the first column gives the number for the night. Note that towards the end of the sequence there are few images per night – that is because we were realizing the necessity to observe ‘stacks’ rather than sequences of single images which results in fewer co-added images.

A histogram of these data is shown above.

Identifying the best of the data

Data reduction issues Posted on Nov 21, 2012 23:35

We have been using the apparent magnitude of the moon compared to the JPL model of its expected magnitude (V band), to check which images suffer from gross problems, such as an incorrect filter in the beam, or light losses from thin cloud, or other reasons. This follows up on previous posts (,

This plot shows the difference between the measured apparent V magnitude of the moon and the JPL magnitude (i.e the difference between instrumental and true magnitude) as a function of airmass. There are about 5000 frames, mainly of averaged data (i.e. up to 100 images averaged together) in the total data set (as of November 2012), of which about 1000 are in the V band, and plotted here.

(a pdf version is also available).

Blue represents one side of the moon (the side with Mare Crisium), green the other side. Opposite sides of the moon have differing distributions of uplands and mare, leading to about a 0.15 magnitude offset between the reflected light and the JPL model (which implicitly assumes a uniform disc). The dotted lines show the expected trend with airmass for an adopted 0.1 magnitudes/airmass extinction — quite close to what we are seeing in the data. There is a lot of structure in this diagram!

To clarify which side of the moon is which: the next plot shows the Crisium side illuminated by sunlight:

This plot shows the Grimaldi side being illuminated.

Now we concentrate on the crescent light on the Crisium Mare side (blue points in the plot):

Concentrating on just this side of the moon, we see a very tight (extinction) relationship with airmass, with notable outliers. They are all on the dim side of the relationship — indicating loss of light relative to the JPL magnitude). The reason for most of these has been identified — most often because of thin cirrus or haze (this affects the scattered light around the bright side of the moon, which we model with a power law with slope “alpha” — these are indicated on the plot with that symbol. The moon went behind the tower and cables on the western horizon in one case “cables” and one image has low S/N. Two have “?” — no reason was found for these outliers, so they are simply dropped. This leaves us with a nice looking sequence of data for this side of the moon. I’ll look into these further as there may still be structure in there as a function of lunar phase.

The other side of the moon, the side with Grimaldi and extensive Mare, shows much less clear trends with airmass than the other side.

The dotted lines are the same as in the previous figure — they bracket the good data on the Crisium side. Strange things are happening at low airmass (z<1.2) – it turned out that most of these peculiar data are from a single night (JD2455912) as the next figure shows (background of these points shown in blue). Note the airmass scale has changed to make these points easier to see.

The entire night of JD2455912 seems to be a dud — as similar strange happenings show up in the other filters — VE1, VE2, IRCUT (but not clearly in B, as only one good B image was taken on that particular night).

Analysis of the remaining (grossly) discrepant points shows that they are often all on a single night. The next plot shows some of these problem nights:

Blue: Totally discrepant data from JD246691 : drop these.

Red: Near full moon!

Purple: data from JD2456061,2,3 (three nights). For almost all, the moon is too bright by 0.3 mag — possibly bias issues? Will look at this.

The bulk of the green points lie along the expected extinction line, with an offset relative to the other side of the moon (Crisium) of about 0.15 mag. We can bracket these data and use them for the ongoing analysis of earthshine. We look at this in the remainder of this post.

The next plot shows both sides of the moon in green and blue crosses, but for all filters (instrumental magnitude – JPL V-magnitude) as a function of airmass. No data have been dropped yet — it’s all the nights we have.

The colour of the moon in the various bands shows up nicely as an offset relative to JPL V-magnitude (e.g. B-V is of order 0.9, second panel from bottom). The Crisium side of the moon shows much cleaner results than the Grimaldi side, in all colours. The dotted lines bracket the good Grimaldi data, and are the same in all panels. Data which appears between these two lines in other bands then V, mean that the V filter was actually in the beam, not the nominated filter. There may be similar effects for all filter combinations, but for the moment we can only pick out misfirings which resulted in the V-filter being in the beam. Only a small fraction of the data seem to be so affected, which is very good.

Extinction depends on the wavelength of the filter : this is clearly seen as differing slopes in the zeropoint (ie. the difference between instrumental and true (JPL) magnitude) with airmass.

Once the grossly erroneous observations/nights are cleared out of the data, the sequence of good data as a function of airmass can be picked off the plot (we call this tunneling).

Firstly, opposite sides of the moon are reduced to the same scale, by adding 0.05 mag to the magnitudes (in all bands, B, V, VE1, VE2, IRCUT) on the Crisium side, and subtracting 0.15 mag from the Grimaldi side photometry. This brings both sides onto more or less the same sequence with airmass. We then estimate (by eye) the extinction as a function of airmass in each filter. The results are:

B 0.15 mag/airmass
V 0.10 mag/airmass
VE1 0.08 mag/airmass
VE2 0.06 mag/airmass
IRCUT 0.12 mag/airmass

Interestingly, the extinctions in IRCUT and VE1 are somewhat different, which we wouldn’t expect as the filters are ostensibly rather similar. We should look into this further. In any case, these extinctions allow us to define lines about 0.3 mag wide around each sequence for each filter and pick off the best data. The result is this diagram:

A list of these frames is here:

There are 3162 frames on this list — cut down from about 5000, so about 60% of the frames have survived the “tunneling” process. This is still a quite conservative list because I haven’t

1) checked for the Lunar centering,

2) checked for “alpha”, the scattered light parameter, except for those very obvious cases which showed up by eye.

Final plot: using just the good data, the extinction corrected apparent magnitudes at MLO are compared to the JPL model as a function of Lunar phase. The scatter is small in some of the filters — we are down to absolute photometry errors of just a ~ 0.1 mag relative to JPL, which is pretty good considering we cover a range of airmass mainly from 1 to 5, but with some of the data out to airmasses of 15!

There is a hint we can do better still, because there is still structure as a function of phase — for example, the B band magnitudes at negative phase are a bit brighter still than the positive phases — indicating we didn’t get the correction for each side of the moon quite right. There are also trends in the reflected light with phase, but this probably reflects the fact that phase and airmass are somewhat correlated (small absolute phases tend to get observed at higher airmass, because the moon starts the night closer to the horizon).


Data reduction issues Posted on Nov 04, 2012 13:29

While waiting for the CCD to be fixed, I am summarizing the discussions we have had on this blog. The material should be used largely in section 6 of the manuscript we submitted, leaving the other parts of that paper – or splitting it into tywo papers. Here is the submitted manuscript as well as the start of a summary the summary is really just a collection of links to this blog.

The Daily Progression

Data reduction issues Posted on Sep 28, 2012 15:25

In this post, we noted that there is a daily progression in the obs/mod ratio of the DS/BS ratio itself (let’s call this thing ‘the relative albedo’ from now on!). We speculated as to its origin.

We have now calculated the “relative albedo” with two sets of terrestrial models – one a uniform Lambert-reflectance sphere and one a non-uniform sphere with ocean-land contrast but with a Lambert reflectance, while keeping the set of observations constant. We compare these two sets of results:

The upper panel shows the “relative albedo” for a non-uniform Earth and the lower one for a uniform Earth. The linear regression slope is printed on each line.

We note that the ‘progression’ for V has changed sign in going from a non-uniform Earth to a uniform earth. The magnitude of the slope on B has changed by 50 %.

We conclude from this that the ‘daily progression’ is due to the terrestrial modelling of surface albedo on Earth.

Or: we can state that we are seeing actual surface albedo variations during a nightly sequence of observations!

During this set of observations (night of JD 2456046) The B-band albedo of Earth decreased slightly, while the V-band albedo rose. This is the albedo of the half-Earth shown in the figures in the post here.

For fun, we plot here the B and V relative albedos for tow nights with overlapping phase intervals:

The two nights are 2456075 and 2456016. For the B progressions it looks like there is a difference, at phase -96 degrees of about 0.15 in the relative albedo, and given the scatter of the points we might be able to confidently see a difference 1/10th of that, I guesstimate. That is 1.5% at a given point in phase (that is, time).

BBSO discuss ‘1-2% per night’.

If we take the above results at face value we seem to be able to track albedo changes through a night. Averaging sucha progression, in order to get the ‘nightly mean value’ would, of course, yield an average value, and an average can be affixed with a ‘standard deviation of the mean’ which can become very small if you have enough points – and if it gives any meaning to average the points. It does not give that much meaning to average the above progressions since they change on scales longer than our observing sequences.

But anyway, we need to start understanding what BBSO does in their data-reduction. Notably – do they assume Earth is a Lambert sphere? If so, we have perhaps arrived at their “A*” quantity? If, even more, so – we now just need to observe for a few decades!

Given only 1.5 years of data we should for now focus on how we present the results, and how we get something ‘meteorological’ out of the data – as opposed to ‘climatological’. For instance, does the level-difference between nights 2456016 and 2456075 correspond to a change in albedo in the relevant areas of Earth that we can determine from Earth observation data from satellites? More clouds, perhaps? 2456016 is the upper sequence.

Let’s see! (more to follow…)

Lunar phase and azimuth

Data reduction issues Posted on Sep 27, 2012 12:49

For all our observations, we plot the lunar phase (in the convention we use) against the azimuth (degrees East of North) of that observations.

We see that most positive phases are observations East of the meridian; most negative phases are West. Only rarely has the Moon been followed past the meridian.

A look at some first results: B and V

Data reduction issues Posted on Sep 27, 2012 11:38

We consider now the DS/BS ‘ratio of ratios’ – that is, the DS/BS ratio in observations divided by the DS/BS ratio in the model images. We have used the Clementine albedo map scaled to the Wildey levels, and Hapke 1963 reflectances in the model images.

We look at the data from several nights in the phase-interval from -130 to -50 degrees (Full Moon is at 0 degrees; negative phases implies we are observing the Moon in the West (setting)):

What we see are short sequences of points – each sequence covering one night. We note how many of these nightly sequences ‘dip’ during the progression of one night. We note the smoothness of each nightly sequnce – this bodes well for our original idea that albedo could be extracted with high precision. Jumps from night to night are more troublesome! But we need to first understand the ‘progression’.

Inspection of the time-stamps on the data in a ‘progression’ shows that time moves to the right – that is the ‘ratio’ plotted starts out high at a larger negative phase and ends as a lower ratio at a less negative phase (moves left to right). Since the observations are in the West we are looking at a setting Moon, or increasing airmass along a progression.

On a single night several things happen:
1) the Moon is seen through different airmasses, and
2) the Earth turns.
[recall that the synthetic models contain all the correct distances as functions of time, so effects due to that should divide out].

1) If the effect is due to the increasing (or decreasing) air mass we need to understand how a ratio of intensities taken from two areas on the Moon, obtained from the same image, can show a dependence on airmass. Differential extinction across the lunar disc? really? We are looking at 10-30% changes in a single progression!

2) Since these observations are all at negative phase they may all represent the same configuration of Earth and Moon – i.e. we could be seeing the effect of the same parts of the Earth rotating into view through the night. From a fixed observatory there is a tendency for the same areas on Earth to reflect light to the Moon on subsequent nights – in our case either the American sector and oceans on either side, or the Austral-Asian sector.

We should not forget 1 – but let us turn to 2 for now. We can extract ‘scenarios’ for the Earth at the moments of observation corresponding to the data-points above. Inspectingh the observing scenarios could perhaps teach us what is going on.

We select the progression near -91 degrees. That is a set of observations from JD 2456047.75 to 2456046.87, or about 3 hours (the date is April 30 2012). At the start the Earth was in this configuration (crosses on the globe show the part illuminated by the Sun):

While this is the configuration at the end islike this:

That is, the Earth has rotated slightly during the observations and
more of the Asian landmass has swung into view, meaning that more of the
earthshine is coming from the Asian landsmass, than at the start.

Can we estimate the expected change in Earth’s reflectivity during this sequence? Yes. We have quite elaborate code for that, know as “Emma’s code” due to a student who worked on it in Lund. It shows the Earth from any angle at any time, with superimposed albedo maps of clouds as well as continents. A reflectance law is imposed. I think there is even some Rayleigh scattering added to give the ‘bright edge’ of the Earth (see any Google Earth image at startup, as you zoom in – the edge is bright, and that may even be realistic).

That code is quite complex and is not run in a jiffy. [We need more students!] But we can make a perspective model of the Earth using simple IDL routines. We can wrap surface albedo and cloud maps onto a sphere and view it from the angle the Moon is at, for any given time. There is no reflectance imposed – just the albedo maps and the perspective view. The two situations above look like this with cloud maps for the moment of observations (taken from the NCEP reanalysis product TCDC):

The simple average of the brightness of the two images above are 42.6 for the former and 40.6 for the latter – so it is getting dimmer by about 4% as time passes – perhaps because bright Australia (cloud-free at this time) is entering dusk. The cloud map images are available at 6 hourly intervals only, so there has not been time for the cloud image to change – it is simply rotated and other bits have become hidden by the day-night sector advancing.

Brightness of opposite sides of the Moon

Data reduction issues Posted on Sep 19, 2012 11:10

We have isolated the very best images in V band, in which the Moon is well centered and well exposed.

We measure the apparent magnitude of the moon in these images by simply measuring the total flux, and applying the standard photometry relations (i.e.

We correct for extinction of 0.10 mag/airmass (from

We then compare the apparent magnitude to the expected apparent V magnitude from the JPL ephemeris for the Moon ( This uses the relation quoted in Allen “Astrophysical Quantities”, which actually comes from Eqn 8 of this paper:

The plot shows the difference in the apparent magnitude as a function of phase (new moon = 0). Blue and green show opposite sides of the moon.

There is a bit of scatter in these data, but there are two clear sequences around phases 50 to 100 showing that opposites sides of the moon differ in luminosity by about 0.1 mag. This is quite a lot less than we expected to see, viz. this post:

Lunar apparent magnitude with phase

Data reduction issues Posted on Sep 11, 2012 05:22

The plot shows the apparent magnitude of the moon in V and B as a function of lunar phase (phase=0 is new moon).

We measured the flux in images in which the filter was reliably V or B and used the transformations determined from NGC6633 (i.e. to get the apparent magnitude.

These data have been obtained for a large range of airmass (z) — from z = 1 to 10, with most of the data in the range z = 1 to 3. We derived extinctions of 0.10 mag/airmass for V and 0.17 mag/airmass for B, by comparing to the apparent V magnitude from the JPL ephemeris for the Moon ( (more about this below). The solid black line shows the V apparent magnitude as a function of
phase after extinction correction, and adjusting the zeropoint by 0.2
mag in V to fit.

Note that B is ~ 1.0 mag fainter (i.e. B-V ~ 1.0, as we’d expect).

The airmass fits are shown above : the plot is the difference between the apparent magnitude from the JPF ephermeris and our transformed instrumental V band (or B band) magnitudes, shown as a function of airmass. The two lines show 0.10 mag/airmass (V band) and 0.17 mag/airmass (B band) — they are not fits. There are some bad outliers, especially in the V band, which are probably due to the incorrect filter being in the beam.

CCD camera rotation

Data reduction issues Posted on Sep 05, 2012 09:09

As discussed here, the CCD camera became twisted in its thread on the telescope at one point. The problem was fixed, but this means that some of our images have a slight rotation about the frame mid-point. This influences the success of subseqeunt data reduction steps: especially the steps that depend on extracting flux from specific areas on the lunar surface.

We therefore tested for the presence of a rotation angle by correlating a synthetic image made for the observing moment with each and every observed image, rotating the synthetic image until the correlation was maximum – in 1 degree steps.

We plot the detected best rotation angle as a function of image sequence number and date:

Top frame: detected image rotation angle vs sequence number, Bottom panel: angle vs observing day since start.

It certainly seems that almost all images up to number 1800 or so has a rotation angle of some -7 to -8 degrees. That seems to correspond to just a few nights near night 40-50. The detection of rotation is a bit spotty so there are also other episodes where a rotation angle other than 0 is detected – such as images 2000-2500. That more intermittent episode corresponds to a few nights near night 180, but there are a few more examples near night 220.

The CCD twist was correctd by Ben on JD2455991, and this datum is shown as a vertical dashed line in the plot above. Since this is not consistent with the angles measured we have to say that the test so far has been inconclusive!

Added later:

Actually, it was not impossible to inspect the relationship between model images and observed images visually and to confirm when an obvious image rotation was present. Partial results (note: more points than above) look like this:

Here color coding indicates in red the images that so far obviously have a rotation problem, and in blue images that show no obvious problem.

The presence of blue symbols at large rotation angle must be due to failure of the algorithm for detecting rotation! It is not an easy problem to solve – at New Moon there is precious little to correlate images on, unless the DS is used – but the presence of the halo gives problems, so that histogram equalization is not an obvious remedy.

A fixed derotation for the detected nights could be implemented – this affects some of the early observing nights where single images were taken (not stacks).

Added even later:

By manual inspection and image comparison, the following de-rotation angles were found for the JD in the beginning of our sequence:
2455856.1078938 7.711
2455857.0817247 6.881
2455858.0931277 6.881
2455859.1269613 0.000
2455864.7037639 5.352
2455865.7157216 6.116
2455886.0356274 0.000
2455905.9722493 0.000
2455912.0991783 0.000
2455917.1285750 0.000
2455923.7124300 0.000
2455924.7257543 0.000

Apart from 2455859 all angles before 2455886 were clear to find. 2455859 was hard to inspect as it is very near New Moon and almost no features were detectable.

As a working hypothesis, let us assume that all images before 2455886 must be rotated by something like 6 or 7 degrees, to bring them into good alignment with their synthetic models.

First results

Data reduction issues Posted on Aug 31, 2012 14:00

Here are some of the first results from applying the BBSO linear method – the reductions are slow so more will be on hand later! Sorry about the low quality in the image – there must be a way to do it better, but …

Rows 1-5 give results for each filter. Ignore columns 1,2 and 3 for now – they are diagnostic. Column 4 shows the ratio of observed terrestrial albedo to modelled terrestrial albedo as a function of lunar phase.

We pick out the DS as either Grimaldi or Crisium depending on which is in the DS. We then calculate the observed DS brightness divided by the brightness of the whole disk. We then do the same for the synthetic model, and plot the ratio of the two.

This may seem a strange quantity to plot, but consider that in the (unlickely case) that we both had perfect observations and the model was correct in all aspects, then we would see a ratio of 1.0.

If the model is somehow wrong – for instance if the phase function it is based on is unrealistic then the ratio would have a phase dependence.

As it is, we do not have perfect observations and we see a fair bit of scatter. The scatter comes about for several reasons – first of all the observations have Poisson noise – we are extracting a small 4×4 pixel patch on what is the dark side where counts are probably on the average of 5 or so. Additionally we have noise from the alignment between the actual Moon and the coordinates we have calculated from which to extract information – there is a missmatch of up to several pixels here, so for a small area like Grimaldi a few pixels error in placement brings you into the bright surrounding areas. For Crisium, which is larger, this is less of a problem.

Finally there is a still un-solved problem with synthetic models and observations apparently being off by some small amount in terms of a small rotation. This may be from the days when the CCD camera was actually physically twisted by a few degrees in its placement on the telescope. In such cases the intensities of the pixels extracted in the synethtic model and observed image are even more different.

So, some things to work on are:

1) Use larger DS patches, so that the Poisson statistics are not as much of an issue.

2) Put the patches in uniform areas on the Moon so that missalignments do not cause acquisition of contrasting areas. Inside large, even Mares or on the brighter highlands.

3) Use better estimates of disc centre and radius.

4) Figure out a way of aligning the synthetic model and the observed image better.

PS: More data are available all the time so the figure above will update now and then.

Added later:

Have implemented 1 and 2, by enlarging the area inside Crisium that is used, and using a rectangle in the highlands south of Grimaldi instead of Grimaldi itself. Also trying 3.

Differential refraction

Data reduction issues Posted on Aug 29, 2012 15:16

We have to determine the radius and centre of the lunar disc in order to reduce observations.

In doing that we must be aware that differential refraction causes the Moon to appear non-circular as it comes closer to the horizon. Using a formula from the Nautical Almanac Explanatory Supplement we generate the following table for 600 mmHg, 10 degrees C and 30% relative humidity:

Z d_refr am
27 1″ 1.1
69 6″ 2.8
75 12″ 3.9

where z is the zenith distance in degrees, d_refr is the differential refraction in arc seconds and am is the airmass. The differential refraction is calculated over a 1-degree distance centred on the given zenith distance. Our FOV is about 1 degree wide and one pixel covers about 7″.

We thus see that the Moon is differentially refracted by less than one pixel up to about 2.6 airmasses. Two pixels are reached about 26 degrees (4 airmasses) above the horizon.

Some of our observations are certainly close to the horizon, as we have tried to observe when the lunar phase is near or less than 30 degrees (at Newish Moon).

At 30 degrees the uncertainty in our determination of disc radius and disc centre (based on a circular assumption) is thus starting to be challenged by the differential refraction.

Driving on the Moon

Data reduction issues Posted on Aug 29, 2012 14:14

In order to extract fluxes from particular areas of the Moon we need to take lunar libration into account. Based on the synthetic model code Hans wrote we can do that now.

On the left is the synthetic model. On the right is the observed image. On it are some dots – they are supposed to be inside Mare Crisium and Mare Nectaris and two other locations on the darkened side. We see that we nail Nectaris and Crisium. Doing this for every observed image we will be able to extract feature fluxes and compile DS/BS statistics.

Fluxes vs phases

Data reduction issues Posted on Aug 27, 2012 16:22

Since the setting of the colour filter (as well as the shutter) was unreliable we must find a way to detect which images are taken through which filters.

Here we plot the raw fluxes (counts divided by nominal exposure time):

Total lunar Fluxes plotted as magnitudes against the lunar phase (New Moon is at 0). Bright is up, faint is down.

The data for each image (black symbols), here expressed in magnitudes, are overplotted with a phase law (red) inspired by that in Allen “Astrophysical Quantities”, except we modify the coefficients in that and use instead:


Notably the coefficient on the linear term is about half of what Allen specifies.

Particularly in VE1 we note the presence of two sequences of data. We have ‘fitted’ (by eye) the sequence that is represented on both sides of the new moon to the Allen phase law. The orphan sequence is below the adopted data suggesting that a filter with less transmission was obtained here when VE1 was requested.

For B it seems that intermittency causes some fluxes to be higher, but they do not fall in a delineated sequence so cannot be identified with a filter.

The V filter seems to have the same problem, although less so.

VE2 is also somewhat ‘broad’ in its distribution.

IRCUT shows two sequences.

The VE1 and IRCUT filters are extremely similar in transmission properties and are quite similar in the plot above, including the presence of the ‘second sequence’.

The ‘second sequence’ is quite similar in flux to the V observations, and it is consistent to say that when IRCUT and VE1 failed to be chosen V was obtained instead.

This seems to also apply to some of the outliers in VE2 and B.

We thus suggest that a working hypothesis for the failure of the filter wheel is that when it failed the V filter was selected instead.

We next proceed to eliminate the outliers in each filter so that a dataset can be defined which will allow identification of the extinction laws in each filter. With that in hand it may be possible to ‘tighten up’ the data and move towards a ‘golden dataset’ from which also DS fluxes are worth extracting.

The presence mainly on the left side of the diagrams of a ‘second sequence’ of data implies that something like a mechanical problem is behind the FW failure – because the telescope is flipped over the meridian to observe mainly in the East or the West depending on whether the Moon is rising or setting (before or after New Moon).