Some things to note.
A quick look at the full data set Tom sent me.
The Orion data set (part); the set analyzed by Michael Richmond.
Spatial correction.
Standard deviations.
The full Orion data set.
Details of methods.
I am still riding the PSF fitting hobby horse and am not even set up to run the standard analysis pipeline. This investigation is thus a largely independent look at the data as compared with that presented by Michael Richmond in TN #89.
I processed the same set of data, provided by Tom Droege, as that used by Michael Richmond. I used my latest PSF fitting routines: these fit "Type A" distributions - that is, Gaussians multiplied by polynomials.
Tom provided me with 3 sets of data each consisting of 49 images each in V-band and I-band.
The following picture shows my fitted PSF widths in X1 (Declination) and X2(Right
Ascension) for the 147 successive images.
Some things to note.
1) These are the first MK IV images that I have seen with the I images sharper than V. The V stars are pretty much circular in the centre of the image but show strong variation north-south across the image. This probably indicates tilt of the lens or of some of its components. The much sharper I stars are elongated in X2 or Right Ascension suggesting, perhaps, a small degree of trail. Gaussian widths barely greater than one pixel are marginal for aperture photometry.
2) The plot is a bit spiky. Some of the positive spikes are the result of mysterious double images - Tom blames these on Right Ascension drive problems but the splitting of the stars is not E-W as one would expect. Some of the variation probably indicate clouds.
Next picture is the preliminary photometric calibration against Tycho2 magnitudes modified to TASS V-band and I-band in the usual way. Magnitude = Base Magnitude - 2.5 Log10(total adu).
I have flagged the measurements for the detected double-star images. These account for some of the negative going excursions but not for all of them. One would suppose that the remaining negative going excursions indicate clouds.
A seven image periodicity is particularly clear at the start of the I-band observations but is much less convincing later. Images were taken at seven successive declinations for each right ascension so this appears to indicate a declination or zenith angle dependence. Transparency improves further from the horizon, which is hardly surprising, but the size of the effect seen (0.05 magnitudes or so over only four degrees) is rather large? Maybe a bank of haze!
The middle set of the three, in the constellation Orion, has been analyzed by Michael Richmond in TN #89: here is the data replotted for comparison with his.
This plot is quite similar to Michael Richmond's but with some obvious differences. The negative excursion for the V-band calibration of image #21 (Michaels's numbering), on the strength of which he rejected that image, is not reproduced. On the other hand, I get a much larger negative excursion for #41. The photometric calibration uses a large number of Tycho2 stars and it is hard to see how the uncertainty introduced by using a slightly different set of Tycho2 stars could account for more than a fraction of the observed scatter: The formal error in the photometric calibration is well below 0.01 magnitudes. It is just not obvious how the difference between the two analyses can be so large!
Finally the estimated position accuracy.
The quantity plotted is the rms residual of the final astrometric fit. Again there is
a hint of periodicity. Again there are spikes, positive this time, some of which
correlate with the double-star images.
And a new feature: the position errors are larger for the 3rd (right hand)
data set. The star density is higher for this data set so the additional variability
is presumably the result of confusion.
Note that the smaller PSF width of the I-band images combines with the greater brightness of the sources in I-band to give much better positional accuracy.
This is the middle of the three sets previously shown and is the same set analyzed by Michael Richmond. For this preliminary analysis, I have selected an area of approximately one degree square in the centre of the area covered. This reduced the volume of data to something I could manage in a spreadsheet. The downside of this is that the stars are seen in various positions on the different images and do not provide good spatial coverage on any one image. The full Orion dat aset is analysed further below.
I used three different Flat files for the processing.
First a null flat that applies no correction; second the flat Tom provided with the data set and third a light box flat made from data Tom sent me a long time ago.
Here is V-band data for stars brighter than magnitude 12.0 for the 3 cases.First the null flat:
Um! Not too similar to Michael Richmond's version (unflattened data: new version of TN #89). The results should differ from his only because he used an aperture sum and I used PSF fitting. Ha! His plot is labelled as having West to the right while I think my plot has West to the left. Nope! Flipping doesn't help much. PSF fitting is different from the aperture sum.
Again, not too similar to each other or to the first plot in Michael Richmond's TN89 which, I think, uses the same Tom Droege flat. Remember, when comparing with TN89, you have to flip the plots to get East and West to match but this does not do much good: the pattern of residuals is quite different! For one thing, the north-south gradient is much more prominent in my results than in Michael Richmond's. At first, I found it hard to believe the size of the difference between the two analyses. However, my fitted PSF also shows a strong north-south variation in the FWHM; the lens or some of its components must be tilted. It is not too unreasonable that fixed sized aperture photometry should give different results compared to PSF fitting in this case. The size of the difference is much larger than I had expected.
There are several other things to notice about both plots:
There are many outliers; note the sprinkling of the larger symbols of both colours (signs) across the middle of the plot where most stars give small residuals plotted as black dots. I have already pruned all the really big residuals, greater than 0.2 magnitudes, before preparing these plots There are, I think, far more outliers than in Michael Richmond's corresponding plots, which suggests that my PSF fitting routines are still rather temperamental! Theoretically, PSF fitting should, of course, give less scatter. That's why I keep trying to get it to work!
Some large positive residuals are the result of my program trying to split bright sources into two components. Even after pruning the very large residuals, some of the remaining large residuals probably arise in this way. This effect can give errors of any size. An error of two magnitudes was noted where a source had been split and the weaker component managed to get selected. Small errors occur when the second source (real or imagined) is faint or not very close. All analysis programs suffer from this problem to some extent. It was the major problem restricting the usefulness of the MK III data and will remain so with the MK IV. I get the impression that my PSF fitting program is even more liable to this problem than the standard aperture sum program.
Many of the larger magnitude residuals of either sign are for sources showing larger fitting errors than the same source on other images. There is thus some prospect of flagging and removing these measurements. This has not been done in this case. My program does flag large fitting errors but uses a fixed threshold that currently fails to catch more than a small fraction of cases. Why the fitting errors are large is an open question.
I kept for plotting only stars that turned up on all three analyses of an image. A surprisingly large number of cases, even of bright non-saturated stars, turned up only once or twice out of the three analyses. This is a mystery that I have not yet had time to explore. Many more cases had large residuals ( > 0.2 magnitudes) in one or two analyses and so got rejected. Does this sensitivity to rather small changes happen with the standard pipeline or is it just a "feature" of the Bennett PSF fitting routines?
All the plots show large scale structure but to differing degrees: the first plot, with the null flat correction, shows the least while the other two, with flat correction, show more. Qualitatively, the patterns are similar but with considerable quantitative differences. With my browser, I can blink the plots using Page-up and Page-down. For example, Tom's flat gives larger areas of consistently large residuals: negative (too bright) in the N-W and N-E corners and positive (too faint) in the South centre.
Now the corresponding results for I-band:
These I-band plots show the same features:
The patterns are different.
The null correction shows the least structure.
There are lots of outliers. They do not all come up in the same places on all three plots.
I fitted spatial polynomials to each of the six sets of residuals plotted above and computed new sets of corrected residuals by subtracting the fitted polynomial. As the order of the spatial polynomial was increased, the visual patterns became less visible and at fourth order (15 coefficients) the residual plot appears random. Sorry - I didn't keep the intermediate plots and am too lazy to go back and recreate them. And the final fourth order plots for this data set are completely dominated by random noise so I am not bothering to show them either. However, my spreadsheet has provided me with plots of the rms scatter of individual sources, which I show below. Note that these estimates of rms scatter are seriously affected by the large residuals - especially the ones I accidentally left in. I should really be using other methods which behave better when afflicted by large residuals. But life is short and I will do what my spreadsheet allows me to do easily!
The plots are a bit confusing, so first I had better tell you what you are supposed to see in them! I have plotted only two out of the three different flat analyses: the light box flat gave results similar to those for Tom's flat. The plots show, or so I claim, that without spatial correction, the null flat processing is superior both for V-band and for I-band. Spatial correction not only reduces the scatter for both flats but reduces it to similar values. Look closely at the V-band data and you will find cases where Tom's flat comes out best after spatial correction and a more or less equal number where the null flat does best.
Results are more clear-cut in I-band: here Tom's flat, while clearly the worst before spatial correction is clearly the best after. Well - most of the time. You are to ignore the sporadic high residuals which are to be blamed on those mysterious outliers.
Both plots show a noise floor at a depressingly high level of scatter - perhaps 0.016 magnitudes in V-band (discounting all the high outliers) and around 0.012 in I-band. Spatial analysis of the flats and of the background level in individual images was carried to 5th order spatial coefficients: 21 coefficients, most of which proved to be significantly different from zero. Indeed, the magnitudes of the coefficients did not appear to be converging to zero even at 5th order. I carried the current investigation only to 4th order because this was the lowest order fit at which the residuals did not show obvious structure and increasing the order of fit is somewhat laborious and of doubtful utility in view of the problem of the large residuals.
To investigate the noise floor further, I selected the brightest sources, brighter than magnitude 10.5, and did the whole thing over, first getting rid of measurements with residuals greater than 0.05 magnitudes. The residual plots are just like those shown above except with far fewer points and far less scatter. The raw residuals show the same patterns as those I have already shown you so I have not plotted them. However, with the smaller scatter, one can now see spatial structure in the corrected data even with the 4th order fit. Here is the V-band data for the null flat and for Tom's flat.
These plots show almost identical spatial structure apart from a few of the apparently inevitable large residuals. The light box flat (not shown) is also almost identical. I initially found this entirely unexpected since it appears to imply that the small scale structure with the 3 flats - the structure not removed by the polynomial fit - is essentially identical even though the large scale structure is not. A little reflection shows that we are actually seeing mostly the time dependent part of the correction. Some of the structure is removed by subtracting the mean residual for each image. I also tried subtracting linear trends estimated for each image and this gave a further slight improvement. This is followed up using the full data set below.
To give some idea of the magnitudes, I have collected various standard deviations computed in the above analysis. I must emphasize again that errors are dominated by the large errors such as splitting sources on one image while joining them on another and that these sources of large residuals have been eliminated manually with greater or less success. Most of the figures quoted are actual standard deviations but where it really matters I have done the extra work to compute probable errors so as to minimize the effect of the large residuals. I have used the conversion
Standard deviation = probable error/0.6745
After all corrections are applied, the residuals come out almost identical, with a correlation sometimes greater than 0.9. But the correlation is not perfect at least in part because of the differences in the flat fielding. Tom's flat and my light box flat should both compensate for pixel-to-pixel variations in sensitivity whereas the null flat obviously does not. This shows in the statistics of the residuals and of the differences between the different sets. Specifically, the null flat results should show higher standard deviation because of the uncorrected pixel-to-pixel errors and should have lower correlation with either of the other sets than they have with each other. Estimating the standard deviation from the probable error of the differences largely eliminates the effect of the large residuals. I also show the standard deviation calculated directly so you can see how well (or otherwise) I have taken out the large residuals.
Null-Tom Tom-Box Box-Null
V, mag < 12.0
r 0.923 0.910 0.881
sd from pe 0.0084 0.0073 0.0089 magnitudes
sd direct 0.0136 0.0149 0.0173
V, mag <10.5
r 0.894 0.924 0.915
sd from pe 0.0049 0.0037 0.0048
sd direct 0.0067 0.0057 0.0058
I, mag < 12.0
r 0.926 0.984 0.876
sd from pe 0.0077 0.0060 0.0084
sd direct 0.0132 0.0102 0.0134
I, mag < 10.5
r 0.797 0.864 0.790
sd from pe 0.0051 0.0034 0.0053
sd direct 0.0067 0.0055 0.0070
All four sets show the expected pattern in the standard deviation calculated from probable error. Directly calculated quadratic quantities, standard deviation and the correlation, r, are more badly affected by the remaining large residuals.
Assuming these residual errors are uncorrelated, one can estimate the individual contributions. For example, that last set (bright I-band sources):
Null-Tom Tom-Box Box-Null
sd from pe 0.0051 0.0034 0.0053
gives
Null Tom Box
sd 0.0046 0.0022 0.0026
V data using Tom's flat
M < 12.0, gross outliers removed s.d. = 0.0553 magnitudes
minus spatial terms to 4th order 0.0376
difference gives spatial terms as 0.0408
M < 10.5; more outliers removed s.d. = 0.0469 magnitudes
minus spatial terms to 4th order 0.0181
difference gives spatial terms as 0.0432
minus image mean 0.0149
difference gives image mean s.d. 0.0103
minus image linear trends 0.0120
difference gives trend s.d. 0.0084
I data using Tom's flat
M < 12.0, gross outliers removed s.d. = 0.0463 magnitudes
minus spatial terms to 4th order 0.0270
difference gives spatial terms as 0.0376
M < 10.5; more outliers removed s.d. = 0.0398 magnitudes
minus spatial terms to 4th order 0.0123
difference gives spatial terms as 0.0378
minus image mean 0.0101
difference gives image mean scatter 0.0070
minus image linear trends 0.0087
difference gives trend scatter 0.0050
That final residual standard deviation, 0.0120 magnitudes for V-band and 0.0087 magnitudes
for I-band, then includes
Photon counting and electronic noise: around 0.0080 in V-band and 0.0050 in I-band for
magnitude near 10
The flat fielding noise estimated above as, for example, 0.0022 magnitudes for bright
sources in I-band if using Tom's flat.
Residual constant spatial structure not fixed by the 4th order polynomial.
Residual image-to-image changes not included in the mean and linear gradients.
Anything else we haven't thought of yet.
And finally, real variation of the brightness of the stars.
Gee! It's beginning to look as though we don't need to put in a huge amount under "Anything else we haven't thought of yet" in order to get the books to balance.
The full set is too big for me to attempt the kind of manual analysis I did on the partial set. I therefore wrote a set of programs to automate what I tried to do manually - data cleaning consists of a rather arbitrary pruning of big residuals. I have tried to err on the side of leaving most of the data in. The results are therefore affected to an unknown degree by all those big residuals.
I have corrected the data by subtracting off:
A pth order 2-D polynomial fit (p = 4) to the whole data set. This is pretty much what I did manually above and gives similar results.
A qth order 2-D polynomial fit (q = 1, 2) to each image. As will be seen, q = 2 gives a substantial improvement, especially in V-band, so one would like to try a higher order fit. I have not done this yet and it is not obvious that there is enough data to justify estimating all the extra parameters.
The fits are achieved iteratively. At each stage, a correction to the fits is computed from the current residuals. This process appears to be stable and is quite quick - about a minute for 8 iterations on a 1GHz Pentium, rereading 300,000 star measurements each iteration. Less than half that time if reprogrammed to read the star measurements only once.
The success can be judged from the rms scatter in the magnitude measurements.The manual analysis above which is the case p=4, q=0, showed a noise floor above 0.01 magnitudes. The results here, p=4, q=1, are better and the lower plot, p=4, q=2, is better still.
Results for the I-band data are similar. The manual analysis noise floor is around 0.01 magnitudes. The upper plot, p=4, q=1, shows a substantial improvement and the lower plot, p-4, q=2, an additional small improvement.
Notice, however, the large number of stars showing substantially worse scatter: the measured scatter for these hardly changes with the improved correction. The larger scatter, in the vast majority of cases, is not an indication of variability but rather of the large errors at least some of which I have traced to the splitting/joining problem. There are far more high scatter stars in I-band than in V-band. I presume this is because the relative sensitivities and noise levels are such that there are far more sources above the analysis threshold in I-band and hence a greater likelihood of finding detectable stars close together.
I have also calculated Welch-Stetson variability parameters. In order to get reasonable looking results, I found it necessary to boost the error estimates above the ones calculated by my PSF fitting program. These are similar to the error estimates computed by the standard pipeline and both underestimate the errors for bright sources. I therefore added a constant amount to the estimated variances so as to make the general scatter similar for bright and for faint sources. The upper plot is for p=4, q=1 and the lower for p=4, q=2.
The scatter is reduced considerably in the q=2 case. The outlying points, both positive and negative, are largely the result of the large errors that I attribute to the splitting/joining problem. None of the outliers that I have checked look remotely like plausible variable stars. In the second case that I checked, both stars of the split pair were retained in the final collected data set so the cause of the problem was quite clear in that case.
I show TASS V-I plotted against Tycho2 B-V. Note that the zero point of TASS I magnitudes is not at all precisely known so this plot may be slid up or down relative to the vertical axis if you want to do so.
Um. Quite a mess. Dare I suggest that the horizontal spreading of the plot at around TASS V-I = 0.6 indicates that Tycho2 has a problem for faint stars? I append the plots for brighter stars:
Many TASS stars have magnitudes very different from the Tycho2 magnitudes. At least some of these really are long period variables. To try to get a handle on these, I have plotted TASS I - Tycho2 "I" versus TASS V - Tycho2 V. I have put Tycho2 "I" in scare quotes to emphasise that it is not an I-band measurement but rather an extrapolation from what was actually measured.
What I expected to see was a big clump around (0, 0) with outliers either way along a line with slope around 0.8 - the typical ratio of I-band change to V-band change for the types of variable that I have come across so far.
At first sight, the plot appears to show that almost all the stars are variable, typically by say a quarter of a magnitude in V-band and approaching half a magnitude in I-band. A ratio nearer to 2 than to 0.8. I don't believe this! So I then plotted only the bright stars.
Then, for these bright stars, I separated the colour ranges B-V < 0.8, which I called "blue" and 0.8 < B-V < 2.0, which I called "red", and threw away the very red stars which scatter all over the place.
Both "blue" and "red" plots show what I would have thought was rather large scatter, forming an ellipse with major axis along dI/dV = 2. This is very much what one would see if a major source of the scatter were in the Tycho2 V-band measurements. Ouch! I thought Tycho2 was better than that.
In TN #90 Tom Droege points out that his flats suffer from north-south gradients due to sky light that vary with declination/zenith angle. He proposes to generate different flats for different declinations. I don't think this is the right thing to do. I suggest it would be better to use one flat file and subsequently apply empirical spatial corrections within each image to deal with the observed magnitude errors. The larger the amount of data averaged to make the flat file, the more reliable it should be and the less noise it introduces into the final image. So making one flat rather than many should be more reliable and give lower final noise.
My preference would be to use only light box flats. Light box images have better signal to noise ratio than sky flats so a flat file can be generated from a very small number of images. Only one light box image is needed as far as precision is concerned, giving around 0.002 magnitudes which is adequate in view of our bright skies. Three images is the minimum number for median processing to eliminate cosmic rays but this doesn't mean that one needs to take three every night.
In view of the observed large variation in the large scale structures of the various flats and their lack of correlation with the required corrections, it would probably be advisable to remove the large scale structure in the process of forming the flat. This is relatively simple to do, if one is content with subtracting a polynomial, and could be applied equally well to sky flats as to light box flats ... don't use this as an excuse not to take light box flats!
So my proposal is to take a small number of darks and light box flats on a regular basis, sufficient to follow the evolution of the system and in sufficient number to allow median processing to eliminate cosmic rays. This should be done often enough to give sufficiently current dark/flat processing since we know that the pattern does evolve, if only by moving the dust around. I think this should be nightly. Tom, with his experience of actually doing these things, obviously prefers less frequent flats! Certainly, the improvement available from frequent flats is small: as I report below, the increased noise level from the use of a very elderly dark/flat combination was rather small.
Taking regular flats, of course, requires a convenient light box that can be positioned reproducibly.
The work reported above shows that varying spatial correction must be applied on an image-by-image basis. This requires a survey pattern that provides some overlap so that at least the linear gradients can be determined from the stars seen at opposite edges of images. Only a small overlap would be needed for this. Maybe 10% each way; about 20% overall. However, the quadratic terms, which I showed above gave an additional improvement, require much greater overlap. Stepping by half the width of the image in both coordinates provides full coverage for more complicated corrections. I think this would be a good idea on a routine basis. This overlap does not directly reduce the amount of surveying that can be done, provided the overlap pairs are taken over a short enough time to neglect changes. The basic compromise of a large area survey between short repeat time and total area covered is largely unaffected by the requirement to space repeats by half the image width in one or both coordinates.
My latest PSF fitting routines fit "Type A" distributions - that is, Gaussians multiplied by polynomials. I suspect that, like me, Type A distributions are a little old-fashioned. They go back to Gauss. I found them in A. C. Aitken, Statistical Mathematics, Oliver and Boyd, 7th edition, 1952. "It can be shown that" any general distribution may be expressed as a sum of Type A contributions ... with some restrictions, such as the moments existing and being bounded. The extension to two dimensional distributions such as the PSF is fairly obvious. The expansion is, however, of rather restricted use in the real world: in most of the cases I have tried, convergence as one increases the number of terms is very slow indeed and, in spite of the theorem that says one can represent "any" (well - almost any) distribution, divergence sometimes occurs in the sense that the fit gets worse as one increases the number of terms. A major drawback is that the expansion, truncated to a finite number of terms, takes on negative values. Ouch: not nice for non-negative things like PSFs. I did not try this type of approximation for the earlier MK IV images because of this problem: the finite Type A expansion does not do well with the core/halo type of distribution. However, the images are now very much better and my multi-component ellipse model had such serious problems that I have been driven to the Type A.
The major drawback of the multi-component ellipse model was that it gave equally good fits with the ellipses permuted (a sixfold redundancy with three ellipses) or reflected across their axes (eightfold redundancy with three ellipses: a total of 48-fold!). Fits could settle onto any of these equivalent representations after the random search. This made it very difficult to compare the fits for different images. An even more pernicious problem was the possibility that the fit could lock on to different equivalent representations in different parts of the image, leading to "solutions" with completely incorrect spatial variation. Think of it as trying to pot a billiard ball blindfold on a table with 48 identical pockets and a very rough playing surface with a large number of pits (local minima) which look just like the real pockets except for being shallower. Then imagine, if you can, that the playing surface is, say, 52-dimensional.
The current version includes quadratic variation across the image for the 2nd order (Gaussian) terms, linear variation for the 3rd order Type A terms and constant 4th order Type A terms. This requires the fitting of 35 coefficients and gives a moderately good fit to the actual variation of the PSF - not as good as I was able to obtain (some of the time) with my previous multi-component ellipses but a great deal easier to fit. The new program does not appear to suffer from the problem of multiple solutions. Coefficients show a regular progression from image to image, apart from images with clouds.
The actual fitting uses the Simplex or Amoeba method. Each image after the first uses the previous fit as a starting point. The initial step size is chosen to be rather large. Returning to the pitted billiard table analogy, one initially thrashes around violently in the hope that even if one accidentally starts out in one of the false pits, one will succeed, equally accidentally, in finding the pocket ... hopefully now one pocket rather than 48! Large step sizes are retained until the solution ceases to change rapidly. This gives rapid convergence if the starting estimate is good balanced with a high probability of climbing out of a local minimum (false solution).
Data is read in, Dark and Flat correction applied and noisy pixels and identified bad areas of the image are flagged. A star list is created by hunting for local maxima. An initial astrometric fit to Tycho2 is made using my version of Michael Richmond's matching routines. The PSF is fitted to a selection of bright but non-saturated Tycho2 stars distributed over the whole image. All stars are then measured using the final PSF. Final astrometry is calculated using the fitted (X1, X2) coordinates, final photometry is set using the median of the Tycho2 stars and the results are added to the output files.
I used three different pairs of Dark and Flat files. I did my initial analysis using my old Data Set 20 files which included a Flat file derived from Tom's lightbox flats. Then I overcame my laziness and converted Tom's master files to my formats and used those. Yes - I use my own private personal formats. The results are rather different. I also used null flat files in which the factor is set to unity for all data points.
As part of the routine processing, my programs compute an estimated standard deviation for the sky background both before and after Dark/Flat processing. Typical figures in adu for the central image of the Orion set:
Raw Null Tom Box
V 36.33 32.51 29.53 30.35
I 45.04 39.47 35.84 37.65
The noise improves with processing, as would be expected since processing takes
out some of the variation due to differences between pixel sensitivities and
offsets. The improvement is greater using Tom's current files instead of my old
ones. I think this tells us that things do change with time, if only by
adding or loosing dust. This supports the common practice of taking lots of
Dark and Flat images. In case you think that the improvement of one or two
adu is pretty small to be worth bothering about, let me remind you that dropping
the digitizer resolution from 16 bits to 12 bits would have the same effect as
increasing an rms noise of 29.53 adu to 29.89 adu. 16 bits was considered essential
and there were those who wanted more ... going to 18 bits would not quite reduce
an rms noise of 29.530 adu to 29.528 adu. But I digress.
I used weighted least squares fitting, the weightings being determined in the usual way from the estimated measurement standard deviations. The slight inconvenience of weighted fitting is vastly outweighed by the advantage of being able to throw in the entire data set without worrying that the many noisy weak stars will dominate the fit at the expense of the few bright stars. On the other hand, one does have to worry a bit when counting effective degrees of freedom. With an unweighted fit, as long as there are far more observations than data points, one doesn't have to worry but with a weighted fit to an inhomogeneous data set the effective number of degrees of freedom may be very much less than the number of observations. What does all this pedantry have to do with the problem at hand? Well, there's no problem fitting the 15 4th order spatial coefficients to say 25,000 observations but fitting the individual images, with maybe 500 points, most of them faint, is another matter. One should really do the sums. I haven't, yet. Meanwhile, I am stopping at 2nd order (6 parameters) for individual images.
The quantities fitted are the residuals from the previous fit and the fits are iterated until they converge. A residuals is the difference between an observation and the mean of all observations for that star: this does not depend on the Tycho2 magnitudes - Tycho2 is used only to select the stars used in the fit. One would be better off using a catalog generated from the data set as the reference catalog, in view of the uncertainties in Tycho2, especially when extrapolated to I-band. I did not generate such a catalog because I judged Tycho2 to be good enough for the purpose. First, all the residuals are fitted to give the change in the higher order overall spatial correction. Then the individual images are fitted without recomputing for this change. The change in the individual image fits is then the fitted amount less the lower order part of the overall change.
With the corrections initially set to zero, the first set of residuals is substantially different from an ideal set which could be fitted to give the final fit in one go. Hence the need for iteration. Changes in the individual image fits change the overall, higher order spatial correction. Experimentally, the iteration does converge in about 8 cycles. I am not sure if it "always" converges. I think it does.