The Estimation of Stellar Image Centers from
                 ▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀
                               Raster Scan Data
                               ▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀

      1. Model of stellar photographic image

         In case that the density of the emulsion does not reaches the
         limit of maximum density the Gaussian two dimensional Normal
         Distribution Function known from error theory is a good
         representation of the photographic density distribution of an
         well exposed stellar image.

         Let's assume an area which contains the image of a star. The
         scanned area might look like as outlined:

              -5  -4  -3  -2  -1   0   1   2   3   4   5
             ┌───┬───┬───┬───┬───┬───┬───┬───┬───┬───┬───┐
           5 │  1│  2│ . │ . │   │   │   │   │ . │ 10│ 11│ 5
             ├───┼───┼───┼───┼───┼───┼───┼───┼───┼───┼───┤
           4 │ 12│ 13│ . │   │   │   │   │   │   │ . │ 22│ 4
             ├───┼───┼───┼───┼───┼───┼───┼───┼───┼───┼───┤
           3 │   │   │   │   │   │   │   │   │   │   │   │ 3
             ├───┼───┼───┼───┼───┼───┼───┼───┼░░─┼───┼───┤
           2 │   │   │   │   │   │   │ ░░░░░░░░░░░   │   │ 2 /│\
             ├───┼───┼───┼───┼───┼──░░░░░░░░░░░░░░░░─┼───┤    │
           1 │   │   │   │   │   │░░░░▒▒▒▒▒▒▒▒▒░░░░░ │   │ 1  │
             ├───┼───┼───┼───┼──░░░░▒▒▒▓▓▓▓▓▓▒▒▒░░░░─┼───┤    │
           0 │   │   │   │   │ ░░░▒▒▒▓▓▓███▓▓▓▒▒░░░  │   │ 0  │ y-axis
             ├───┼───┼───┼───┼░░░▒▒▓▓▓█████▓▓▓▒▒░░░──┼───┤    │
          -1 │   │   │   │   │░░░▒▒▓▓▓███▓▓▓▒▒▒░░░   │   │-1  │
             ├───┼───┼───┼───░░░░▒▒▒▓▓▓▓▓▓▒▒▒░░░░┼───┼───┤    │
          -2 │   │   │   │   ░░░░░▒▒▒▒▒▒▒▒▒░░░░  │   │   │-2  │
             ├───┼───┼───┼───░░░░░░░░░░░░░░░░┼───┼───┼───┤
          -3 │   │   │   │   │ ░░░░░░░░░░░   │   │   │   │-3
             ├───┼───┼───┼───┼───░░──┼───┼───┼───┼───┼───┤
          -4 │   │   │   │   │   │   │   │   │   │   │   │-4
             ├───┼───┼───┼───┼───┼───┼───┼───┼───┼───┼───┤
          -5 │111│ . │   │   │   │   │   │   │   │ . │121│-5
             └───┴───┴───┴───┴───┴───┴───┴───┴───┴───┴───┘
              -5  -4  -3  -2  -1   0   1   2   3   4   5
                        ────────────────────────>
                                x-axis

         With regard to the fog of the platen we may write by
         introducing a Cartesian coordinate system centered at the
         center of the stellar image (a5,a7):


                   x-a5         y-a7
               X = ────  ,  Y = ────  :
                    a6           a8


                                                         ┌                                       ┐
                                                         │   ┌                               ┐a10│
                                  x-a5      y-a7         │ 1 │   1     ┌                    ┐│   │
                 d(x,y) = a1 + a2∙──── + a3∙──── + a4∙exp│-─∙│────────∙│ X² - 2∙a9∙X∙Y + Y² ││   │
                                   a6        a8          │ 2 │ 1 - a9² └                    ┘│   │
                                                         │   └                               ┘   │.
                                                         └                                       ┘

         Here are  a1 : fog of platen
                   a2 : x-slope of fog in units of a6
                   a3 : y-slope of fog in units of a8
                   a4 : center density above fog
                   a5 : center x-coordinate of stellar image
                   a6 : measure of width in x-coordinate
                   a7 : center y-coordinate of stellar image
                   a8 : measure of width in y-coordinate
                   a9 : 'correlation' between x and y, │a9│ < 1
                   a10: accounts for flattened center, a10 > 0.

         a6, a8, and a9 alone will describe any kind of elliptical
         shape.  However, under normal conditions we may expect that:

            1. a6 = a8 and a9 = 0, i.e. a circular stellar image,
            2. a10 ≈ 1, i.e. not strongly saturated stellar image,
            3. a1 ≈ 0.25, i.e. no unusually high fog on plate.

         In case of guiding errors the images may be somewhat elliptic.
         Due to the non linearity of photographic emulsions also a
         strong magnitude dependence may be introduced additionally which
         has to be accounted for.

         If there are bright stars also to be measured the coefficient
         a10 will be significantly greater then 1 (up to 3).  However,
         this has no influence on the accuracy of the stellar position.

         However, in case of linear photosensitive photon recording, i.e.
         CCD's, no flattened stellar image is expected and hence a10 = 1.


      2. Data Reduction by Least Squares

         Let's assume a consequent numbering of the pixels.  With
         pixel size DX,DY and pixel numbers NX,NY in x- and y-direction
         respectively we assume our pixels (xi,yi, i=1 to NX∙NY) to be
         defined with respect to the center coordinates xc,yc of the
         rectangular data area as follows:

               xi = xc - DX∙( ½∙(NX-1) + (i-1) mod NX )

               yi = yc + DY∙( ½∙(NY-1) - (i-1) / NX )   ,

               i = 1 to NX∙NY .

         To solve the equations

               D(xi,yi) ≡ Di = d(xi,yi; a1, ..., a10) + vi ,
                                        _
               D(xi,yi)      = d(xi,yi; a) + vi ,  i = 1 to NX∙NY,

         we first have to linearize d(xi,yi; a1, ..., a10) with respect
         to initial coefficients a01, a02, ..., a010. Therefore we need
         the partial derivations with respect to a1, ..., a10 at the
         point a01, ..., a010. With the abbreviations:

               X(x,â) = (x-a5)/a6 ,  Y(y,â) = (y-a7)/a8 ,


                   ┌    ┐
               _   │ a1 │
               a = │ .. │  ,
                   │ a10│
                   └    ┘

                      1     ┌                    ┐
               A = ────────∙│ X² - 2∙a9∙X∙Y + Y² │   ,
                    1 - a9² └                    ┘

                     _          ┌     a10 ┐
               u(x,y;a) = a4∙exp│ -½∙A    │
                                └         ┘
               we derive:

                dd
               ──── = 1  ,
                da1

                dd
               ──── = X  ,
                da2

                dd
               ──── = Y  ,
                da3

                dd      u
               ──── = ────   ,
                da4    a2

                dd    -1    a10   1      a10-1
               ──── = ── + ────∙─────∙u∙A     ∙(X - a9∙Y)   ,
                da5   a6    a6  1-a9²

                dd       dd
               ──── = X∙────   ,
                da6      da5

                dd    -1    a10   1      a10-1
               ──── = ── + ────∙─────∙u∙A     ∙(Y - a9∙X)   ,
                da7   a8    a8  1-a9²

                dd       dd
               ──── = Y∙────   ,
                da8      da7

                dd          1      a10-1
               ──── = a10∙─────∙u∙A     ∙(X∙Y - a9∙A)   ,
                da9       1-a9²

                dd      1     a10
               ──── = -───∙u∙A   ∙ln(A)   .
               da10     2

      The integrated intensity of this 2-dimensional modified
      Gaußfunction is given by:


              Γ(1/a10)             1/a10        ½
          I = ────────∙π∙a4∙a6∙a8∙2     ∙(1-a9²)
                 a10


                 1/a10
                2                               ½
            = π∙──────∙Γ(1/a10)∙a4∙a6∙a8∙(1-a9²)   .
                 a10

      For estimation of the mean error σI it's necessary to linearize
      I = I(a4,a6,a8,a9,a10) with respect to the parameters a4,.. ,a10.


                   1/a10
          dI      2                               ½    I
          ─── = π∙──────∙Γ(1/a10)∙a4∙a6∙a8∙(1-a9²)  = ────
          da4      a10                                 a4

          dI                                           I
          ─── =                                       ────
          da6                                          a6

          dI                                           I
          ─── =                                       ────
          da8                                          a8

                         -½
          dI    ½∙(1-a9²)  ∙(-2a9)      -a9
          ─── = ──────────────────∙I = ─────∙I
          da9              ½           1-a9²
                    (1-a9²)

                                         ┌               ┐
                                         │ 1/a10         │
           dI                      ½  d  │2              │
          ──── = π∙a4∙a6∙a8∙(1-a9²) ∙────│──────∙Γ(1/a10)│
          da10                       da10│ a10           │
                                         └               ┘

      The term in brackets (u) is:

               1/a10               ┌ 1/a10┐
              2                  ln│2     │-ln(a10)
          u = ──────∙Γ(1/a10) = e  └      ┘        ∙Γ(1/a10)
               a10

                                 1/a10∙ln(2)-ln(a10)
                              = e                   ∙Γ(1/a10)

      It follows for the derivation of u with respect to a10:


           du      ┌ -1           1 ┐      u     dΓ(1/a10)
          ──── = u∙│─── ∙ln(2) - ───│ + ─────── ∙─────────
          da10     └a10²         a10┘   Γ(1/a10)   da10

                     ┌                                  ┐
                  u  │ -ln(2)          a10    dΓ(1/a10) │
               = ───∙│ ────── - 1 +  ────────∙───────── │
                 a10 │  a10          Γ(1/a10)   da10    │ ,
                     └                                  ┘

                      ┌                ,        ┐
           du     -u  │               Γ (1/a10) │
          ──── = ────∙│ ln(2) + a10 + ───────── │
          da10   a10² │               Γ(1/a10)  │  .
                      └                         ┘