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 }