public software.sextractor

[/] [trunk/] [src/] [fitswcs.c] - Blame information for rev 298

Details | Compare with Previous | View Log

Line No. Rev Author Line
1 95 bertin 
/*
2 233 bertin 
*                               fitswcs.c
3 95 bertin 
*
4 233 bertin 
* Manage World Coordinate System data.
5 95 bertin 
*
6 233 bertin 
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
7 95 bertin 
*
8 233 bertin 
*       This file part of:      AstrOmatic software
9 95 bertin 
*
10 284 bertin 
*       Copyright:              (C) 1993-2012 Emmanuel Bertin -- IAP/CNRS/UPMC
11 95 bertin 
*
12 233 bertin 
*       License:                GNU General Public License
13 
*
14 
*       AstrOmatic software is free software: you can redistribute it and/or
15 
*       modify it under the terms of the GNU General Public License as
16 
*       published by the Free Software Foundation, either version 3 of the
17 
*       License, or (at your option) any later version.
18 
*       AstrOmatic software is distributed in the hope that it will be useful,
19 
*       but WITHOUT ANY WARRANTY; without even the implied warranty of
20 
*       MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
21 
*       GNU General Public License for more details.
22 
*       You should have received a copy of the GNU General Public License
23 
*       along with AstrOmatic software.
24 
*       If not, see <http://www.gnu.org/licenses/>.
25 
*
26 298 bertin 
*       Last modified:          13/07/2012
27 233 bertin 
*
28 
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
29 95 bertin 
 
30 
#ifdef HAVE_CONFIG_H
31 
#include        "config.h"
32 
#endif
33 
 
34 
#ifdef HAVE_MATHIMF_H
35 
#include <mathimf.h>
36 
#else
37 
#include <math.h>
38 
#endif
39 
#include        <stdio.h>
40 
#include        <stdlib.h>
41 
#include        <string.h>
42 
 
43 
#include        "fits/fitscat_defs.h"
44 
#include        "fits/fitscat.h"
45 
#include        "fitswcs.h"
46 
#include        "wcscelsys.h"
47 
#include        "wcs/wcs.h"
48 
#include        "wcs/lin.h"
49 
#include        "wcs/tnx.h"
50 
#include        "wcs/poly.h"
51 
 
52 
/******* copy_wcs ************************************************************
53 
PROTO   wcsstruct *copy_wcs(wcsstruct *wcsin)
54 
PURPOSE Copy a WCS (World Coordinate System) structure.
55 
INPUT   WCS structure to be copied.
56 
OUTPUT  pointer to a copy of the input structure.
57 
NOTES   Actually, only FITS parameters are copied. Lower-level structures
58 
        such as those created by the WCS or TNX libraries are generated.
59 
AUTHOR  E. Bertin (IAP)
60 
VERSION 31/08/2002
61 
 ***/
62 
wcsstruct       *copy_wcs(wcsstruct *wcsin)
63 
 
64 
  {
65 
   wcsstruct    *wcs;
66 
 
67 
/* Copy the basic stuff */
68 
  QMEMCPY(wcsin, wcs, wcsstruct, 1);
69 
/* The PROJP WCS parameters */
70 
  QMEMCPY(wcsin->projp, wcs->projp, double, wcs->naxis*100);
71 
 
72 
/* Set other structure pointers to NULL (they'll have to be reallocated) */
73 
  wcs->wcsprm = NULL;
74 
  wcs->lin = NULL;
75 
  wcs->cel = NULL;
76 
  wcs->prj = NULL;
77 
  wcs->tnx_lngcor = copy_tnxaxis(wcsin->tnx_lngcor);
78 
  wcs->tnx_latcor = copy_tnxaxis(wcsin->tnx_latcor);
79 
  wcs->inv_x = wcs->inv_y = NULL;
80 
 
81 
  QCALLOC(wcs->wcsprm, struct wcsprm, 1);
82 
/* Test if the WCS is recognized and a celestial pair is found */
83 
  wcsset(wcs->naxis,(const char(*)[9])wcs->ctype, wcs->wcsprm);
84 
 
85 
/* Initialize other WCS structures */
86 
  init_wcs(wcs);
87 
/* Invert projection corrections */
88 
  invert_wcs(wcs);
89 
/* Find the range of coordinates */
90 
  range_wcs(wcs);
91 
 
92 
  return wcs;
93 
  }
94 
 
95 
 
96 
/******* create_wcs ***********************************************************
97 
PROTO   wcsstruct *create_wcs(char **ctype, double *crval, double *crpix,
98 
                        double *cdelt, int *naxisn, int naxis)
99 
PURPOSE Generate a simple WCS (World Coordinate System) structure.
100 
INPUT   Pointer to an array of char strings with WCS projection on each axis,
101 
        pointer to an array of center coordinates (double),
102 
        pointer to an array of device coordinates (double),
103 
        pointer to an array of pixel scales (double),
104 
        pointer to an array of image dimensions (int),
105 
        number of dimensions.
106 
OUTPUT  pointer to a WCS structure.
107 
NOTES   If a pointer is set to null, the corresponding variables are set to
108 
        default values.
109 
AUTHOR  E. Bertin (IAP)
110 
VERSION 09/08/2006
111 
 ***/
112 
wcsstruct       *create_wcs(char **ctype, double *crval, double *crpix,
113 
                        double *cdelt, int *naxisn, int naxis)
114 
 
115 
  {
116 
   wcsstruct    *wcs;
117 
   int          l;
118 
 
119 
  QCALLOC(wcs, wcsstruct, 1);
120 
  wcs->naxis = naxis;
121 
  QCALLOC(wcs->projp, double, naxis*100);
122 
  wcs->nprojp = 0;
123 
 
124 
  wcs->longpole = wcs->latpole = 999.0;
125 
  for (l=0; l<naxis; l++)
126 
    {
127 
    wcs->naxisn[l] = naxisn? naxisn[l] : 360.0;
128 
/*-- The default WCS projection system is an all-sky Aitoff projection */
129 
    if (ctype)
130 
      strncpy(wcs->ctype[l], ctype[l], 8);
131 
    else if (l==0)
132 
      strncpy(wcs->ctype[l], "RA---AIT", 8);
133 
    else if (l==1)
134 
      strncpy(wcs->ctype[l], "DEC--AIT", 8);
135 
    wcs->crval[l] = crval? crval[l]: 0.0;
136 
    wcs->crpix[l] = crpix? crpix[l]: 0.0;
137 
    wcs->cdelt[l] = 1.0;
138 
    wcs->cd[l*(naxis+1)] = cdelt? cdelt[l] : 1.0;
139 
    }
140 
 
141 
  wcs->epoch = wcs->equinox = 2000.0;
142 
  QCALLOC(wcs->wcsprm, struct wcsprm, 1);
143 
 
144 
/* Test if the WCS is recognized and a celestial pair is found */
145 
  wcsset(wcs->naxis,(const char(*)[9])wcs->ctype, wcs->wcsprm);
146 
 
147 
/* Initialize other WCS structures */
148 
  init_wcs(wcs);
149 
/* Invert projection corrections */
150 
  invert_wcs(wcs);
151 
/* Find the range of coordinates */
152 
  range_wcs(wcs);
153 
 
154 
  return wcs;
155 
  }
156 
 
157 
 
158 
/******* init_wcs ************************************************************
159 
PROTO   void init_wcs(wcsstruct *wcs)
160 
PURPOSE Initialize astrometry and WCS (World Coordinate System) structures.
161 
INPUT   WCS structure.
162 
OUTPUT  -.
163 
NOTES   -.
164 
AUTHOR  E. Bertin (IAP)
165 105 bertin 
VERSION 17/05/2007
166 95 bertin 
 ***/
167 
void    init_wcs(wcsstruct *wcs)
168 
 
169 
  {
170 
   int          l,n,lng,lat,naxis;
171 
 
172 
  naxis = wcs->naxis;
173 
  if (wcs->lin)
174 
    {
175 
    free(wcs->lin->cdelt);
176 
    free(wcs->lin->crpix);
177 
    free(wcs->lin->pc);
178 
    free(wcs->lin->piximg);
179 
    free(wcs->lin->imgpix);
180 
    free(wcs->lin);
181 
    }
182 
  QCALLOC(wcs->lin, struct linprm, 1);
183 
  QCALLOC(wcs->lin->cdelt, double, naxis);
184 
  QCALLOC(wcs->lin->crpix, double, naxis);
185 
  QCALLOC(wcs->lin->pc, double, naxis*naxis);
186 
 
187 
 
188 
  if (wcs->cel)
189 
    free(wcs->cel);
190 
  QCALLOC(wcs->cel, struct celprm, 1);
191 
 
192 
  if (wcs->prj)
193 
    free(wcs->prj);
194 
  QCALLOC(wcs->prj, struct prjprm, 1);
195 
 
196 
  if (wcs->inv_x)
197 
    {
198 
    poly_end(wcs->inv_x);
199 
    wcs->inv_x = NULL;
200 
    }
201 
  if (wcs->inv_y)
202 
    {
203 
    poly_end(wcs->inv_y);
204 
    wcs->inv_y = NULL;
205 
    }
206 
 
207 
/* Set WCS flags to 0: structures will be reinitialized by the WCS library */
208 
  wcs->lin->flag = wcs->cel->flag = wcs->prj->flag = 0;
209 
  wcs->lin->naxis = naxis;
210 
 
211 
/* wcsprm structure */
212 
  lng = wcs->lng = wcs->wcsprm->lng;
213 
  lat = wcs->lat = wcs->wcsprm->lat;
214 
 
215 
/* linprm structure */
216 
  for (l=0; l<naxis; l++)
217 
    {
218 
    wcs->lin->crpix[l] = wcs->crpix[l];
219 
    wcs->lin->cdelt[l] = 1.0;
220 
    }
221 
 
222 
  for (l=0; l<naxis*naxis; l++)
223 
    wcs->lin->pc[l] = wcs->cd[l];
224 
 
225 
/* celprm structure */
226 
  if (lng>=0)
227 
    {
228 
    wcs->cel->ref[0] = wcs->crval[lng];
229 
    wcs->cel->ref[1] = wcs->crval[lat];
230 
    }
231 
  else
232 
    {
233 
    wcs->cel->ref[0] = wcs->crval[0];
234 
    wcs->cel->ref[1] = wcs->crval[1];
235 
    }
236 
  wcs->cel->ref[2] = wcs->longpole;
237 
  wcs->cel->ref[3] = wcs->latpole;
238 
 
239 
/* prjprm structure */
240 
  wcs->prj->r0 = wcs->r0;
241 
  wcs->prj->tnx_lngcor = wcs->tnx_lngcor;
242 
  wcs->prj->tnx_latcor = wcs->tnx_latcor;
243 
  if (lng>=0)
244 
    {
245 
    n = 0;
246 
    for (l=100; l--;)
247 
      {
248 105 bertin 
      wcs->prj->p[l] = wcs->projp[l+lat*100];   /* lat comes first for ... */
249 
      wcs->prj->p[l+100] = wcs->projp[l+lng*100];/* ... compatibility reasons */
250 95 bertin 
      if (!n && (wcs->prj->p[l] || wcs->prj->p[l+100]))
251 
        n = l+1;
252 
      }
253 
    wcs->nprojp = n;
254 
    }
255 
 
256 
/* Check-out chirality */
257 
  wcs->chirality = wcs_chirality(wcs);
258 
 
259 
/* Initialize Equatorial <=> Celestial coordinate system transforms */
260 
  init_wcscelsys(wcs);
261 
 
262 
  return;
263 
  }
264 
 
265 
 
266 
/******* init_wcscelsys *******************************************************
267 
PROTO   void init_wcscelsys(wcsstruct *wcs)
268 
PURPOSE Initialize Equatorial <=> Celestial coordinate system transforms.
269 
INPUT   WCS structure.
270 
OUTPUT  -.
271 
NOTES   -.
272 
AUTHOR  E. Bertin (IAP)
273 
VERSION 18/07/2006
274 
 ***/
275 
void    init_wcscelsys(wcsstruct *wcs)
276 
 
