fitswcs.c

/*
*				fitswcs.c
*
* Manage World Coordinate System data.
*
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
*
*	This file part of:	AstrOmatic software
*
*	Copyright:		(C) 1993-2016 Emmanuel Bertin -- IAP/CNRS/UPMC
*
*	License:		GNU General Public License
*
*	AstrOmatic software is free software: you can redistribute it and/or
*	modify it under the terms of the GNU General Public License as
*	published by the Free Software Foundation, either version 3 of the
*	License, or (at your option) any later version.
*	AstrOmatic software is distributed in the hope that it will be useful,
*	but WITHOUT ANY WARRANTY; without even the implied warranty of
*	MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
*	GNU General Public License for more details.
*	You should have received a copy of the GNU General Public License
*	along with AstrOmatic software.
*	If not, see <http://www.gnu.org/licenses/>.
*
*	Last modified:		08/03/2016
*
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/

#ifdef HAVE_CONFIG_H
#include	"config.h"
#endif

#ifdef HAVE_MATHIMF_H
#include <mathimf.h>
#else
#include <math.h>
#endif
#include	<stdio.h>
#include	<stdlib.h>
#include	<string.h>

#include	"fits/fitscat_defs.h"
#include	"fits/fitscat.h"
#include	"fitswcs.h"
#include	"wcscelsys.h"
#include	"wcs/wcs.h"
#include	"wcs/lin.h"
#include	"wcs/tnx.h"
#include	"wcs/poly.h"

/******* copy_wcs ************************************************************
PROTO	wcsstruct *copy_wcs(wcsstruct *wcsin)
PURPOSE	Copy a WCS (World Coordinate System) structure.
INPUT	WCS structure to be copied.
OUTPUT	pointer to a copy of the input structure.
NOTES	Actually, only FITS parameters are copied. Lower-level structures
	such as those created by the WCS or TNX libraries are generated.
AUTHOR	E. Bertin (IAP)
VERSION	31/08/2002
 ***/
wcsstruct	*copy_wcs(wcsstruct *wcsin)

  {
   wcsstruct	*wcs;

/* Copy the basic stuff */
  QMEMCPY(wcsin, wcs, wcsstruct, 1);
/* The PROJP WCS parameters */
  QMEMCPY(wcsin->projp, wcs->projp, double, wcs->naxis*100);

/* Set other structure pointers to NULL (they'll have to be reallocated) */
  wcs->wcsprm = NULL;
  wcs->lin = NULL;
  wcs->cel = NULL;
  wcs->prj = NULL;
  wcs->tnx_lngcor = copy_tnxaxis(wcsin->tnx_lngcor);
  wcs->tnx_latcor = copy_tnxaxis(wcsin->tnx_latcor);
  wcs->inv_x = wcs->inv_y = NULL;

  QCALLOC(wcs->wcsprm, struct wcsprm, 1);
/* Test if the WCS is recognized and a celestial pair is found */
  wcsset(wcs->naxis,(const char(*)[9])wcs->ctype, wcs->wcsprm);

/* Initialize other WCS structures */
  init_wcs(wcs);
/* Invert projection corrections */
  invert_wcs(wcs);
/* Find the range of coordinates */
  range_wcs(wcs);  

  return wcs;
  }


/******* create_wcs ***********************************************************
PROTO	wcsstruct *create_wcs(char **ctype, double *crval, double *crpix,
			double *cdelt, int *naxisn, int naxis)
PURPOSE	Generate a simple WCS (World Coordinate System) structure.
INPUT	Pointer to an array of char strings with WCS projection on each axis,
	pointer to an array of center coordinates (double),
	pointer to an array of device coordinates (double),
	pointer to an array of pixel scales (double),
	pointer to an array of image dimensions (int),
	number of dimensions.
OUTPUT	pointer to a WCS structure.
NOTES	If a pointer is set to null, the corresponding variables are set to
	default values.
AUTHOR	E. Bertin (IAP)
VERSION	09/08/2006
 ***/
wcsstruct	*create_wcs(char **ctype, double *crval, double *crpix,
			double *cdelt, int *naxisn, int naxis)

  {
   wcsstruct	*wcs;
   int		l;

  QCALLOC(wcs, wcsstruct, 1);
  wcs->naxis = naxis;
  QCALLOC(wcs->projp, double, naxis*100);
  wcs->nprojp = 0;

  wcs->longpole = wcs->latpole = 999.0;
  for (l=0; l<naxis; l++)
    {
    wcs->naxisn[l] = naxisn? naxisn[l] : 360.0;
/*-- The default WCS projection system is an all-sky Aitoff projection */
    if (ctype)
      strncpy(wcs->ctype[l], ctype[l], 8);
    else if (l==0)
      strncpy(wcs->ctype[l], "RA---AIT", 8);
    else if (l==1)
      strncpy(wcs->ctype[l], "DEC--AIT", 8);
    wcs->crval[l] = crval? crval[l]: 0.0;
    wcs->crpix[l] = crpix? crpix[l]: 0.0;
    wcs->cdelt[l] = 1.0;
    wcs->cd[l*(naxis+1)] = cdelt? cdelt[l] : 1.0;
    }

  wcs->epoch = wcs->equinox = 2000.0;
  QCALLOC(wcs->wcsprm, struct wcsprm, 1);

/* Test if the WCS is recognized and a celestial pair is found */
  wcsset(wcs->naxis,(const char(*)[9])wcs->ctype, wcs->wcsprm);

/* Initialize other WCS structures */
  init_wcs(wcs);
/* Invert projection corrections */
  invert_wcs(wcs);
/* Find the range of coordinates */
  range_wcs(wcs);  

  return wcs;
  }


