sys/src/cmd/scat/posn.c

   1 #include        <u.h>
   2 #include        <libc.h>
   3 #include        <bio.h>
   4 #include        "sky.h"
   5
   6 void
   7 amdinv(Header *h, Angle ra, Angle dec, float mag, float col)
   8 {
   9         int i, max_iterations;
  10         float tolerance;
  11         float object_x, object_y, delta_x, delta_y, f, fx, fy, g, gx, gy;
  12         float x1, x2, x3, x4;
  13         float y1, y2, y3, y4;
  14
  15         /*
  16          *  Initialize
  17          */
  18         max_iterations = 50;
  19         tolerance = 0.0000005;
  20
  21         /*
  22          *  Convert RA and Dec to St.coords
  23          */
  24         traneqstd(h, ra, dec);
  25
  26         /*
  27          *  Set initial value for x,y
  28          */
  29         object_x = h->xi/h->param[Ppltscale];
  30         object_y = h->eta/h->param[Ppltscale];
  31
  32         /*
  33          *  Iterate by Newtons method
  34          */
  35         for(i = 0; i < max_iterations; i++) {
  36                 /*
  37                  *  X plate model
  38                  */
  39                 x1 = object_x;
  40                 x2 = x1 * object_x;
  41                 x3 = x1 * object_x;
  42                 x4 = x1 * object_x;
  43
  44                 y1 = object_y;
  45                 y2 = y1 * object_y;
  46                 y3 = y1 * object_y;
  47                 y4 = y1 * object_y;
  48
  49                 f =
  50                         h->param[Pamdx1] * x1 +
  51                         h->param[Pamdx2] * y1 +
  52                         h->param[Pamdx3] +
  53                         h->param[Pamdx4] * x2 +
  54                         h->param[Pamdx5] * x1*y1 +
  55                         h->param[Pamdx6] * y2 +
  56                         h->param[Pamdx7] * (x2+y2) +
  57                         h->param[Pamdx8] * x2*x1 +
  58                         h->param[Pamdx9] * x2*y1 +
  59                         h->param[Pamdx10] * x1*y2 +
  60                         h->param[Pamdx11] * y3 +
  61                         h->param[Pamdx12] * x1* (x2+y2) +
  62                         h->param[Pamdx13] * x1 * (x2+y2) * (x2+y2) +
  63                         h->param[Pamdx14] * mag +
  64                         h->param[Pamdx15] * mag*mag +
  65                         h->param[Pamdx16] * mag*mag*mag +
  66                         h->param[Pamdx17] * mag*x1 +
  67                         h->param[Pamdx18] * mag * (x2+y2) +
  68                         h->param[Pamdx19] * mag*x1 * (x2+y2) +
  69                         h->param[Pamdx20] * col;
  70
  71                 /*
  72                  *  Derivative of X model wrt x
  73                  */
  74                 fx =
  75                         h->param[Pamdx1] +
  76                         h->param[Pamdx4] * 2*x1 +
  77                         h->param[Pamdx5] * y1 +
  78                         h->param[Pamdx7] * 2*x1 +
  79                         h->param[Pamdx8] * 3*x2 +
  80                         h->param[Pamdx9] * 2*x1*y1 +
  81                         h->param[Pamdx10] * y2 +
  82                         h->param[Pamdx12] * (3*x2+y2) +
  83                         h->param[Pamdx13] * (5*x4 + 6*x2*y2 + y4) +
  84                         h->param[Pamdx17] * mag +
  85                         h->param[Pamdx18] * mag*2*x1 +
  86                         h->param[Pamdx19] * mag*(3*x2+y2);
  87
  88                 /*
  89                  *  Derivative of X model wrt y
  90                  */
  91                 fy =
  92                         h->param[Pamdx2] +
  93                         h->param[Pamdx5] * x1 +
  94                         h->param[Pamdx6] * 2*y1 +
  95                         h->param[Pamdx7] * 2*y1 +
  96                         h->param[Pamdx9] * x2 +
  97                         h->param[Pamdx10] * x1*2*y1 +
  98                         h->param[Pamdx11] * 3*y2 +
  99                         h->param[Pamdx12] * 2*x1*y1 +
 100                         h->param[Pamdx13] * 4*x1*y1* (x2+y2) +
 101                         h->param[Pamdx18] * mag*2*y1 +
 102                         h->param[Pamdx19] * mag*2*x1*y1;
 103                 /*
 104                  *  Y plate model
 105                  */
 106                 g =
 107                         h->param[Pamdy1] * y1 +
 108                         h->param[Pamdy2] * x1 +
 109                         h->param[Pamdy3] +
 110                         h->param[Pamdy4] * y2 +
 111                         h->param[Pamdy5] * y1*x1 +
 112                         h->param[Pamdy6] * x2 +
 113                         h->param[Pamdy7] * (x2+y2) +
 114                         h->param[Pamdy8] * y3 +
 115                         h->param[Pamdy9] * y2*x1 +
 116                         h->param[Pamdy10] * y1*x3 +
 117                         h->param[Pamdy11] * x3 +
 118                         h->param[Pamdy12] * y1 * (x2+y2) +
 119                         h->param[Pamdy13] * y1 * (x2+y2) * (x2+y2) +
 120                         h->param[Pamdy14] * mag +
 121                         h->param[Pamdy15] * mag*mag +
 122                         h->param[Pamdy16] * mag*mag*mag +
 123                         h->param[Pamdy17] * mag*y1 +
 124                         h->param[Pamdy18] * mag * (x2+y2) +
 125                         h->param[Pamdy19] * mag*y1 * (x2+y2) +
 126                         h->param[Pamdy20] * col;
 127
 128                 /*
 129                  *  Derivative of Y model wrt x
 130                  */
 131                 gx =
 132                         h->param[Pamdy2] +
 133                         h->param[Pamdy5] * y1 +
 134                         h->param[Pamdy6] * 2*x1 +
 135                         h->param[Pamdy7] * 2*x1 +
 136                         h->param[Pamdy9] * y2 +
 137                         h->param[Pamdy10] * y1*2*x1 +
 138                         h->param[Pamdy11] * 3*x2 +
 139                         h->param[Pamdy12] * 2*x1*y1 +
 140                         h->param[Pamdy13] * 4*x1*y1 * (x2+y2) +
 141                         h->param[Pamdy18] * mag*2*x1 +
 142                         h->param[Pamdy19] * mag*y1*2*x1;
 143
 144                 /*
 145                  *  Derivative of Y model wrt y
 146                  */
 147                 gy =
 148                         h->param[Pamdy1] +
 149                         h->param[Pamdy4] * 2*y1 +
 150                         h->param[Pamdy5] * x1 +
 151                         h->param[Pamdy7] * 2*y1 +
 152                         h->param[Pamdy8] * 3*y2 +
 153                         h->param[Pamdy9] * 2*y1*x1 +
 154                         h->param[Pamdy10] * x2 +
 155                         h->param[Pamdy12] * 3*y2 +
 156                         h->param[Pamdy13] * (5*y4 + 6*x2*y2 + x4) +
 157                         h->param[Pamdy17] * mag +
 158                         h->param[Pamdy18] * mag*2*y1 +
 159                         h->param[Pamdy19] * mag*(x2 + 3*y2);
 160
 161                 f = f - h->xi;
 162                 g = g - h->eta;
 163                 delta_x = (-f*gy+g*fy) / (fx*gy-fy*gx);
 164                 delta_y = (-g*fx+f*gx) / (fx*gy-fy*gx);
 165                 object_x = object_x + delta_x;
 166                 object_y = object_y + delta_y;
 167                 if((fabs(delta_x) < tolerance) && (fabs(delta_y) < tolerance))
 168                         break;
 169         }
 170
 171         /*
 172          *  Convert mm from plate center to pixels
 173          */
 174         h->x = (h->param[Pppo3]-object_x*1000.0)/h->param[Pxpixelsz];
 175         h->y = (h->param[Pppo6]+object_y*1000.0)/h->param[Pypixelsz];
 176 }
 177
 178 void
 179 ppoinv(Header *h, Angle ra, Angle dec)
 180 {
 181
 182         /*
 183          *  Convert RA and Dec to standard coords.
 184          */
 185         traneqstd(h, ra, dec);
 186
 187         /*
 188          *  Convert st.coords from arcsec to radians
 189          */
 190         h->xi  /= ARCSECONDS_PER_RADIAN;
 191         h->eta /= ARCSECONDS_PER_RADIAN;
 192
 193         /*
 194          *  Compute PDS coordinates from solution
 195          */
 196         h->x =
 197                 h->param[Pppo1]*h->xi +
 198                 h->param[Pppo2]*h->eta +
 199                 h->param[Pppo3];
 200         h->y =
 201                 h->param[Pppo4]*h->xi +
 202                 h->param[Pppo5]*h->eta +
 203                 h->param[Pppo6];
 204         /*
 205          * Convert x,y from microns to pixels
 206          */
 207         h->x /= h->param[Pxpixelsz];
 208         h->y /= h->param[Pypixelsz];
 209
 210 }
 211
 212 void
 213 traneqstd(Header *h, Angle object_ra, Angle object_dec)
 214 {
 215         float div;
 216
 217         /*
 218          *  Find divisor
 219          */
 220         div =
 221                 (sin(object_dec)*sin(h->param[Ppltdec]) +
 222                 cos(object_dec)*cos(h->param[Ppltdec]) *
 223                 cos(object_ra - h->param[Ppltra]));
 224
 225         /*
 226          *  Compute standard coords and convert to arcsec
 227          */
 228         h->xi = cos(object_dec) *
 229                 sin(object_ra - h->param[Ppltra]) *
 230                 ARCSECONDS_PER_RADIAN/div;
 231
 232         h->eta = (sin(object_dec)*cos(h->param[Ppltdec])-
 233                 cos(object_dec)*sin(h->param[Ppltdec])*
 234                 cos(object_ra - h->param[Ppltra]))*
 235                 ARCSECONDS_PER_RADIAN/div;
 236
 237 }
 238
 239 void
 240 xypos(Header *h, Angle ra, Angle dec, float mag, float col)
 241 {
 242         if (h->amdflag) {
 243                 amdinv(h, ra, dec, mag, col);
 244         } else {
 245                 ppoinv(h, ra, dec);
 246         }
 247 }