277 
  {
278 
  double        *mat,
279 
                a0,d0,ap,dp,ap2,y;
280 
  int           s,lng,lat;
281 
 
282 
  lng = wcs->wcsprm->lng;
283 
  lat = wcs->wcsprm->lat;
284 
/* Is it a celestial system? If not, exit! */
285 
  if (lng==lat)
286 
    {
287 
    wcs->celsysconvflag = 0;
288 
    return;
289 
    }
290 
/* Find the celestial system */
291 
  for (s=0; *celsysname[s][0] && strncmp(wcs->ctype[lng], celsysname[s][0], 4);
292 
        s++);
293 
/* Is it a known, non-equatorial system? If not, exit! */
294 
  if (!s || !*celsysname[s][0])
295 
    {
296 
    wcs->celsysconvflag = 0;
297 
    return;
298 
    }
299 
  wcs->celsys = (celsysenum)s;
300 
/* Some shortcuts */
301 
  a0 = celsysorig[s][0]*DEG;
302 
  d0 = celsysorig[s][1]*DEG;
303 
  ap = celsyspole[s][0]*DEG;
304 
  dp = celsyspole[s][1]*DEG;
305 
/* First compute in the output referential the longitude of the south pole */
306 
  y = sin(ap - a0);
307 
/*
308 
  x = cos(d0)*(cos(d0)*sin(dp)*cos(ap-a0)-sin(d0)*cos(dp));
309 
  ap2 = atan2(y,x);
310 
*/
311 
  ap2 = asin(cos(d0)*y) ;
312 
/* Equatorial <=> Celestial System transformation parameters */
313 
  mat = wcs->celsysmat;
314 
  mat[0] = ap;
315 
  mat[1] = ap2;
316 
  mat[2] = cos(dp);
317 
  mat[3] = sin(dp);
318 
 
319 
  wcs->celsysconvflag = 1;
320 
  return;
321 
  }
322 
 
323 
 
324 
/******* read_wcs *************************************************************
325 
PROTO   wcsstruct *read_wcs(tabstruct *tab)
326 
PURPOSE Read WCS (World Coordinate System) info in the FITS header.
327 
INPUT   tab structure.
328 
OUTPUT  -.
329 
NOTES   -.
330 
AUTHOR  E. Bertin (IAP)
331 293 bertin 
VERSION 18/06/2012
332 95 bertin 
 ***/
333 
wcsstruct       *read_wcs(tabstruct *tab)
334 
 
335 
  {
336 
#define FITSREADF(buf, k, val, def) \
337 
                {if (fitsread(buf,k, &val, H_FLOAT,T_DOUBLE) != RETURN_OK) \
338 
                   val = def; \
339 
                }
340 
 
341 
#define FITSREADI(buf, k, val, def) \
342 
                {if (fitsread(buf,k, &val, H_INT,T_LONG) != RETURN_OK) \
343 
                   val = def; \
344 
                }
345 
 
346 
#define FITSREADS(buf, k, str, def) \
347 
                {if (fitsread(buf,k,str, H_STRING,T_STRING) != RETURN_OK) \
348 
                   strcpy(str, (def)); \
349 
                }
350 
   char         str[MAXCHARS];
351 
   char         wstr1[TNX_MAXCHARS], wstr2[TNX_MAXCHARS];
352 
 
353 
   wcsstruct    *wcs;
354 
   double       drota;
355 
   int          j, l, naxis;
356 
   char         name[16],
357 
                *buf, *filename, *ptr;
358 
 
359 
  buf = tab->headbuf;
360 
  filename = (tab->cat? tab->cat->filename : strcpy(name, "internal header"));
361 
 
362 
  FITSREADS(buf, "OBJECT  ", str, "Unnamed");
363 
 
364 
  QCALLOC(wcs, wcsstruct, 1);
365 105 bertin 
  if (tab->naxis > NAXIS)
366 
    {
367 
    warning("Maximum number of dimensions supported by this version of the ",
368 
        "software exceeded\n");
369 
    tab->naxis = 2;
370 
    }
371 
 
372 95 bertin 
  wcs->naxis = naxis = tab->naxis;
373 
  QCALLOC(wcs->projp, double, naxis*100);
374 
 
375 
  for (l=0; l<naxis; l++)
376 
    {
377 
    wcs->naxisn[l] = tab->naxisn[l];
378 
    sprintf(str, "CTYPE%-3d", l+1);
379 
    FITSREADS(buf, str, str, "");
380 
    strncpy(wcs->ctype[l], str, 8);
381 
    sprintf(str, "CUNIT%-3d", l+1);
382 
    FITSREADS(buf, str, str, "deg");
383 
    strncpy(wcs->cunit[l], str, 32);
384 
    sprintf(str, "CRVAL%-3d", l+1);
385 
    FITSREADF(buf, str, wcs->crval[l], 0.0);
386 
    sprintf(str, "CRPIX%-3d", l+1);
387 
    FITSREADF(buf, str, wcs->crpix[l], 1.0);
388 
    sprintf(str, "CDELT%-3d", l+1);
389 
    FITSREADF(buf, str, wcs->cdelt[l], 1.0);
390 
    sprintf(str, "CRDER%-3d", l+1);
391 
    FITSREADF(buf, str, wcs->crder[l], 0.0);
392 
    sprintf(str, "CSYER%-3d", l+1);
393 
    FITSREADF(buf, str, wcs->csyer[l], 0.0);
394 
    if (fabs(wcs->cdelt[l]) < 1e-30)
395 
      error(EXIT_FAILURE, "*Error*: CDELT parameters out of range in ",
396 
        filename);
397 
    }
398 
 
399 
  if (fitsfind(buf, "CD?_????")!=RETURN_ERROR)
400 
    {
401 
/*-- If CD keywords exist, use them for the linear mapping terms... */
402 
    for (l=0; l<naxis; l++)
403 
      for (j=0; j<naxis; j++)
404 
        {
405 
        sprintf(str, "CD%d_%d", l+1, j+1);
406 
        FITSREADF(buf, str, wcs->cd[l*naxis+j], l==j?1.0:0.0)
407 
        }
408 
    }
409 
  else if (fitsfind(buf, "PC00?00?")!=RETURN_ERROR)
410 
/*-- ...If PC keywords exist, use them for the linear mapping terms... */
411 
    for (l=0; l<naxis; l++)
412 
      for (j=0; j<naxis; j++)
413 
        {
414 
        sprintf(str, "PC%03d%03d", l+1, j+1);
415 
        FITSREADF(buf, str, wcs->cd[l*naxis+j], l==j?1.0:0.0)
416 
        wcs->cd[l*naxis+j] *= wcs->cdelt[l];
417 
        }
418 
  else
419 
    {
420 
/*-- ...otherwise take the obsolete CROTA2 parameter */
421 
    FITSREADF(buf, "CROTA2  ", drota, 0.0)
422 
    wcs->cd[3] = wcs->cd[0] = cos(drota*DEG);
423 
    wcs->cd[1] = -(wcs->cd[2] = sin(drota*DEG));
424 
    wcs->cd[0] *= wcs->cdelt[0];
425 
    wcs->cd[2] *= wcs->cdelt[0];
426 
    wcs->cd[1] *= wcs->cdelt[1];
427 
    wcs->cd[3] *= wcs->cdelt[1];
428 
    }
429 
  QCALLOC(wcs->wcsprm, struct wcsprm, 1);
430 
 
431 
/* Test if the WCS is recognized and a celestial pair is found */
432 
  if (!wcsset(wcs->naxis,(const char(*)[9])wcs->ctype, wcs->wcsprm)
433 
        && wcs->wcsprm->flag<999)
434 
    {
435 
     char       *pstr;
436 
     double     date;
437 
     int        biss, dpar[3];
438 
 
439 
/*-- Coordinate reference frame */
440 
/*-- Search for an observation date expressed in Julian days */
441 
    FITSREADF(buf, "MJD-OBS ", date, -1.0);
442 293 bertin 
    if (date<0.0)
443 
      FITSREADF(buf, "MJDSTART", date, -1.0);
444 95 bertin 
/*-- Precession date (defined from Ephemerides du Bureau des Longitudes) */
445 
/*-- in Julian years from 2000.0 */
446 
    if (date>0.0)
447 
      wcs->obsdate = 2000.0 - (MJD2000 - date)/365.25;
448 
    else
449 
      {
450 
/*---- Search for an observation date expressed in "civilian" format */
451 225 bertin 
      FITSREADS(buf, "DATE-OBS", str, "");
452 95 bertin 
      if (*str)
453 
        {
454 
/*------ Decode DATE-OBS format: DD/MM/YY or YYYY-MM-DD */
455 
        for (l=0; l<3 && (pstr = strtok_r(l?NULL:str,"/- ", &ptr)); l++)
456 
          dpar[l] = atoi(pstr);
457 
        if (l<3 || !dpar[0] || !dpar[1] || !dpar[2])
458 
          {
459 
/*-------- If DATE-OBS value corrupted or incomplete, assume 2000-1-1 */
460 
          warning("Invalid DATE-OBS value in header: ", str);
461 
          dpar[0] = 2000; dpar[1] = 1; dpar[2] = 1;
462 
          }
463 
        else if (strchr(str, '/') && dpar[0]<32 && dpar[2]<100)
464 
          {
465 
          j = dpar[0];
466 
          dpar[0] = dpar[2]+1900;
467 
          dpar[2] = j;
468 
          }
469 
 
470 
        biss = (dpar[0]%4)?0:1;
471 
/*------ Convert date to MJD */
472 
        date = -678956 + (365*dpar[0]+dpar[0]/4) - biss
473 
                        + ((dpar[1]>2?((int)((dpar[1]+1)*30.6)-63+biss)
474 
                :((dpar[1]-1)*(63+biss))/2) + dpar[2]);
475 
        wcs->obsdate = 2000.0 - (MJD2000 - date)/365.25;
476 
        }
477 
      else
478 
/*------ Well if really no date is found */
479 
        wcs->obsdate = 0.0;
480 
      }
481 
 
482 
    FITSREADF(buf, "EPOCH", wcs->epoch, 2000.0);
483 
    FITSREADF(buf, "EQUINOX", wcs->equinox, wcs->epoch);
484 290 bertin 
    if (fitsread(buf, "RADESYS", str, H_STRING,T_STRING) != RETURN_OK)
485 
      FITSREADS(buf, "RADECSYS", str,
486 95 bertin 
        wcs->equinox >= 2000.0? "ICRS" : (wcs->equinox<1984.0? "FK4" : "FK5"));
487 
    if (!strcmp(str, "ICRS"))
488 
      wcs->radecsys = RDSYS_ICRS;
489 
    else if (!strcmp(str, "FK5"))
490 
      wcs->radecsys = RDSYS_FK5;
491 
    else if (!strcmp(str, "FK4"))
492 
      {
493 
      if (wcs->equinox == 2000.0)
494 
        {
495 
        FITSREADF(buf, "EPOCH  ", wcs->equinox, 1950.0);
496 
        FITSREADF(buf, "EQUINOX", wcs->equinox, wcs->equinox);
497 
        }
498 
      wcs->radecsys = RDSYS_FK4;
499 
      warning("FK4 precession formulae not yet implemented:\n",
500 
                "            Astrometry may be slightly inaccurate");
501 
      }
502 
    else if (!strcmp(str, "FK4-NO-E"))
503 
      {
504 
      if (wcs->equinox == 2000.0)
505 
        {
506 
        FITSREADF(buf, "EPOCH", wcs->equinox, 1950.0);
507 
        FITSREADF(buf, "EQUINOX", wcs->equinox, wcs->equinox);
508 
        }
509 
      wcs->radecsys = RDSYS_FK4_NO_E;
510 
      warning("FK4 precession formulae not yet implemented:\n",
511 
                "            Astrometry may be slightly inaccurate");
512 
      }
513 
    else if (!strcmp(str, "GAPPT"))
514 
      {
515 
      wcs->radecsys = RDSYS_GAPPT;
516 
      warning("GAPPT reference frame not yet implemented:\n",
517 
                "            Astrometry may be slightly inaccurate");
518 
      }
519 
    else
520 
      {
521 
      warning("Using ICRS instead of unknown astrometric reference frame: ",
522 
                str);
523 
      wcs->radecsys = RDSYS_ICRS;
524 
      }
525 
 
526 
/*-- Projection parameters */
527 
    if (!strcmp(wcs->wcsprm->pcode, "TNX"))
528 
      {
529 
/*---- IRAF's TNX projection: decode these #$!?@#!! WAT parameters */
530 
      if (fitsfind(buf, "WAT?????") != RETURN_ERROR)
531 
        {
532 
/*------ First we need to concatenate strings */
533 
        pstr = wstr1;
534 
        sprintf(str, "WAT1_001");
535 
        for (j=2; fitsread(buf,str,pstr,H_STRINGS,T_STRING)==RETURN_OK; j++)
536 
          {
537 
          sprintf(str, "WAT1_%03d", j);
538 
          pstr += strlen(pstr);
539 
          }
540 
        pstr = wstr2;
541 
        sprintf(str, "WAT2_001");
542 
        for (j=2; fitsread(buf,str,pstr,H_STRINGS,T_STRING)==RETURN_OK; j++)
543 
          {
544 
          sprintf(str, "WAT2_%03d", j);
545 
          pstr += strlen(pstr);
546 
          }
547 
/*------ LONGPOLE defaulted to 180 deg if not found */
548 
        if ((pstr = strstr(wstr1, "longpole"))
549 
                || (pstr = strstr(wstr2, "longpole")))
550 
          pstr = strpbrk(pstr, "1234567890-+.");
551 
        wcs->longpole = pstr? atof(pstr) : 999.0;
552 
        wcs->latpole = 999.0;
553 
/*------ RO defaulted to 180/PI if not found */
554 
        if ((pstr = strstr(wstr1, "ro"))
555 
                || (pstr = strstr(wstr2, "ro")))
556 
          pstr = strpbrk(pstr, "1234567890-+.");
557 
        wcs->r0 = pstr? atof(pstr) : 0.0;
558 
/*------ Read the remaining TNX parameters */
559 
        if ((pstr = strstr(wstr1, "lngcor"))
560 
                || (pstr = strstr(wstr2, "lngcor")))
561 
          wcs->tnx_lngcor = read_tnxaxis(pstr);
562 
        if (!wcs->tnx_lngcor)
563 
          error(EXIT_FAILURE, "*Error*: incorrect TNX parameters in ",
564 
                        filename);
565 
        if ((pstr = strstr(wstr1, "latcor"))
566 
                || (pstr = strstr(wstr2, "latcor")))
567 
          wcs->tnx_latcor = read_tnxaxis(pstr);
568 
        if (!wcs->tnx_latcor)
569 
          error(EXIT_FAILURE, "*Error*: incorrect TNX parameters in ",
570 
                        filename);
571 
        }
572 
      }
573 
    else
574 
      {
575 
      FITSREADF(buf, "LONGPOLE", wcs->longpole, 999.0);
576 
      FITSREADF(buf, "LATPOLE ", wcs->latpole, 999.0);
577 
/*---- Old convention */
578 
      if (fitsfind(buf, "PROJP???") != RETURN_ERROR)
579 
        for (j=0; j<10; j++)
580 
          {
581 
          sprintf(str, "PROJP%-3d", j);
582 
          FITSREADF(buf, str, wcs->projp[j], 0.0);
583 
          }
584 
/*---- New convention */
585 
      if (fitsfind(buf, "PV?_????") != RETURN_ERROR)
586 
        for (l=0; l<naxis; l++)
587 
          for (j=0; j<100; j++)
588 
            {
589 177 bertin 
            sprintf(str, "PV%d_%d ", l+1, j);
590 95 bertin 
            FITSREADF(buf, str, wcs->projp[j+l*100], 0.0);
591 
            }
592 
      }
593 
    }
594 
 
595 
/* Initialize other WCS structures */
596 
  init_wcs(wcs);
597 105 bertin 
 
598 95 bertin 
/* Find the range of coordinates */
599 
  range_wcs(wcs);
600 
/* Invert projection corrections */
601 
  invert_wcs(wcs);
602 
 
603 
#undef FITSREADF
604 
#undef FITSREADI
605 
#undef FITSREADS
606 
 
607 
  return wcs;
608 
  }