/******* init_wcs ************************************************************
PROTO	void init_wcs(wcsstruct *wcs)
PURPOSE	Initialize astrometry and WCS (World Coordinate System) structures.
INPUT	WCS structure.
OUTPUT	-.
NOTES	-.
AUTHOR	E. Bertin (IAP)
VERSION	17/05/2007
 ***/
void	init_wcs(wcsstruct *wcs)

  {
   int		l,n,lng,lat,naxis;

  naxis = wcs->naxis;
  if (wcs->lin)
    {
    free(wcs->lin->cdelt);
    free(wcs->lin->crpix);
    free(wcs->lin->pc);
    free(wcs->lin->piximg);
    free(wcs->lin->imgpix);
    free(wcs->lin);
    }
  QCALLOC(wcs->lin, struct linprm, 1);
  QCALLOC(wcs->lin->cdelt, double, naxis);
  QCALLOC(wcs->lin->crpix, double, naxis);
  QCALLOC(wcs->lin->pc, double, naxis*naxis);


  if (wcs->cel)
    free(wcs->cel);
  QCALLOC(wcs->cel, struct celprm, 1);

  if (wcs->prj)
    free(wcs->prj);
  QCALLOC(wcs->prj, struct prjprm, 1);

  if (wcs->inv_x)
    {
    poly_end(wcs->inv_x);
    wcs->inv_x = NULL;
    }
  if (wcs->inv_y)
    {
    poly_end(wcs->inv_y);
    wcs->inv_y = NULL;
    }

/* Set WCS flags to 0: structures will be reinitialized by the WCS library */
  wcs->lin->flag = wcs->cel->flag = wcs->prj->flag = 0;
  wcs->lin->naxis = naxis;

/* wcsprm structure */
  lng = wcs->lng = wcs->wcsprm->lng;
  lat = wcs->lat = wcs->wcsprm->lat;

/* linprm structure */
  for (l=0; l<naxis; l++)
    {
    wcs->lin->crpix[l] = wcs->crpix[l];
    wcs->lin->cdelt[l] = 1.0;
    }

  for (l=0; l<naxis*naxis; l++)
    wcs->lin->pc[l] = wcs->cd[l];

/* celprm structure */
  if (lng>=0)
    {
    wcs->cel->ref[0] = wcs->crval[lng];
    wcs->cel->ref[1] = wcs->crval[lat];
    }
  else
    {
    wcs->cel->ref[0] = wcs->crval[0];
    wcs->cel->ref[1] = wcs->crval[1];
    }
  wcs->cel->ref[2] = wcs->longpole;
  wcs->cel->ref[3] = wcs->latpole;

/* prjprm structure */
  wcs->prj->r0 = wcs->r0;
  wcs->prj->tnx_lngcor = wcs->tnx_lngcor;
  wcs->prj->tnx_latcor = wcs->tnx_latcor;
  if (lng>=0)
    {
    n = 0;
    for (l=100; l--;)
      {
      wcs->prj->p[l] = wcs->projp[l+lat*100];	/* lat comes first for ... */
      wcs->prj->p[l+100] = wcs->projp[l+lng*100];/* ... compatibility reasons */
      if (!n && (wcs->prj->p[l] || wcs->prj->p[l+100]))
        n = l+1;
      }
    wcs->nprojp = n;
    }

/* Check-out chirality */
  wcs->chirality = wcs_chirality(wcs);

/* Initialize Equatorial <=> Celestial coordinate system transforms */
  init_wcscelsys(wcs);

  return;
  }


/******* init_wcscelsys *******************************************************
PROTO	void init_wcscelsys(wcsstruct *wcs)
PURPOSE	Initialize Equatorial <=> Celestial coordinate system transforms.
INPUT	WCS structure.
OUTPUT	-.
NOTES	-.
AUTHOR	E. Bertin (IAP)
VERSION	18/07/2006
 ***/
