1 8a3b2ceb 2004-04-24 devnull #include <u.h>
2 8a3b2ceb 2004-04-24 devnull #include <libc.h>
3 8a3b2ceb 2004-04-24 devnull #include <bio.h>
4 8a3b2ceb 2004-04-24 devnull #include "sky.h"
7 8a3b2ceb 2004-04-24 devnull amdinv(Header *h, Angle ra, Angle dec, float mag, float col)
9 8a3b2ceb 2004-04-24 devnull int i, max_iterations;
10 8a3b2ceb 2004-04-24 devnull float tolerance;
11 8a3b2ceb 2004-04-24 devnull float object_x, object_y, delta_x, delta_y, f, fx, fy, g, gx, gy;
12 8a3b2ceb 2004-04-24 devnull float x1, x2, x3, x4;
13 8a3b2ceb 2004-04-24 devnull float y1, y2, y3, y4;
16 8a3b2ceb 2004-04-24 devnull * Initialize
18 8a3b2ceb 2004-04-24 devnull max_iterations = 50;
19 8a3b2ceb 2004-04-24 devnull tolerance = 0.0000005;
22 8a3b2ceb 2004-04-24 devnull * Convert RA and Dec to St.coords
24 8a3b2ceb 2004-04-24 devnull traneqstd(h, ra, dec);
27 8a3b2ceb 2004-04-24 devnull * Set initial value for x,y
29 8a3b2ceb 2004-04-24 devnull object_x = h->xi/h->param[Ppltscale];
30 8a3b2ceb 2004-04-24 devnull object_y = h->eta/h->param[Ppltscale];
33 8a3b2ceb 2004-04-24 devnull * Iterate by Newtons method
35 8a3b2ceb 2004-04-24 devnull for(i = 0; i < max_iterations; i++) {
37 8a3b2ceb 2004-04-24 devnull * X plate model
39 8a3b2ceb 2004-04-24 devnull x1 = object_x;
40 8a3b2ceb 2004-04-24 devnull x2 = x1 * object_x;
41 8a3b2ceb 2004-04-24 devnull x3 = x1 * object_x;
42 8a3b2ceb 2004-04-24 devnull x4 = x1 * object_x;
44 8a3b2ceb 2004-04-24 devnull y1 = object_y;
45 8a3b2ceb 2004-04-24 devnull y2 = y1 * object_y;
46 8a3b2ceb 2004-04-24 devnull y3 = y1 * object_y;
47 8a3b2ceb 2004-04-24 devnull y4 = y1 * object_y;
50 8a3b2ceb 2004-04-24 devnull h->param[Pamdx1] * x1 +
51 8a3b2ceb 2004-04-24 devnull h->param[Pamdx2] * y1 +
52 8a3b2ceb 2004-04-24 devnull h->param[Pamdx3] +
53 8a3b2ceb 2004-04-24 devnull h->param[Pamdx4] * x2 +
54 8a3b2ceb 2004-04-24 devnull h->param[Pamdx5] * x1*y1 +
55 8a3b2ceb 2004-04-24 devnull h->param[Pamdx6] * y2 +
56 8a3b2ceb 2004-04-24 devnull h->param[Pamdx7] * (x2+y2) +
57 8a3b2ceb 2004-04-24 devnull h->param[Pamdx8] * x2*x1 +
58 8a3b2ceb 2004-04-24 devnull h->param[Pamdx9] * x2*y1 +
59 8a3b2ceb 2004-04-24 devnull h->param[Pamdx10] * x1*y2 +
60 8a3b2ceb 2004-04-24 devnull h->param[Pamdx11] * y3 +
61 8a3b2ceb 2004-04-24 devnull h->param[Pamdx12] * x1* (x2+y2) +
62 8a3b2ceb 2004-04-24 devnull h->param[Pamdx13] * x1 * (x2+y2) * (x2+y2) +
63 8a3b2ceb 2004-04-24 devnull h->param[Pamdx14] * mag +
64 8a3b2ceb 2004-04-24 devnull h->param[Pamdx15] * mag*mag +
65 8a3b2ceb 2004-04-24 devnull h->param[Pamdx16] * mag*mag*mag +
66 8a3b2ceb 2004-04-24 devnull h->param[Pamdx17] * mag*x1 +
67 8a3b2ceb 2004-04-24 devnull h->param[Pamdx18] * mag * (x2+y2) +
68 8a3b2ceb 2004-04-24 devnull h->param[Pamdx19] * mag*x1 * (x2+y2) +
69 8a3b2ceb 2004-04-24 devnull h->param[Pamdx20] * col;
72 8a3b2ceb 2004-04-24 devnull * Derivative of X model wrt x
75 8a3b2ceb 2004-04-24 devnull h->param[Pamdx1] +
76 8a3b2ceb 2004-04-24 devnull h->param[Pamdx4] * 2*x1 +
77 8a3b2ceb 2004-04-24 devnull h->param[Pamdx5] * y1 +
78 8a3b2ceb 2004-04-24 devnull h->param[Pamdx7] * 2*x1 +
79 8a3b2ceb 2004-04-24 devnull h->param[Pamdx8] * 3*x2 +
80 8a3b2ceb 2004-04-24 devnull h->param[Pamdx9] * 2*x1*y1 +
81 8a3b2ceb 2004-04-24 devnull h->param[Pamdx10] * y2 +
82 8a3b2ceb 2004-04-24 devnull h->param[Pamdx12] * (3*x2+y2) +
83 8a3b2ceb 2004-04-24 devnull h->param[Pamdx13] * (5*x4 + 6*x2*y2 + y4) +
84 8a3b2ceb 2004-04-24 devnull h->param[Pamdx17] * mag +
85 8a3b2ceb 2004-04-24 devnull h->param[Pamdx18] * mag*2*x1 +
86 8a3b2ceb 2004-04-24 devnull h->param[Pamdx19] * mag*(3*x2+y2);
89 8a3b2ceb 2004-04-24 devnull * Derivative of X model wrt y
92 8a3b2ceb 2004-04-24 devnull h->param[Pamdx2] +
93 8a3b2ceb 2004-04-24 devnull h->param[Pamdx5] * x1 +
94 8a3b2ceb 2004-04-24 devnull h->param[Pamdx6] * 2*y1 +
95 8a3b2ceb 2004-04-24 devnull h->param[Pamdx7] * 2*y1 +
96 8a3b2ceb 2004-04-24 devnull h->param[Pamdx9] * x2 +
97 8a3b2ceb 2004-04-24 devnull h->param[Pamdx10] * x1*2*y1 +
98 8a3b2ceb 2004-04-24 devnull h->param[Pamdx11] * 3*y2 +
99 8a3b2ceb 2004-04-24 devnull h->param[Pamdx12] * 2*x1*y1 +
100 8a3b2ceb 2004-04-24 devnull h->param[Pamdx13] * 4*x1*y1* (x2+y2) +
101 8a3b2ceb 2004-04-24 devnull h->param[Pamdx18] * mag*2*y1 +
102 8a3b2ceb 2004-04-24 devnull h->param[Pamdx19] * mag*2*x1*y1;
104 8a3b2ceb 2004-04-24 devnull * Y plate model
107 8a3b2ceb 2004-04-24 devnull h->param[Pamdy1] * y1 +
108 8a3b2ceb 2004-04-24 devnull h->param[Pamdy2] * x1 +
109 8a3b2ceb 2004-04-24 devnull h->param[Pamdy3] +
110 8a3b2ceb 2004-04-24 devnull h->param[Pamdy4] * y2 +
111 8a3b2ceb 2004-04-24 devnull h->param[Pamdy5] * y1*x1 +
112 8a3b2ceb 2004-04-24 devnull h->param[Pamdy6] * x2 +
113 8a3b2ceb 2004-04-24 devnull h->param[Pamdy7] * (x2+y2) +
114 8a3b2ceb 2004-04-24 devnull h->param[Pamdy8] * y3 +
115 8a3b2ceb 2004-04-24 devnull h->param[Pamdy9] * y2*x1 +
116 8a3b2ceb 2004-04-24 devnull h->param[Pamdy10] * y1*x3 +
117 8a3b2ceb 2004-04-24 devnull h->param[Pamdy11] * x3 +
118 8a3b2ceb 2004-04-24 devnull h->param[Pamdy12] * y1 * (x2+y2) +
119 8a3b2ceb 2004-04-24 devnull h->param[Pamdy13] * y1 * (x2+y2) * (x2+y2) +
120 8a3b2ceb 2004-04-24 devnull h->param[Pamdy14] * mag +
121 8a3b2ceb 2004-04-24 devnull h->param[Pamdy15] * mag*mag +
122 8a3b2ceb 2004-04-24 devnull h->param[Pamdy16] * mag*mag*mag +
123 8a3b2ceb 2004-04-24 devnull h->param[Pamdy17] * mag*y1 +
124 8a3b2ceb 2004-04-24 devnull h->param[Pamdy18] * mag * (x2+y2) +
125 8a3b2ceb 2004-04-24 devnull h->param[Pamdy19] * mag*y1 * (x2+y2) +
126 8a3b2ceb 2004-04-24 devnull h->param[Pamdy20] * col;
129 8a3b2ceb 2004-04-24 devnull * Derivative of Y model wrt x
132 8a3b2ceb 2004-04-24 devnull h->param[Pamdy2] +
133 8a3b2ceb 2004-04-24 devnull h->param[Pamdy5] * y1 +
134 8a3b2ceb 2004-04-24 devnull h->param[Pamdy6] * 2*x1 +
135 8a3b2ceb 2004-04-24 devnull h->param[Pamdy7] * 2*x1 +
136 8a3b2ceb 2004-04-24 devnull h->param[Pamdy9] * y2 +
137 8a3b2ceb 2004-04-24 devnull h->param[Pamdy10] * y1*2*x1 +
138 8a3b2ceb 2004-04-24 devnull h->param[Pamdy11] * 3*x2 +
139 8a3b2ceb 2004-04-24 devnull h->param[Pamdy12] * 2*x1*y1 +
140 8a3b2ceb 2004-04-24 devnull h->param[Pamdy13] * 4*x1*y1 * (x2+y2) +
141 8a3b2ceb 2004-04-24 devnull h->param[Pamdy18] * mag*2*x1 +
142 8a3b2ceb 2004-04-24 devnull h->param[Pamdy19] * mag*y1*2*x1;
145 8a3b2ceb 2004-04-24 devnull * Derivative of Y model wrt y
148 8a3b2ceb 2004-04-24 devnull h->param[Pamdy1] +
149 8a3b2ceb 2004-04-24 devnull h->param[Pamdy4] * 2*y1 +
150 8a3b2ceb 2004-04-24 devnull h->param[Pamdy5] * x1 +
151 8a3b2ceb 2004-04-24 devnull h->param[Pamdy7] * 2*y1 +
152 8a3b2ceb 2004-04-24 devnull h->param[Pamdy8] * 3*y2 +
153 8a3b2ceb 2004-04-24 devnull h->param[Pamdy9] * 2*y1*x1 +
154 8a3b2ceb 2004-04-24 devnull h->param[Pamdy10] * x2 +
155 8a3b2ceb 2004-04-24 devnull h->param[Pamdy12] * 3*y2 +
156 8a3b2ceb 2004-04-24 devnull h->param[Pamdy13] * (5*y4 + 6*x2*y2 + x4) +
157 8a3b2ceb 2004-04-24 devnull h->param[Pamdy17] * mag +
158 8a3b2ceb 2004-04-24 devnull h->param[Pamdy18] * mag*2*y1 +
159 8a3b2ceb 2004-04-24 devnull h->param[Pamdy19] * mag*(x2 + 3*y2);
161 8a3b2ceb 2004-04-24 devnull f = f - h->xi;
162 8a3b2ceb 2004-04-24 devnull g = g - h->eta;
163 8a3b2ceb 2004-04-24 devnull delta_x = (-f*gy+g*fy) / (fx*gy-fy*gx);
164 8a3b2ceb 2004-04-24 devnull delta_y = (-g*fx+f*gx) / (fx*gy-fy*gx);
165 8a3b2ceb 2004-04-24 devnull object_x = object_x + delta_x;
166 8a3b2ceb 2004-04-24 devnull object_y = object_y + delta_y;
167 8a3b2ceb 2004-04-24 devnull if((fabs(delta_x) < tolerance) && (fabs(delta_y) < tolerance))
172 8a3b2ceb 2004-04-24 devnull * Convert mm from plate center to pixels
174 8a3b2ceb 2004-04-24 devnull h->x = (h->param[Pppo3]-object_x*1000.0)/h->param[Pxpixelsz];
175 8a3b2ceb 2004-04-24 devnull h->y = (h->param[Pppo6]+object_y*1000.0)/h->param[Pypixelsz];
179 8a3b2ceb 2004-04-24 devnull ppoinv(Header *h, Angle ra, Angle dec)
183 8a3b2ceb 2004-04-24 devnull * Convert RA and Dec to standard coords.
185 8a3b2ceb 2004-04-24 devnull traneqstd(h, ra, dec);
188 8a3b2ceb 2004-04-24 devnull * Convert st.coords from arcsec to radians
190 8a3b2ceb 2004-04-24 devnull h->xi /= ARCSECONDS_PER_RADIAN;
191 8a3b2ceb 2004-04-24 devnull h->eta /= ARCSECONDS_PER_RADIAN;
194 8a3b2ceb 2004-04-24 devnull * Compute PDS coordinates from solution
197 8a3b2ceb 2004-04-24 devnull h->param[Pppo1]*h->xi +
198 8a3b2ceb 2004-04-24 devnull h->param[Pppo2]*h->eta +
199 8a3b2ceb 2004-04-24 devnull h->param[Pppo3];
201 8a3b2ceb 2004-04-24 devnull h->param[Pppo4]*h->xi +
202 8a3b2ceb 2004-04-24 devnull h->param[Pppo5]*h->eta +
203 8a3b2ceb 2004-04-24 devnull h->param[Pppo6];
205 8a3b2ceb 2004-04-24 devnull * Convert x,y from microns to pixels
207 8a3b2ceb 2004-04-24 devnull h->x /= h->param[Pxpixelsz];
208 8a3b2ceb 2004-04-24 devnull h->y /= h->param[Pypixelsz];
213 8a3b2ceb 2004-04-24 devnull traneqstd(Header *h, Angle object_ra, Angle object_dec)
215 8a3b2ceb 2004-04-24 devnull float div;
218 8a3b2ceb 2004-04-24 devnull * Find divisor
221 8a3b2ceb 2004-04-24 devnull (sin(object_dec)*sin(h->param[Ppltdec]) +
222 8a3b2ceb 2004-04-24 devnull cos(object_dec)*cos(h->param[Ppltdec]) *
223 8a3b2ceb 2004-04-24 devnull cos(object_ra - h->param[Ppltra]));
226 8a3b2ceb 2004-04-24 devnull * Compute standard coords and convert to arcsec
228 8a3b2ceb 2004-04-24 devnull h->xi = cos(object_dec) *
229 8a3b2ceb 2004-04-24 devnull sin(object_ra - h->param[Ppltra]) *
230 8a3b2ceb 2004-04-24 devnull ARCSECONDS_PER_RADIAN/div;
232 8a3b2ceb 2004-04-24 devnull h->eta = (sin(object_dec)*cos(h->param[Ppltdec])-
233 8a3b2ceb 2004-04-24 devnull cos(object_dec)*sin(h->param[Ppltdec])*
234 8a3b2ceb 2004-04-24 devnull cos(object_ra - h->param[Ppltra]))*
235 8a3b2ceb 2004-04-24 devnull ARCSECONDS_PER_RADIAN/div;
240 8a3b2ceb 2004-04-24 devnull xypos(Header *h, Angle ra, Angle dec, float mag, float col)
242 8a3b2ceb 2004-04-24 devnull if (h->amdflag) {
243 8a3b2ceb 2004-04-24 devnull amdinv(h, ra, dec, mag, col);
244 8a3b2ceb 2004-04-24 devnull } else {
245 8a3b2ceb 2004-04-24 devnull ppoinv(h, ra, dec);