609 
 
610 
 
611 
/******* write_wcs ***********************************************************
612 
PROTO   void write_wcs(tabstruct *tab, wcsstruct *wcs)
613 
PURPOSE Write WCS (World Coordinate System) info in the FITS header.
614 
INPUT   tab structure,
615 
        WCS structure.
616 
OUTPUT  -.
617 
NOTES   -.
618 
AUTHOR  E. Bertin (IAP)
619 233 bertin 
VERSION 01/09/2010
620 95 bertin 
 ***/
621 
void    write_wcs(tabstruct *tab, wcsstruct *wcs)
622 
 
623 
  {
624 233 bertin 
   double       mjd;
625 95 bertin 
   char         str[MAXCHARS];
626 
   int          j, l, naxis;
627 
 
628 
  naxis = wcs->naxis;
629 
  addkeywordto_head(tab, "BITPIX  ", "Bits per pixel");
630 
  fitswrite(tab->headbuf, "BITPIX  ", &tab->bitpix, H_INT, T_LONG);
631 
  addkeywordto_head(tab, "NAXIS   ", "Number of axes");
632 
  fitswrite(tab->headbuf, "NAXIS   ", &wcs->naxis, H_INT, T_LONG);
633 
  for (l=0; l<naxis; l++)
634 
    {
635 
    sprintf(str, "NAXIS%-3d", l+1);
636 
    addkeywordto_head(tab, str, "Number of pixels along this axis");
637 
    fitswrite(tab->headbuf, str, &wcs->naxisn[l], H_INT, T_LONG);
638 
    }
639 
  addkeywordto_head(tab, "EQUINOX ", "Mean equinox");
640 
  fitswrite(tab->headbuf, "EQUINOX ", &wcs->equinox, H_FLOAT, T_DOUBLE);
641 233 bertin 
  if (wcs->obsdate!=0.0)
642 
    {
643 
    mjd = (wcs->obsdate-2000.0)*365.25 + MJD2000;
644 
    addkeywordto_head(tab, "MJD-OBS ", "Modified Julian date at start");
645 
    fitswrite(tab->headbuf, "MJD-OBS ", &mjd, H_EXPO,T_DOUBLE);
646 
    }
647 95 bertin 
  addkeywordto_head(tab, "RADECSYS", "Astrometric system");
648 
  switch(wcs->radecsys)
649 
    {
650 
    case RDSYS_ICRS:
651 
      fitswrite(tab->headbuf, "RADECSYS", "ICRS", H_STRING, T_STRING);
652 
      break;
653 
    case RDSYS_FK5:
654 
      fitswrite(tab->headbuf, "RADECSYS", "FK5", H_STRING, T_STRING);
655 
      break;
656 
    case RDSYS_FK4:
657 
      fitswrite(tab->headbuf, "RADECSYS", "FK4", H_STRING, T_STRING);
658 
      break;
659 
    case RDSYS_FK4_NO_E:
660 
      fitswrite(tab->headbuf, "RADECSYS", "FK4-NO-E", H_STRING, T_STRING);
661 
      break;
662 
    case RDSYS_GAPPT:
663 
      fitswrite(tab->headbuf, "RADECSYS", "GAPPT", H_STRING, T_STRING);
664 
      break;
665 
    default:
666 
      error(EXIT_FAILURE, "*Error*: unknown RADECSYS type in write_wcs()", "");
667 
    }
668 
  for (l=0; l<naxis; l++)
669 
    {
670 
    sprintf(str, "CTYPE%-3d", l+1);
671 
    addkeywordto_head(tab, str, "WCS projection type for this axis");
672 
    fitswrite(tab->headbuf, str, wcs->ctype[l], H_STRING, T_STRING);
673 
    sprintf(str, "CUNIT%-3d", l+1);
674 
    addkeywordto_head(tab, str, "Axis unit");
675 
    fitswrite(tab->headbuf, str, wcs->cunit[l], H_STRING, T_STRING);
676 
    sprintf(str, "CRVAL%-3d", l+1);
677 
    addkeywordto_head(tab, str, "World coordinate on this axis");
678 
    fitswrite(tab->headbuf, str, &wcs->crval[l], H_EXPO, T_DOUBLE);
679 
    sprintf(str, "CRPIX%-3d", l+1);
680 
    addkeywordto_head(tab, str, "Reference pixel on this axis");
681 
    fitswrite(tab->headbuf, str, &wcs->crpix[l], H_EXPO, T_DOUBLE);
682 
    for (j=0; j<naxis; j++)
683 
      {
684 
      sprintf(str, "CD%d_%d", l+1, j+1);
685 
      addkeywordto_head(tab, str, "Linear projection matrix");
686 
      fitswrite(tab->headbuf, str, &wcs->cd[l*naxis+j], H_EXPO, T_DOUBLE);
687 
      }
688 
    for (j=0; j<100; j++)
689 
      if (wcs->projp[j+100*l] != 0.0)
690 
        {
691 
        sprintf(str, "PV%d_%d", l+1, j);
692 
        addkeywordto_head(tab, str, "Projection distortion parameter");
693 
        fitswrite(tab->headbuf, str, &wcs->projp[j+100*l], H_EXPO, T_DOUBLE);
694 
        }
695 
    }
696 
 
697 
/* Update the tab data */
698 
  readbasic_head(tab);
699 
 
700 
  return;
701 
  }
702 
 
703 
 
704 
/******* end_wcs **************************************************************
705 
PROTO   void end_wcs(wcsstruct *wcs)
706 
PURPOSE Free WCS (World Coordinate System) infos.
707 
INPUT   WCS structure.
708 
OUTPUT  -.
709 
NOTES   .
710 
AUTHOR  E. Bertin (IAP)
711 
VERSION 24/05/2000
712 
 ***/
713 
void    end_wcs(wcsstruct *wcs)
714 
 
715 
  {
716 
  if (wcs)
717 
    {
718 
    if (wcs->lin)
719 
      {
720 
      free(wcs->lin->cdelt);
721 
      free(wcs->lin->crpix);
722 
      free(wcs->lin->pc);
723 
      free(wcs->lin->piximg);
724 
      free(wcs->lin->imgpix);
725 
      free(wcs->lin);
726 
      }
727 
    free(wcs->cel);
728 
    free(wcs->prj);
729 
    free(wcs->wcsprm);
730 
    free_tnxaxis(wcs->tnx_lngcor);
731 
    free_tnxaxis(wcs->tnx_latcor);
732 
    poly_end(wcs->inv_x);
733 
    poly_end(wcs->inv_y);
734 
    free(wcs->projp);
735 
    free(wcs);
736 
    }
737 
 
738 
  return;
739 
  }
740 
 
741 
 
742 
/******* wcs_supproj *********************************************************
743 
PROTO   int wcs_supproj(char *name)
744 
PURPOSE Tell if a projection system is supported or not.
745 
INPUT   Proposed projection code name.
746 
OUTPUT  RETURN_OK if projection is supported, RETURN_ERROR otherwise.
747 
NOTES   -.
748 
AUTHOR  E. Bertin (IAP)
749 293 bertin 
VERSION 14/06/2012
750 95 bertin 
 ***/
751 
int     wcs_supproj(char *name)
752 
 
753 
  {
754 284 bertin 
   char projcode[28][5] =
755 293 bertin 
        {"TAN", "TPV", "AZP", "SIN", "STG", "ARC", "ZPN", "ZEA", "AIR", "CYP",
756 
        "CAR", "MER", "CEA", "COP", "COD", "COE", "COO", "BON", "PCO", "GLS",
757 
        "PAR", "AIT", "MOL", "CSC", "QSC", "TSC", "TNX", "NONE"};
758 
 
759 95 bertin 
   int  i;
760 
 
761 284 bertin 
  for (i=0; i<28; i++)
762 95 bertin 
    if (!strcmp(name, projcode[i]))
763 
      return RETURN_OK;
764 
 
765 
  return RETURN_ERROR;
766 
  }
767 
 
768 
 
769 
/******* invert_wcs ***********************************************************
770 
PROTO   void invert_wcs(wcsstruct *wcs)
771 
PURPOSE Invert WCS projection mapping (using a polynomial).
772 
INPUT   WCS structure.
773 
OUTPUT  -.
774 
NOTES   .
775 
AUTHOR  E. Bertin (IAP)
776 298 bertin 
VERSION 13/07/2012
777 95 bertin 
 ***/
778 
void    invert_wcs(wcsstruct *wcs)
779 
 
