posn.c (5532B)
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 }