void	init_wcscelsys(wcsstruct *wcs)

  {
  double	*mat,
		a0,d0,ap,dp,ap2,y;
  int		s,lng,lat;

  lng = wcs->wcsprm->lng;
  lat = wcs->wcsprm->lat;
/* Is it a celestial system? If not, exit! */
  if (lng==lat)
    {
    wcs->celsysconvflag = 0;
    return;
    }
/* Find the celestial system */
  for (s=0; *celsysname[s][0] && strncmp(wcs->ctype[lng], celsysname[s][0], 4);
	s++);
/* Is it a known, non-equatorial system? If not, exit! */
  if (!s || !*celsysname[s][0])
    {
    wcs->celsysconvflag = 0;
    return;
    }
  wcs->celsys = (celsysenum)s;
/* Some shortcuts */
  a0 = celsysorig[s][0]*DEG;
  d0 = celsysorig[s][1]*DEG;
  ap = celsyspole[s][0]*DEG;
  dp = celsyspole[s][1]*DEG;
/* First compute in the output referential the longitude of the south pole */
  y = sin(ap - a0);
/*
  x = cos(d0)*(cos(d0)*sin(dp)*cos(ap-a0)-sin(d0)*cos(dp));
  ap2 = atan2(y,x);
*/
  ap2 = asin(cos(d0)*y) ;
/* Equatorial <=> Celestial System transformation parameters */
  mat = wcs->celsysmat;
  mat[0] = ap;
  mat[1] = ap2;
  mat[2] = cos(dp);
  mat[3] = sin(dp);

  wcs->celsysconvflag = 1;
  return;
  }


/******* read_wcs *************************************************************
PROTO	wcsstruct *read_wcs(tabstruct *tab)
PURPOSE	Read WCS (World Coordinate System) info in the FITS header.
INPUT	tab structure.
OUTPUT	-.
NOTES	-.
AUTHOR	E. Bertin (IAP)
VERSION	15/11/2013
 ***/