780 
  {
781 
   polystruct           *poly;
782 
   double               pixin[NAXIS],raw[NAXIS],rawmin[NAXIS];
783 
   double               *outpos,*outpost, *lngpos,*lngpost,
784 
                        *latpos,*latpost,
785 
                        lngstep,latstep, rawsize, epsilon;
786 
   int                  group[] = {1,1};
787 
                                /* Don't ask, this is needed by poly_init()! */
788 
   int          i,j,lng,lat,deg, tnxflag, maxflag;
789 
 
790 
/* Check first that inversion is not straightforward */
791 
  lng = wcs->wcsprm->lng;
792 
  lat = wcs->wcsprm->lat;
793 
  if (!strcmp(wcs->wcsprm->pcode, "TNX"))
794 
    tnxflag = 1;
795 298 bertin 
  else if ((!strcmp(wcs->wcsprm->pcode, "TAN")
796 293 bertin 
        || !strcmp(wcs->wcsprm->pcode, "TPV"))
797 95 bertin 
                && (wcs->projp[1+lng*100] || wcs->projp[1+lat*100]))
798 
    tnxflag = 0;
799 
  else
800 
    return;
801 
 
802 
/* We define x as "longitude" and y as "latitude" projections */
803 
/* We assume that PCxx cross-terms with additional dimensions are small */
804 
/* Sample the whole image with a regular grid */
805 
  lngstep = wcs->naxisn[lng]/(WCS_NGRIDPOINTS-1.0);
806 
  latstep = wcs->naxisn[lat]/(WCS_NGRIDPOINTS-1.0);
807 
  QMALLOC(outpos, double, 2*WCS_NGRIDPOINTS2);
808 
  QMALLOC(lngpos, double, WCS_NGRIDPOINTS2);
809 
  QMALLOC(latpos, double, WCS_NGRIDPOINTS2);
810 
  for (i=0; i<wcs->naxis; i++)
811 
    raw[i] = rawmin[i] = 0.5;
812 
  outpost = outpos;
813 
  lngpost = lngpos;
814 
  latpost = latpos;
815 
  for (j=WCS_NGRIDPOINTS; j--; raw[lat]+=latstep)
816 
    {
817 
    raw[lng] = rawmin[lng];
818 
    for (i=WCS_NGRIDPOINTS; i--; raw[lng]+=lngstep)
819 
      {
820 
      if (linrev(raw, wcs->lin, pixin))
821 
        error(EXIT_FAILURE, "*Error*: incorrect linear conversion in ",
822 
                wcs->wcsprm->pcode);
823 
      *(lngpost++) = pixin[lng];
824 
      *(latpost++) = pixin[lat];
825 
      if (tnxflag)
826 
        {
827 
        *(outpost++) = pixin[lng]
828 
                        +raw_to_tnxaxis(wcs->tnx_lngcor,pixin[lng],pixin[lat]);
829 
        *(outpost++) = pixin[lat]
830 
                        +raw_to_tnxaxis(wcs->tnx_latcor,pixin[lng],pixin[lat]);
831 
        }
832 
      else
833 
        {
834 
        raw_to_pv(wcs->prj,pixin[lng],pixin[lat], outpost, outpost+1);
835 
        outpost += 2;
836 
        }
837 
      }
838 
    }
839 
 
840 
/* Invert "longitude" */
841 
/* Compute the extent of the pixel in reduced projected coordinates */
842 
  linrev(rawmin, wcs->lin, pixin);
843 
  pixin[lng] += ARCSEC/DEG;
844 
  linfwd(pixin, wcs->lin, raw);
845 
  rawsize = sqrt((raw[lng]-rawmin[lng])*(raw[lng]-rawmin[lng])
846 
                +(raw[lat]-rawmin[lat])*(raw[lat]-rawmin[lat]))*DEG/ARCSEC;
847 
  if (!rawsize)
848 
    error(EXIT_FAILURE, "*Error*: incorrect linear conversion in ",
849 
                wcs->wcsprm->pcode);
850 
  epsilon = WCS_INVACCURACY/rawsize;
851 
/* Find the lowest degree polynom */
852 
  poly = NULL;  /* to avoid gcc -Wall warnings */
853 
  maxflag = 1;
854 
  for (deg=1; deg<=WCS_INVMAXDEG && maxflag; deg++)
855 
    {
856 
    if (deg>1)
857 
      poly_end(poly);
858 
    poly = poly_init(group, 2, &deg, 1);
859 
    poly_fit(poly, outpos, lngpos, NULL, WCS_NGRIDPOINTS2, NULL);
860 
    maxflag = 0;
861 
    outpost = outpos;
862 
    lngpost = lngpos;
863 
    for (i=WCS_NGRIDPOINTS2; i--; outpost+=2)
864 
      if (fabs(poly_func(poly, outpost)-*(lngpost++))>epsilon)
865 
        {
866 
        maxflag = 1;
867 
        break;
868 
        }
869 
    }
870 
  if (maxflag)
871 
    warning("Significant inaccuracy likely to occur in projection","");
872 
/* Now link the created structure */
873 
  wcs->prj->inv_x = wcs->inv_x = poly;
874 
 
875 
/* Invert "latitude" */
876 
/* Compute the extent of the pixel in reduced projected coordinates */
877 
  linrev(rawmin, wcs->lin, pixin);
878 
  pixin[lat] += ARCSEC/DEG;
879 
  linfwd(pixin, wcs->lin, raw);
880 
  rawsize = sqrt((raw[lng]-rawmin[lng])*(raw[lng]-rawmin[lng])
881 
                +(raw[lat]-rawmin[lat])*(raw[lat]-rawmin[lat]))*DEG/ARCSEC;
882 
  if (!rawsize)
883 
    error(EXIT_FAILURE, "*Error*: incorrect linear conversion in ",
884 
                wcs->wcsprm->pcode);
885 
  epsilon = WCS_INVACCURACY/rawsize;
886 
/* Find the lowest degree polynom */
887 
  maxflag = 1;
888 
  for (deg=1; deg<=WCS_INVMAXDEG && maxflag; deg++)
889 
    {
890 
    if (deg>1)
891 
      poly_end(poly);
892 
    poly = poly_init(group, 2, &deg, 1);
893 
    poly_fit(poly, outpos, latpos, NULL, WCS_NGRIDPOINTS2, NULL);
894 
    maxflag = 0;
895 
    outpost = outpos;
896 
    latpost = latpos;
897 
    for (i=WCS_NGRIDPOINTS2; i--; outpost+=2)
898 
      if (fabs(poly_func(poly, outpost)-*(latpost++))>epsilon)
899 
        {
900 
        maxflag = 1;
901 
        break;
902 
        }
903 
    }
904 
  if (maxflag)
905 
    warning("Significant inaccuracy likely to occur in projection","");
906 
/* Now link the created structure */
907 
  wcs->prj->inv_y = wcs->inv_y = poly;
908 
 
909 
/* Free memory */
910 
  free(outpos);
911 
  free(lngpos);
912 
  free(latpos);
913 
 
914 
  return;
915 
  }
916 
 
917 
 
918 
/******* range_wcs ***********************************************************
919 
PROTO   void range_wcs(wcsstruct *wcs)
920 
PURPOSE Find roughly the range of WCS coordinates on all axes,
921 
        and typical pixel scales.
922 
INPUT   WCS structure.
923 
OUTPUT  -.
924 
NOTES   .
925 
AUTHOR  E. Bertin (IAP)
926 233 bertin 
VERSION 24/08/2010
927 95 bertin 
 ***/
928 
void    range_wcs(wcsstruct *wcs)
929 
 
930 
  {
931 
   double               step[NAXIS], raw[NAXIS], rawmin[NAXIS],
932 
                        world[NAXIS], world2[NAXIS];
933 
   double               *worldmin, *worldmax, *scale, *worldc,
934 
                        rad, radmax, lc;
935 
   int                  linecount[NAXIS];
936 
   int                  i,j, naxis, npoints, lng,lat;
937 
 
938 
  naxis = wcs->naxis;
939 
 
940 
/* World range */
941 
  npoints = 1;
942 
  worldmin = wcs->wcsmin;
943 
  worldmax = wcs->wcsmax;
944 
/* First, find the center and use it as a reference point for lng */
945 
  lng = wcs->lng;
946 
  lat = wcs->lat;
947 
  for (i=0; i<naxis; i++)
948 
    raw[i] = (wcs->naxisn[i]+1.0)/2.0;
949 
  if (raw_to_wcs(wcs, raw, world))
950 
    {
951 
/*-- Oops no mapping there! So explore the image in an increasingly large */
952 
/*-- domain to find  a better "center" (now we know there must be angular */
953 
/*-- coordinates) */
954 
    for (j=0; j<100; j++)
955 
      {
956 
      for (i=0; i<naxis; i++)
957 
        raw[i] += wcs->naxisn[i]/100.0*(0.5-(double)rand()/RAND_MAX);
958 
      if (!raw_to_wcs(wcs, raw, world))
959 
        break;
960 
      }
961 
    }
962 
 
963 
  if (lng!=lat)
964 233 bertin 
    lc = world[lng];
965 95 bertin 
  else
966 
    {
967 
    lc = 0.0;   /* to avoid gcc -Wall warnings */
968 
    lng = -1;
969 
    }
970 
 
971 
/* Pixel scales at image center */
972 
  scale = wcs->wcsscale;
973 
  for (i=0; i<naxis; i++)
974 
    {
975 
    if ((i==lng || i==lat) && lng!=lat)
976 
      wcs->pixscale = scale[i] = sqrt(wcs_scale(wcs, raw));
977 
    else
978 
      {
979 
      raw[i] += 1.0;
980 
      raw_to_wcs(wcs, raw, world2);
981 
      scale[i] = fabs(world2[i] - world[i]);
982 
      raw[i] -= 1.0;
983 
      if (lng==lat)
984 
        wcs->pixscale = scale[i];
985 
      }
986 
    wcs->wcsscalepos[i] = world[i];
987 
    }
988 
 
989 
 
990 
/* Find "World limits" */
991 
  for (i=0; i<naxis; i++)
992 
    {
993 
    raw[i] = rawmin[i] = 0.5;
994 
    step[i] = wcs->naxisn[i]/(WCS_NRANGEPOINTS-1.0);
995 
    npoints *= WCS_NRANGEPOINTS;
996 
    worldmax[i] = -(worldmin[i] = 1e31);
997 
    linecount[i] = 0;
998 
    }
999 
 
1000 
  radmax = 0.0;
1001 
  worldc = wcs->wcsscalepos;
1002 
 
1003 
  for (j=npoints; j--;)
1004 
    {
1005 
    raw_to_wcs(wcs, raw, world);
1006 233 bertin 
/*-- Compute maximum distance to center */
1007 
    if ((rad=wcs_dist(wcs, world, worldc)) > radmax)
1008 
      radmax = rad;
1009 95 bertin 
    for (i=0; i<naxis; i++)
1010 
      {
1011 
/*---- Handle longitudes around 0 */
1012 233 bertin 
      if (i==lng)
1013 
        {
1014 
        world[i] -= lc;
1015 
        if (world[i]>180.0)
1016 
          world[i] -= 360.0;
1017 
        else if (world[i] <= -180.0)
1018 
          world[i] += 360.0;
1019 
        }
1020 95 bertin 
      if (world[i]<worldmin[i])
1021 
        worldmin[i] = world[i];
1022 
      if (world[i]>worldmax[i])
1023 
        worldmax[i] = world[i];
1024 
      }
1025 
 
1026 
 
1027 
    for (i=0; i<naxis; i++)
1028 
      {
1029 
      raw[i] += step[i];
1030 
      if (++linecount[i]<WCS_NRANGEPOINTS)
1031 
        break;
1032 
      else
1033 
        {
1034 
        linecount[i] = 0;       /* No need to initialize it to 0! */
1035 
        raw[i] = rawmin[i];
1036 
        }
1037 
      }
1038 
    }
1039 
 
1040 
  wcs->wcsmaxradius = radmax;
1041 
 
1042 
  if (lng!=lat)
1043 
    {
1044 233 bertin 
    worldmin[lng] = fmod_0_p360(worldmin[lng]+lc);
1045 
    worldmax[lng] = fmod_0_p360(worldmax[lng]+lc);
1046 95 bertin 
    if (worldmax[lat]<-90.0)
1047 
      worldmax[lat] = -90.0;
1048 
    if (worldmax[lat]>90.0)
1049 
      worldmax[lat] = 90.0;
1050 
    }
1051 
 
1052 
  return;
1053 
  }
1054 
 
1055 
 
1056 
/******* frame_wcs ***********************************************************
1057 284 bertin 
PROTO   int frame_wcs(wcsstruct *wcsin, wcsstruct *wcsout)
1058 95 bertin 
PURPOSE Find the x and y limits of an input frame in an output image.
1059 
INPUT   WCS structure of the input frame,
1060 
        WCS structure of the output frame.
1061 284 bertin 
OUTPUT  1 if frames overlap, 0 otherwise.
1062 
NOTES   -.
1063 95 bertin 
AUTHOR  E. Bertin (IAP)
1064 284 bertin 
VERSION 26/03/2012
1065 95 bertin 
 ***/
