1 #include <math.h>
2 #include "fitsio2.h"
3 #define D2R 0.01745329252
4 #define TWOPI 6.28318530717959
5
6 /*--------------------------------------------------------------------------*/
ffwldp(double xpix,double ypix,double xref,double yref,double xrefpix,double yrefpix,double xinc,double yinc,double rot,char * type,double * xpos,double * ypos,int * status)7 int ffwldp(double xpix, double ypix, double xref, double yref,
8 double xrefpix, double yrefpix, double xinc, double yinc, double rot,
9 char *type, double *xpos, double *ypos, int *status)
10
11 /* This routine is based on the classic AIPS WCS routine.
12
13 It converts from pixel location to RA,Dec for 9 projective geometries:
14 "-CAR", "-SIN", "-TAN", "-ARC", "-NCP", "-GLS", "-MER", "-AIT" and "-STG".
15 */
16
17 /*-----------------------------------------------------------------------*/
18 /* routine to determine accurate position for pixel coordinates */
19 /* returns 0 if successful otherwise: */
20 /* 501 = angle too large for projection; */
21 /* does: -CAR, -SIN, -TAN, -ARC, -NCP, -GLS, -MER, -AIT -STG projections*/
22 /* Input: */
23 /* f xpix x pixel number (RA or long without rotation) */
24 /* f ypiy y pixel number (dec or lat without rotation) */
25 /* d xref x reference coordinate value (deg) */
26 /* d yref y reference coordinate value (deg) */
27 /* f xrefpix x reference pixel */
28 /* f yrefpix y reference pixel */
29 /* f xinc x coordinate increment (deg) */
30 /* f yinc y coordinate increment (deg) */
31 /* f rot rotation (deg) (from N through E) */
32 /* c *type projection type code e.g. "-SIN"; */
33 /* Output: */
34 /* d *xpos x (RA) coordinate (deg) */
35 /* d *ypos y (dec) coordinate (deg) */
36 /*-----------------------------------------------------------------------*/
37 {double cosr, sinr, dx, dy, dz, temp, x, y, z;
38 double sins, coss, dect, rat, dt, l, m, mg, da, dd, cos0, sin0;
39 double dec0, ra0;
40 double geo1, geo2, geo3;
41 double deps = 1.0e-5;
42 char *cptr;
43
44 if (*status > 0)
45 return(*status);
46
47 /* Offset from ref pixel */
48 dx = (xpix-xrefpix) * xinc;
49 dy = (ypix-yrefpix) * yinc;
50
51 /* Take out rotation */
52 cosr = cos(rot * D2R);
53 sinr = sin(rot * D2R);
54 if (rot != 0.0) {
55 temp = dx * cosr - dy * sinr;
56 dy = dy * cosr + dx * sinr;
57 dx = temp;
58 }
59
60 /* convert to radians */
61 ra0 = xref * D2R;
62 dec0 = yref * D2R;
63
64 l = dx * D2R;
65 m = dy * D2R;
66 sins = l*l + m*m;
67 cos0 = cos(dec0);
68 sin0 = sin(dec0);
69
70 if (*type != '-') { /* unrecognized projection code */
71 return(*status = 504);
72 }
73
74 cptr = type + 1;
75
76 if (*cptr == 'C') { /* linear -CAR */
77 if (*(cptr + 1) != 'A' || *(cptr + 2) != 'R') {
78 return(*status = 504);
79 }
80 rat = ra0 + l;
81 dect = dec0 + m;
82
83 } else if (*cptr == 'T') { /* -TAN */
84 if (*(cptr + 1) != 'A' || *(cptr + 2) != 'N') {
85 return(*status = 504);
86 }
87 x = cos0*cos(ra0) - l*sin(ra0) - m*cos(ra0)*sin0;
88 y = cos0*sin(ra0) + l*cos(ra0) - m*sin(ra0)*sin0;
89 z = sin0 + m* cos0;
90 rat = atan2( y, x );
91 dect = atan ( z / sqrt(x*x+y*y) );
92
93 } else if (*cptr == 'S') {
94
95 if (*(cptr + 1) == 'I' && *(cptr + 2) == 'N') { /* -SIN */
96 if (sins>1.0)
97 return(*status = 501);
98 coss = sqrt (1.0 - sins);
99 dt = sin0 * coss + cos0 * m;
100 if ((dt>1.0) || (dt<-1.0))
101 return(*status = 501);
102 dect = asin (dt);
103 rat = cos0 * coss - sin0 * m;
104 if ((rat==0.0) && (l==0.0))
105 return(*status = 501);
106 rat = atan2 (l, rat) + ra0;
107
108 } else if (*(cptr + 1) == 'T' && *(cptr + 2) == 'G') { /* -STG Sterographic*/
109 dz = (4.0 - sins) / (4.0 + sins);
110 if (fabs(dz)>1.0)
111 return(*status = 501);
112 dect = dz * sin0 + m * cos0 * (1.0+dz) / 2.0;
113 if (fabs(dect)>1.0)
114 return(*status = 501);
115 dect = asin (dect);
116 rat = cos(dect);
117 if (fabs(rat)<deps)
118 return(*status = 501);
119 rat = l * (1.0+dz) / (2.0 * rat);
120 if (fabs(rat)>1.0)
121 return(*status = 501);
122 rat = asin (rat);
123 mg = 1.0 + sin(dect) * sin0 + cos(dect) * cos0 * cos(rat);
124 if (fabs(mg)<deps)
125 return(*status = 501);
126 mg = 2.0 * (sin(dect) * cos0 - cos(dect) * sin0 * cos(rat)) / mg;
127 if (fabs(mg-m)>deps)
128 rat = TWOPI /2.0 - rat;
129 rat = ra0 + rat;
130 } else {
131 return(*status = 504);
132 }
133
134 } else if (*cptr == 'A') {
135
136 if (*(cptr + 1) == 'R' && *(cptr + 2) == 'C') { /* ARC */
137 if (sins>=TWOPI*TWOPI/4.0)
138 return(*status = 501);
139 sins = sqrt(sins);
140 coss = cos (sins);
141 if (sins!=0.0)
142 sins = sin (sins) / sins;
143 else
144 sins = 1.0;
145 dt = m * cos0 * sins + sin0 * coss;
146 if ((dt>1.0) || (dt<-1.0))
147 return(*status = 501);
148 dect = asin (dt);
149 da = coss - dt * sin0;
150 dt = l * sins * cos0;
151 if ((da==0.0) && (dt==0.0))
152 return(*status = 501);
153 rat = ra0 + atan2 (dt, da);
154
155 } else if (*(cptr + 1) == 'I' && *(cptr + 2) == 'T') { /* -AIT Aitoff */
156 dt = yinc*cosr + xinc*sinr;
157 if (dt==0.0)
158 dt = 1.0;
159 dt = dt * D2R;
160 dy = yref * D2R;
161 dx = sin(dy+dt)/sqrt((1.0+cos(dy+dt))/2.0) -
162 sin(dy)/sqrt((1.0+cos(dy))/2.0);
163 if (dx==0.0)
164 dx = 1.0;
165 geo2 = dt / dx;
166 dt = xinc*cosr - yinc* sinr;
167 if (dt==0.0)
168 dt = 1.0;
169 dt = dt * D2R;
170 dx = 2.0 * cos(dy) * sin(dt/2.0);
171 if (dx==0.0) dx = 1.0;
172 geo1 = dt * sqrt((1.0+cos(dy)*cos(dt/2.0))/2.0) / dx;
173 geo3 = geo2 * sin(dy) / sqrt((1.0+cos(dy))/2.0);
174 rat = ra0;
175 dect = dec0;
176 if ((l != 0.0) || (m != 0.0)) {
177 dz = 4.0 - l*l/(4.0*geo1*geo1) - ((m+geo3)/geo2)*((m+geo3)/geo2) ;
178 if ((dz>4.0) || (dz<2.0)) return(*status = 501);
179 dz = 0.5 * sqrt (dz);
180 dd = (m+geo3) * dz / geo2;
181 if (fabs(dd)>1.0) return(*status = 501);
182 dd = asin (dd);
183 if (fabs(cos(dd))<deps) return(*status = 501);
184 da = l * dz / (2.0 * geo1 * cos(dd));
185 if (fabs(da)>1.0) return(*status = 501);
186 da = asin (da);
187 rat = ra0 + 2.0 * da;
188 dect = dd;
189 }
190 } else {
191 return(*status = 504);
192 }
193
194 } else if (*cptr == 'N') { /* -NCP North celestial pole*/
195 if (*(cptr + 1) != 'C' || *(cptr + 2) != 'P') {
196 return(*status = 504);
197 }
198 dect = cos0 - m * sin0;
199 if (dect==0.0)
200 return(*status = 501);
201 rat = ra0 + atan2 (l, dect);
202 dt = cos (rat-ra0);
203 if (dt==0.0)
204 return(*status = 501);
205 dect = dect / dt;
206 if ((dect>1.0) || (dect<-1.0))
207 return(*status = 501);
208 dect = acos (dect);
209 if (dec0<0.0) dect = -dect;
210
211 } else if (*cptr == 'G') { /* -GLS global sinusoid */
212 if (*(cptr + 1) != 'L' || *(cptr + 2) != 'S') {
213 return(*status = 504);
214 }
215 dect = dec0 + m;
216 if (fabs(dect)>TWOPI/4.0)
217 return(*status = 501);
218 coss = cos (dect);
219 if (fabs(l)>TWOPI*coss/2.0)
220 return(*status = 501);
221 rat = ra0;
222 if (coss>deps) rat = rat + l / coss;
223
224 } else if (*cptr == 'M') { /* -MER mercator*/
225 if (*(cptr + 1) != 'E' || *(cptr + 2) != 'R') {
226 return(*status = 504);
227 }
228 dt = yinc * cosr + xinc * sinr;
229 if (dt==0.0) dt = 1.0;
230 dy = (yref/2.0 + 45.0) * D2R;
231 dx = dy + dt / 2.0 * D2R;
232 dy = log (tan (dy));
233 dx = log (tan (dx));
234 geo2 = dt * D2R / (dx - dy);
235 geo3 = geo2 * dy;
236 geo1 = cos (yref*D2R);
237 if (geo1<=0.0) geo1 = 1.0;
238 rat = l / geo1 + ra0;
239 if (fabs(rat - ra0) > TWOPI)
240 return(*status = 501);
241 dt = 0.0;
242 if (geo2!=0.0) dt = (m + geo3) / geo2;
243 dt = exp (dt);
244 dect = 2.0 * atan (dt) - TWOPI / 4.0;
245
246 } else {
247 return(*status = 504);
248 }
249
250 /* correct for RA rollover */
251 if (rat-ra0>TWOPI/2.0) rat = rat - TWOPI;
252 if (rat-ra0<-TWOPI/2.0) rat = rat + TWOPI;
253 if (rat < 0.0) rat += TWOPI;
254
255 /* convert to degrees */
256 *xpos = rat / D2R;
257 *ypos = dect / D2R;
258 return(*status);
259 }
260 /*--------------------------------------------------------------------------*/
ffxypx(double xpos,double ypos,double xref,double yref,double xrefpix,double yrefpix,double xinc,double yinc,double rot,char * type,double * xpix,double * ypix,int * status)261 int ffxypx(double xpos, double ypos, double xref, double yref,
262 double xrefpix, double yrefpix, double xinc, double yinc, double rot,
263 char *type, double *xpix, double *ypix, int *status)
264
265 /* This routine is based on the classic AIPS WCS routine.
266
267 It converts from RA,Dec to pixel location to for 9 projective geometries:
268 "-CAR", "-SIN", "-TAN", "-ARC", "-NCP", "-GLS", "-MER", "-AIT" and "-STG".
269 */
270 /*-----------------------------------------------------------------------*/
271 /* routine to determine accurate pixel coordinates for an RA and Dec */
272 /* returns 0 if successful otherwise: */
273 /* 501 = angle too large for projection; */
274 /* 502 = bad values */
275 /* does: -SIN, -TAN, -ARC, -NCP, -GLS, -MER, -AIT projections */
276 /* anything else is linear */
277 /* Input: */
278 /* d xpos x (RA) coordinate (deg) */
279 /* d ypos y (dec) coordinate (deg) */
280 /* d xref x reference coordinate value (deg) */
281 /* d yref y reference coordinate value (deg) */
282 /* f xrefpix x reference pixel */
283 /* f yrefpix y reference pixel */
284 /* f xinc x coordinate increment (deg) */
285 /* f yinc y coordinate increment (deg) */
286 /* f rot rotation (deg) (from N through E) */
287 /* c *type projection type code e.g. "-SIN"; */
288 /* Output: */
289 /* f *xpix x pixel number (RA or long without rotation) */
290 /* f *ypiy y pixel number (dec or lat without rotation) */
291 /*-----------------------------------------------------------------------*/
292 {
293 double dx, dy, dz, r, ra0, dec0, ra, dec, coss, sins, dt, da, dd, sint;
294 double l, m, geo1, geo2, geo3, sinr, cosr, cos0, sin0;
295 double deps=1.0e-5;
296 char *cptr;
297
298 if (*type != '-') { /* unrecognized projection code */
299 return(*status = 504);
300 }
301
302 cptr = type + 1;
303
304 dt = (xpos - xref);
305 if (dt > 180) xpos -= 360;
306 if (dt < -180) xpos += 360;
307 /* NOTE: changing input argument xpos is OK (call-by-value in C!) */
308
309 /* default values - linear */
310 dx = xpos - xref;
311 dy = ypos - yref;
312
313 /* Correct for rotation */
314 r = rot * D2R;
315 cosr = cos (r);
316 sinr = sin (r);
317 dz = dx*cosr + dy*sinr;
318 dy = dy*cosr - dx*sinr;
319 dx = dz;
320
321 /* check axis increments - bail out if either 0 */
322 if ((xinc==0.0) || (yinc==0.0)) {*xpix=0.0; *ypix=0.0;
323 return(*status = 502);}
324
325 /* convert to pixels */
326 *xpix = dx / xinc + xrefpix;
327 *ypix = dy / yinc + yrefpix;
328
329 if (*cptr == 'C') { /* linear -CAR */
330 if (*(cptr + 1) != 'A' || *(cptr + 2) != 'R') {
331 return(*status = 504);
332 }
333
334 return(*status); /* done if linear */
335 }
336
337 /* Non linear position */
338 ra0 = xref * D2R;
339 dec0 = yref * D2R;
340 ra = xpos * D2R;
341 dec = ypos * D2R;
342
343 /* compute direction cosine */
344 coss = cos (dec);
345 sins = sin (dec);
346 cos0 = cos (dec0);
347 sin0 = sin (dec0);
348 l = sin(ra-ra0) * coss;
349 sint = sins * sin0 + coss * cos0 * cos(ra-ra0);
350
351 /* process by case */
352 if (*cptr == 'T') { /* -TAN tan */
353 if (*(cptr + 1) != 'A' || *(cptr + 2) != 'N') {
354 return(*status = 504);
355 }
356
357 if (sint<=0.0)
358 return(*status = 501);
359 if( cos0<0.001 ) {
360 /* Do a first order expansion around pole */
361 m = (coss * cos(ra-ra0)) / (sins * sin0);
362 m = (-m + cos0 * (1.0 + m*m)) / sin0;
363 } else {
364 m = ( sins/sint - sin0 ) / cos0;
365 }
366 if( fabs(sin(ra0)) < 0.3 ) {
367 l = coss*sin(ra)/sint - cos0*sin(ra0) + m*sin(ra0)*sin0;
368 l /= cos(ra0);
369 } else {
370 l = coss*cos(ra)/sint - cos0*cos(ra0) + m*cos(ra0)*sin0;
371 l /= -sin(ra0);
372 }
373
374 } else if (*cptr == 'S') {
375
376 if (*(cptr + 1) == 'I' && *(cptr + 2) == 'N') { /* -SIN */
377 if (sint<0.0)
378 return(*status = 501);
379 m = sins * cos(dec0) - coss * sin(dec0) * cos(ra-ra0);
380
381 } else if (*(cptr + 1) == 'T' && *(cptr + 2) == 'G') { /* -STG Sterographic*/
382 da = ra - ra0;
383 if (fabs(dec)>TWOPI/4.0)
384 return(*status = 501);
385 dd = 1.0 + sins * sin(dec0) + coss * cos(dec0) * cos(da);
386 if (fabs(dd)<deps)
387 return(*status = 501);
388 dd = 2.0 / dd;
389 l = l * dd;
390 m = dd * (sins * cos(dec0) - coss * sin(dec0) * cos(da));
391
392 } else {
393 return(*status = 504);
394 }
395
396 } else if (*cptr == 'A') {
397
398 if (*(cptr + 1) == 'R' && *(cptr + 2) == 'C') { /* ARC */
399 m = sins * sin(dec0) + coss * cos(dec0) * cos(ra-ra0);
400 if (m<-1.0) m = -1.0;
401 if (m>1.0) m = 1.0;
402 m = acos (m);
403 if (m!=0)
404 m = m / sin(m);
405 else
406 m = 1.0;
407 l = l * m;
408 m = (sins * cos(dec0) - coss * sin(dec0) * cos(ra-ra0)) * m;
409
410 } else if (*(cptr + 1) == 'I' && *(cptr + 2) == 'T') { /* -AIT Aitoff */
411 da = (ra - ra0) / 2.0;
412 if (fabs(da)>TWOPI/4.0)
413 return(*status = 501);
414 dt = yinc*cosr + xinc*sinr;
415 if (dt==0.0) dt = 1.0;
416 dt = dt * D2R;
417 dy = yref * D2R;
418 dx = sin(dy+dt)/sqrt((1.0+cos(dy+dt))/2.0) -
419 sin(dy)/sqrt((1.0+cos(dy))/2.0);
420 if (dx==0.0) dx = 1.0;
421 geo2 = dt / dx;
422 dt = xinc*cosr - yinc* sinr;
423 if (dt==0.0) dt = 1.0;
424 dt = dt * D2R;
425 dx = 2.0 * cos(dy) * sin(dt/2.0);
426 if (dx==0.0) dx = 1.0;
427 geo1 = dt * sqrt((1.0+cos(dy)*cos(dt/2.0))/2.0) / dx;
428 geo3 = geo2 * sin(dy) / sqrt((1.0+cos(dy))/2.0);
429 dt = sqrt ((1.0 + cos(dec) * cos(da))/2.0);
430 if (fabs(dt)<deps)
431 return(*status = 503);
432 l = 2.0 * geo1 * cos(dec) * sin(da) / dt;
433 m = geo2 * sin(dec) / dt - geo3;
434
435 } else {
436 return(*status = 504);
437 }
438
439 } else if (*cptr == 'N') { /* -NCP North celestial pole*/
440 if (*(cptr + 1) != 'C' || *(cptr + 2) != 'P') {
441 return(*status = 504);
442 }
443
444 if (dec0==0.0)
445 return(*status = 501); /* can't stand the equator */
446 else
447 m = (cos(dec0) - coss * cos(ra-ra0)) / sin(dec0);
448
449 } else if (*cptr == 'G') { /* -GLS global sinusoid */
450 if (*(cptr + 1) != 'L' || *(cptr + 2) != 'S') {
451 return(*status = 504);
452 }
453
454 dt = ra - ra0;
455 if (fabs(dec)>TWOPI/4.0)
456 return(*status = 501);
457 if (fabs(dec0)>TWOPI/4.0)
458 return(*status = 501);
459 m = dec - dec0;
460 l = dt * coss;
461
462 } else if (*cptr == 'M') { /* -MER mercator*/
463 if (*(cptr + 1) != 'E' || *(cptr + 2) != 'R') {
464 return(*status = 504);
465 }
466
467 dt = yinc * cosr + xinc * sinr;
468 if (dt==0.0) dt = 1.0;
469 dy = (yref/2.0 + 45.0) * D2R;
470 dx = dy + dt / 2.0 * D2R;
471 dy = log (tan (dy));
472 dx = log (tan (dx));
473 geo2 = dt * D2R / (dx - dy);
474 geo3 = geo2 * dy;
475 geo1 = cos (yref*D2R);
476 if (geo1<=0.0) geo1 = 1.0;
477 dt = ra - ra0;
478 l = geo1 * dt;
479 dt = dec / 2.0 + TWOPI / 8.0;
480 dt = tan (dt);
481 if (dt<deps)
482 return(*status = 502);
483 m = geo2 * log (dt) - geo3;
484
485 } else {
486 return(*status = 504);
487 }
488
489 /* convert to degrees */
490 dx = l / D2R;
491 dy = m / D2R;
492
493 /* Correct for rotation */
494 dz = dx*cosr + dy*sinr;
495 dy = dy*cosr - dx*sinr;
496 dx = dz;
497
498 /* convert to pixels */
499 *xpix = dx / xinc + xrefpix;
500 *ypix = dy / yinc + yrefpix;
501 return(*status);
502 }
503