wcsstruct	*read_wcs(tabstruct *tab)

  {
#define	FITSREADF(buf, k, val, def) \
		{if (fitsread(buf,k, &val, H_FLOAT,T_DOUBLE) != RETURN_OK) \
		   val = def; \
		}

#define	FITSREADI(buf, k, val, def) \
		{if (fitsread(buf,k, &val, H_INT,T_LONG) != RETURN_OK) \
		   val = def; \
		}

#define	FITSREADS(buf, k, str, def) \
		{if (fitsread(buf,k,str, H_STRING,T_STRING) != RETURN_OK) \
		   strcpy(str, (def)); \
		}
   char		str[MAXCHARS];
   char		wstr1[TNX_MAXCHARS], wstr2[TNX_MAXCHARS];

   wcsstruct	*wcs;
   double	drota;
   int		j, l, naxis;
   char		name[16],
		*buf, *filename, *ptr;

  buf = tab->headbuf;
  filename = (tab->cat? tab->cat->filename : strcpy(name, "internal header"));

  FITSREADS(buf, "OBJECT  ", str, "Unnamed");

  QCALLOC(wcs, wcsstruct, 1);
  if (tab->naxis > NAXIS)
    {
    warning("Maximum number of dimensions supported by this version of the ",
	"software exceeded\n");
    tab->naxis = 2;
    }

  wcs->naxis = naxis = tab->naxis;
  QCALLOC(wcs->projp, double, naxis*100);

  for (l=0; l<naxis; l++)
    {
    wcs->naxisn[l] = tab->naxisn[l];
    sprintf(str, "CTYPE%-3d", l+1);
    FITSREADS(buf, str, str, "");
    strncpy(wcs->ctype[l], str, 8);
    sprintf(str, "CUNIT%-3d", l+1);
    FITSREADS(buf, str, str, "deg");
    strncpy(wcs->cunit[l], str, 32);
    sprintf(str, "CRVAL%-3d", l+1);
    FITSREADF(buf, str, wcs->crval[l], 0.0);
    sprintf(str, "CRPIX%-3d", l+1);
    FITSREADF(buf, str, wcs->crpix[l], 1.0);
    sprintf(str, "CDELT%-3d", l+1);
    FITSREADF(buf, str, wcs->cdelt[l], 1.0);
    sprintf(str, "CRDER%-3d", l+1);
    FITSREADF(buf, str, wcs->crder[l], 0.0);
    sprintf(str, "CSYER%-3d", l+1);
    FITSREADF(buf, str, wcs->csyer[l], 0.0);
    if (fabs(wcs->cdelt[l]) < 1e-30)
      error(EXIT_FAILURE, "*Error*: CDELT parameters out of range in ",
	filename);
    }

  if (fitsfind(buf, "CD?_????")!=RETURN_ERROR)
/*-- If CD keywords exist, use them for the linear mapping terms... */
    for (l=0; l<naxis; l++)
      for (j=0; j<naxis; j++)
        {
        sprintf(str, "CD%d_%d", l+1, j+1);
        FITSREADF(buf, str, wcs->cd[l*naxis+j], l==j?1.0:0.0)
        }
  else if (fitsfind(buf, "PC?_????")!=RETURN_ERROR)
/*-- ...If PC keywords exist, use them for the linear mapping terms... */
    for (l=0; l<naxis; l++)
      for (j=0; j<naxis; j++)
        {
        sprintf(str, "PC%d_%d", l+1, j+1);
        FITSREADF(buf, str, wcs->cd[l*naxis+j], l==j?1.0:0.0)
        wcs->cd[l*naxis+j] *= wcs->cdelt[l];
        }
  else if (fitsfind(buf, "PC0??0??")!=RETURN_ERROR)
/*-- ...If PC keywords exist, use them for the linear mapping terms... */
    for (l=0; l<naxis; l++)
      for (j=0; j<naxis; j++)
        {
        sprintf(str, "PC%03d%03d", l+1, j+1);
        FITSREADF(buf, str, wcs->cd[l*naxis+j], l==j?1.0:0.0)
        wcs->cd[l*naxis+j] *= wcs->cdelt[l];
        }
  else
    {
/*-- ...otherwise take the obsolete CROTA2 parameter */
    FITSREADF(buf, "CROTA2  ", drota, 0.0)
    wcs->cd[3] = wcs->cd[0] = cos(drota*DEG);
    wcs->cd[1] = -(wcs->cd[2] = sin(drota*DEG));
    wcs->cd[0] *= wcs->cdelt[0];
    wcs->cd[2] *= wcs->cdelt[0];
    wcs->cd[1] *= wcs->cdelt[1];
    wcs->cd[3] *= wcs->cdelt[1];
    }
  QCALLOC(wcs->wcsprm, struct wcsprm, 1);

/* Test if the WCS is recognized and a celestial pair is found */
  if (!wcsset(wcs->naxis,(const char(*)[9])wcs->ctype, wcs->wcsprm)
	&& wcs->wcsprm->flag<999)
    {
     char	*pstr;
     double	date;
     int	biss, dpar[3];

/*-- Coordinate reference frame */
/*-- Search for an observation date expressed in Julian days */
    FITSREADF(buf, "MJD-OBS ", date, -1.0);
    if (date<0.0)
      FITSREADF(buf, "MJDSTART", date, -1.0);
/*-- Precession date (defined from Ephemerides du Bureau des Longitudes) */
/*-- in Julian years from 2000.0 */
    if (date>0.0)
      wcs->obsdate = 2000.0 - (MJD2000 - date)/365.25;
    else
      {
/*---- Search for an observation date expressed in "civilian" format */
      FITSREADS(buf, "DATE-OBS", str, "");
      if (*str)
        {
/*------ Decode DATE-OBS format: DD/MM/YY or YYYY-MM-DD */
        for (l=0; l<3 && (pstr = strtok_r(l?NULL:str,"/- ", &ptr)); l++)
          dpar[l] = atoi(pstr);
        if (l<3 || !dpar[0] || !dpar[1] || !dpar[2])
          {
/*-------- If DATE-OBS value corrupted or incomplete, assume 2000-1-1 */
          warning("Invalid DATE-OBS value in header: ", str);
          dpar[0] = 2000; dpar[1] = 1; dpar[2] = 1;
          }
        else if (strchr(str, '/') && dpar[0]<32 && dpar[2]<100)
          {
          j = dpar[0];
          dpar[0] = dpar[2]+1900;
          dpar[2] = j;
          }

        biss = (dpar[0]%4)?0:1;
/*------ Convert date to MJD */
        date = -678956 + (365*dpar[0]+dpar[0]/4) - biss
			+ ((dpar[1]>2?((int)((dpar[1]+1)*30.6)-63+biss)
		:((dpar[1]-1)*(63+biss))/2) + dpar[2]);
        wcs->obsdate = 2000.0 - (MJD2000 - date)/365.25;
        }
      else
/*------ Well if really no date is found */
        wcs->obsdate = 0.0;
      }

    FITSREADF(buf, "EPOCH", wcs->epoch, 2000.0);
    FITSREADF(buf, "EQUINOX", wcs->equinox, wcs->epoch);
    if (fitsread(buf, "RADESYS", str, H_STRING,T_STRING) != RETURN_OK)
      FITSREADS(buf, "RADECSYS", str,
	wcs->equinox >= 2000.0? "ICRS" : (wcs->equinox<1984.0? "FK4" : "FK5"));
    if (!strcmp(str, "ICRS"))
      wcs->radecsys = RDSYS_ICRS;
    else if (!strcmp(str, "FK5"))
      wcs->radecsys = RDSYS_FK5;
    else if (!strcmp(str, "FK4"))
      {
      if (wcs->equinox == 2000.0)
        {
        FITSREADF(buf, "EPOCH  ", wcs->equinox, 1950.0);
        FITSREADF(buf, "EQUINOX", wcs->equinox, wcs->equinox);
        }
      wcs->radecsys = RDSYS_FK4;
      warning("FK4 precession formulae not yet implemented:\n",
		"            Astrometry may be slightly inaccurate");
      }
    else if (!strcmp(str, "FK4-NO-E"))
      {
      if (wcs->equinox == 2000.0)
        {
        FITSREADF(buf, "EPOCH", wcs->equinox, 1950.0);
        FITSREADF(buf, "EQUINOX", wcs->equinox, wcs->equinox);
        }
      wcs->radecsys = RDSYS_FK4_NO_E;
      warning("FK4 precession formulae not yet implemented:\n",
		"            Astrometry may be slightly inaccurate");
      }
    else if (!strcmp(str, "GAPPT"))
      {
      wcs->radecsys = RDSYS_GAPPT;
      warning("GAPPT reference frame not yet implemented:\n",
		"            Astrometry may be slightly inaccurate");
      }
    else
      {
      warning("Using ICRS instead of unknown astrometric reference frame: ",
		str);
      wcs->radecsys = RDSYS_ICRS;
      }

/*-- Projection parameters */
    if (!strcmp(wcs->wcsprm->pcode, "TNX"))
      {
/*---- IRAF's TNX projection: decode these #$!?@#!! WAT parameters */
      if (fitsfind(buf, "WAT?????") != RETURN_ERROR)
        {
/*------ First we need to concatenate strings */
        pstr = wstr1;
        sprintf(str, "WAT1_001");
        for (j=2; fitsread(buf,str,pstr,H_STRINGS,T_STRING)==RETURN_OK; j++)
	  {
          sprintf(str, "WAT1_%03d", j);
          pstr += strlen(pstr);
	  }
        pstr = wstr2;
        sprintf(str, "WAT2_001");
        for (j=2; fitsread(buf,str,pstr,H_STRINGS,T_STRING)==RETURN_OK; j++)
	  {
          sprintf(str, "WAT2_%03d", j);
          pstr += strlen(pstr);
	  }
/*------ LONGPOLE defaulted to 180 deg if not found */
        if ((pstr = strstr(wstr1, "longpole"))
		|| (pstr = strstr(wstr2, "longpole")))
          pstr = strpbrk(pstr, "1234567890-+.");
        wcs->longpole = pstr? atof(pstr) : 999.0;
        wcs->latpole = 999.0;
/*------ RO defaulted to 180/PI if not found */
        if ((pstr = strstr(wstr1, "ro"))
		|| (pstr = strstr(wstr2, "ro")))
          pstr = strpbrk(pstr, "1234567890-+.");
        wcs->r0 = pstr? atof(pstr) : 0.0;
/*------ Read the remaining TNX parameters */
        if ((pstr = strstr(wstr1, "lngcor"))
		|| (pstr = strstr(wstr2, "lngcor")))
          wcs->tnx_lngcor = read_tnxaxis(pstr);
        if (!wcs->tnx_lngcor)
          error(EXIT_FAILURE, "*Error*: incorrect TNX parameters in ",
			filename);
        if ((pstr = strstr(wstr1, "latcor"))
		|| (pstr = strstr(wstr2, "latcor")))
          wcs->tnx_latcor = read_tnxaxis(pstr);
        if (!wcs->tnx_latcor)
          error(EXIT_FAILURE, "*Error*: incorrect TNX parameters in ",
			filename);
        }
      }
    else
      {
      if (fitsread(buf, "LONPOLE",&wcs->longpole,H_FLOAT,T_DOUBLE) != RETURN_OK)
        FITSREADF(buf, "LONGPOLE", wcs->longpole, 999.0);
      FITSREADF(buf, "LATPOLE ", wcs->latpole, 999.0);
/*---- Old convention */
      if (fitsfind(buf, "PROJP???") != RETURN_ERROR)
        for (j=0; j<10; j++)
          {
          sprintf(str, "PROJP%-3d", j);
          FITSREADF(buf, str, wcs->projp[j], 0.0);
          }
/*---- New convention */
      if (fitsfind(buf, "PV?_????") != RETURN_ERROR)
        for (l=0; l<naxis; l++)
          for (j=0; j<100; j++)
            {
            sprintf(str, "PV%d_%d ", l+1, j);
            FITSREADF(buf, str, wcs->projp[j+l*100], 0.0);
            }
      }
    }

/* Initialize other WCS structures */
  init_wcs(wcs);

/* Find the range of coordinates */
  range_wcs(wcs);
/* Invert projection corrections */
  invert_wcs(wcs);

#undef FITSREADF
#undef FITSREADI
#undef FITSREADS

  return wcs;
  }