1066 284 bertin 
int     frame_wcs(wcsstruct *wcsin, wcsstruct *wcsout)
1067 95 bertin 
 
1068 
  {
1069 
   double               rawin[NAXIS], rawout[NAXIS], world[NAXIS];
1070 
   int                  linecount[NAXIS];
1071 
   double               worldc;
1072 
   int                  *min, *max,
1073 284 bertin 
                        i,j, naxis, npoints, out, swapflag, overlapflag;
1074 95 bertin 
 
1075 
  naxis = wcsin->naxis;
1076 
 
1077 
/* World range */
1078 
  npoints = 1;
1079 
  min = wcsin->outmin;
1080 
  max = wcsin->outmax;
1081 
  for (i=0; i<naxis; i++)
1082 
    {
1083 
    rawin[i] = 0.5;     /* Lower pixel limits */
1084 
    npoints *= WCS_NRANGEPOINTS;
1085 
    max[i] = -(min[i] = 1<<30);
1086 
    linecount[i] = 0;
1087 
    }
1088 
 
1089 
/* Check if lng and lat are swapped between in and out wcs (vicious idea!) */
1090 
  swapflag = (((wcsin->lng != wcsout->lng) || (wcsin->lat != wcsout->lat))
1091 
        && (wcsin->lng != wcsin->lat) && (wcsout->lng != wcsout->lat));
1092 
 
1093 
  for (j=npoints; j--;)
1094 
    {
1095 
    if (!raw_to_wcs(wcsin, rawin, world))
1096 
      {
1097 
      if (swapflag)
1098 
        {
1099 
        worldc = world[wcsout->lat];
1100 
        world[wcsout->lat] = world[wcsin->lat];
1101 
        world[wcsin->lat] = worldc;
1102 
        }
1103 
      if (!wcs_to_raw(wcsout, world, rawout))
1104 
        for (i=0; i<naxis; i++)
1105 
          {
1106 
          if ((out=(int)(rawout[i]+0.499))<min[i])
1107 
            min[i] = out;
1108 
          if (out>max[i])
1109 
            max[i] = out;
1110 
          }
1111 
      }
1112 
    for (i=0; i<naxis; i++)
1113 
      {
1114 
      rawin[i] = 0.5 + 0.5*wcsin->naxisn[i]
1115 
        *(1-cos(PI*(linecount[i]+1.0)/(WCS_NRANGEPOINTS-1)));
1116 
      if (++linecount[i]<WCS_NRANGEPOINTS)
1117 
        break;
1118 
      else
1119 
        {
1120 
        linecount[i] = 0;       /* No need to initialize it to 0! */
1121 
        rawin[i] =  0.5;
1122 
        }
1123 
      }
1124 
    }
1125 
 
1126 284 bertin 
/* Add a little margin, just in case... */
1127 95 bertin 
  for (i=0; i<naxis; i++)
1128 
    {
1129 
    if (min[i]>-2147483647)
1130 
      min[i] -= 2;
1131 
    if (max[i]>2147483647)
1132 
      max[i] += 2;
1133 
    }
1134 
 
1135 284 bertin 
/* Check overlap */
1136 
  overlapflag = 1;
1137 
  for (i=0; i<naxis; i++)
1138 
    if (min[i]>wcsout->naxisn[i] || max[i]<0)
1139 
      {
1140 
      overlapflag = 0;
1141 
      break;
1142 
      }
1143 
 
1144 
  return overlapflag;
1145 95 bertin 
  }
1146 
 
1147 
 
1148 
/******* reaxe_wcs ***********************************************************
1149 
PROTO   int reaxe_wcs(wcsstruct *wcs, int lng, int lat)
1150 
PURPOSE Reformulate a wcs structure to match lng and lat axis indices
1151 
INPUT   WCS structure,
1152 
        longitude index,
1153 
        latitude index.
1154 
OUTPUT  RETURN_OK if successful, RETURN_ERROR otherwise.
1155 
NOTES   .
1156 
AUTHOR  E. Bertin (IAP)
1157 
VERSION 20/11/2003
1158 
 ***/
1159 
int     reaxe_wcs(wcsstruct *wcs, int lng, int lat)
1160 
 
1161 
  {
1162 
   char         strlng[80], strlat[80];
1163 
   double       dlng,dlat,dlng2,dlat2;
1164 
   int          l, ilng,ilat,olng,olat, naxis;
1165 
 
1166 
  olng = wcs->lng;
1167 
  olat = wcs->lat;
1168 
  if (lng<0 || lat<0 || olng<0 || olat<0)
1169 
    return RETURN_ERROR;
1170 
 
1171 
  ilng = wcs->naxisn[olng];
1172 
  ilat = wcs->naxisn[olat];
1173 
  wcs->naxisn[lng] = ilng;
1174 
  wcs->naxisn[lat] = ilat;
1175 
  strcpy(strlng, wcs->ctype[olng]);
1176 
  strcpy(strlat, wcs->ctype[olat]);
1177 
  strcpy(wcs->ctype[lng], strlng);
1178 
  strcpy(wcs->ctype[lat], strlat);
1179 
  dlng = wcs->crval[olng];
1180 
  dlat = wcs->crval[olat];
1181 
  wcs->crval[lng] = dlng;
1182 
  wcs->crval[lat] = dlat;
1183 
  naxis = wcs->naxis;
1184 
  dlng =  wcs->cd[olng+olng*naxis];
1185 
  dlng2 = wcs->cd[olng+olat*naxis];
1186 
  dlat =  wcs->cd[olat+olat*naxis];
1187 
  dlat2 = wcs->cd[olat+olng*naxis];
1188 
  wcs->cd[lng+lng*naxis] = dlng2;
1189 
  wcs->cd[lng+lat*naxis] = dlng;
1190 
  wcs->cd[lat+lat*naxis] = dlat2;
1191 
  wcs->cd[lat+lng*naxis] = dlat;
1192 
  for (l=0; l<100; l++)
1193 
    {
1194 
    dlng = wcs->projp[l+olng*100];
1195 
    dlat = wcs->projp[l+olat*100];
1196 
    wcs->projp[l+lng*100] = dlng;
1197 
    wcs->projp[l+lat*100] = dlat;
1198 
    }
1199 
 
1200 
/*-- Reinitialize wcs */
1201 
    wcsset(wcs->naxis,(const char(*)[9])wcs->ctype, wcs->wcsprm);
1202 
 
1203 
/*-- Initialize other WCS structures */
1204 
    init_wcs(wcs);
1205 
/*-- Find the range of coordinates */
1206 
    range_wcs(wcs);
1207 
 
1208 
  return RETURN_OK;
1209 
  }
1210 
 
1211 
 
1212 105 bertin 
/******* celsys_to_eq *********************************************************
1213 
PROTO   int celsys_to_eq(wcsstruct *wcs, double *wcspos)
1214 
PURPOSE Convert arbitrary celestial coordinates to equatorial.
1215 
INPUT   WCS structure,
1216 
        Coordinate vector.
1217 
OUTPUT  RETURN_OK if mapping successful, RETURN_ERROR otherwise.
1218 
NOTES   -.
1219 
AUTHOR  E. Bertin (IAP)
1220 
VERSION 08/02/2007
1221 
 ***/
1222 
int     celsys_to_eq(wcsstruct *wcs, double *wcspos)
1223 
 
1224 
  {
1225 
   double       *mat,
1226 
                a2,d2,sd2,cd2cp,sd,x,y;
1227 
   int          lng, lat;
1228 
 
1229 
  mat = wcs->celsysmat;
1230 
  a2 = wcspos[lng = wcs->wcsprm->lng]*DEG - mat[1];
1231 
  d2 = wcspos[lat = wcs->wcsprm->lat]*DEG;
1232 
/* A bit of spherical trigonometry... */
1233 
/* Compute the latitude... */
1234 
  sd2 = sin(d2);
1235 
  cd2cp = cos(d2)*mat[2];
1236 
  sd = sd2*mat[3]-cd2cp*cos(a2);
1237 
/* ...and the longitude */
1238 
  y = cd2cp*sin(a2);
1239 
  x = sd2 - sd*mat[3];
1240 
  wcspos[lng] = fmod((atan2(y,x) + mat[0])/DEG+360.0, 360.0);
1241 
  wcspos[lat] = asin(sd)/DEG;
1242 
 
1243 
  return RETURN_OK;
1244 
  }
1245 
 
1246 
 
1247 
/******* eq_to_celsys *********************************************************
1248 
PROTO   int eq_to_celsys(wcsstruct *wcs, double *wcspos)
1249 
PURPOSE Convert equatorial to arbitrary celestial coordinates.
1250 
INPUT   WCS structure,
1251 
        Coordinate vector.
1252 
OUTPUT  RETURN_OK if mapping successful, RETURN_ERROR otherwise.
1253 
NOTES   -.
1254 
AUTHOR  E. Bertin (IAP)
1255 
VERSION 08/02/2007
1256 
 ***/
1257 
int     eq_to_celsys(wcsstruct *wcs, double *wcspos)
1258 
 
1259 
  {
1260 
   double       *mat,
1261 
                a,d,sd2,cdcp,sd,x,y;
1262 
   int          lng, lat;
1263 
 
1264 
  mat = wcs->celsysmat;
1265 
  a = wcspos[lng = wcs->wcsprm->lng]*DEG - mat[0];
1266 
  d = wcspos[lat = wcs->wcsprm->lat]*DEG;
1267 
/* A bit of spherical trigonometry... */
1268 
/* Compute the latitude... */
1269 
  sd = sin(d);
1270 
  cdcp = cos(d)*mat[2];
1271 
  sd2 = sd*mat[3]+cdcp*cos(a);
1272 
/* ...and the longitude */
1273 
  y = cdcp*sin(a);
1274 
  x = sd2*mat[3]-sd;
1275 
  wcspos[lng] = fmod((atan2(y,x) + mat[1])/DEG+360.0, 360.0);
1276 
  wcspos[lat] = asin(sd2)/DEG;
1277 
 
1278 
  return RETURN_OK;
1279 
  }
1280 
 
1281 
 
1282 95 bertin 
/******* raw_to_wcs ***********************************************************
1283 
PROTO   int raw_to_wcs(wcsstruct *, double *, double *)
1284 
PURPOSE Convert raw (pixel) coordinates to WCS (World Coordinate System).
1285 
INPUT   WCS structure,
1286 
        Pointer to the array of input coordinates,
1287 
        Pointer to the array of output coordinates.
1288 
OUTPUT  RETURN_OK if mapping successful, RETURN_ERROR otherwise.
1289 
NOTES   -.
1290 
AUTHOR  E. Bertin (IAP)
1291 105 bertin 
VERSION 08/02/2007
1292 95 bertin 
 ***/
1293 
int     raw_to_wcs(wcsstruct *wcs, double *pixpos, double *wcspos)
1294 
 
1295 
  {
1296 105 bertin 
   double       imgcrd[NAXIS],
1297 
                phi,theta;
1298 
   int          i;
1299 95 bertin 
 
1300 
  if (wcsrev((const char(*)[9])wcs->ctype, wcs->wcsprm, pixpos,
1301 
        wcs->lin,imgcrd, wcs->prj, &phi, &theta, wcs->crval, wcs->cel, wcspos))
1302 
    {
1303 
    for (i=0; i<wcs->naxis; i++)
1304 
      wcspos[i] = WCS_NOCOORD;
1305 
    return RETURN_ERROR;
1306 
    }
1307 
 
1308 
/* If needed, convert from a different coordinate system to equatorial */
1309 
  if (wcs->celsysconvflag)
1310 105 bertin 
    celsys_to_eq(wcs, wcspos);
1311 95 bertin 
 
1312 
  return RETURN_OK;
1313 
  }
1314 
 
1315 
 
1316 
/******* wcs_to_raw ***********************************************************
1317 
PROTO   int wcs_to_raw(wcsstruct *, double *, double *)
1318 
PURPOSE Convert WCS (World Coordinate System) coords to raw (pixel) coords.
1319 
INPUT   WCS structure,
1320 
        Pointer to the array of input coordinates,
1321 
        Pointer to the array of output coordinates.
1322 
OUTPUT  RETURN_OK if mapping successful, RETURN_ERROR otherwise.
1323 
NOTES   -.
1324 
AUTHOR  E. Bertin (IAP)
1325 105 bertin 
VERSION 08/02/2007
1326 95 bertin 
 ***/
1327 
int     wcs_to_raw(wcsstruct *wcs, double *wcspos, double *pixpos)
1328 
 
