"Fossies" - the Fresh Open Source Software Archive

Member "cfitsio-4.0.0/wcssub.c" (20 May 2021, 35267 Bytes) of package /linux/misc/cfitsio-4.0.0.tar.gz:


As a special service "Fossies" has tried to format the requested source page into HTML format using (guessed) C and C++ source code syntax highlighting (style: standard) with prefixed line numbers and code folding option. Alternatively you can here view or download the uninterpreted source code file. For more information about "wcssub.c" see the Fossies "Dox" file reference documentation.

    1 #include <stdlib.h>
    2 #include <math.h>
    3 #include <string.h>
    4 #include "fitsio2.h"
    5 
    6 /*--------------------------------------------------------------------------*/
    7 int fits_read_wcstab(
    8    fitsfile   *fptr, /* I - FITS file pointer           */
    9    int  nwtb,        /* Number of arrays to be read from the binary table(s) */
   10    wtbarr *wtb,      /* Address of the first element of an array of wtbarr
   11                          typedefs.  This wtbarr typedef is defined below to
   12                          match the wtbarr struct defined in WCSLIB.  An array
   13                          of such structs returned by the WCSLIB function
   14                          wcstab(). */
   15    int  *status)
   16 
   17 /*
   18 *   Author: Mark Calabretta, Australia Telescope National Facility
   19 *   http://www.atnf.csiro.au/~mcalabre/index.html
   20 *
   21 *   fits_read_wcstab() extracts arrays from a binary table required in
   22 *   constructing -TAB coordinates.  This helper routine is intended for
   23 *   use by routines in the WCSLIB library when dealing with the -TAB table
   24 *   look up WCS convention.
   25 */
   26 
   27 {
   28    int  anynul, colnum, hdunum, iwtb, m, naxis, nostat;
   29    long *naxes = 0, nelem;
   30    wtbarr *wtbp;
   31 
   32 
   33    if (*status) return *status;
   34 
   35    if (fptr == 0) {
   36       return (*status = NULL_INPUT_PTR);
   37    }
   38 
   39    if (nwtb == 0) return 0;
   40 
   41    /* Zero the array pointers. */
   42    wtbp = wtb;
   43    for (iwtb = 0; iwtb < nwtb; iwtb++, wtbp++) {
   44      *wtbp->arrayp = 0x0;
   45    }
   46 
   47    /* Save HDU number so that we can move back to it later. */
   48    fits_get_hdu_num(fptr, &hdunum);
   49 
   50    wtbp = wtb;
   51    for (iwtb = 0; iwtb < nwtb; iwtb++, wtbp++) {
   52       /* Move to the required binary table extension. */
   53       if (fits_movnam_hdu(fptr, BINARY_TBL, (char *)(wtbp->extnam),
   54           wtbp->extver, status)) {
   55          goto cleanup;
   56       }
   57 
   58       /* Locate the table column. */
   59       if (fits_get_colnum(fptr, CASEINSEN, (char *)(wtbp->ttype), &colnum,
   60           status)) {
   61          goto cleanup;
   62       }
   63 
   64       /* Get the array dimensions and check for consistency. */
   65       if (wtbp->ndim < 1) {
   66          *status = NEG_AXIS;
   67          goto cleanup;
   68       }
   69 
   70       if (!(naxes = calloc(wtbp->ndim, sizeof(long)))) {
   71          *status = MEMORY_ALLOCATION;
   72          goto cleanup;
   73       }
   74 
   75       if (fits_read_tdim(fptr, colnum, wtbp->ndim, &naxis, naxes, status)) {
   76          goto cleanup;
   77       }
   78 
   79       if (naxis != wtbp->ndim) {
   80          if (wtbp->kind == 'c' && wtbp->ndim == 2) {
   81             /* Allow TDIMn to be omitted for degenerate coordinate arrays. */
   82             naxis = 2;
   83             naxes[1] = naxes[0];
   84             naxes[0] = 1;
   85          } else {
   86             *status = BAD_TDIM;
   87             goto cleanup;
   88          }
   89       }
   90 
   91       if (wtbp->kind == 'c') {
   92          /* Coordinate array; calculate the array size. */
   93          nelem = naxes[0];
   94          for (m = 0; m < naxis-1; m++) {
   95             *(wtbp->dimlen + m) = naxes[m+1];
   96             nelem *= naxes[m+1];
   97          }
   98       } else {
   99          /* Index vector; check length. */
  100          if ((nelem = naxes[0]) != *(wtbp->dimlen)) {
  101             /* N.B. coordinate array precedes the index vectors. */
  102             *status = BAD_TDIM;
  103             goto cleanup;
  104          }
  105       }
  106 
  107       free(naxes);
  108       naxes = 0;
  109 
  110       /* Allocate memory for the array. */
  111       if (!(*wtbp->arrayp = calloc((size_t)nelem, sizeof(double)))) {
  112          *status = MEMORY_ALLOCATION;
  113          goto cleanup;
  114       }
  115 
  116       /* Read the array from the table. */
  117       if (fits_read_col_dbl(fptr, colnum, wtbp->row, 1L, nelem, 0.0,
  118           *wtbp->arrayp, &anynul, status)) {
  119          goto cleanup;
  120       }
  121    }
  122 
  123 cleanup:
  124    /* Move back to the starting HDU. */
  125    nostat = 0;
  126    fits_movabs_hdu(fptr, hdunum, 0, &nostat);
  127 
  128    /* Release allocated memory. */
  129    if (naxes) free(naxes);
  130    if (*status) {
  131       wtbp = wtb;
  132       for (iwtb = 0; iwtb < nwtb; iwtb++, wtbp++) {
  133          if (*wtbp->arrayp) free(*wtbp->arrayp);
  134       }
  135    }
  136 
  137    return *status;
  138 }
  139 /*--------------------------------------------------------------------------*/
  140 int ffgiwcs(fitsfile *fptr,  /* I - FITS file pointer                    */
  141            char **header,   /* O - pointer to the WCS related keywords  */
  142            int *status)     /* IO - error status                        */
  143 /*
  144   int fits_get_image_wcs_keys 
  145   return a string containing all the image WCS header keywords.
  146   This string is then used as input to the wcsinit WCSlib routine.
  147   
  148   THIS ROUTINE IS DEPRECATED. USE fits_hdr2str INSTEAD
  149 */
  150 {
  151     int hdutype;
  152 
  153     if (*status > 0)
  154         return(*status);
  155 
  156     fits_get_hdu_type(fptr, &hdutype, status);
  157     if (hdutype != IMAGE_HDU)
  158     {
  159       ffpmsg(
  160      "Error in ffgiwcs. This HDU is not an image. Can't read WCS keywords");
  161       return(*status = NOT_IMAGE);
  162     }
  163 
  164     /* read header keywords into a long string of chars */
  165     if (ffh2st(fptr, header, status) > 0)
  166     {
  167         ffpmsg("error creating string of image WCS keywords (ffgiwcs)");
  168         return(*status);
  169     }
  170 
  171     return(*status);
  172 }
  173 
  174 /*--------------------------------------------------------------------------*/
  175 int ffgics(fitsfile *fptr,    /* I - FITS file pointer           */
  176            double *xrval,     /* O - X reference value           */
  177            double *yrval,     /* O - Y reference value           */
  178            double *xrpix,     /* O - X reference pixel           */
  179            double *yrpix,     /* O - Y reference pixel           */
  180            double *xinc,      /* O - X increment per pixel       */
  181            double *yinc,      /* O - Y increment per pixel       */
  182            double *rot,       /* O - rotation angle (degrees)    */
  183            char *type,        /* O - type of projection ('-tan') */
  184            int *status)       /* IO - error status               */
  185 /*
  186        read the values of the celestial coordinate system keywords.
  187        These values may be used as input to the subroutines that
  188        calculate celestial coordinates. (ffxypx, ffwldp)
  189 
  190        Modified in Nov 1999 to convert the CD matrix keywords back
  191        to the old CDELTn form, and to swap the axes if the dec-like
  192        axis is given first, and to assume default values if any of the
  193        keywords are not present.
  194 */
  195 {
  196     int tstat = 0, cd_exists = 0, pc_exists = 0;
  197     char ctype[FLEN_VALUE];
  198     double cd11 = 0.0, cd21 = 0.0, cd22 = 0.0, cd12 = 0.0;
  199     double pc11 = 1.0, pc21 = 0.0, pc22 = 1.0, pc12 = 0.0;
  200     double pi =  3.1415926535897932;
  201     double phia, phib, temp;
  202     double toler = .0002;  /* tolerance for angles to agree (radians) */
  203                            /*   (= approximately 0.01 degrees) */
  204 
  205     if (*status > 0)
  206        return(*status);
  207 
  208     tstat = 0;
  209     if (ffgkyd(fptr, "CRVAL1", xrval, NULL, &tstat))
  210        *xrval = 0.;
  211 
  212     tstat = 0;
  213     if (ffgkyd(fptr, "CRVAL2", yrval, NULL, &tstat))
  214        *yrval = 0.;
  215 
  216     tstat = 0;
  217     if (ffgkyd(fptr, "CRPIX1", xrpix, NULL, &tstat))
  218         *xrpix = 0.;
  219 
  220     tstat = 0;
  221     if (ffgkyd(fptr, "CRPIX2", yrpix, NULL, &tstat))
  222         *yrpix = 0.;
  223 
  224     /* look for CDELTn first, then CDi_j keywords */
  225     tstat = 0;
  226     if (ffgkyd(fptr, "CDELT1", xinc, NULL, &tstat))
  227     {
  228         /* CASE 1: no CDELTn keyword, so look for the CD matrix */
  229         tstat = 0;
  230         if (ffgkyd(fptr, "CD1_1", &cd11, NULL, &tstat))
  231             tstat = 0;  /* reset keyword not found error */
  232         else
  233             cd_exists = 1;  /* found at least 1 CD_ keyword */
  234 
  235         if (ffgkyd(fptr, "CD2_1", &cd21, NULL, &tstat))
  236             tstat = 0;  /* reset keyword not found error */
  237         else
  238             cd_exists = 1;  /* found at least 1 CD_ keyword */
  239 
  240         if (ffgkyd(fptr, "CD1_2", &cd12, NULL, &tstat))
  241             tstat = 0;  /* reset keyword not found error */
  242         else
  243             cd_exists = 1;  /* found at least 1 CD_ keyword */
  244 
  245         if (ffgkyd(fptr, "CD2_2", &cd22, NULL, &tstat))
  246             tstat = 0;  /* reset keyword not found error */
  247         else
  248             cd_exists = 1;  /* found at least 1 CD_ keyword */
  249 
  250         if (cd_exists)  /* convert CDi_j back to CDELTn */
  251         {
  252             /* there are 2 ways to compute the angle: */
  253             phia = atan2( cd21, cd11);
  254             phib = atan2(-cd12, cd22);
  255 
  256             /* ensure that phia <= phib */
  257             temp = minvalue(phia, phib);
  258             phib = maxvalue(phia, phib);
  259             phia = temp;
  260 
  261             /* there is a possible 180 degree ambiguity in the angles */
  262             /* so add 180 degress to the smaller value if the values  */
  263             /* differ by more than 90 degrees = pi/2 radians.         */
  264             /* (Later, we may decide to take the other solution by    */
  265             /* subtracting 180 degrees from the larger value).        */
  266 
  267             if ((phib - phia) > (pi / 2.))
  268                phia += pi;
  269 
  270             if (fabs(phia - phib) > toler) 
  271             {
  272                /* angles don't agree, so looks like there is some skewness */
  273                /* between the axes.  Return with an error to be safe. */
  274                *status = APPROX_WCS_KEY;
  275             }
  276       
  277             phia = (phia + phib) /2.;  /* use the average of the 2 values */
  278             *xinc = cd11 / cos(phia);
  279             *yinc = cd22 / cos(phia);
  280             *rot = phia * 180. / pi;
  281 
  282             /* common usage is to have a positive yinc value.  If it is */
  283             /* negative, then subtract 180 degrees from rot and negate  */
  284             /* both xinc and yinc.  */
  285 
  286             if (*yinc < 0)
  287             {
  288                 *xinc = -(*xinc);
  289                 *yinc = -(*yinc);
  290                 *rot = *rot - 180.;
  291             }
  292         }
  293         else   /* no CD matrix keywords either */
  294         {
  295             *xinc = 1.;
  296 
  297             /* there was no CDELT1 keyword, but check for CDELT2 just in case */
  298             tstat = 0;
  299             if (ffgkyd(fptr, "CDELT2", yinc, NULL, &tstat))
  300                 *yinc = 1.;
  301 
  302             tstat = 0;
  303             if (ffgkyd(fptr, "CROTA2", rot, NULL, &tstat))
  304                 *rot=0.;
  305         }
  306     }
  307     else  /* Case 2: CDELTn + optional PC matrix */
  308     {
  309         if (ffgkyd(fptr, "CDELT2", yinc, NULL, &tstat))
  310             *yinc = 1.;
  311 
  312         tstat = 0;
  313         if (ffgkyd(fptr, "CROTA2", rot, NULL, &tstat))
  314         {
  315             *rot=0.;
  316 
  317             /* no CROTA2 keyword, so look for the PC matrix */
  318             tstat = 0;
  319             if (ffgkyd(fptr, "PC1_1", &pc11, NULL, &tstat))
  320                 tstat = 0;  /* reset keyword not found error */
  321             else
  322                 pc_exists = 1;  /* found at least 1 PC_ keyword */
  323 
  324             if (ffgkyd(fptr, "PC2_1", &pc21, NULL, &tstat))
  325                 tstat = 0;  /* reset keyword not found error */
  326             else
  327                 pc_exists = 1;  /* found at least 1 PC_ keyword */
  328 
  329             if (ffgkyd(fptr, "PC1_2", &pc12, NULL, &tstat))
  330                 tstat = 0;  /* reset keyword not found error */
  331             else
  332                 pc_exists = 1;  /* found at least 1 PC_ keyword */
  333 
  334             if (ffgkyd(fptr, "PC2_2", &pc22, NULL, &tstat))
  335                 tstat = 0;  /* reset keyword not found error */
  336             else
  337                 pc_exists = 1;  /* found at least 1 PC_ keyword */
  338 
  339             if (pc_exists)  /* convert PCi_j back to CDELTn */
  340             {
  341                 /* there are 2 ways to compute the angle: */
  342                 phia = atan2( pc21, pc11);
  343                 phib = atan2(-pc12, pc22);
  344 
  345                 /* ensure that phia <= phib */
  346                 temp = minvalue(phia, phib);
  347                 phib = maxvalue(phia, phib);
  348                 phia = temp;
  349 
  350                 /* there is a possible 180 degree ambiguity in the angles */
  351                 /* so add 180 degress to the smaller value if the values  */
  352                 /* differ by more than 90 degrees = pi/2 radians.         */
  353                 /* (Later, we may decide to take the other solution by    */
  354                 /* subtracting 180 degrees from the larger value).        */
  355 
  356                 if ((phib - phia) > (pi / 2.))
  357                    phia += pi;
  358 
  359                 if (fabs(phia - phib) > toler) 
  360                 {
  361                   /* angles don't agree, so looks like there is some skewness */
  362                   /* between the axes.  Return with an error to be safe. */
  363                   *status = APPROX_WCS_KEY;
  364                 }
  365       
  366                 phia = (phia + phib) /2.;  /* use the average of the 2 values */
  367                 *rot = phia * 180. / pi;
  368             }
  369         }
  370     }
  371 
  372     /* get the type of projection, if any */
  373     tstat = 0;
  374     if (ffgkys(fptr, "CTYPE1", ctype, NULL, &tstat))
  375          type[0] = '\0';
  376     else
  377     {
  378         /* copy the projection type string */
  379         strncpy(type, &ctype[4], 4);
  380         type[4] = '\0';
  381 
  382         /* check if RA and DEC are inverted */
  383         if (!strncmp(ctype, "DEC-", 4) || !strncmp(ctype+1, "LAT", 3))
  384         {
  385             /* the latitudinal axis is given first, so swap them */
  386 
  387 /*
  388  this case was removed on 12/9.  Apparently not correct.
  389 
  390             if ((*xinc / *yinc) < 0. )  
  391                 *rot = -90. - (*rot);
  392             else
  393 */
  394             *rot = 90. - (*rot);
  395 
  396             /* Empirical tests with ds9 show the y-axis sign must be negated */
  397             /* and the xinc and yinc values must NOT be swapped. */
  398             *yinc = -(*yinc);
  399 
  400             temp = *xrval;
  401             *xrval = *yrval;
  402             *yrval = temp;
  403         }   
  404     }
  405 
  406     return(*status);
  407 }
  408 /*--------------------------------------------------------------------------*/
  409 int ffgicsa(fitsfile *fptr,    /* I - FITS file pointer           */
  410            char version,      /* I - character code of desired version */
  411                           /*     A - Z or blank */
  412            double *xrval,     /* O - X reference value           */
  413            double *yrval,     /* O - Y reference value           */
  414            double *xrpix,     /* O - X reference pixel           */
  415            double *yrpix,     /* O - Y reference pixel           */
  416            double *xinc,      /* O - X increment per pixel       */
  417            double *yinc,      /* O - Y increment per pixel       */
  418            double *rot,       /* O - rotation angle (degrees)    */
  419            char *type,        /* O - type of projection ('-tan') */
  420            int *status)       /* IO - error status               */
  421 /*
  422        read the values of the celestial coordinate system keywords.
  423        These values may be used as input to the subroutines that
  424        calculate celestial coordinates. (ffxypx, ffwldp)
  425 
  426        Modified in Nov 1999 to convert the CD matrix keywords back
  427        to the old CDELTn form, and to swap the axes if the dec-like
  428        axis is given first, and to assume default values if any of the
  429        keywords are not present.
  430 */
  431 {
  432     int tstat = 0, cd_exists = 0, pc_exists = 0;
  433     char ctype[FLEN_VALUE], keyname[FLEN_VALUE], alt[2];
  434     double cd11 = 0.0, cd21 = 0.0, cd22 = 0.0, cd12 = 0.0;
  435     double pc11 = 1.0, pc21 = 0.0, pc22 = 1.0, pc12 = 0.0;
  436     double pi =  3.1415926535897932;
  437     double phia, phib, temp;
  438     double toler = .0002;  /* tolerance for angles to agree (radians) */
  439                            /*   (= approximately 0.01 degrees) */
  440 
  441     if (*status > 0)
  442        return(*status);
  443 
  444     if (version == ' ') {
  445       ffgics(fptr, xrval, yrval, xrpix, yrpix, xinc, yinc, rot, type, status);
  446       return (*status);
  447     }
  448 
  449     if (version > 'Z' || version < 'A') {
  450       ffpmsg("ffgicsa: illegal WCS version code (must be A - Z or blank)");
  451       return(*status = WCS_ERROR);
  452     }
  453 
  454     alt[0] = version;
  455     alt[1] = '\0';
  456     
  457     tstat = 0;
  458     strcpy(keyname, "CRVAL1");
  459     strcat(keyname, alt);
  460     if (ffgkyd(fptr, keyname, xrval, NULL, &tstat))
  461        *xrval = 0.;
  462 
  463     tstat = 0;
  464     strcpy(keyname, "CRVAL2");
  465     strcat(keyname, alt);
  466     if (ffgkyd(fptr, keyname, yrval, NULL, &tstat))
  467        *yrval = 0.;
  468 
  469     tstat = 0;
  470     strcpy(keyname, "CRPIX1");
  471     strcat(keyname, alt);
  472     if (ffgkyd(fptr, keyname, xrpix, NULL, &tstat))
  473         *xrpix = 0.;
  474 
  475     tstat = 0;
  476     strcpy(keyname, "CRPIX2");
  477     strcat(keyname, alt);
  478      if (ffgkyd(fptr, keyname, yrpix, NULL, &tstat))
  479         *yrpix = 0.;
  480 
  481     /* look for CDELTn first, then CDi_j keywords */
  482     tstat = 0;
  483     strcpy(keyname, "CDELT1");
  484     strcat(keyname, alt);
  485     if (ffgkyd(fptr, keyname, xinc, NULL, &tstat))
  486     {
  487         /* CASE 1: no CDELTn keyword, so look for the CD matrix */
  488         tstat = 0;
  489         strcpy(keyname, "CD1_1");
  490         strcat(keyname, alt);
  491         if (ffgkyd(fptr, keyname, &cd11, NULL, &tstat))
  492             tstat = 0;  /* reset keyword not found error */
  493         else
  494             cd_exists = 1;  /* found at least 1 CD_ keyword */
  495 
  496         strcpy(keyname, "CD2_1");
  497         strcat(keyname, alt);
  498         if (ffgkyd(fptr, keyname, &cd21, NULL, &tstat))
  499             tstat = 0;  /* reset keyword not found error */
  500         else
  501             cd_exists = 1;  /* found at least 1 CD_ keyword */
  502 
  503         strcpy(keyname, "CD1_2");
  504         strcat(keyname, alt);
  505         if (ffgkyd(fptr, keyname, &cd12, NULL, &tstat))
  506             tstat = 0;  /* reset keyword not found error */
  507         else
  508             cd_exists = 1;  /* found at least 1 CD_ keyword */
  509 
  510         strcpy(keyname, "CD2_2");
  511         strcat(keyname, alt);
  512         if (ffgkyd(fptr, keyname, &cd22, NULL, &tstat))
  513             tstat = 0;  /* reset keyword not found error */
  514         else
  515             cd_exists = 1;  /* found at least 1 CD_ keyword */
  516 
  517         if (cd_exists)  /* convert CDi_j back to CDELTn */
  518         {
  519             /* there are 2 ways to compute the angle: */
  520             phia = atan2( cd21, cd11);
  521             phib = atan2(-cd12, cd22);
  522 
  523             /* ensure that phia <= phib */
  524             temp = minvalue(phia, phib);
  525             phib = maxvalue(phia, phib);
  526             phia = temp;
  527 
  528             /* there is a possible 180 degree ambiguity in the angles */
  529             /* so add 180 degress to the smaller value if the values  */
  530             /* differ by more than 90 degrees = pi/2 radians.         */
  531             /* (Later, we may decide to take the other solution by    */
  532             /* subtracting 180 degrees from the larger value).        */
  533 
  534             if ((phib - phia) > (pi / 2.))
  535                phia += pi;
  536 
  537             if (fabs(phia - phib) > toler) 
  538             {
  539                /* angles don't agree, so looks like there is some skewness */
  540                /* between the axes.  Return with an error to be safe. */
  541                *status = APPROX_WCS_KEY;
  542             }
  543       
  544             phia = (phia + phib) /2.;  /* use the average of the 2 values */
  545             *xinc = cd11 / cos(phia);
  546             *yinc = cd22 / cos(phia);
  547             *rot = phia * 180. / pi;
  548 
  549             /* common usage is to have a positive yinc value.  If it is */
  550             /* negative, then subtract 180 degrees from rot and negate  */
  551             /* both xinc and yinc.  */
  552 
  553             if (*yinc < 0)
  554             {
  555                 *xinc = -(*xinc);
  556                 *yinc = -(*yinc);
  557                 *rot = *rot - 180.;
  558             }
  559         }
  560         else   /* no CD matrix keywords either */
  561         {
  562             *xinc = 1.;
  563 
  564             /* there was no CDELT1 keyword, but check for CDELT2 just in case */
  565             tstat = 0;
  566             strcpy(keyname, "CDELT2");
  567             strcat(keyname, alt);
  568             if (ffgkyd(fptr, keyname, yinc, NULL, &tstat))
  569                 *yinc = 1.;
  570 
  571             tstat = 0;
  572             strcpy(keyname, "CROTA2");
  573             strcat(keyname, alt);
  574             if (ffgkyd(fptr, keyname, rot, NULL, &tstat))
  575                 *rot=0.;
  576         }
  577     }
  578     else  /* Case 2: CDELTn + optional PC matrix */
  579     {
  580         strcpy(keyname, "CDELT2");
  581         strcat(keyname, alt);
  582         if (ffgkyd(fptr, keyname, yinc, NULL, &tstat))
  583             *yinc = 1.;
  584 
  585         tstat = 0;
  586         strcpy(keyname, "CROTA2");
  587         strcat(keyname, alt);
  588         if (ffgkyd(fptr, keyname, rot, NULL, &tstat))
  589         {
  590             *rot=0.;
  591 
  592             /* no CROTA2 keyword, so look for the PC matrix */
  593             tstat = 0;
  594             strcpy(keyname, "PC1_1");
  595             strcat(keyname, alt);
  596             if (ffgkyd(fptr, keyname, &pc11, NULL, &tstat))
  597                 tstat = 0;  /* reset keyword not found error */
  598             else
  599                 pc_exists = 1;  /* found at least 1 PC_ keyword */
  600 
  601             strcpy(keyname, "PC2_1");
  602             strcat(keyname, alt);
  603             if (ffgkyd(fptr, keyname, &pc21, NULL, &tstat))
  604                 tstat = 0;  /* reset keyword not found error */
  605             else
  606                 pc_exists = 1;  /* found at least 1 PC_ keyword */
  607 
  608             strcpy(keyname, "PC1_2");
  609             strcat(keyname, alt);
  610             if (ffgkyd(fptr, keyname, &pc12, NULL, &tstat))
  611                 tstat = 0;  /* reset keyword not found error */
  612             else
  613                 pc_exists = 1;  /* found at least 1 PC_ keyword */
  614 
  615             strcpy(keyname, "PC2_2");
  616             strcat(keyname, alt);
  617             if (ffgkyd(fptr, keyname, &pc22, NULL, &tstat))
  618                 tstat = 0;  /* reset keyword not found error */
  619             else
  620                 pc_exists = 1;  /* found at least 1 PC_ keyword */
  621 
  622             if (pc_exists)  /* convert PCi_j back to CDELTn */
  623             {
  624                 /* there are 2 ways to compute the angle: */
  625                 phia = atan2( pc21, pc11);
  626                 phib = atan2(-pc12, pc22);
  627 
  628                 /* ensure that phia <= phib */
  629                 temp = minvalue(phia, phib);
  630                 phib = maxvalue(phia, phib);
  631                 phia = temp;
  632 
  633                 /* there is a possible 180 degree ambiguity in the angles */
  634                 /* so add 180 degress to the smaller value if the values  */
  635                 /* differ by more than 90 degrees = pi/2 radians.         */
  636                 /* (Later, we may decide to take the other solution by    */
  637                 /* subtracting 180 degrees from the larger value).        */
  638 
  639                 if ((phib - phia) > (pi / 2.))
  640                    phia += pi;
  641 
  642                 if (fabs(phia - phib) > toler) 
  643                 {
  644                   /* angles don't agree, so looks like there is some skewness */
  645                   /* between the axes.  Return with an error to be safe. */
  646                   *status = APPROX_WCS_KEY;
  647                 }
  648       
  649                 phia = (phia + phib) /2.;  /* use the average of the 2 values */
  650                 *rot = phia * 180. / pi;
  651             }
  652         }
  653     }
  654 
  655     /* get the type of projection, if any */
  656     tstat = 0;
  657     strcpy(keyname, "CTYPE1");
  658     strcat(keyname, alt);
  659     if (ffgkys(fptr, keyname, ctype, NULL, &tstat))
  660          type[0] = '\0';
  661     else
  662     {
  663         /* copy the projection type string */
  664         strncpy(type, &ctype[4], 4);
  665         type[4] = '\0';
  666 
  667         /* check if RA and DEC are inverted */
  668         if (!strncmp(ctype, "DEC-", 4) || !strncmp(ctype+1, "LAT", 3))
  669         {
  670             /* the latitudinal axis is given first, so swap them */
  671 
  672             *rot = 90. - (*rot);
  673 
  674             /* Empirical tests with ds9 show the y-axis sign must be negated */
  675             /* and the xinc and yinc values must NOT be swapped. */
  676             *yinc = -(*yinc);
  677 
  678             temp = *xrval;
  679             *xrval = *yrval;
  680             *yrval = temp;
  681         }   
  682     }
  683 
  684     return(*status);
  685 }
  686 /*--------------------------------------------------------------------------*/
  687 int ffgtcs(fitsfile *fptr,    /* I - FITS file pointer           */
  688            int xcol,          /* I - column containing the RA coordinate  */
  689            int ycol,          /* I - column containing the DEC coordinate */
  690            double *xrval,     /* O - X reference value           */
  691            double *yrval,     /* O - Y reference value           */
  692            double *xrpix,     /* O - X reference pixel           */
  693            double *yrpix,     /* O - Y reference pixel           */
  694            double *xinc,      /* O - X increment per pixel       */
  695            double *yinc,      /* O - Y increment per pixel       */
  696            double *rot,       /* O - rotation angle (degrees)    */
  697            char *type,        /* O - type of projection ('-sin') */
  698            int *status)       /* IO - error status               */
  699 /*
  700        read the values of the celestial coordinate system keywords
  701        from a FITS table where the X and Y or RA and DEC coordinates
  702        are stored in separate column.  Do this by converting the
  703        table to a temporary FITS image, then reading the keywords
  704        from the image file.
  705        These values may be used as input to the subroutines that
  706        calculate celestial coordinates. (ffxypx, ffwldp)
  707 */
  708 {
  709     int colnum[2];
  710     long naxes[2];
  711     fitsfile *tptr;
  712 
  713     if (*status > 0)
  714        return(*status);
  715 
  716     colnum[0] = xcol;
  717     colnum[1] = ycol;
  718     naxes[0] = 10;
  719     naxes[1] = 10;
  720 
  721     /* create temporary  FITS file, in memory */
  722     ffinit(&tptr, "mem://", status);
  723     
  724     /* create a temporary image; the datatype and size are not important */
  725     ffcrim(tptr, 32, 2, naxes, status);
  726     
  727     /* now copy the relevant keywords from the table to the image */
  728     fits_copy_pixlist2image(fptr, tptr, 9, 2, colnum, status);
  729 
  730     /* write default WCS keywords, if they are not present */
  731     fits_write_keys_histo(fptr, tptr, 2, colnum, status);
  732 
  733     if (*status > 0)
  734        return(*status);
  735          
  736     /* read the WCS keyword values from the temporary image */
  737     ffgics(tptr, xrval, yrval, xrpix, yrpix, xinc, yinc, rot, type, status); 
  738 
  739     if (*status > 0)
  740     {
  741       ffpmsg
  742       ("ffgtcs could not find all the celestial coordinate keywords");
  743       return(*status = NO_WCS_KEY); 
  744     }
  745 
  746     /* delete the temporary file */
  747     fits_delete_file(tptr, status);
  748     
  749     return(*status);
  750 }
  751 /*--------------------------------------------------------------------------*/
  752 int ffgtwcs(fitsfile *fptr,  /* I - FITS file pointer              */
  753            int xcol,        /* I - column number for the X column  */
  754            int ycol,        /* I - column number for the Y column  */
  755            char **header,   /* O - string of all the WCS keywords  */
  756            int *status)     /* IO - error status                   */
  757 /*
  758   int fits_get_table_wcs_keys
  759   Return string containing all the WCS keywords appropriate for the 
  760   pair of X and Y columns containing the coordinate
  761   of each event in an event list table.  This string may then be passed
  762   to Doug Mink's WCS library wcsinit routine, to create and initialize the
  763   WCS structure.  The calling routine must free the header character string
  764   when it is no longer needed. 
  765 
  766   THIS ROUTINE IS DEPRECATED. USE fits_hdr2str INSTEAD
  767 */
  768 {
  769     int hdutype, ncols, tstatus, length;
  770     int naxis1 = 1, naxis2 = 1;
  771     long tlmin, tlmax;
  772     char keyname[FLEN_KEYWORD];
  773     char valstring[FLEN_VALUE];
  774     char comm[2];
  775     char *cptr;
  776     /*  construct a string of 80 blanks, for adding fill to the keywords */
  777                  /*  12345678901234567890123456789012345678901234567890123456789012345678901234567890 */
  778     char blanks[] = "                                                                                ";
  779 
  780     if (*status > 0)
  781         return(*status);
  782 
  783     fits_get_hdu_type(fptr, &hdutype, status);
  784     if (hdutype == IMAGE_HDU)
  785     {
  786         ffpmsg("Can't read table WSC keywords. This HDU is not a table");
  787         return(*status = NOT_TABLE);
  788     }
  789 
  790     fits_get_num_cols(fptr, &ncols, status);
  791     
  792     if (xcol < 1 || xcol > ncols)
  793     {
  794         ffpmsg("illegal X axis column number in fftwcs");
  795         return(*status = BAD_COL_NUM);
  796     }
  797 
  798     if (ycol < 1 || ycol > ncols)
  799     {
  800         ffpmsg("illegal Y axis column number in fftwcs");
  801         return(*status = BAD_COL_NUM);
  802     }
  803 
  804     /* allocate character string for all the WCS keywords */
  805     *header = calloc(1, 2401);  /* room for up to 30 keywords */
  806     if (*header == 0)
  807     {
  808         ffpmsg("error allocating memory for WCS header keywords (fftwcs)");
  809         return(*status = MEMORY_ALLOCATION);
  810     }
  811 
  812     cptr = *header;
  813     comm[0] = '\0';
  814     
  815     tstatus = 0;
  816     ffkeyn("TLMIN",xcol,keyname,status);
  817     ffgkyj(fptr,keyname, &tlmin,NULL,&tstatus);
  818 
  819     if (!tstatus)
  820     {
  821         ffkeyn("TLMAX",xcol,keyname,status);
  822         ffgkyj(fptr,keyname, &tlmax,NULL,&tstatus);
  823     }
  824 
  825     if (!tstatus)
  826     {
  827         naxis1 = tlmax - tlmin + 1;
  828     }
  829 
  830     tstatus = 0;
  831     ffkeyn("TLMIN",ycol,keyname,status);
  832     ffgkyj(fptr,keyname, &tlmin,NULL,&tstatus);
  833 
  834     if (!tstatus)
  835     {
  836         ffkeyn("TLMAX",ycol,keyname,status);
  837         ffgkyj(fptr,keyname, &tlmax,NULL,&tstatus);
  838     }
  839 
  840     if (!tstatus)
  841     {
  842         naxis2 = tlmax - tlmin + 1;
  843     }
  844 
  845     /*            123456789012345678901234567890    */
  846     strcat(cptr, "NAXIS   =                    2");
  847     strncat(cptr, blanks, 50);
  848     cptr += 80;
  849 
  850     ffi2c(naxis1, valstring, status);   /* convert to formatted string */
  851     ffmkky("NAXIS1", valstring, comm, cptr, status);  /* construct the keyword*/
  852     strncat(cptr, blanks, 50);  /* pad with blanks */
  853     cptr += 80;
  854 
  855     strcpy(keyname, "NAXIS2");
  856     ffi2c(naxis2, valstring, status);   /* convert to formatted string */
  857     ffmkky(keyname, valstring, comm, cptr, status);  /* construct the keyword*/
  858     strncat(cptr, blanks, 50);  /* pad with blanks */
  859     cptr += 80;
  860 
  861     /* read the required header keywords (use defaults if not found) */
  862 
  863     /*  CTYPE1 keyword */
  864     tstatus = 0;
  865     ffkeyn("TCTYP",xcol,keyname,status);
  866     if (ffgkey(fptr, keyname, valstring, NULL, &tstatus) )
  867        valstring[0] =  '\0';
  868     ffmkky("CTYPE1", valstring, comm, cptr, status);  /* construct the keyword*/
  869     length = strlen(cptr);
  870     strncat(cptr, blanks, 80 - length);  /* pad with blanks */
  871     cptr += 80;
  872 
  873     /*  CTYPE2 keyword */
  874     tstatus = 0;
  875     ffkeyn("TCTYP",ycol,keyname,status);
  876     if (ffgkey(fptr, keyname, valstring, NULL, &tstatus) )
  877        valstring[0] =  '\0';
  878     ffmkky("CTYPE2", valstring, comm, cptr, status);  /* construct the keyword*/
  879     length = strlen(cptr);
  880     strncat(cptr, blanks, 80 - length);  /* pad with blanks */
  881     cptr += 80;
  882 
  883     /*  CRPIX1 keyword */
  884     tstatus = 0;
  885     ffkeyn("TCRPX",xcol,keyname,status);
  886     if (ffgkey(fptr, keyname, valstring, NULL, &tstatus) )
  887        strcpy(valstring, "1");
  888     ffmkky("CRPIX1", valstring, comm, cptr, status);  /* construct the keyword*/
  889     strncat(cptr, blanks, 50);  /* pad with blanks */
  890     cptr += 80;
  891 
  892     /*  CRPIX2 keyword */
  893     tstatus = 0;
  894     ffkeyn("TCRPX",ycol,keyname,status);
  895     if (ffgkey(fptr, keyname, valstring, NULL, &tstatus) )
  896        strcpy(valstring, "1");
  897     ffmkky("CRPIX2", valstring, comm, cptr, status);  /* construct the keyword*/
  898     strncat(cptr, blanks, 50);  /* pad with blanks */
  899     cptr += 80;
  900 
  901     /*  CRVAL1 keyword */
  902     tstatus = 0;
  903     ffkeyn("TCRVL",xcol,keyname,status);
  904     if (ffgkey(fptr, keyname, valstring, NULL, &tstatus) )
  905        strcpy(valstring, "1");
  906     ffmkky("CRVAL1", valstring, comm, cptr, status);  /* construct the keyword*/
  907     strncat(cptr, blanks, 50);  /* pad with blanks */
  908     cptr += 80;
  909 
  910     /*  CRVAL2 keyword */
  911     tstatus = 0;
  912     ffkeyn("TCRVL",ycol,keyname,status);
  913     if (ffgkey(fptr, keyname, valstring, NULL, &tstatus) )
  914        strcpy(valstring, "1");
  915     ffmkky("CRVAL2", valstring, comm, cptr, status);  /* construct the keyword*/
  916     strncat(cptr, blanks, 50);  /* pad with blanks */
  917     cptr += 80;
  918 
  919     /*  CDELT1 keyword */
  920     tstatus = 0;
  921     ffkeyn("TCDLT",xcol,keyname,status);
  922     if (ffgkey(fptr, keyname, valstring, NULL, &tstatus) )
  923        strcpy(valstring, "1");
  924     ffmkky("CDELT1", valstring, comm, cptr, status);  /* construct the keyword*/
  925     strncat(cptr, blanks, 50);  /* pad with blanks */
  926     cptr += 80;
  927 
  928     /*  CDELT2 keyword */
  929     tstatus = 0;
  930     ffkeyn("TCDLT",ycol,keyname,status);
  931     if (ffgkey(fptr, keyname, valstring, NULL, &tstatus) )
  932        strcpy(valstring, "1");
  933     ffmkky("CDELT2", valstring, comm, cptr, status);  /* construct the keyword*/
  934     strncat(cptr, blanks, 50);  /* pad with blanks */
  935     cptr += 80;
  936 
  937     /* the following keywords may not exist */
  938 
  939     /*  CROTA2 keyword */
  940     tstatus = 0;
  941     ffkeyn("TCROT",ycol,keyname,status);
  942     if (ffgkey(fptr, keyname, valstring, NULL, &tstatus) == 0 )
  943     {
  944         ffmkky("CROTA2", valstring, comm, cptr, status);  /* construct keyword*/
  945         strncat(cptr, blanks, 50);  /* pad with blanks */
  946         cptr += 80;
  947     }
  948 
  949     /*  EPOCH keyword */
  950     tstatus = 0;
  951     if (ffgkey(fptr, "EPOCH", valstring, NULL, &tstatus) == 0 )
  952     {
  953         ffmkky("EPOCH", valstring, comm, cptr, status);  /* construct keyword*/
  954         length = strlen(cptr);
  955         strncat(cptr, blanks, 80 - length);  /* pad with blanks */
  956         cptr += 80;
  957     }
  958 
  959     /*  EQUINOX keyword */
  960     tstatus = 0;
  961     if (ffgkey(fptr, "EQUINOX", valstring, NULL, &tstatus) == 0 )
  962     {
  963         ffmkky("EQUINOX", valstring, comm, cptr, status); /* construct keyword*/
  964         length = strlen(cptr);
  965         strncat(cptr, blanks, 80 - length);  /* pad with blanks */
  966         cptr += 80;
  967     }
  968 
  969     /*  RADECSYS keyword */
  970     tstatus = 0;
  971     if (ffgkey(fptr, "RADECSYS", valstring, NULL, &tstatus) == 0 )
  972     {
  973         ffmkky("RADECSYS", valstring, comm, cptr, status); /*construct keyword*/
  974         length = strlen(cptr);
  975         strncat(cptr, blanks, 80 - length);  /* pad with blanks */
  976         cptr += 80;
  977     }
  978 
  979     /*  TELESCOPE keyword */
  980     tstatus = 0;
  981     if (ffgkey(fptr, "TELESCOP", valstring, NULL, &tstatus) == 0 )
  982     {
  983         ffmkky("TELESCOP", valstring, comm, cptr, status); 
  984         length = strlen(cptr);
  985         strncat(cptr, blanks, 80 - length);  /* pad with blanks */
  986         cptr += 80;
  987     }
  988 
  989     /*  INSTRUME keyword */
  990     tstatus = 0;
  991     if (ffgkey(fptr, "INSTRUME", valstring, NULL, &tstatus) == 0 )
  992     {
  993         ffmkky("INSTRUME", valstring, comm, cptr, status);  
  994         length = strlen(cptr);
  995         strncat(cptr, blanks, 80 - length);  /* pad with blanks */
  996         cptr += 80;
  997     }
  998 
  999     /*  DETECTOR keyword */
 1000     tstatus = 0;
 1001     if (ffgkey(fptr, "DETECTOR", valstring, NULL, &tstatus) == 0 )
 1002     {
 1003         ffmkky("DETECTOR", valstring, comm, cptr, status);  
 1004         length = strlen(cptr);
 1005         strncat(cptr, blanks, 80 - length);  /* pad with blanks */
 1006         cptr += 80;
 1007     }
 1008 
 1009     /*  MJD-OBS keyword */
 1010     tstatus = 0;
 1011     if (ffgkey(fptr, "MJD-OBS", valstring, NULL, &tstatus) == 0 )
 1012     {
 1013         ffmkky("MJD-OBS", valstring, comm, cptr, status);  
 1014         length = strlen(cptr);
 1015         strncat(cptr, blanks, 80 - length);  /* pad with blanks */
 1016         cptr += 80;
 1017     }
 1018 
 1019     /*  DATE-OBS keyword */
 1020     tstatus = 0;
 1021     if (ffgkey(fptr, "DATE-OBS", valstring, NULL, &tstatus) == 0 )
 1022     {
 1023         ffmkky("DATE-OBS", valstring, comm, cptr, status);  
 1024         length = strlen(cptr);
 1025         strncat(cptr, blanks, 80 - length);  /* pad with blanks */
 1026         cptr += 80;
 1027     }
 1028 
 1029     /*  DATE keyword */
 1030     tstatus = 0;
 1031     if (ffgkey(fptr, "DATE", valstring, NULL, &tstatus) == 0 )
 1032     {
 1033         ffmkky("DATE", valstring, comm, cptr, status);  
 1034         length = strlen(cptr);
 1035         strncat(cptr, blanks, 80 - length);  /* pad with blanks */
 1036         cptr += 80;
 1037     }
 1038 
 1039     strcat(cptr, "END");
 1040     strncat(cptr, blanks, 77);
 1041 
 1042     return(*status);
 1043 }