1 #include <u.h>
2 #include <libc.h>
3 #include <bio.h>
4 #include "sky.h"
5
6 void
amdinv(Header * h,Angle ra,Angle dec,float mag,float col)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
ppoinv(Header * h,Angle ra,Angle dec)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
traneqstd(Header * h,Angle object_ra,Angle object_dec)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
xypos(Header * h,Angle ra,Angle dec,float mag,float col)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 }
248