1329 
  {
1330 105 bertin 
   double       imgcrd[NAXIS],
1331 
                phi,theta;
1332 
   int          i;
1333 95 bertin 
 
1334 
/* If needed, convert to a coordinate system different from equatorial */
1335 
  if (wcs->celsysconvflag)
1336 105 bertin 
    eq_to_celsys(wcs, wcspos);
1337 95 bertin 
 
1338 
  if (wcsfwd((const char(*)[9])wcs->ctype,wcs->wcsprm,wcspos,
1339 
        wcs->crval, wcs->cel,&phi,&theta,wcs->prj, imgcrd,wcs->lin,pixpos))
1340 
    {
1341 
    for (i=0; i<wcs->naxis; i++)
1342 
      pixpos[i] = WCS_NOCOORD;
1343 
    return RETURN_ERROR;
1344 
    }
1345 
 
1346 
  return RETURN_OK;
1347 
  }
1348 
 
1349 
 
1350 
/******* red_to_raw **********************************************************
1351 
PROTO   int red_to_raw(wcsstruct *, double *, double *)
1352 
PURPOSE Convert reduced (World Coordinate System) coords to raw (pixel)
1353 
        coords.
1354 
INPUT   WCS structure,
1355 
        Pointer to the array of input (reduced) coordinates,
1356 
        Pointer to the array of output (pixel) coordinates.
1357 
OUTPUT  RETURN_OK if mapping successful, RETURN_ERROR otherwise.
1358 
NOTES   -.
1359 
AUTHOR  E. Bertin (IAP)
1360 
VERSION 23/10/2003
1361 
 ***/
1362 
int     red_to_raw(wcsstruct *wcs, double *redpos, double *pixpos)
1363 
 
1364 
  {
1365 
   struct wcsprm        *wcsprm;
1366 
   double               offset;
1367 
 
1368 
  wcsprm = wcs->wcsprm;
1369 
/* Initialize if required */
1370 
  if (wcsprm && wcsprm->flag != WCSSET)
1371 
    {
1372 
    if (wcsset(wcs->naxis, (const char(*)[9])wcs->ctype, wcsprm))
1373 
      return RETURN_ERROR;
1374 
    }
1375 
 
1376 
  if (wcsprm && wcsprm->flag != 999 && wcsprm->cubeface != -1)
1377 
    {
1378 
/*-- Separation between faces */
1379 
    offset = (wcs->prj->r0 == 0.0 ? 90.0 : wcs->prj->r0*PI/2.0);
1380 
/*-- Stack faces in a cube */
1381 
    if (redpos[wcs->lat] < -0.5*offset)
1382 
      {
1383 
      redpos[wcs->lat] += offset;
1384 
      redpos[wcsprm->cubeface] = 5.0;
1385 
      }
1386 
    else if (redpos[wcs->lat] > 0.5*offset)
1387 
      {
1388 
      redpos[wcs->lat] -= offset;
1389 
      redpos[wcsprm->cubeface] = 0.0;
1390 
      }
1391 
    else if (redpos[wcs->lng] > 2.5*offset)
1392 
      {
1393 
      redpos[wcs->lng] -= 3.0*offset;
1394 
      redpos[wcsprm->cubeface] = 4.0;
1395 
      }
1396 
    else if (redpos[wcs->lng] > 1.5*offset)
1397 
      {
1398 
      redpos[wcs->lng] -= 2.0*offset;
1399 
      redpos[wcsprm->cubeface] = 3.0;
1400 
      }
1401 
    else if (redpos[wcs->lng] > 0.5*offset)
1402 
      {
1403 
      redpos[wcs->lng] -= offset;
1404 
      redpos[wcsprm->cubeface] = 2.0;
1405 
      }
1406 
    else
1407 
      redpos[wcsprm->cubeface] = 1.0;
1408 
    }
1409 
 
1410 
/* Apply forward linear transformation */
1411 
  if (linfwd(redpos, wcs->lin, pixpos))
1412 
      return RETURN_ERROR;
1413 
 
1414 
  return RETURN_OK;
1415 
  }
1416 
 
1417 
 
1418 
/******* raw_to_red **********************************************************
1419 
PROTO   int raw_to_red(wcsstruct *, double *, double *)
1420 
PURPOSE Convert raw (pixel) coordinates to reduced WCS coordinates.
1421 
INPUT   WCS structure,
1422 
        Pointer to the array of input (pixel) coordinates,
1423 
        Pointer to the array of output (reduced) coordinates.
1424 
OUTPUT  RETURN_OK if mapping successful, RETURN_ERROR otherwise.
1425 
NOTES   -.
1426 
AUTHOR  E. Bertin (IAP)
1427 
VERSION 23/10/2003
1428 
 ***/
1429 
int     raw_to_red(wcsstruct *wcs, double *pixpos, double *redpos)
1430 
 
1431 
  {
1432 
   struct wcsprm        *wcsprm;
1433 
   double               offset;
1434 
   int                  face;
1435 
 
1436 
  wcsprm = wcs->wcsprm;
1437 
/* Initialize if required */
1438 
  if (wcsprm && wcsprm->flag != WCSSET)
1439 
    {
1440 
    if (wcsset(wcs->naxis, (const char(*)[9])wcs->ctype, wcsprm))
1441 
      return RETURN_ERROR;
1442 
    }
1443 
 
1444 
/* Apply reverse linear transformation */
1445 
   if (linrev(pixpos, wcs->lin, redpos))
1446 
     return RETURN_ERROR;
1447 
 
1448 
  if (wcsprm && wcsprm->flag != 999 && wcsprm->cubeface != -1)
1449 
    {
1450 
/*-- Do we have a CUBEFACE axis? */
1451 
    face = (int)(redpos[wcsprm->cubeface] + 0.5);
1452 
    if (fabs(redpos[wcsprm->cubeface]-face) > 1e-10)
1453 
      return RETURN_ERROR;
1454 
 
1455 
/*-- Separation between faces. */
1456 
    offset = (wcs->prj->r0 == 0.0 ? 90.0 : wcs->prj->r0*PI/2.0);
1457 
/*-- Lay out faces in a plane. */
1458 
    switch (face)
1459 
      {
1460 
      case 0:
1461 
        redpos[wcs->lat] += offset;
1462 
        break;
1463 
      case 1:
1464 
        break;
1465 
      case 2:
1466 
        redpos[wcs->lng] += offset;
1467 
        break;
1468 
      case 3:
1469 
        redpos[wcs->lng] += offset*2;
1470 
        break;
1471 
      case 4:
1472 
        redpos[wcs->lng] += offset*3;
1473 
        break;
1474 
      case 5:
1475 
        redpos[wcs->lat] -= offset;
1476 
        break;
1477 
      default:
1478 
        return RETURN_ERROR;
1479 
      }
1480 
    }
1481 
 
1482 
  return RETURN_OK;
1483 
  }
1484 
 
1485 
 
1486 
/******* wcs_dist ***********************************************************
1487 
PROTO   double wcs_dist(wcsstruct *wcs, double *wcspos1, double *wcspos2)
1488 
PURPOSE Compute the angular distance between 2 points on the sky.
1489 
INPUT   WCS structure,
1490 
        Pointer to the first array of world coordinates,
1491 
        Pointer to the second array of world coordinates.
1492 
OUTPUT  Angular distance (in degrees) between points.
1493 
NOTES   -.
1494 
AUTHOR  E. Bertin (IAP)
1495 
VERSION 24/07/2002
1496 
 ***/
1497 
double  wcs_dist(wcsstruct *wcs, double *wcspos1, double *wcspos2)
1498 
 
1499 
  {
1500 
  double        d, dp;
1501 
  int           i, lng, lat;
1502 
 
1503 
  lng = wcs->lng;
1504 
  lat = wcs->lat;
1505 
  if (lat!=lng)
1506 
    {
1507 
/*-- We are operating in angular coordinates */
1508 
    d = sin(wcspos1[lat]*DEG)*sin(wcspos2[lat]*DEG)
1509 
        + cos(wcspos1[lat]*DEG)*cos(wcspos2[lat]*DEG)
1510 
                *cos((wcspos1[lng]-wcspos2[lng])*DEG);
1511 
    return d>-1.0? (d<1.0 ? acos(d)/DEG : 0.0) : 180.0;
1512 
    }
1513 
  else
1514 
    {
1515 
    d = 0.0;
1516 
    for (i=0; i<wcs->naxis; i++)
1517 
      {
1518 
      dp = wcspos1[i] - wcspos2[i];
1519 
      d += dp*dp;
1520 
      }
1521 
    return sqrt(d);
1522 
    }
1523 
  }
1524 
 
1525 
 
1526 
/******* wcs_scale ***********************************************************
1527 
PROTO   double wcs_scale(wcsstruct *wcs, double *pixpos)
1528 
PURPOSE Compute the sky area equivalent to a local pixel.
1529 
INPUT   WCS structure,
1530 
        Pointer to the array of local raw coordinates,
1531 
OUTPUT  -.
1532 
NOTES   -.
1533 
AUTHOR  E. Bertin (IAP)
1534 106 bertin 
VERSION 03/01/2008
1535 95 bertin 
 ***/
1536 
double  wcs_scale(wcsstruct *wcs, double *pixpos)
1537 
 
1538 
  {
1539 
   double       wcspos[NAXIS], wcspos1[NAXIS], wcspos2[NAXIS], pixpos2[NAXIS];
1540 
   double       dpos1,dpos2;
1541 
   int          lng, lat;
1542 
 
1543 
  if (raw_to_wcs(wcs, pixpos, wcspos))
1544 
    return 0.0;
1545 
 
1546 106 bertin 
  lng = wcs->lng;
1547 
  lat = wcs->lat;
1548 
  if (lng == lat)
1549 
    {
1550 
    lng = 0;
1551 
    lat = 1;
1552 
    }
1553 
 
1554 95 bertin 
/* Compute pixel scale */
1555 
  pixpos2[lng] = pixpos[lng] + 1.0;
1556 
  pixpos2[lat] = pixpos[lat];
1557 
  if (raw_to_wcs(wcs, pixpos2, wcspos1))
1558 
    return 0.0;
1559 
  pixpos2[lng] -= 1.0;
1560 
  pixpos2[lat] += 1.0;
1561 
  if (raw_to_wcs(wcs, pixpos2, wcspos2))
1562 
    return 0.0;
1563 
  dpos1 = wcspos1[lng]-wcspos[lng];
1564 
  dpos2 = wcspos2[lng]-wcspos[lng];
1565 106 bertin 
  if (wcs->lng!=wcs->lat)
1566 
    {
1567 
    if (dpos1>180.0)
1568 
      dpos1 -= 360.0;
1569 
    else if (dpos1<-180.0)
1570 
      dpos1 += 360.0;
1571 
    if (dpos2>180.0)
1572 
      dpos2 -= 360.0;
1573 
    else if (dpos2<-180.0)
1574 
      dpos2 += 360.0;
1575 
    return fabs((dpos1*(wcspos2[lat]-wcspos[lat])
1576 
                -(wcspos1[lat]-wcspos[lat])*dpos2)*cos(wcspos[lat]*DEG));
1577 
    }
1578 
  else
1579 
    return fabs((dpos1*(wcspos2[lat]-wcspos[lat])
1580 
                -(wcspos1[lat]-wcspos[lat])*dpos2));
1581 95 bertin 
  }
1582 
 
1583 
 
1584 
/****** wcs jacobian *********************************************************
1585 
PROTO   double wcs_jacobian(wcsstruct *wcs, double *pixpos, double *jacob)
1586 284 bertin 
PURPOSE Compute the local Jacobian matrix of the astrometric deprojection.
1587 95 bertin 
INPUT   WCS structure,
1588 
        Pointer to the array of local raw coordinates,
1589 
        Pointer to the jacobian array (output).
1590 
OUTPUT  Determinant over spatial coordinates (=pixel area), or -1.0 if mapping
1591 
        was unsuccesful.
1592 
NOTES   Memory must have been allocated (naxis*naxis*sizeof(double)) for the
1593 
        Jacobian array.
1594 
AUTHOR  E. Bertin (IAP)
1595 
VERSION 11/10/2007
1596 
 ***/
