| | 1 | (started by dgm 9 Feb 07 after conversations with Gene) |
| | 2 | |
| | 3 | Big Picture: |
| | 4 | |
| | 5 | Due to the effects of seeing (mostly) and other stuff (a little), one get better astrometric accuracy by solving |
| | 6 | for frame-to-frame mapping coefficients and then star parameters (relative position, motion, parallax, etc.) |
| | 7 | in a small chunk of sky than one gets in a big chunk of sky. The division between "small" and "big" comes from |
| | 8 | being able to locally characterize the seeing and other errors by a low order polynomial (or similar) as disctinct |
| | 9 | from actually understanding all of the effects and taking them out on all (or at least larger) scales. You |
| | 10 | get better accuracy, but the price is doing differential measurements. You get relative position, motion, |
| | 11 | parallax, etc., with respect to the mean of all stars and galaxies in the soltuion. For bright stars, this |
| | 12 | is a big problem because the correction from relative to absolute parallax, for example, can be a few milliarcseconds |
| | 13 | at 15th magnitude. |
| | 14 | |
| | 15 | Brief Description of Current Algorithm: |
| | 16 | |
| | 17 | Outer loop: All portions of the visible sky. At the moment, I have adopted the HEALPix tessellation using |
| | 18 | NSide=512 which gives NPixel=3,145,728 and a squarish box of side about 6.9 arcminutes. This may be too |
| | 19 | big, but it is a starting point. |
| | 20 | |
| | 21 | Task: Consult database and get the list of all stars that lie in this HEALPix. Life is a lot easier if |
| | 22 | all stars have unique IDs and even more so if the database has a column containing the HEALPix. |
| | 23 | |
| | 24 | Task: Given the list of star IDs, extract all measures of all stars from all observations. We use the |
| | 25 | chip-based observations (x and y in pixels, brightnes in counts), and we rely on the database |
| | 26 | to provide metadata such exposure data as epoch, bore sight (RA,Dec), camera rotation, filter, etc. |
| | 27 | |
| | 28 | Task: Go to the reference catalog and extract "first guess" data such as position (RA,Dec), magnitude, |
| | 29 | color, star-or-galaxy, parallax, proper motion, etc. |
| | 30 | |
| | 31 | Task: Identify each combination of Observation+CCD that contains enough HEALPix stars, and compute the |
| | 32 | expected position for each reference catalog star. |
| | 33 | |
| | 34 | (At the moment, and possibly for not well understood reasons, I have chosen to project each |
| | 35 | reference star to the tangent plane of an idealized exposure of the nominal field center |
| | 36 | at the instant of the observation based on the SLALIB DS2TP() routine. Since my simulation |
| | 37 | starts with no knowledge of proper motions or parallaxes, the starting guess is Apparent |
| | 38 | places, as distinguished from SLALIB's Catalog and Mean places, at the instant of 2000.0. |
| | 39 | Maybe this needs to be fixed.) |
| | 40 | |
| | 41 | Inner Loop: Since neither the places nor transformations are known, a relaxation is needed. Since most |
| | 42 | do not move and have no parallax, the catalog positions should be reasonably accurate. Hence, the |
| | 43 | first step of the relaxation is to compute the transformations, and the second is to compute the |
| | 44 | mean positions. |
| | 45 | |
| | 46 | Task: Given the set of actual positions (x and y in pixels) and the true positions (xi, eta; starting at |
| | 47 | the catalog values but moving as the iterations proceed), use least squares to compute the |
| | 48 | transformation between (x,y) and (xi,eta). |
| | 49 | |
| | 50 | (At the moment, I use a cubic polynomial in each of X and Y. Many experiments have shown that |
| | 51 | the seeing does all sorts of things, but that a cubic gives a decent approximation without inserting |
| | 52 | steep gradients, etc.) |
| | 53 | |
| | 54 | (Numerical experiments have shown two areas where this solution has problems. The first is where |
| | 55 | the stars used in the fit occupy only a small portion of the HEALPix. This means that the solution |
| | 56 | will be extrapolated for stars in the HEALPix, not included in the soltuion, but whose parameters |
| | 57 | need to be estimated. It is much better to insert these with a low weight than to omit them. |
| | 58 | The other is for the first solution of the first iteration. Noise is bad enough, but stars with |
| | 59 | significant motion or parallax tend to distort the solution. This distortion can persist through |
| | 60 | many iterations. It is much better to restrict this first solution to galaxies or faint stars |
| | 61 | whose motion is zero or negligible.) |
| | 62 | |
| | 63 | (By using (x,y) and (xi,eta), the scale (arcsec/pixel for example) will appear in the coefficients. |
| | 64 | This isn't an extra free parameter because the linear term includes effects from seeing and other |
| | 65 | things, but it makes these terms nowhere near unity. If you convert (xi,eta) into idealized pixels, |
| | 66 | the constant terms can get really big. At some level, this choice is driven by personal preference.) |
| | 67 | |
| | 68 | |
| | 69 | (Remember that this fit is mostly to remove the seeing, residual errors in the telescope model, and |
| | 70 | similar. Since a single Observation+CCD contains many HEALPix, the same pairing will have different |
| | 71 | transformation coefficients in solutions for different HEALPix. One should avoid detailed interpretation |
| | 72 | of these transformation coefficients and attribution of their value to causes other than seeing.) |
| | 73 | |
| | 74 | Task: Given the transforamtion coeffients for each Observation+CCD, transform the measures for each star into |
| | 75 | the mean coordinate frame, and then fit these for position, proper motion, parallax, refraction, |
| | 76 | perturbations from unseen compaions, and anything else you can think of. |
| | 77 | |
| | 78 | (This step needs lots of metadata. A useful choice is to retain times in Julian years counted from |
| | 79 | 2000.0. This makes the constant term the one you need for the catalog and gives you annual proper |
| | 80 | motion. My preference is to use SLALIB to compute the parallax factors and refraction coefficients by |
| | 81 | |
| | 82 | pi = 4.0D00*ATAN(1.0D00) |
| | 83 | halfpi = pi/2.0D00 |
| | 84 | twopi = pi*2.0D00 |
| | 85 | z = SIN(obs_ha)*SIN(halfpi-obs_lat)/SIN(obs_zd) |
| | 86 | obs_pa = ASIN(z) |
| | 87 | z = TAN(obs_zd) |
| | 88 | obs_rra = z*SIN(obs_pa) |
| | 89 | obs_rdec = z*COS(obs_pa) |
| | 90 | CALL sla_MAPQK(r0,d0, ZERO,ZERO,1.0D00,ZERO, |
| | 91 | * amprms, ara,adec) |
| | 92 | CALL sla_MAPQK(r0,d0, ZERO,ZERO,ZERO,ZERO, |
| | 93 | * amprms, bra,bdec) |
| | 94 | obs_pfx = (ara-bra)*3600.0D00*radian*COS(d0) |
| | 95 | obs_pfy = (adec-bdec)*3600.0D00*radian |
| | 96 | |
| | 97 | Sorry for FORTRAN but the C would be very similar. These follow the assumption of the idealized projection |
| | 98 | to a tangent plane axes oriented along RA and Dec. The parallax factors are in units of arcseconds, and the |
| | 99 | refraction vectors are dimensionless. The least squares fit for each star is |
| | 100 | |
| | 101 | X = a + b*t + c*obs_pfx + d*obs_rra |
| | 102 | Y = A + B*t + C*obs_pfy + D*obs_rdec |
| | 103 | |
| | 104 | where X and Y are computed for the star from the appropriate transformation, and obs_pfx, obs_pfy, obs_rra, |
| | 105 | and obs_rry are taken from each observation's metadata.) |
| | 106 | |
| | 107 | (One can argue that there is only one parallax for the star, and that the soltuions for X and Y should be |
| | 108 | coupled into a single solution. I really don't like to do this, but will omit all the reasons here.) |
| | 109 | |
| | 110 | Task: Using the solution coefficients, the predicted place for each star on each CCD of each exposure is updated. |
| | 111 | |
| | 112 | End of Inner Loop: Three iterations is usually enough for convergence. |
| | 113 | |
| | 114 | Task: If you have a list of known galaxies, then their mean motion and parallax should be computed and removed |
| | 115 | from the soltuion. |
| | 116 | |
| | 117 | Task: Update the database with the values for proper motion and parallax. It would be useful to save the |
| | 118 | coefficients of refraction for the first guess when you do this again after new observations become available. |
| | 119 | You probably shouldn't use the constants to update the RA and Dec of your reference star catalog. Too many |
| | 120 | other things use this catalog, and it shouldn't be changed unless there is a really good reason to do so. |
| | 121 | |
| | 122 | End of Outer Loop: The sky is done. Start writing letters to the ApJ or Nature. |
| | 123 | |
| | 124 | Comments, brickbats, and flames. |
| | 125 | |
| | 126 | (dgm) This is a big and messy operation. Simulations based on USNO-B and the LSST focal plane showed that it would take |
| | 127 | several months on a single machine, and was limited by disk access. Further experiments highlighted the need to |
| | 128 | balance disk I/O with the limitations of physical memory. For LSST, the simulation was coded using three different |
| | 129 | disk storage models. The first was one file per observation, the second was one file per raft (9 CCDs), and the |
| | 130 | third was one file per CCD. Just as Goldilocks found out, the middle was pretty good. Saving by observation put |
| | 131 | too many measures into memory, and virtual memory could only hold a year or so. Saving by CCD was amazingly slow |
| | 132 | because of the need to load (and reload) so many little chunks for each HEALPix. The raft-based storage needed |
| | 133 | a few hundred megabytes of memory, and this seems pretty reasonable. |
| | 134 | |
| | 135 | (dgm) The HEALPix tessellation isn't great, but it seems useful. The C language library gives a subroutine that returns |
| | 136 | a unique HEALPix ID for any (RA,Dec), and adding this to the structure for each star assures that stars do not |
| | 137 | move between tiles. It also allows indexing so that the extraction goes a lot faster. The logic for boxes of |
| | 138 | (RA,Dec) gets messy, particularly near the poles, and it really need only be done once per star. |