/******* write_wcs ***********************************************************
PROTO	void write_wcs(tabstruct *tab, wcsstruct *wcs)
PURPOSE	Write WCS (World Coordinate System) info in the FITS header.
INPUT	tab structure,
	WCS structure.
OUTPUT	-.
NOTES	-.
AUTHOR	E. Bertin (IAP)
VERSION	27/11/2013
 ***/
void	write_wcs(tabstruct *tab, wcsstruct *wcs)

  {
   double	mjd;
   char		str[MAXCHARS];
   int		j, l, naxis;

  naxis = wcs->naxis;
  addkeywordto_head(tab, "BITPIX  ", "Bits per pixel");
  fitswrite(tab->headbuf, "BITPIX  ", &tab->bitpix, H_INT, T_LONG);
  addkeywordto_head(tab, "NAXIS   ", "Number of axes");
  fitswrite(tab->headbuf, "NAXIS   ", &wcs->naxis, H_INT, T_LONG);
  for (l=0; l<naxis; l++)
    {
    sprintf(str, "NAXIS%-3d", l+1);
    addkeywordto_head(tab, str, "Number of pixels along this axis");
    fitswrite(tab->headbuf, str, &wcs->naxisn[l], H_INT, T_LONG);
    }
  addkeywordto_head(tab, "EQUINOX ", "Mean equinox");
  fitswrite(tab->headbuf, "EQUINOX ", &wcs->equinox, H_FLOAT, T_DOUBLE);
  if (wcs->obsdate!=0.0)
    {
    mjd = (wcs->obsdate-2000.0)*365.25 + MJD2000;
    addkeywordto_head(tab, "MJD-OBS ", "Modified Julian date at start");
    fitswrite(tab->headbuf, "MJD-OBS ", &mjd, H_EXPO,T_DOUBLE);
    }
  addkeywordto_head(tab, "RADESYS ", "Astrometric system");
  switch(wcs->radecsys)
    {
    case RDSYS_ICRS:
      fitswrite(tab->headbuf, "RADESYS ", "ICRS", H_STRING, T_STRING);
      break;
    case RDSYS_FK5:
      fitswrite(tab->headbuf, "RADESYS ", "FK5", H_STRING, T_STRING);
      break;
    case RDSYS_FK4:
      fitswrite(tab->headbuf, "RADESYS ", "FK4", H_STRING, T_STRING);
      break;
    case RDSYS_FK4_NO_E:
      fitswrite(tab->headbuf, "RADESYS ", "FK4-NO-E", H_STRING, T_STRING);
      break;
    case RDSYS_GAPPT:
      fitswrite(tab->headbuf, "RADESYS ", "GAPPT", H_STRING, T_STRING);
      break;
    default:
      error(EXIT_FAILURE, "*Error*: unknown RADESYS type in write_wcs()", "");
    }
  for (l=0; l<naxis; l++)
    {
    sprintf(str, "CTYPE%-3d", l+1);
    addkeywordto_head(tab, str, "WCS projection type for this axis");
    fitswrite(tab->headbuf, str, wcs->ctype[l], H_STRING, T_STRING);
    sprintf(str, "CUNIT%-3d", l+1);
    addkeywordto_head(tab, str, "Axis unit");
    fitswrite(tab->headbuf, str, wcs->cunit[l], H_STRING, T_STRING);
    sprintf(str, "CRVAL%-3d", l+1);
    addkeywordto_head(tab, str, "World coordinate on this axis");
    fitswrite(tab->headbuf, str, &wcs->crval[l], H_EXPO, T_DOUBLE);
    sprintf(str, "CRPIX%-3d", l+1);
    addkeywordto_head(tab, str, "Reference pixel on this axis");
    fitswrite(tab->headbuf, str, &wcs->crpix[l], H_EXPO, T_DOUBLE);
    for (j=0; j<naxis; j++)
      {
      sprintf(str, "CD%d_%d", l+1, j+1);
      addkeywordto_head(tab, str, "Linear projection matrix");
      fitswrite(tab->headbuf, str, &wcs->cd[l*naxis+j], H_EXPO, T_DOUBLE);
      }
    for (j=0; j<100; j++)
      if (fabs(wcs->projp[j+100*l]) > TINY)
        {
        sprintf(str, "PV%d_%d", l+1, j);
        addkeywordto_head(tab, str, "Projection distortion parameter");
        fitswrite(tab->headbuf, str, &wcs->projp[j+100*l], H_EXPO, T_DOUBLE);
        }
    }

/* Update the tab data */
  readbasic_head(tab);

  return;
  }


