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