1597 
double  wcs_jacobian(wcsstruct *wcs, double *pixpos, double *jacob)
1598 
  {
1599 
   double       pixpos0[NAXIS], wcspos0[NAXIS], wcspos[NAXIS],
1600 
                dpos;
1601 
   int          i,j, lng,lat,naxis;
1602 
 
1603 
  lng = wcs->lng;
1604 
  lat = wcs->lat;
1605 
  naxis = wcs->naxis;
1606 
  for (i=0; i<naxis; i++)
1607 
    pixpos0[i] = pixpos[i];
1608 
  if (raw_to_wcs(wcs, pixpos0, wcspos0) == RETURN_ERROR)
1609 
    return -1.0;
1610 
  for (i=0; i<naxis; i++)
1611 
    {
1612 
    pixpos0[i] += 1.0;
1613 
    if (raw_to_wcs(wcs, pixpos0, wcspos) == RETURN_ERROR)
1614 
      return -1.0;
1615 
    pixpos0[i] -= 1.0;
1616 
    for (j=0; j<naxis; j++)
1617 
      {
1618 
      dpos = wcspos[j]-wcspos0[j];
1619 
      if (lng!=lat && j==lng)
1620 
        {
1621 
        if (dpos>180.0)
1622 
          dpos -= 360.0;
1623 
        else if (dpos<-180.0)
1624 
          dpos += 360.0;
1625 
        dpos *= cos(wcspos0[lat]*DEG);
1626 
        }
1627 
      jacob[j*naxis+i] = dpos;
1628 
      }
1629 
    }
1630 
 
1631 
  if (lng==lat)
1632 
    {
1633 
    lng = 0;
1634 
    lat = 1;
1635 
    }
1636 
 
1637 
  return fabs(jacob[lng+naxis*lng]*jacob[lat+naxis*lat]
1638 
                - jacob[lat+naxis*lng]*jacob[lng+naxis*lat]);
1639 
  }
1640 
 
1641 
 
1642 284 bertin 
/****** wcs rawtoraw *********************************************************
1643 
PROTO   double wcs_rawtoraw(wcsstruct *wcsin, wcsstruct *wcsout,
1644 
                double *pixpos, double *jacob)
1645 
PURPOSE Convert raw coordinates from input wcs structure to raw coordinates
1646 
        from output wcs structure, and the local Jacobian matrix of the
1647 
        reprojection.
1648 
INPUT   WCS input structure,
1649 
        WCS output structure,
1650 
        pointer to the array of local input raw coordinates,
1651 
        pointer to the array of local output raw coordinates (output),
1652 
        pointer to the jacobian array (output).
1653 
OUTPUT  Determinant over spatial coordinates (ratio of pixel areas),
1654 
        0.0 if jacob is NULL (Jacobian not computed in that case),
1655 
        or -1.0 if mapping was unsuccesful.
1656 
NOTES   Memory must have been allocated (naxis*naxis*sizeof(double)) for the
1657 
        Jacobian array.
1658 
AUTHOR  E. Bertin (IAP)
1659 293 bertin 
VERSION 12/06/2012
1660 284 bertin 
 ***/
1661 
double  wcs_rawtoraw(wcsstruct *wcsin, wcsstruct *wcsout,
1662 
                double *pixposin, double *pixposout, double *jacob)
1663 
  {
1664 293 bertin 
   double       pixpos0[NAXIS], pixpos2[NAXIS], wcspos[NAXIS];
1665 284 bertin 
   int          i,j, lng,lat,naxis;
1666 
 
1667 
  naxis = wcsin->naxis;
1668 
  for (i=0; i<naxis; i++)
1669 
    pixpos0[i] = pixposin[i];
1670 
 
1671 
/* Coordinate transformation */
1672 
  if (raw_to_wcs(wcsin, pixposin, wcspos) == RETURN_ERROR)
1673 
    return -1.0;
1674 
  if (wcs_to_raw(wcsout, wcspos, pixposout) == RETURN_ERROR)
1675 
    return -1.0;
1676 
 
1677 
/* Jacobian */
1678 
  if (!jacob)
1679 
    return 0.0;
1680 
 
1681 
  lng = wcsin->lng;
1682 
  lat = wcsin->lat;
1683 
  for (i=0; i<naxis; i++)
1684 
    {
1685 
    pixpos0[i] += 1.0;
1686 
    if (raw_to_wcs(wcsin, pixpos0, wcspos) == RETURN_ERROR)
1687 
      return -1.0;
1688 
    if (wcs_to_raw(wcsout, wcspos, pixpos2) == RETURN_ERROR)
1689 
      return -1.0;
1690 
    pixpos0[i] -= 1.0;
1691 
    for (j=0; j<naxis; j++)
1692 
      jacob[j*naxis+i] = pixpos2[j] - pixposout[j];
1693 
    }
1694 
 
1695 
  if (lng==lat)
1696 
    {
1697 
    lng = 0;
1698 
    lat = 1;
1699 
    }
1700 
 
1701 
  return fabs(jacob[lng+naxis*lng]*jacob[lat+naxis*lat]
1702 
                - jacob[lat+naxis*lng]*jacob[lng+naxis*lat]);
1703 
  }
1704 
 
1705 
 
1706 95 bertin 
/******* wcs_chirality *******************************************************
1707 
PROTO   int wcs_chirality(wcsstruct *wcs)
1708 
PURPOSE Compute the chirality of a WCS projection.
1709 
INPUT   WCS structure.
1710 
OUTPUT  +1 if determinant of matrix is positive, -1 if negative, 0 if null.
1711 
NOTES   -.
1712 
AUTHOR  E. Bertin (IAP)
1713 
VERSION 26/09/2006
1714 
 ***/
1715 
int     wcs_chirality(wcsstruct *wcs)
1716 
 
1717 
  {
1718 
   double       a;
1719 
   int          lng,lat, naxis;
1720 
 
1721 
  lng = wcs->lng;
1722 
  lat = wcs->lat;
1723 
  naxis = wcs->naxis;
1724 
  if (lng==lat && naxis>=2)
1725 
    {
1726 
    lng = 0;
1727 
    lat = 1;
1728 
    }
1729 
 
1730 
  a = wcs->cd[lng*naxis+lng]*wcs->cd[lat*naxis+lat]
1731 
        - wcs->cd[lng*naxis+lat]*wcs->cd[lat*naxis+lng];
1732 
  return a>TINY? 1 : (a<-TINY? -1 : 0);
1733 
  }
1734 
 
1735 
 
1736 
/****** precess_wcs **********************************************************
1737 
PROTO   void precess_wcs(wcsstruct *wcs, double yearin, double yearout)
1738 
PURPOSE Precess the content of a WCS structure according to the equinox.
1739 
INPUT   WCS structure,
1740 
        Input year,
1741 
        Output year.
1742 
OUTPUT  -.
1743 
NOTES   Epoch for coordinates should be J2000 (FK5 system).
1744 
AUTHOR  E. Bertin (IAP)
1745 106 bertin 
VERSION 04/01/2008
1746 95 bertin 
 ***/
1747 
void    precess_wcs(wcsstruct *wcs, double yearin, double yearout)
1748 
 
1749 
  {
1750 
   double       crval[NAXIS],a[NAXIS*NAXIS],b[NAXIS*NAXIS],
1751 
                *c,*at,
1752 
                val, cas, sas, angle, dalpha;
1753 
   int          i,j,k, lng,lat, naxis;
1754 
 
1755 
  lng = wcs->lng;
1756 
  lat = wcs->lat;
1757 
  if (lat==lng || yearin==yearout)
1758 
    return;
1759 
  naxis = wcs->naxis;
1760 
/* Precess to year out */
1761 
  precess(yearin, wcs->crval[lng], wcs->crval[lat], yearout,
1762 
        &crval[lng], &crval[lat]);
1763 
 
1764 
  dalpha = (crval[lng] - wcs->crval[lng])*DEG;
1765 
 
1766 
/* Compute difference angle with the north axis between start and end */
1767 
  angle = (dalpha!=0.0 && (crval[lat] - wcs->crval[lat])*DEG != 0.0) ?
1768 
        180.0 - (atan2(sin(dalpha),
1769 
        cos(crval[lat]*DEG)*tan(wcs->crval[lat]*DEG)
1770 
        - sin(crval[lat]*DEG)*cos(dalpha))
1771 
        + atan2(sin(dalpha),
1772 
        cos(wcs->crval[lat]*DEG)*tan(crval[lat]*DEG)
1773 
        - sin(wcs->crval[lat]*DEG)*cos(dalpha)))/DEG
1774 
        : 0.0;
1775 
 
1776 106 bertin 
/* A = C*B */
1777 95 bertin 
  c = wcs->cd;
1778 
/* The B matrix is made of 2 numbers */
1779 
 
1780 
  cas = cos(angle*DEG);
1781 105 bertin 
  sas = sin(-angle*DEG);
1782 95 bertin 
  for (i=0; i<naxis; i++)
1783 
    b[i+i*naxis] = 1.0;
1784 
  b[lng+lng*naxis] = cas;
1785 
  b[lat+lng*naxis] = -sas;
1786 
  b[lng+lat*naxis] = sas;
1787 
  b[lat+lat*naxis] = cas;
1788 
  at = a;
1789 
  for (j=0; j<naxis; j++)
1790 
    for (i=0; i<naxis; i++)
1791 
      {
1792 
      val = 0.0;
1793 
      for (k=0; k<naxis; k++)
1794 106 bertin 
        val += c[k+j*naxis]*b[i+k*naxis];
1795 95 bertin 
      *(at++) = val;
1796 
      }
1797 
 
1798 
  at = a;
1799 
 
1800 
  for (i=0; i<naxis*naxis; i++)
1801 
    *(c++) = *(at++);
1802 
 
1803 
  wcs->crval[lng] = crval[lng];
1804 
  wcs->crval[lat] = crval[lat];
1805 
  wcs->equinox = yearout;
1806 
 
1807 
/* Initialize other WCS structures */
1808 
  init_wcs(wcs);
1809 
/* Find the range of coordinates */
1810 
  range_wcs(wcs);
1811 
/* Invert projection corrections */
1812 
  invert_wcs(wcs);
1813 
 
1814 
  return;
1815 
  }
1816 
 
1817 
 
1818 
/********************************* precess ***********************************/
1819 
/*
1820 
precess equatorial coordinates according to the equinox (from Ephemerides du
1821 
Bureau des Longitudes 1992). Epoch for coordinates should be J2000
1822 
(FK5 system).
1823 
*/
1824 
void    precess(double yearin, double alphain, double deltain,
1825 
                double yearout, double *alphaout, double *deltaout)
1826 
 
1827 
  {
1828 
   double       dzeta,theta,z, t1,t1t1, t2,t2t2,t2t2t2,
1829 
                cddsadz, cddcadz, cdd, sdd, adz, cdin,sdin,ct,st,caindz;
1830 
 
1831 
  alphain *= DEG;
1832 
  deltain *= DEG;
1833 
 
1834 
  t1 = (yearin - 2000.0)/1000.0;
1835 
  t2 = (yearout - yearin)/1000.0;
1836 
  t1t1 = t1*t1;
1837 
  t2t2t2 = (t2t2 = t2*t2)*t2;
1838 
  theta = (97171.735e-06 - 413.691e-06*t1 - 1.052e-06 * t1t1) * t2
1839 
        + (-206.846e-06 - 1.052e-06*t1) * t2t2 - 202.812e-06 * t2t2t2;
1840 
  dzeta = (111808.609e-06 + 677.071e-06*t1 - 0.674e-06 * t1t1) * t2
1841 
        + (146.356e-06 - 1.673e-06*t1) * t2t2 + 87.257e-06 * t2t2t2;
1842 
  z = (111808.609e-06 +677.071e-06*t1 - 0.674e-06 * t1t1) * t2
1843 
        + (530.716e-06 + 0.320e-06*t1) * t2t2 + 88.251e-06 * t2t2t2;
1844 
  cddsadz = (cdin=cos(deltain)) * sin(alphain+dzeta);
1845 
  cddcadz = -(sdin=sin(deltain))*(st=sin(theta))
1846 
        +cdin*(ct=cos(theta))*(caindz=cos(alphain+dzeta));
1847 
  sdd = sdin*ct + cdin*st*caindz;
1848 
  cdd = cos(*deltaout = asin(sdd));
1849 
  adz = asin(cddsadz/cdd);
1850 
  if (cddcadz<0.0)
1851 
    adz = PI - adz;
1852 
  if (adz<0.0)
1853 
    adz += 2.0*PI;
1854 
  adz += z;
1855 
  *alphaout = adz/DEG;
1856 
  *deltaout /= DEG;
1857 
 
1858 
  return;
1859 
  }
1860 
 
1861 
 