/******* wipe_wcs ***********************************************************
PROTO	void wipe_wcs(tabstruct *tab)
PURPOSE	Remove all WCS (World Coordinate System) info in a FITS header.
INPUT	tab structure.
OUTPUT	-.
NOTES	-.
AUTHOR	E. Bertin (IAP)
VERSION	27/11/2013
 ***/
void	wipe_wcs(tabstruct *tab)

  {
  removekeywordfrom_head(tab, "CRVAL???");
  removekeywordfrom_head(tab, "CTYPE???");
  removekeywordfrom_head(tab, "CUNIT???");
  removekeywordfrom_head(tab, "CRPIX???");
  removekeywordfrom_head(tab, "CRDER???");
  removekeywordfrom_head(tab, "CSYER???");
  removekeywordfrom_head(tab, "CDELT???");
  removekeywordfrom_head(tab, "CROTA???");
  removekeywordfrom_head(tab, "CD?_????");
  removekeywordfrom_head(tab, "PROJP_??");
  removekeywordfrom_head(tab, "PV?_????");
  removekeywordfrom_head(tab, "PC?_????");
  removekeywordfrom_head(tab, "PC0??0??");
  removekeywordfrom_head(tab, "EQUINOX?");
  removekeywordfrom_head(tab, "RADESYS?");
  removekeywordfrom_head(tab, "RADECSYS");
  removekeywordfrom_head(tab, "LONPOLE?");
  removekeywordfrom_head(tab, "LONGPOLE");
  removekeywordfrom_head(tab, "LATPOLE?");
  removekeywordfrom_head(tab, "WAT?????");

  return;
  }


