1 /*** File saoimage/wcslib/platepos.c
2 *** February 29, 2000
3 *** By Jessica Mink, jmink@cfa.harvard.edu
4 *** Harvard-Smithsonian Center for Astrophysics
5 *** Copyright (C) 1998-2002
6 *** Smithsonian Astrophysical Observatory, Cambridge, MA, USA
7
8 This library is free software; you can redistribute it and/or
9 modify it under the terms of the GNU Lesser General Public
10 License as published by the Free Software Foundation; either
11 version 2 of the License, or (at your option) any later version.
12
13 This library is distributed in the hope that it will be useful,
14 but WITHOUT ANY WARRANTY; without even the implied warranty of
15 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
16 Lesser General Public License for more details.
17
18 You should have received a copy of the GNU Lesser General Public
19 License along with this library; if not, write to the Free Software
20 Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
21
22 Correspondence concerning WCSTools should be addressed as follows:
23 Internet email: jmink@cfa.harvard.edu
24 Postal address: Jessica Mink
25 Smithsonian Astrophysical Observatory
26 60 Garden St.
27 Cambridge, MA 02138 USA
28
29 * Module: platepos.c (Plate solution WCS conversion
30 * Purpose: Compute WCS from plate fit
31 * Subroutine: platepos() converts from pixel location to RA,Dec
32 * Subroutine: platepix() converts from RA,Dec to pixel location
33
34 These functions are based on the astrmcal.c portion of GETIMAGE by
35 J. Doggett and the documentation distributed with the Digital Sky Survey.
36
37 */
38
39 #include <math.h>
40 #include <string.h>
41 #include <stdio.h>
42 #include "wcs.h"
43
44 int
platepos(xpix,ypix,wcs,xpos,ypos)45 platepos (xpix, ypix, wcs, xpos, ypos)
46
47 /* Routine to determine accurate position for pixel coordinates */
48 /* returns 0 if successful otherwise 1 = angle too large for projection; */
49 /* based on amdpos() from getimage */
50
51 /* Input: */
52 double xpix; /* x pixel number (RA or long without rotation) */
53 double ypix; /* y pixel number (dec or lat without rotation) */
54 struct WorldCoor *wcs; /* WCS parameter structure */
55
56 /* Output: */
57 double *xpos; /* Right ascension or longitude in degrees */
58 double *ypos; /* Declination or latitude in degrees */
59
60 {
61 double x, y, x2, y2, x3, y3, r2;
62 double xi, xir, eta, etar, raoff, ra, dec, ra0, dec0;
63 double twopi = 6.28318530717959;
64 double ctan, ccos;
65 int ncoeff1 = wcs->ncoeff1;
66 int ncoeff2 = wcs->ncoeff2;
67
68 /* Ignore magnitude and color terms
69 double mag = 0.0;
70 double color = 0.0; */
71
72 /* Convert from pixels to millimeters */
73 x = xpix - wcs->crpix[0];
74 y = ypix - wcs->crpix[1];
75 x2 = x * x;
76 y2 = y * y;
77 x3 = x * x2;
78 y3 = y * y2;
79 r2 = x2 + y2;
80
81 /* Compute xi,eta coordinates in degrees from x,y and plate model */
82 xi = wcs->x_coeff[ 0] + wcs->x_coeff[ 1]*x +
83 wcs->x_coeff[ 2]*y + wcs->x_coeff[ 3]*x2 +
84 wcs->x_coeff[ 4]*y2 + wcs->x_coeff[ 5]*x*y;
85
86 if (ncoeff1 > 6)
87 xi = xi + wcs->x_coeff[ 6]*x3 + wcs->x_coeff[ 7]*y3;
88
89 if (ncoeff1 > 8) {
90 xi = xi + wcs->x_coeff[ 8]*x2*y + wcs->x_coeff[ 9]*x*y2 +
91 wcs->x_coeff[10]*(r2) + wcs->x_coeff[11]*x*r2 +
92 wcs->x_coeff[12]*y*r2;
93 }
94
95 eta = wcs->y_coeff[ 0] + wcs->y_coeff[ 1]*x +
96 wcs->y_coeff[ 2]*y + wcs->y_coeff[ 3]*x2 +
97 wcs->y_coeff[ 4]*y2 + wcs->y_coeff[ 5]*x*y;
98
99 if (ncoeff2 > 6)
100 eta = eta + wcs->y_coeff[ 6]*x3 + wcs->y_coeff[ 7]*y3;
101
102 if (ncoeff2 > 8) {
103 eta = eta + wcs->y_coeff[ 8]*x2*y + wcs->y_coeff[ 9]*y2*x +
104 wcs->y_coeff[10]*r2 + wcs->y_coeff[11]*x*r2 +
105 wcs->y_coeff[12]*y*r2;
106 }
107
108 /* Convert to radians */
109 xir = degrad (xi);
110 etar = degrad (eta);
111
112 /* Convert to RA and Dec */
113 ra0 = degrad (wcs->crval[0]);
114 dec0 = degrad (wcs->crval[1]);
115 ctan = tan (dec0);
116 ccos = cos (dec0);
117 raoff = atan2 (xir / ccos, 1.0 - etar * ctan);
118 ra = raoff + ra0;
119 if (ra < 0.0) ra = ra + twopi;
120 *xpos = raddeg (ra);
121
122 dec = atan (cos (raoff) / ((1.0 - (etar * ctan)) / (etar + ctan)));
123 *ypos = raddeg (dec);
124 return 0;
125 }
126
127
128 int
platepix(xpos,ypos,wcs,xpix,ypix)129 platepix (xpos, ypos, wcs, xpix, ypix)
130
131 /* Routine to determine pixel coordinates for sky position */
132 /* returns 0 if successful otherwise 1 = angle too large for projection; */
133 /* based on amdinv() from getimage */
134
135 /* Input: */
136 double xpos; /* Right ascension or longitude in degrees */
137 double ypos; /* Declination or latitude in degrees */
138 struct WorldCoor *wcs; /* WCS parameter structure */
139
140 /* Output: */
141 double *xpix; /* x pixel number (RA or long without rotation) */
142 double *ypix; /* y pixel number (dec or lat without rotation) */
143
144 {
145 double xi,eta,x,y,xy,x2,y2,x2y,y2x,x3,y3,r2,dx,dy;
146 double tdec,ctan,ccos,traoff, craoff, etar, xir;
147 double f,fx,fy,g,gx,gy;
148 double ra0, dec0, ra, dec;
149 double tolerance = 0.0000005;
150 int max_iterations = 50;
151 int i;
152 int ncoeff1 = wcs->ncoeff1;
153 int ncoeff2 = wcs->ncoeff2;
154
155 /* Convert RA and Dec in radians to standard coordinates on a plate */
156 ra = degrad (xpos);
157 dec = degrad (ypos);
158 tdec = tan (dec);
159 ra0 = degrad (wcs->crval[0]);
160 dec0 = degrad (wcs->crval[1]);
161 ctan = tan (dec0);
162 ccos = cos (dec0);
163 traoff = tan (ra - ra0);
164 craoff = cos (ra - ra0);
165 etar = (1.0 - ctan * craoff / tdec) / (ctan + (craoff / tdec));
166 xir = traoff * ccos * (1.0 - (etar * ctan));
167 xi = raddeg (xir);
168 eta = raddeg (etar);
169
170 /* Set initial value for x,y */
171 x = xi * wcs->dc[0] + eta * wcs->dc[1];
172 y = xi * wcs->dc[2] + eta * wcs->dc[3];
173
174 /* if (wcs->x_coeff[1] == 0.0)
175 x = xi - wcs->x_coeff[0];
176 else
177 x = (xi - wcs->x_coeff[0]) / wcs->x_coeff[1];
178 if (wcs->y_coeff[2] == 0.0)
179 y = eta - wcs->y_coeff[0];
180 else
181 y = (eta - wcs->y_coeff[0]) / wcs->y_coeff[2]; */
182
183 /* Iterate by Newton's method */
184 for (i = 0; i < max_iterations; i++) {
185
186 /* X plate model */
187 xy = x * y;
188 x2 = x * x;
189 y2 = y * y;
190 x3 = x2 * x;
191 y3 = y2 * y;
192 x2y = x2 * y;
193 y2x = y2 * x;
194 r2 = x2 + y2;
195
196 f = wcs->x_coeff[0] + wcs->x_coeff[1]*x +
197 wcs->x_coeff[2]*y + wcs->x_coeff[3]*x2 +
198 wcs->x_coeff[4]*y2 + wcs->x_coeff[5]*xy;
199
200 /* Derivative of X model wrt x */
201 fx = wcs->x_coeff[1] + wcs->x_coeff[3]*2.0*x +
202 wcs->x_coeff[5]*y;
203
204 /* Derivative of X model wrt y */
205 fy = wcs->x_coeff[2] + wcs->x_coeff[4]*2.0*y +
206 wcs->x_coeff[5]*x;
207
208 if (ncoeff1 > 6) {
209 f = f + wcs->x_coeff[6]*x3 + wcs->x_coeff[7]*y3;
210 fx = fx + wcs->x_coeff[6]*3.0*x2;
211 fy = fy + wcs->x_coeff[7]*3.0*y2;
212 }
213
214 if (ncoeff1 > 8) {
215 f = f +
216 wcs->x_coeff[8]*x2y + wcs->x_coeff[9]*y2x +
217 wcs->x_coeff[10]*r2 + wcs->x_coeff[11]*x*r2 +
218 wcs->x_coeff[12]*y*r2;
219
220 fx = fx + wcs->x_coeff[8]*2.0*xy +
221 wcs->x_coeff[9]*y2 +
222 wcs->x_coeff[10]*2.0*x +
223 wcs->x_coeff[11]*(3.0*x2+y2) +
224 wcs->x_coeff[12]*2.0*xy;
225
226 fy = fy + wcs->x_coeff[8]*x2 +
227 wcs->x_coeff[9]*2.0*xy +
228 wcs->x_coeff[10]*2.0*y +
229 wcs->x_coeff[11]*2.0*xy +
230 wcs->x_coeff[12]*(3.0*y2+x2);
231 }
232
233 /* Y plate model */
234 g = wcs->y_coeff[0] + wcs->y_coeff[1]*x +
235 wcs->y_coeff[2]*y + wcs->y_coeff[3]*x2 +
236 wcs->y_coeff[4]*y2 + wcs->y_coeff[5]*xy;
237
238 /* Derivative of Y model wrt x */
239 gx = wcs->y_coeff[1] + wcs->y_coeff[3]*2.0*x +
240 wcs->y_coeff[5]*y;
241
242 /* Derivative of Y model wrt y */
243 gy = wcs->y_coeff[2] + wcs->y_coeff[4]*2.0*y +
244 wcs->y_coeff[5]*x;
245
246 if (ncoeff2 > 6) {
247 g = g + wcs->y_coeff[6]*x3 + wcs->y_coeff[7]*y3;
248 gx = gx + wcs->y_coeff[6]*3.0*x2;
249 gy = gy + wcs->y_coeff[7]*3.0*y2;
250 }
251
252 if (ncoeff2 > 8) {
253 g = g +
254 wcs->y_coeff[8]*x2y + wcs->y_coeff[9]*y2x +
255 wcs->y_coeff[10]*r2 + wcs->y_coeff[11]*x*r2 +
256 wcs->y_coeff[12]*y*r2;
257
258 gx = gx + wcs->y_coeff[8]*2.0*xy +
259 wcs->y_coeff[9]*y2 +
260 wcs->y_coeff[10]*2.0*x +
261 wcs->y_coeff[11]*(3.0*x2+y2) +
262 wcs->y_coeff[12]*2.0*xy;
263
264 gy = gy + wcs->y_coeff[8]*x2 +
265 wcs->y_coeff[9]*2.0*xy +
266 wcs->y_coeff[10]*2.0*y +
267 wcs->y_coeff[11]*2.0*xy +
268 wcs->y_coeff[12]*(3.0*y2+x2);
269 }
270
271 f = f - xi;
272 g = g - eta;
273 dx = ((-f * gy) + (g * fy)) / ((fx * gy) - (fy * gx));
274 dy = ((-g * fx) + (f * gx)) / ((fx * gy) - (fy * gx));
275 x = x + dx;
276 y = y + dy;
277 if ((fabs(dx) < tolerance) && (fabs(dy) < tolerance)) break;
278 }
279
280 /* Convert from plate pixels to image pixels */
281 *xpix = x + wcs->crpix[0];
282 *ypix = y + wcs->crpix[1];
283
284 /* If position is off of the image, return offscale code */
285 if (*xpix < 0.5 || *xpix > wcs->nxpix+0.5)
286 return -1;
287 if (*ypix < 0.5 || *ypix > wcs->nypix+0.5)
288 return -1;
289
290 return 0;
291 }
292
293
294 /* Set plate fit coefficients in structure from arguments */
295 int
SetPlate(wcs,ncoeff1,ncoeff2,coeff)296 SetPlate (wcs, ncoeff1, ncoeff2, coeff)
297
298 struct WorldCoor *wcs; /* World coordinate system structure */
299 int ncoeff1; /* Number of coefficients for x */
300 int ncoeff2; /* Number of coefficients for y */
301 double *coeff; /* Plate fit coefficients */
302
303 {
304 int i;
305
306 if (nowcs (wcs) || (ncoeff1 < 1 && ncoeff2 < 1))
307 return 1;
308
309 wcs->ncoeff1 = ncoeff1;
310 wcs->ncoeff2 = ncoeff2;
311 wcs->prjcode = WCS_PLT;
312
313 for (i = 0; i < 20; i++) {
314 if (i < ncoeff1)
315 wcs->x_coeff[i] = coeff[i];
316 else
317 wcs->x_coeff[i] = 0.0;
318 }
319
320 for (i = 0; i < 20; i++) {
321 if (i < ncoeff2)
322 wcs->y_coeff[i] = coeff[ncoeff1+i];
323 else
324 wcs->y_coeff[i] = 0.0;
325 }
326 return 0;
327 }
328
329
330 /* Return plate fit coefficients from structure in arguments */
331 int
GetPlate(wcs,ncoeff1,ncoeff2,coeff)332 GetPlate (wcs, ncoeff1, ncoeff2, coeff)
333
334 struct WorldCoor *wcs; /* World coordinate system structure */
335 int *ncoeff1; /* Number of coefficients for x */
336 int *ncoeff2; /* Number of coefficients for y) */
337 double *coeff; /* Plate fit coefficients */
338
339 {
340 int i;
341
342 if (nowcs (wcs))
343 return 1;
344
345 *ncoeff1 = wcs->ncoeff1;
346 *ncoeff2 = wcs->ncoeff2;
347
348 for (i = 0; i < *ncoeff1; i++)
349 coeff[i] = wcs->x_coeff[i];
350
351 for (i = 0; i < *ncoeff2; i++)
352 coeff[*ncoeff1+i] = wcs->y_coeff[i];
353
354 return 0;
355 }
356
357
358 /* Set FITS header plate fit coefficients from structure */
359 void
SetFITSPlate(header,wcs)360 SetFITSPlate (header, wcs)
361
362 char *header; /* Image FITS header */
363 struct WorldCoor *wcs; /* WCS structure */
364
365 {
366 char keyword[16];
367 int i;
368
369 for (i = 0; i < wcs->ncoeff1; i++) {
370 sprintf (keyword,"CO1_%d",i+1);
371 hputnr8 (header, keyword, -15, wcs->x_coeff[i]);
372 }
373 for (i = 0; i < wcs->ncoeff2; i++) {
374 sprintf (keyword,"CO2_%d",i+1);
375 hputnr8 (header, keyword, -15, wcs->y_coeff[i]);
376 }
377 return;
378 }
379
380 /* Mar 27 1998 New subroutines for direct image pixel <-> sky polynomials
381 * Apr 10 1998 Make terms identical for both x and y polynomials
382 * Apr 10 1998 Allow different numbers of coefficients for x and y
383 * Apr 16 1998 Drom NCOEFF header parameter
384 * Apr 28 1998 Change projection flags to WCS_*
385 * Sep 10 1998 Check for xc1 and yc2 divide by zero after Allen Harris, SAO
386 *
387 * Oct 21 1999 Drop unused variables after lint
388 *
389 * Feb 29 2000 Use inverse CD matrix to get initial X,Y in platepix()
390 * as suggested by Paolo Montegriffo from Bologna Ast. Obs.
391 */
392