1 /*** File saoimage/wcslib/dsspos.c
2 *** October 21, 1999
3 *** By Jessica Mink, jmink@cfa.harvard.edu
4 *** Harvard-Smithsonian Center for Astrophysics
5 *** Copyright (C) 1995-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: dsspos.c (Plate solution WCS conversion)
30 * Purpose: Compute WCS from Digital Sky Survey plate fit
31 * Subroutine: dsspos() converts from pixel location to RA,Dec
32 * Subroutine: dsspix() 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
dsspos(xpix,ypix,wcs,xpos,ypos)45 dsspos (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, xmm, ymm, xmm2, ymm2, xmm3, ymm3, x2y2;
62 double xi, xir, eta, etar, raoff, ra, dec;
63 double cond2r = 1.745329252e-2;
64 double cons2r = 206264.8062470964;
65 double twopi = 6.28318530717959;
66 double ctan, ccos;
67
68 /* Ignore magnitude and color terms
69 double mag = 0.0;
70 double color = 0.0; */
71
72 /* Convert from image pixels to plate pixels */
73 x = xpix + wcs->x_pixel_offset - 1.0 + 0.5;
74 y = ypix + wcs->y_pixel_offset - 1.0 + 0.5;
75
76 /* Convert from pixels to millimeters */
77 xmm = (wcs->ppo_coeff[2] - x * wcs->x_pixel_size) / 1000.0;
78 ymm = (y * wcs->y_pixel_size - wcs->ppo_coeff[5]) / 1000.0;
79 xmm2 = xmm * xmm;
80 ymm2 = ymm * ymm;
81 xmm3 = xmm * xmm2;
82 ymm3 = ymm * ymm2;
83 x2y2 = xmm2 + ymm2;
84
85 /* Compute coordinates from x,y and plate model */
86
87 xi = wcs->x_coeff[ 0]*xmm + wcs->x_coeff[ 1]*ymm +
88 wcs->x_coeff[ 2] + wcs->x_coeff[ 3]*xmm2 +
89 wcs->x_coeff[ 4]*xmm*ymm + wcs->x_coeff[ 5]*ymm2 +
90 wcs->x_coeff[ 6]*(x2y2) + wcs->x_coeff[ 7]*xmm3 +
91 wcs->x_coeff[ 8]*xmm2*ymm + wcs->x_coeff[ 9]*xmm*ymm2 +
92 wcs->x_coeff[10]*ymm3 + wcs->x_coeff[11]*xmm*(x2y2) +
93 wcs->x_coeff[12]*xmm*x2y2*x2y2;
94
95 /* Ignore magnitude and color terms
96 + wcs->x_coeff[13]*mag + wcs->x_coeff[14]*mag*mag +
97 wcs->x_coeff[15]*mag*mag*mag + wcs->x_coeff[16]*mag*xmm +
98 wcs->x_coeff[17]*mag*x2y2 + wcs->x_coeff[18]*mag*xmm*x2y2 +
99 wcs->x_coeff[19]*color; */
100
101 eta = wcs->y_coeff[ 0]*ymm + wcs->y_coeff[ 1]*xmm +
102 wcs->y_coeff[ 2] + wcs->y_coeff[ 3]*ymm2 +
103 wcs->y_coeff[ 4]*xmm*ymm + wcs->y_coeff[ 5]*xmm2 +
104 wcs->y_coeff[ 6]*(x2y2) + wcs->y_coeff[ 7]*ymm3 +
105 wcs->y_coeff[ 8]*ymm2*xmm + wcs->y_coeff[ 9]*ymm*xmm2 +
106 wcs->y_coeff[10]*xmm3 + wcs->y_coeff[11]*ymm*(x2y2) +
107 wcs->y_coeff[12]*ymm*x2y2*x2y2;
108
109 /* Ignore magnitude and color terms
110 + wcs->y_coeff[13]*mag + wcs->y_coeff[14]*mag*mag +
111 wcs->y_coeff[15]*mag*mag*mag + wcs->y_coeff[16]*mag*ymm +
112 wcs->y_coeff[17]*mag*x2y2) + wcs->y_coeff[18]*mag*ymm*x2y2 +
113 wcs->y_coeff[19]*color; */
114
115 /* Convert to radians */
116
117 xir = xi / cons2r;
118 etar = eta / cons2r;
119
120 /* Convert to RA and Dec */
121
122 ctan = tan (wcs->plate_dec);
123 ccos = cos (wcs->plate_dec);
124 raoff = atan2 (xir / ccos, 1.0 - etar * ctan);
125 ra = raoff + wcs->plate_ra;
126 if (ra < 0.0) ra = ra + twopi;
127 *xpos = ra / cond2r;
128
129 dec = atan (cos (raoff) * ((etar + ctan) / (1.0 - (etar * ctan))));
130 *ypos = dec / cond2r;
131 return 0;
132 }
133
134
135 int
dsspix(xpos,ypos,wcs,xpix,ypix)136 dsspix (xpos, ypos, wcs, xpix, ypix)
137
138 /* Routine to determine pixel coordinates for sky position */
139 /* returns 0 if successful otherwise 1 = angle too large for projection; */
140 /* based on amdinv() from getimage */
141
142 /* Input: */
143 double xpos; /* Right ascension or longitude in degrees */
144 double ypos; /* Declination or latitude in degrees */
145 struct WorldCoor *wcs; /* WCS parameter structure */
146
147 /* Output: */
148 double *xpix; /* x pixel number (RA or long without rotation) */
149 double *ypix; /* y pixel number (dec or lat without rotation) */
150
151 {
152 double div,xi,eta,x,y,xy,x2,y2,x2y,y2x,x3,y3,x4,y4,x2y2,cjunk,dx,dy;
153 double sypos,cypos,syplate,cyplate,sxdiff,cxdiff;
154 double f,fx,fy,g,gx,gy, xmm, ymm;
155 double conr2s = 206264.8062470964;
156 double tolerance = 0.0000005;
157 int max_iterations = 50;
158 int i;
159 double xr, yr; /* position in radians */
160
161 *xpix = 0.0;
162 *ypix = 0.0;
163
164 /* Convert RA and Dec in radians to standard coordinates on a plate */
165 xr = degrad (xpos);
166 yr = degrad (ypos);
167 sypos = sin (yr);
168 cypos = cos (yr);
169 if (wcs->plate_dec == 0.0)
170 wcs->plate_dec = degrad (wcs->yref);
171 syplate = sin (wcs->plate_dec);
172 cyplate = cos (wcs->plate_dec);
173 if (wcs->plate_ra == 0.0)
174 wcs->plate_ra = degrad (wcs->yref);
175 sxdiff = sin (xr - wcs->plate_ra);
176 cxdiff = cos (xr - wcs->plate_ra);
177 div = (sypos * syplate) + (cypos * cyplate * cxdiff);
178 if (div == 0.0)
179 return (1);
180 xi = cypos * sxdiff * conr2s / div;
181 eta = ((sypos * cyplate) - (cypos * syplate * cxdiff)) * conr2s / div;
182
183 /* Set initial value for x,y */
184 if (wcs->plate_scale == 0.0)
185 return (1);
186 xmm = xi / wcs->plate_scale;
187 ymm = eta / wcs->plate_scale;
188
189 /* Iterate by Newton's method */
190 for (i = 0; i < max_iterations; i++) {
191
192 /* X plate model */
193 xy = xmm * ymm;
194 x2 = xmm * xmm;
195 y2 = ymm * ymm;
196 x2y = x2 * ymm;
197 y2x = y2 * xmm;
198 x2y2 = x2 + y2;
199 cjunk = x2y2 * x2y2;
200 x3 = x2 * xmm;
201 y3 = y2 * ymm;
202 x4 = x2 * x2;
203 y4 = y2 * y2;
204 f = wcs->x_coeff[0]*xmm + wcs->x_coeff[1]*ymm +
205 wcs->x_coeff[2] + wcs->x_coeff[3]*x2 +
206 wcs->x_coeff[4]*xy + wcs->x_coeff[5]*y2 +
207 wcs->x_coeff[6]*x2y2 + wcs->x_coeff[7]*x3 +
208 wcs->x_coeff[8]*x2y + wcs->x_coeff[9]*y2x +
209 wcs->x_coeff[10]*y3 + wcs->x_coeff[11]*xmm*x2y2 +
210 wcs->x_coeff[12]*xmm*cjunk;
211 /* magnitude and color terms ignored
212 + wcs->x_coeff[13]*mag +
213 wcs->x_coeff[14]*mag*mag + wcs->x_coeff[15]*mag*mag*mag +
214 wcs->x_coeff[16]*mag*xmm + wcs->x_coeff[17]*mag*(x2+y2) +
215 wcs->x_coeff[18]*mag*xmm*(x2+y2) + wcs->x_coeff[19]*color;
216 */
217
218 /* Derivative of X model wrt x */
219 fx = wcs->x_coeff[0] + wcs->x_coeff[3]*2.0*xmm +
220 wcs->x_coeff[4]*ymm + wcs->x_coeff[6]*2.0*xmm +
221 wcs->x_coeff[7]*3.0*x2 + wcs->x_coeff[8]*2.0*xy +
222 wcs->x_coeff[9]*y2 + wcs->x_coeff[11]*(3.0*x2+y2) +
223 wcs->x_coeff[12]*(5.0*x4 +6.0*x2*y2+y4);
224 /* magnitude and color terms ignored
225 wcs->x_coeff[16]*mag + wcs->x_coeff[17]*mag*2.0*xmm +
226 wcs->x_coeff[18]*mag*(3.0*x2+y2);
227 */
228
229 /* Derivative of X model wrt y */
230 fy = wcs->x_coeff[1] + wcs->x_coeff[4]*xmm +
231 wcs->x_coeff[5]*2.0*ymm + wcs->x_coeff[6]*2.0*ymm +
232 wcs->x_coeff[8]*x2 + wcs->x_coeff[9]*2.0*xy +
233 wcs->x_coeff[10]*3.0*y2 + wcs->x_coeff[11]*2.0*xy +
234 wcs->x_coeff[12]*4.0*xy*x2y2;
235 /* magnitude and color terms ignored
236 wcs->x_coeff[17]*mag*2.0*ymm +
237 wcs->x_coeff[18]*mag*2.0*xy;
238 */
239
240 /* Y plate model */
241 g = wcs->y_coeff[0]*ymm + wcs->y_coeff[1]*xmm +
242 wcs->y_coeff[2] + wcs->y_coeff[3]*y2 +
243 wcs->y_coeff[4]*xy + wcs->y_coeff[5]*x2 +
244 wcs->y_coeff[6]*x2y2 + wcs->y_coeff[7]*y3 +
245 wcs->y_coeff[8]*y2x + wcs->y_coeff[9]*x2y +
246 wcs->y_coeff[10]*x3 + wcs->y_coeff[11]*ymm*x2y2 +
247 wcs->y_coeff[12]*ymm*cjunk;
248 /* magnitude and color terms ignored
249 wcs->y_coeff[13]*mag + wcs->y_coeff[14]*mag*mag +
250 wcs->y_coeff[15]*mag*mag*mag + wcs->y_coeff[16]*mag*ymm +
251 wcs->y_coeff[17]*mag*x2y2 +
252 wcs->y_coeff[18]*mag*ymm*x2y2 + wcs->y_coeff[19]*color;
253 */
254
255 /* Derivative of Y model wrt x */
256 gx = wcs->y_coeff[1] + wcs->y_coeff[4]*ymm +
257 wcs->y_coeff[5]*2.0*xmm + wcs->y_coeff[6]*2.0*xmm +
258 wcs->y_coeff[8]*y2 + wcs->y_coeff[9]*2.0*xy +
259 wcs->y_coeff[10]*3.0*x2 + wcs->y_coeff[11]*2.0*xy +
260 wcs->y_coeff[12]*4.0*xy*x2y2;
261 /* magnitude and color terms ignored
262 wcs->y_coeff[17]*mag*2.0*xmm +
263 wcs->y_coeff[18]*mag*ymm*2.0*xmm;
264 */
265
266 /* Derivative of Y model wrt y */
267 gy = wcs->y_coeff[0] + wcs->y_coeff[3]*2.0*ymm +
268 wcs->y_coeff[4]*xmm + wcs->y_coeff[6]*2.0*ymm +
269 wcs->y_coeff[7]*3.0*y2 + wcs->y_coeff[8]*2.0*xy +
270 wcs->y_coeff[9]*x2 + wcs->y_coeff[11]*(x2+3.0*y2) +
271 wcs->y_coeff[12]*(5.0*y4 + 6.0*x2*y2 + x4);
272 /* magnitude and color terms ignored
273 wcs->y_coeff[16]*mag + wcs->y_coeff[17]*mag*2.0*ymm +
274 wcs->y_coeff[18]*mag*(x2+3.0*y2);
275 */
276
277 f = f - xi;
278 g = g - eta;
279 dx = ((-f * gy) + (g * fy)) / ((fx * gy) - (fy * gx));
280 dy = ((-g * fx) + (f * gx)) / ((fx * gy) - (fy * gx));
281 xmm = xmm + dx;
282 ymm = ymm + dy;
283 if ((fabs(dx) < tolerance) && (fabs(dy) < tolerance)) break;
284 }
285
286 /* Convert mm from plate center to plate pixels */
287 if (wcs->x_pixel_size == 0.0 || wcs->y_pixel_size == 0.0)
288 return (1);
289 x = (wcs->ppo_coeff[2] - xmm*1000.0) / wcs->x_pixel_size;
290 y = (wcs->ppo_coeff[5] + ymm*1000.0) / wcs->y_pixel_size;
291
292 /* Convert from plate pixels to image pixels */
293 *xpix = x - wcs->x_pixel_offset + 1.0 - 0.5;
294 *ypix = y - wcs->y_pixel_offset + 1.0 - 0.5;
295
296 /* If position is off of the image, return offscale code */
297 if (*xpix < 0.5 || *xpix > wcs->nxpix+0.5)
298 return -1;
299 if (*ypix < 0.5 || *ypix > wcs->nypix+0.5)
300 return -1;
301
302 return 0;
303 }
304 /* Mar 6 1995 Original version of this code
305 * May 4 1995 Fix eta cross terms which were all in y
306 * Jun 21 1995 Add inverse routine
307 * Oct 17 1995 Fix inverse routine (degrees -> radians)
308 * Nov 7 1995 Add half pixel to image coordinates to get astrometric
309 * plate coordinates
310 * Feb 26 1996 Fix plate to image pixel conversion error
311 *
312 * Mar 23 1998 Change names from plate*() to dss*()
313 * Apr 7 1998 Change amd_i_coeff to i_coeff
314 * Sep 4 1998 Fix possible divide by zero in dsspos() from Allen Harris, SAO
315 * Sep 10 1998 Fix possible divide by zero in dsspix() from Allen Harris, SAO
316 *
317 * Oct 21 1999 Drop declaration of cond2r in dsspix()
318 */
319