/******* end_wcs **************************************************************
PROTO	void end_wcs(wcsstruct *wcs)
PURPOSE	Free WCS (World Coordinate System) infos.
INPUT	WCS structure.
OUTPUT	-.
NOTES	.
AUTHOR	E. Bertin (IAP)
VERSION	24/05/2000
 ***/
void	end_wcs(wcsstruct *wcs)

  {
  if (wcs)
    {
    if (wcs->lin)
      {
      free(wcs->lin->cdelt);
      free(wcs->lin->crpix);
      free(wcs->lin->pc);
      free(wcs->lin->piximg);
      free(wcs->lin->imgpix);
      free(wcs->lin);
      }
    free(wcs->cel);
    free(wcs->prj);
    free(wcs->wcsprm);
    free_tnxaxis(wcs->tnx_lngcor);
    free_tnxaxis(wcs->tnx_latcor);
    poly_end(wcs->inv_x);
    poly_end(wcs->inv_y);
    free(wcs->projp);
    free(wcs);
    }

  return;
  }


/******* wcs_supproj *********************************************************
PROTO	int wcs_supproj(char *name)
PURPOSE	Tell if a projection system is supported or not.
INPUT	Proposed projection code name.
OUTPUT	RETURN_OK if projection is supported, RETURN_ERROR otherwise.
NOTES	-.
AUTHOR	E. Bertin (IAP)
VERSION	14/06/2012
 ***/
int	wcs_supproj(char *name)

  {
   char	projcode[28][5] =
	{"TAN", "TPV", "AZP", "SIN", "STG", "ARC", "ZPN", "ZEA", "AIR", "CYP",
	"CAR", "MER", "CEA", "COP", "COD", "COE", "COO", "BON", "PCO", "GLS",
	"PAR", "AIT", "MOL", "CSC", "QSC", "TSC", "TNX", "NONE"};

   int	i;

  for (i=0; i<28; i++)
    if (!strcmp(name, projcode[i]))
      return RETURN_OK;

  return RETURN_ERROR;
  }


/******* invert_wcs ***********************************************************
PROTO	void invert_wcs(wcsstruct *wcs)
PURPOSE	Invert WCS projection mapping (using a polynomial).
INPUT	WCS structure.
OUTPUT	-.
NOTES	.
AUTHOR	E. Bertin (IAP)
VERSION	08/03/2016
 ***/