1862 105 bertin 
/********************************* b2j ***********************************/
1863 
/*
1864 
conver equatorial coordinates from equinox and epoch B1950 to equinox and
1865 
epoch J2000 for extragalactic sources (from Aoki et al. 1983).
1866 
*/
1867 
void    b2j(double yearobs, double alphain, double deltain,
1868 
                double *alphaout, double *deltaout)
1869 
  {
1870 
   int          i,j;
1871 
   double       a[3] = {-1.62557e-6, -0.31919e-6, -0.13843e-6},
1872 
                ap[3] = {1.245e-3, -1.580e-3, -0.659e-3},
1873 
                m[6][6] = {
1874 
  { 0.9999256782,     -0.0111820611,     -0.0048579477,
1875 
    0.00000242395018, -0.00000002710663, -0.00000001177656},
1876 
  { 0.0111820610,      0.9999374784,     -0.0000271765,
1877 
    0.00000002710663,  0.00000242397878, -0.00000000006587},
1878 
  { 0.0048579479,     -0.0000271474,      0.9999881997,
1879 
    0.00000001177656, -0.00000000006582,  0.00000242410173},
1880 
  {-0.000551,        -0.238565,           0.435739,
1881 
    0.99994704,      -0.01118251,        -0.00485767},
1882 
  { 0.238514,        -0.002662,          -0.008541,
1883 
    0.01118251,       0.99995883,        -0.00002718},
1884 
  {-0.435623,         0.012254,           0.002117,
1885 
    0.00485767,      -0.00002714,         1.00000956}},
1886 
                a1[3], r[3], ro[3], r1[3], r2[3], v1[3], v[3];
1887 
   double               cai, sai, cdi, sdi, dotp, rmod, alpha, delta,
1888 
                        t1 = (yearobs - 1950.0)/100.0;
1889 
 
1890 
  alphain *= PI/180.0;
1891 
  deltain *= PI/180.0;
1892 
  cai = cos(alphain);
1893 
  sai = sin(alphain);
1894 
  cdi = cos(deltain);
1895 
  sdi = sin(deltain);
1896 
  ro[0] = cdi*cai;
1897 
  ro[1] = cdi*sai;
1898 
  ro[2] = sdi;
1899 
  dotp = 0.0;
1900 
  for (i=0; i<3; i++)
1901 
    {
1902 
    a1[i] = a[i]+ap[i]*ARCSEC*t1;
1903 
    dotp += a1[i]*ro[i];
1904 
    }
1905 
  for (i=0; i<3; i++)
1906 
    {
1907 
    r1[i] = ro[i] - a1[i] + dotp*ro[i];
1908 
    r[i] = v[i] = v1[i] = 0.0;
1909 
    }
1910 
  for (j=0; j<6; j++)
1911 
    for (i=0; i<6; i++)
1912 
      {
1913 
      if (j<3)
1914 
        r[j] += m[j][i]*(i<3?r1[i]:v1[i-3]);
1915 
      else
1916 
         v[j-3] += m[j][i]*(i<3?r1[i]:v1[i-3]);
1917 
      }
1918 
  rmod = 0.0;
1919 
  for (i=0; i<3; i++)
1920 
    {
1921 
    r2[i] = r[i]+v[i]*ARCSEC*(t1-0.5);
1922 
    rmod += r2[i]*r2[i];
1923 
    }
1924 
  rmod = sqrt(rmod);
1925 
  delta = asin(r2[2]/rmod);
1926 
  alpha = acos(r2[0]/cos(delta)/rmod);
1927 
  if (r2[1]<0)
1928 
    alpha = 2*PI - alpha;
1929 
  *alphaout = alpha*180.0/PI;
1930 
  *deltaout = delta*180.0/PI;
1931 
 
1932 
  return;
1933 
  }
1934 
 
1935 
 
1936 95 bertin 
/*********************************** j2b *************************************/
1937 
/*
1938 
conver equatorial coordinates from equinox and epoch J2000 to equinox and
1939 
epoch B1950 for extragalactic sources (from Aoki et al. 1983, after
1940 
inversion of their matrix and some custom arrangements).
1941 
*/
1942 
void    j2b(double yearobs, double alphain, double deltain,
1943 
        double *alphaout, double *deltaout)
1944 
  {
1945 
   int                  i,j;
1946 
   double               a[3] = {-1.62557e-6, -0.31919e-6, -0.13843e-6},
1947 
                        ap[3] = {1.245e-3, -1.580e-3, -0.659e-3},
1948 
                        m[6][6] = {
1949 
  { 0.9999256794678425,    0.01118148281196562,   0.004859003848996022,
1950 
   -2.423898417033081e-06,-2.710547600126671e-08,-1.177738063266745e-08},
1951 
  {-0.01118148272969232,   0.9999374849247641,   -2.717708936468247e-05,
1952 
    2.710547578707874e-08,-2.423927042585208e-06, 6.588254898401055e-11},
1953 
  {-0.00485900399622881,  -2.715579322970546e-05, 0.999988194643078,
1954 
    1.177738102358923e-08, 6.582788892816657e-11,-2.424049920613325e-06},
1955 
  {-0.0005508458576414713, 0.2384844384742432,   -0.4356144527773499,
1956 
    0.9999043171308133,    0.01118145410120206,   0.004858518651645554},
1957 
  {-0.2385354433560954,   -0.002664266996872802,  0.01225282765749546,
1958 
   -0.01118145417187502,   0.9999161290795875,   -2.717034576263522e-05},
1959 
  { 0.4357269351676567,   -0.008536768476441086,  0.002113420799663768,
1960 
   -0.004858518477064975, -2.715994547222661e-05, 0.9999668385070383}},
1961 
                        a1[3], r[3], ro[3], r1[3], r2[3], v1[3], v[3];
1962 
   double               cai, sai, cdi, sdi, dotp, rmod, alpha, delta, t1;
1963 
 
1964 
/* Convert Julian years from J2000.0 to tropic centuries from B1950.0 */
1965 
  t1 = ((yearobs - 2000.0) + (MJD2000 - MJD1950)/365.25)*JU2TROP/100.0;
1966 
  alphain *= DEG;
1967 
  deltain *= DEG;
1968 
  cai = cos(alphain);
1969 
  sai = sin(alphain);
1970 
  cdi = cos(deltain);
1971 
  sdi = sin(deltain);
1972 
  r[0] = cdi*cai;
1973 
  r[1] = cdi*sai;
1974 
  r[2] = sdi;
1975 
  for (i=0; i<3; i++)
1976 
    v[i] = r2[i] = v1[i] = 0.0;
1977 
  for (j=0; j<6; j++)
1978 
    for (i=0; i<6; i++)
1979 
      if (j<3)
1980 
        r2[j] += m[j][i]*(i<3?r[i]:v[i-3]);
1981 
      else
1982 
        v1[j-3] += m[j][i]*(i<3?r[i]:v[i-3]);
1983 
 
1984 
  for (i=0; i<3; i++)
1985 
    r1[i] = r2[i]+v1[i]*ARCSEC*t1;
1986 
 
1987 
  dotp = 0.0;
1988 
  for (i=0; i<3; i++)
1989 
    {
1990 
    a1[i] = a[i]+ap[i]*ARCSEC*t1;
1991 
    dotp += a1[i]*(r1[i]+a1[i]);
1992 
    }
1993 
  dotp = 2.0/(sqrt(1+4.0*dotp)+1.0);
1994 
  rmod = 0.0;
1995 
  for (i=0; i<3; i++)
1996 
    {
1997 
    ro[i] = dotp*(r1[i]+a1[i]);
1998 
    rmod += ro[i]*ro[i];
1999 
    }
2000 
  rmod = sqrt(rmod);
2001 
  delta = asin(ro[2]/rmod);
2002 
  alpha = acos(ro[0]/cos(delta)/rmod);
2003 
  if (ro[1]<0)
2004 
    alpha = 2.0*PI - alpha;
2005 
  *alphaout = alpha/DEG;
2006 
  *deltaout = delta/DEG;
2007 
 
2008 
  return;
2009 
  }
2010 
 
2011 
 
2012 
/******************************** degtosexal *********************************/
2013 
/*
2014 
Convert degrees to hh mm ss.xx alpha coordinates.
2015 
*/
2016 
char    *degtosexal(double alpha, char *str)
2017 
 
2018 
  {
2019 
   int          hh, mm;
2020 
   double       ss;
2021 
 
2022 
  if (alpha>=0.0 && alpha <360.0)
2023 
    {
2024 
    hh = (int)(alpha/15.0);
2025 
    mm = (int)(60.0*(alpha/15.0 - hh));
2026 
    ss = 60.0*(60.0*(alpha/15.0 - hh) - mm);
2027 
    }
2028 
  else
2029 
    hh = mm = ss = 0.0;
2030 
  sprintf(str,"%02d:%02d:%05.2f", hh, mm, ss);
2031 
 
2032 
  return str;
2033 
  }
2034 
 
2035 
 
2036 
/******************************** degtosexde *********************************/
2037 
/*
2038 
Convert degrees to dd dm ds.x delta coordinates.
2039 
*/
2040 
char    *degtosexde(double delta, char *str)
2041 
 
2042 
  {
2043 
   char         sign;
2044 
   double       ds;
2045 
   int          dd, dm;
2046 
 
2047 
  sign = delta<0.0?'-':'+';
2048 
  delta = fabs(delta);
2049 
  if (delta>=-90.0 && delta <=90.0)
2050 
    {
2051 
    dd = (int)delta;
2052 
    dm = (int)(60.0*(delta - dd));
2053 
    ds = 60.0*fabs(60.0*(delta - dd) - dm);
2054 
    }
2055 
  else
2056 
    dd = dm = ds = 0.0;
2057 
  sprintf(str,"%c%02d:%02d:%04.1f", sign, dd, dm, ds);
2058 
  return str;
2059 
  }
2060 
 
2061 
 
2062 
/******************************** sextodegal *********************************/
2063 
/*
2064 
Convert hh mm ss.xxx alpha coordinates to degrees.
2065 
*/
2066 
double  sextodegal(char *hms)
2067 
  {
2068 
   double       val;
2069 
   char         *ptr;
2070 
 
2071 
  val = atof(strtok_r(hms, ": \t", &ptr))*15.0;         /* Hours */
2072 
  val += atof(strtok_r(NULL, ": \t", &ptr))/4.0;        /* Minutes */
2073 
  val += atof(strtok_r(NULL, ": \t", &ptr))/240.0;      /* Seconds */
2074 
 
2075 
  return val;
2076 
  }
2077 
 
2078 
 
2079 
/******************************** sextodegde *********************************/
2080 
/*
2081 
Convert dd dm ds.xxx delta coordinates to degrees.
2082 
*/
2083 
double  sextodegde(char *dms)
2084 
  {
2085 
   double       val, sgn;
2086 
   char         *str, *ptr;
2087 
 
2088 
  str = strtok_r(dms, ": \t", &ptr);
2089 
  sgn = (strchr(str, '-') ? -1.0:1.0);
2090 
  val = atof(dms);                                      /* Degrees */
2091 
  val += atof(strtok_r(NULL, ": \t", &ptr))*sgn/60.0;   /* Minutes */
2092 
  val += atof(strtok_r(NULL, ": \t", &ptr))*sgn/3600.0; /* Seconds */
2093 
 
2094 
  return val;
2095 
  }
2096 
 
2097 
 
2098 123 bertin 
/******************************** fmod_0_p360 *******************************/
2099 
/*
2100 
Fold input angle in the [0,+360[ domain.
2101 
*/
2102 
double  fmod_0_p360(double angle)
2103 
  {
2104 
  return angle>0.0? fmod(angle,360.0) : fmod(angle,360.0)+360.0;
2105 
  }
2106 95 bertin 
 
2107 123 bertin 
 
2108 
/******************************** fmod_m90_p90 *******************************/
2109 
/*
2110 
Fold input angle in the [-90,+90[ domain.
2111 
*/
2112 
double  fmod_m90_p90(double angle)
2113 
  {
2114 
  return angle>0.0? fmod(angle+90.0,180.0)-90.0 : fmod(angle-90.0,180.0)+90.0;
2115 
  }
2116 
 
2117 
 
2118 233 bertin 
/********************************* fcmp_0_p360 *******************************/
2119 
/*
2120 
Compare angles in the [0,+360[ domain: return 1 if anglep>anglem, 0 otherwise.
2121 
*/
2122 
int  fcmp_0_p360(double anglep, double anglem)
2123 
  {
2124 
   double dval = anglep - anglem;
2125 
 
2126 
  return (int)((dval>0.0 && dval<180.0) || dval<-180.0);
2127 
  }
2128 
 
2129