void	invert_wcs(wcsstruct *wcs)

  {
   polystruct		*poly;
   double		pixin[NAXIS],raw[NAXIS],rawmin[NAXIS];
   double		*outpos,*outpost, *lngpos,*lngpost,
			*latpos,*latpost,
			lngstep,latstep, rawsize, epsilon;
   int			group[] = {1,1};
				/* Don't ask, this is needed by poly_init()! */
   int		i,j,lng,lat,deg, tnxflag, maxflag, maxflagflag;

/* Check first that inversion is not straightforward */
  lng = wcs->wcsprm->lng;
  lat = wcs->wcsprm->lat;
  if (!strcmp(wcs->wcsprm->pcode, "TNX"))
    tnxflag = 1;
  else if ((!strcmp(wcs->wcsprm->pcode, "TAN")
	|| !strcmp(wcs->wcsprm->pcode, "TPV"))
		&& (wcs->projp[1+lng*100] || wcs->projp[1+lat*100]))
    tnxflag = 0;
  else
    return;

/* We define x as "longitude" and y as "latitude" projections */
/* We assume that PCxx cross-terms with additional dimensions are small */
/* Sample the whole image with a regular grid */
  lngstep = wcs->naxisn[lng]/(WCS_NGRIDPOINTS-1.0);
  latstep = wcs->naxisn[lat]/(WCS_NGRIDPOINTS-1.0);
  QMALLOC(outpos, double, 2*WCS_NGRIDPOINTS2);
  QMALLOC(lngpos, double, WCS_NGRIDPOINTS2);
  QMALLOC(latpos, double, WCS_NGRIDPOINTS2);
  for (i=0; i<wcs->naxis; i++)
    raw[i] = rawmin[i] = 0.5;
  outpost = outpos;
  lngpost = lngpos;
  latpost = latpos;
  for (j=WCS_NGRIDPOINTS; j--; raw[lat]+=latstep)
    {
    raw[lng] = rawmin[lng];
    for (i=WCS_NGRIDPOINTS; i--; raw[lng]+=lngstep)
      {
      if (linrev(raw, wcs->lin, pixin))
        error(EXIT_FAILURE, "*Error*: incorrect linear conversion in ",
		wcs->wcsprm->pcode);
      *(lngpost++) = pixin[lng];
      *(latpost++) = pixin[lat];
      if (tnxflag)
        {
        *(outpost++) = pixin[lng]
			+raw_to_tnxaxis(wcs->tnx_lngcor,pixin[lng],pixin[lat]);
        *(outpost++) = pixin[lat]
			+raw_to_tnxaxis(wcs->tnx_latcor,pixin[lng],pixin[lat]);
        }
      else
        {
        raw_to_pv(wcs->prj,pixin[lng],pixin[lat], outpost, outpost+1);
        outpost += 2;
        }
      }
    }

/* Invert "longitude" */
/* Compute the extent of the pixel in reduced projected coordinates */
  linrev(rawmin, wcs->lin, pixin);
  pixin[lng] += ARCSEC/DEG;
  linfwd(pixin, wcs->lin, raw);
  rawsize = sqrt((raw[lng]-rawmin[lng])*(raw[lng]-rawmin[lng])
		+(raw[lat]-rawmin[lat])*(raw[lat]-rawmin[lat]))*DEG/ARCSEC;
  if (!rawsize)
    error(EXIT_FAILURE, "*Error*: incorrect linear conversion in ",
		wcs->wcsprm->pcode);
  epsilon = WCS_INVACCURACY/rawsize;
/* Find the lowest degree polynom */
  poly = NULL;  /* to avoid gcc -Wall warnings */
  maxflag = 1;
  for (deg=1; deg<=WCS_INVMAXDEG && maxflag; deg++)
    {
    if (deg>1)
      poly_end(poly);
    poly = poly_init(group, 2, &deg, 1);
    poly_fit(poly, outpos, lngpos, NULL, WCS_NGRIDPOINTS2, NULL, 0.0);
    maxflag = 0;
    outpost = outpos;
    lngpost = lngpos;
    for (i=WCS_NGRIDPOINTS2; i--; outpost+=2)
      if (fabs(poly_func(poly, outpost)-*(lngpost++))>epsilon)
        {
        maxflag = 1;
        break;
        }
    }
/* Now link the created structure */
  wcs->prj->inv_x = wcs->inv_x = poly;

  maxflagflag = maxflag;

/* Invert "latitude" */
/* Compute the extent of the pixel in reduced projected coordinates */
  linrev(rawmin, wcs->lin, pixin);
  pixin[lat] += ARCSEC/DEG;
  linfwd(pixin, wcs->lin, raw);
  rawsize = sqrt((raw[lng]-rawmin[lng])*(raw[lng]-rawmin[lng])
		+(raw[lat]-rawmin[lat])*(raw[lat]-rawmin[lat]))*DEG/ARCSEC;
  if (!rawsize)
    error(EXIT_FAILURE, "*Error*: incorrect linear conversion in ",
		wcs->wcsprm->pcode);
  epsilon = WCS_INVACCURACY/rawsize;
/* Find the lowest degree polynom */
  maxflag = 1;
  for (deg=1; deg<=WCS_INVMAXDEG && maxflag; deg++)
    {
    if (deg>1)
      poly_end(poly);
    poly = poly_init(group, 2, &deg, 1);
    poly_fit(poly, outpos, latpos, NULL, WCS_NGRIDPOINTS2, NULL, 0.0);
    maxflag = 0;
    outpost = outpos;
    latpost = latpos;
    for (i=WCS_NGRIDPOINTS2; i--; outpost+=2)
      if (fabs(poly_func(poly, outpost)-*(latpost++))>epsilon)
        {
        maxflag = 1;
        break;
        }
    }
/* Now link the created structure */
  wcs->prj->inv_y = wcs->inv_y = poly;

  maxflagflag |= maxflag;

  if (maxflagflag)
    warning("Significant inaccuracy likely to occur in projection","");

/* Free memory */
  free(outpos);
  free(lngpos);
  free(latpos);

  return;
  }


/******* range_wcs ***********************************************************
PROTO	void range_wcs(wcsstruct *wcs)
PURPOSE	Find roughly the range of WCS coordinates on all axes,
	and typical pixel scales.
INPUT	WCS structure.
OUTPUT	-.
NOTES	.
AUTHOR	E. Bertin (IAP)
VERSION	24/08/2010
 ***/
void	range_wcs(wcsstruct *wcs)

  {
   double		step[NAXIS], raw[NAXIS], rawmin[NAXIS],
			world[NAXIS], world2[NAXIS];
   double		*worldmin, *worldmax, *scale, *worldc,
			rad, radmax, lc;
   int			linecount[NAXIS];
   int			i,j, naxis, npoints, lng,lat;

  naxis = wcs->naxis;

/* World range */
  npoints = 1;
  worldmin = wcs->wcsmin;
  worldmax = wcs->wcsmax;
/* First, find the center and use it as a reference point for lng */
  lng = wcs->lng;
  lat = wcs->lat;
  for (i=0; i<naxis; i++)
    raw[i] = (wcs->naxisn[i]+1.0)/2.0;
  if (raw_to_wcs(wcs, raw, world))
    {
/*-- Oops no mapping there! So explore the image in an increasingly large */
/*-- domain to find  a better "center" (now we know there must be angular */