/*
 				fitswcs.c

*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
*
*	Part of:	LDACTools+
*
*	Author:		E.BERTIN (IAP)
*
*	Contents:       Read and write WCS header info.
*
*	Last modify:	04/03/2001
*
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
*/
#include	<math.h>
#include	<stdio.h>
#include	<stdlib.h>
#include	<string.h>

#include	"fitscat_defs.h"
#include	"fitscat.h"
#include	"fitswcs.h"
#include	"wcscelsys.h"
#include	"wcs/wcs.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	18/08/2000
 ***/
wcsstruct	*copy_wcs(wcsstruct *wcsin)

  {
   wcsstruct	*wcs;

/* Copy the basic stuff */
  QMEMCPY(wcsin, wcs, wcsstruct, 1);

  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);

/* The PROJP WCS parameters */
  QMEMCPY(wcsin->projp, wcs->projp, double, wcs->naxis*100);
/* 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	13/04/2000
 ***/
void	init_wcs(wcsstruct *wcs)

  {
   int		l,naxis;

  QCALLOC(wcs->lin, struct linprm, 1);
  QCALLOC(wcs->cel, struct celprm, 1);
  QCALLOC(wcs->prj, struct prjprm, 1);
  QCALLOC(wcs->lin->cdelt, double, 2);
  QCALLOC(wcs->lin->crpix, double, 2);
  QCALLOC(wcs->lin->pc, double, 4);
  wcs->lin->flag = wcs->cel->flag = wcs->prj->flag = 0;
  naxis = wcs->lin->naxis = wcs->naxis;

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

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

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

/* celprm structure */
  wcs->cel->ref[0] = wcs->crval[wcs->wcsprm->lng];
  wcs->cel->ref[1] = wcs->crval[wcs->wcsprm->lat];
  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;
  for (l=0; l<100; l++)
    {
    wcs->prj->p[l] = wcs->projp[l+wcs->wcsprm->lng*100];
    wcs->prj->p[l+100] = wcs->projp[l+wcs->wcsprm->lat*100];
    }

/* 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	14/01/2001
 ***/
void	init_wcscelsys(wcsstruct *wcs)

  {
  double	*mat,
		a0,d0,ap,dp,ap2,x,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 = 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);
/* 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	25/08/2000
 ***/
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	FITSREADS(buf, k, str, def) \
		{if (fitsread(buf,k,str, H_STRING,T_STRING) != RETURN_OK) \
		   strcpy(str, (def)); \
		}
   static char	str[MAXCHARS];
   static char	wstr1[TNX_MAXCHARS], wstr2[TNX_MAXCHARS];

   wcsstruct	*wcs;
   catstruct	*cat;
   double	drota;
   int		j, l, naxis;
   char		*buf;

  buf = tab->headbuf;
  if (!(cat = tab->cat))
    error(EXIT_FAILURE, "*Internal Error*: Table has no parent catalog","!");

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

  QCALLOC(wcs, wcsstruct, 1);
  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 ",
	cat->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, "PC00?00?")!=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);
/*-- Precession date (defined from Ephemerides du Bureau des Longitudes) */
/*-- in Julian years from 2000.0 */
    if (date>0.0)
      wcs->equinox = 2000.0 - (MJD2000 - date)/365.25;
    else
      {
/*---- Search for an observation date expressed in "civil" 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(l?NULL:str,"/- ")); 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->equinox = 2000.0 - (MJD2000 - date)/365.25;
        }
      else
/*------ Well if really no date is found */
        wcs->equinox = 2000.0;
      }

    FITSREADF(buf, "EPOCH", wcs->equinox, wcs->equinox);
    FITSREADF(buf, "EQUINOX", wcs->equinox, wcs->equinox);
    FITSREADS(buf, "RADECSYS", str, wcs->equinox<1984.0?"FK4":"FK5");
    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 FK5 instead of unknown astrometric reference frame: ",
		str);
      wcs->radecsys = RDSYS_FK5;
      }

/*-- 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 ",
			cat->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 ",
			cat->filename);
        }
      }
    else
      {
      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);

  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	08/11/2000
 ***/
void	write_wcs(tabstruct *tab, wcsstruct *wcs)

  {
#define	FITSADD(k, com) \
		{if (!((++n)%36)) \
		   { \
		   tab->headnblock++; \
		   QREALLOC(tab->headbuf, char, tab->headnblock*FBSIZE); \
		   if (fitsadd(tab->headbuf, k, com) == RETURN_OK) \
		     { \
                     n--; \
		     tab->headnblock--; \
		     QREALLOC(tab->headbuf, char, tab->headnblock*FBSIZE); \
		     } \
		   } \
                 else \
                   if (fitsadd(tab->headbuf, k, com) == RETURN_OK) \
                     n--; \
		}

   static char	str[MAXCHARS];
   int		n, l, naxis, dummy;

  n = fitsfind(tab->headbuf, "END     ");
  naxis = wcs->naxis;

  FITSADD("BITPIX  ", "Bits per pixel");
  dummy = -32;
  fitswrite(tab->headbuf, "BITPIX  ", &dummy, H_INT, T_LONG);
  FITSADD("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);
    FITSADD(str, "Number of pixels along this axis");
    fitswrite(tab->headbuf, str, &wcs->naxisn[l], H_INT, T_LONG);
    }
  FITSADD("EQUINOX ", "Mean equinox");
  fitswrite(tab->headbuf, "EQUINOX ", &wcs->equinox, H_FLOAT, T_DOUBLE);
  FITSADD("RADECSYS ", "Astrometric system");
  switch(wcs->radecsys)
    {
    case RDSYS_FK5:
      fitswrite(tab->headbuf, "RADECSYS", "FK5", H_STRING, T_STRING);
      break;
    case RDSYS_FK4:
      fitswrite(tab->headbuf, "RADECSYS", "FK4", H_STRING, T_STRING);
      break;
    case RDSYS_FK4_NO_E:
      fitswrite(tab->headbuf, "RADECSYS", "FK4-NO-E", H_STRING, T_STRING);
      break;
    case RDSYS_GAPPT:
      fitswrite(tab->headbuf, "RADECSYS", "GAPPT", H_STRING, T_STRING);
      break;
    default:
      error(EXIT_FAILURE, "*Error*: unknown RADECSYS type in write_wcs()", "");
    }
  for (l=0; l<naxis; l++)
    {
    sprintf(str, "CTYPE%-3d", l+1);
    FITSADD(str, "WCS projection type for this axis");
    fitswrite(tab->headbuf, str, wcs->ctype[l], H_STRING, T_STRING);
    sprintf(str, "CUNIT%-3d", l+1);
    FITSADD(str, "Axis unit");
    fitswrite(tab->headbuf, str, wcs->cunit[l], H_STRING, T_STRING);
    sprintf(str, "CRVAL%-3d", l+1);
    FITSADD(str, "World coordinate on this axis");
    fitswrite(tab->headbuf, str, &wcs->crval[l], H_EXPO, T_DOUBLE);
    sprintf(str, "CRPIX%-3d", l+1);
    FITSADD(str, "Reference pixel on this axis");
    fitswrite(tab->headbuf, str, &wcs->crpix[l], H_EXPO, T_DOUBLE);
    sprintf(str, "CDELT%-3d", l+1);
    FITSADD(str, "Pixel step along this axis");
    fitswrite(tab->headbuf, str, &wcs->cdelt[l], H_EXPO, T_DOUBLE);
/*
    for (j=0; j<naxis; j++)
      {
      sprintf(str, "CD%d_%d", l+1, j+1);
      FITSADD(str, "Linear projection matrix");
      fitswrite(tab->headbuf, str, &wcs->cd[l*naxis+j], H_EXPO, T_DOUBLE);
      }
*/
    }

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

  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	24/05/2000
 ***/
int	wcs_supproj(char *name)

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

  for (i=0; i<26; 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	01/09/2000
 ***/
void	invert_wcs(wcsstruct *wcs)

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

/* 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")
		&& (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]-1)/(WCS_NGRIDPOINTS-1);
  latstep = (wcs->naxisn[lat]-1)/(WCS_NGRIDPOINTS-1);
  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.0;
  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] += 1.0;
  linfwd(pixin, wcs->lin, raw);
  rawsize = sqrt((raw[lng]-rawmin[lng])*(raw[lng]-rawmin[lng])
		+(raw[lat]-rawmin[lat])*(raw[lat]-rawmin[lat]));
  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, lngpos, NULL, WCS_NGRIDPOINTS2, NULL);
    maxflag = 0;
    outpost = outpos;
    lngpost = lngpos;
    for (i=WCS_NGRIDPOINTS2; i--; outpost+=2)
      if (fabs(poly_func(poly, outpost)-*(lngpost++))>epsilon)
        {
        maxflag = 1;
        break;
        }
    }
  if (maxflag)
    warning("Significant inaccuracy likely to occur in projection","");
/* Now link the created structure */
  wcs->prj->inv_x = wcs->inv_x = poly;

/* Invert "latitude" */
/* Compute the extent of the pixel in reduced projected coordinates */
  linrev(rawmin, wcs->lin, pixin);
  pixin[lat] += 1.0;
  linfwd(pixin, wcs->lin, raw);
  rawsize = sqrt((raw[lng]-rawmin[lng])*(raw[lng]-rawmin[lng])
		+(raw[lat]-rawmin[lat])*(raw[lat]-rawmin[lat]));
  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);
    maxflag = 0;
    outpost = outpos;
    latpost = latpos;
    for (i=WCS_NGRIDPOINTS2; i--; outpost+=2)
      if (fabs(poly_func(poly, outpost)-*(latpost++))>epsilon)
        {
        maxflag = 1;
        break;
        }
    }
  if (maxflag)
    warning("Significant inaccuracy likely to occur in projection","");
/* Now link the created structure */
  wcs->prj->inv_y = wcs->inv_y = poly;

/* 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	04/03/2001
 ***/
void	range_wcs(wcsstruct *wcs)

  {
   static double	step[NAXIS], raw[NAXIS], rawmin[NAXIS],
			world[NAXIS], world2[NAXIS];
   double		*worldmin, *worldmax, *scale,
			lc;
   static 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 */
/*-- coordinates) */
    for (j=0; j<100; j++)
      {
      for (i=0; i<naxis; i++)
        raw[i] += wcs->naxisn[i]/100.0*(0.5-(double)rand()/RAND_MAX);      
      if (!raw_to_wcs(wcs, raw, world))
        break;
      }
    }

  if (lng!=lat)
    lc = fmod(world[lng]+180.0, 360.0);
  else
    lng = -1;

/* Pixel scales at image center */
  scale = wcs->wcsscale;
  for (i=0; i<naxis; i++)
    {
    if ((i==lng || i==lat) && lng!=lat)
      scale[i] = sqrt(wcs_scale(wcs, raw));
    else
      {
      raw[i] += 1.0;
      raw_to_wcs(wcs, raw, world2);
      scale[i] = (world2[i] - world[i]);
      raw[i] -= 1.0;
      }
    wcs->wcsscalepos[i] = world[i];
    }

/* Find "World limits" */
  for (i=0; i<naxis; i++)
    {
    raw[i] = rawmin[i] = 1.0;
    step[i] = (wcs->naxisn[i]-1)/(WCS_NRANGEPOINTS-1);
    npoints *= WCS_NRANGEPOINTS;
    worldmax[i] = -(worldmin[i] = 1e31);
    }

  for (j=npoints; j--;)
    {
    raw_to_wcs(wcs, raw, world);
    for (i=0; i<naxis; i++)
      {
/*---- Handle longitudes around 0 */
      if (i==lng && world[i]>lc)
        world[i] -= 359.9999;
      if (world[i]<worldmin[i])
        worldmin[i] = world[i];
      if (world[i]>worldmax[i])
        worldmax[i] = world[i];
      }

    for (i=0; i<naxis; i++)
      {
      raw[i] += step[i];
      if (++linecount[i]<WCS_NRANGEPOINTS)
        break;
      else
        {
        linecount[i] = 0;       /* No need to initialize it to 0! */
        raw[i] = rawmin[i];
        }
      }
    }
  
  if (lng!=lat)
    {
    if (worldmax[lat]<-90.0)
      worldmax[lat] = -90.0;
    if (worldmax[lat]>90.0)
      worldmax[lat] = 90.0;
    }

  return;
  }


/******* frame_wcs ***********************************************************
PROTO	void frame_wcs(wcsstruct *wcsin, wcsstruct *wcsout)
PURPOSE	Find the x and y limits of an input frame in an output image.
INPUT	WCS structure of the input frame,
	WCS structure of the output frame.
OUTPUT	-.
NOTES	.
AUTHOR	E. Bertin (IAP)
VERSION	11/11/2000
 ***/
void	frame_wcs(wcsstruct *wcsin, wcsstruct *wcsout)

  {
   static double	rawin[NAXIS], rawout[NAXIS], world[NAXIS];
   static int		linecount[NAXIS];
   int			*min, *max,
			i,j, naxis, npoints, out;

  naxis = wcsin->naxis;

/* World range */
  npoints = 1;
  min = wcsin->outmin;
  max = wcsin->outmax;
  for (i=0; i<naxis; i++)
    {
    rawin[i] = 0.5;	/* Lower pixel limits */
    npoints *= WCS_NRANGEPOINTS;
    max[i] = -(min[i] = 1<<30);
    }

  for (j=npoints; j--;)
    {
    if (!raw_to_wcs(wcsin, rawin, world) && !wcs_to_raw(wcsout, world, rawout))
      for (i=0; i<naxis; i++)
        {
        if ((out=(int)(rawout[i]+0.499))<min[i])
          min[i] = out;
        if (out>max[i])
          max[i] = out;
        }
    for (i=0; i<naxis; i++)
      {
      rawin[i] = 0.5 + 0.5*wcsin->naxisn[i]
	*(1-cos(PI*(linecount[i]+1.0)/(WCS_NRANGEPOINTS-1)));
      if (++linecount[i]<WCS_NRANGEPOINTS)
        break;
      else
        {
        linecount[i] = 0;       /* No need to initialize it to 0! */
        rawin[i] =  0.5;
        }
      }
    }

/* Just add a little margin, in case of... */
  for (i=0; i<naxis; i++)
    {
    if ((min[i]-=2) < 0)
      min[i] = 0;
    if ((max[i]+=2) >= wcsout->naxisn[i] - 1)
      max[i] = wcsout->naxisn[i] - 1;
    }

  return;
  }


/******* raw_to_wcs ***********************************************************
PROTO	void raw_to_wcs(wcsstruct *, double *, double *)
PURPOSE	Convert raw (pixel) coordinates to WCS (World Coordinate System).
INPUT	WCS structure,
	Pointer to the array of input coordinates,
	Pointer to the array of output coordinates.
OUTPUT	RETURN_OK if mapping successful, RETURN_ERROR otherwise.
NOTES	-.
AUTHOR	E. Bertin (IAP)
VERSION	14/01/2001
 ***/
int	raw_to_wcs(wcsstruct *wcs, double *pixpos, double *wcspos)

  {
   double	*mat,
		imgcrd[2],
		phi,theta, a2,d2,sd2,cd2cp,sd,x,y;
   int		i,lng,lat;

  if (wcsrev((const char(*)[9])wcs->ctype, wcs->wcsprm, pixpos,
	wcs->lin,imgcrd, wcs->prj, &phi, &theta, wcs->crval, wcs->cel, wcspos))
    {
    for (i=0; i<wcs->naxis; i++)
      wcspos[i] = WCS_NOCOORD;
    return RETURN_ERROR;
    }

/* If needed, convert from a different coordinate system to equatorial */
  if (wcs->celsysconvflag)
    {
    mat = wcs->celsysmat;
    a2 = wcspos[lng = wcs->wcsprm->lng]*DEG - mat[1];
    d2 = wcspos[lat = wcs->wcsprm->lat]*DEG;
/*-- A bit of spherical trigonometry...
/*-- Compute the latitude... */
    sd2 = sin(d2);
    cd2cp = cos(d2)*mat[2];
    sd = sd2*mat[3]-cd2cp*cos(a2);
/*-- ...and the longitude */
    y = cd2cp*sin(a2);
    x = sd2 - sd*mat[3];
    wcspos[lng] = fmod((atan2(y,x) + mat[0])/DEG+360.0, 360.0);
    wcspos[lat] = asin(sd)/DEG;
    }

  return RETURN_OK;
  }


/******* wcs_to_raw ***********************************************************
PROTO	void wcs_to_raw(wcsstruct *, double *, double *)
PURPOSE	Convert WCS (World Coordinate System) coords to raw (pixel) coords.
INPUT	WCS structure,
	Pointer to the array of input coordinates,
	Pointer to the array of output coordinates.
OUTPUT	RETURN_OK if mapping successful, RETURN_ERROR otherwise.
NOTES	-.
AUTHOR	E. Bertin (IAP)
VERSION	13/01/2001
 ***/
int	wcs_to_raw(wcsstruct *wcs, double *wcspos, double *pixpos)

  {
   double	*mat,
		imgcrd[2],
		phi,theta, a,d,sd,cdcp,sd2,x,y;
   int		i,lng,lat;

/* If needed, convert to a coordinate system different from equatorial */
  if (wcs->celsysconvflag)
    {
    mat = wcs->celsysmat;
    a = wcspos[lng = wcs->wcsprm->lng]*DEG - mat[0];
    d = wcspos[lat = wcs->wcsprm->lat]*DEG;
/*-- A bit of spherical trigonometry...
/*-- Compute the latitude... */
    sd = sin(d);
    cdcp = cos(d)*mat[2];
    sd2 = sd*mat[3]+cdcp*cos(a);
/*-- ...and the longitude */
    y = cdcp*sin(a);
    x = sd2*mat[3]-sd;
    wcspos[lng] = fmod((atan2(y,x) + mat[1])/DEG+360.0, 360.0);
    wcspos[lat] = asin(sd2)/DEG;
    }

  if (wcsfwd((const char(*)[9])wcs->ctype,wcs->wcsprm,wcspos,
	wcs->crval, wcs->cel,&phi,&theta,wcs->prj, imgcrd,wcs->lin,pixpos))
    {
    for (i=0; i<wcs->naxis; i++)
      pixpos[i] = WCS_NOCOORD;
    return RETURN_ERROR;
    }

  return RETURN_OK;
  }


/******* wcs_scale ***********************************************************
PROTO	double wcs_scale(wcsstruct *wcs, double *pixpos)
PURPOSE	Compute the sky area equivalent to a local pixel.
INPUT	WCS structure,
	Pointer to the array of local raw coordinates,
OUTPUT	-.
NOTES	-.
AUTHOR	E. Bertin (IAP)
VERSION	10/11/2000
 ***/
double	wcs_scale(wcsstruct *wcs, double *pixpos)

  {
  static double	wcspos[NAXIS], wcspos1[NAXIS], wcspos2[NAXIS];
  double	dpos1,dpos2;
  int		lng, lat;

  lng = wcs->lng;
  lat = wcs->lat;
  if (raw_to_wcs(wcs, pixpos, wcspos))
    return 0.0;
/* Compute pixel scale */
  pixpos[lng] += 1.0;
  if (raw_to_wcs(wcs, pixpos, wcspos1))
    return 0.0;
  pixpos[lng] -= 1.0;
  pixpos[lat] += 1.0;
  if (raw_to_wcs(wcs, pixpos, wcspos2))
    return 0.0;
  dpos1 = wcspos1[lng]-wcspos[lng];
  if (dpos1>180.0)
    dpos1 -= 360.0;
  else if (dpos1<-180.0)
    dpos1 += 360.0;
  dpos2 = wcspos2[lng]-wcspos[lng];
  if (dpos2>180.0)
    dpos2 -= 360.0;
  else if (dpos2<-180.0)
    dpos2 += 360.0;
  return fabs((dpos1*(wcspos2[lat]-wcspos[lat])
		-(wcspos1[lat]-wcspos[lat])*dpos2)
		*cos(wcspos[lat]*DEG));
  }


/********************************* precess ***********************************/
/*
precess equatorial coordinates according to the equinox (from Ephemerides du
Bureau des Longitudes 1992). Epoch for coordinates should be J2000
(FK5 system).
*/
void	precess(double yearin, double alphain, double deltain,
		double yearout, double *alphaout, double *deltaout)

  {
   double	dzeta,theta,z, t1,t1t1, t2,t2t2,t2t2t2,
		cddsadz, cddcadz, cdd, sdd, adz, cdin,sdin,ct,st,caindz;

  alphain *= DEG;
  deltain *= DEG;

  t1 = (yearin - 2000.0)/1000.0;
  t2 = (yearout - yearin)/1000.0;
  t1t1 = t1*t1;
  t2t2t2 = (t2t2 = t2*t2)*t2;
  theta = (97171.735e-06 - 413.691e-06*t1 - 1.052e-06 * t1t1) * t2
	+ (-206.846e-06 - 1.052e-06*t1) * t2t2 - 202.812e-06 * t2t2t2;
  dzeta = (111808.609e-06 + 677.071e-06*t1 - 0.674e-06 * t1t1) * t2
	+ (146.356e-06 - 1.673e-06*t1) * t2t2 + 87.257e-06 * t2t2t2;
  z = (111808.609e-06 +677.071e-06*t1 - 0.674e-06 * t1t1) * t2
	+ (530.716e-06 + 0.320e-06*t1) * t2t2 + 88.251e-06 * t2t2t2;
  cddsadz = (cdin=cos(deltain)) * sin(alphain+dzeta);
  cddcadz = -(sdin=sin(deltain))*(st=sin(theta))
	+cdin*(ct=cos(theta))*(caindz=cos(alphain+dzeta));
  sdd = sdin*ct + cdin*st*caindz;
  cdd = cos(*deltaout = asin(sdd));
  adz = asin(cddsadz/cdd);
  if (cddcadz<0.0)
    adz = PI - adz;
  if (adz<0.0)
    adz += 2.0*PI;
  adz += z;
  *alphaout = adz/DEG;
  *deltaout /= DEG;

  return;
  }


/*********************************** j2b *************************************/
/*
conver equatorial coordinates from equinox and epoch J2000 to equinox and
epoch B1950 for extragalactic sources (from Aoki et al. 1983, after
inversion of their matrix and some custom arrangements).
*/
void    j2b(double yearobs, double alphain, double deltain,
	double *alphaout, double *deltaout)
  {
   int			i,j;
   static double	a[3] = {-1.62557e-6, -0.31919e-6, -0.13843e-6},
                        ap[3] = {1.245e-3, -1.580e-3, -0.659e-3},
                        m[6][6] = {
  { 0.9999256794678425,    0.01118148281196562,   0.004859003848996022,
   -2.423898417033081e-06,-2.710547600126671e-08,-1.177738063266745e-08},
  {-0.01118148272969232,   0.9999374849247641,   -2.717708936468247e-05,
    2.710547578707874e-08,-2.423927042585208e-06, 6.588254898401055e-11},
  {-0.00485900399622881,  -2.715579322970546e-05, 0.999988194643078,
    1.177738102358923e-08, 6.582788892816657e-11,-2.424049920613325e-06},
  {-0.0005508458576414713, 0.2384844384742432,   -0.4356144527773499,
    0.9999043171308133,    0.01118145410120206,   0.004858518651645554},
  {-0.2385354433560954,   -0.002664266996872802,  0.01225282765749546,
   -0.01118145417187502,   0.9999161290795875,   -2.717034576263522e-05},
  { 0.4357269351676567,   -0.008536768476441086,  0.002113420799663768,
   -0.004858518477064975, -2.715994547222661e-05, 0.9999668385070383}},
			a1[3], r[3], ro[3], r1[3], r2[3], v1[3], v[3];
   static double	cai, sai, cdi, sdi, dotp, rmod, alpha, delta, t1;

/* Convert Julian years from J2000.0 to tropic centuries from B1950.0 */
  t1 = ((yearobs - 2000.0) + (MJD2000 - MJD1950)/365.25)*JU2TROP/100.0;
  alphain *= DEG;
  deltain *= DEG;
  cai = cos(alphain);
  sai = sin(alphain);
  cdi = cos(deltain);
  sdi = sin(deltain);
  r[0] = cdi*cai;
  r[1] = cdi*sai;
  r[2] = sdi;
  for (i=0; i<3; i++)
    v[i] = r2[i] = v1[i] = 0.0;
  for (j=0; j<6; j++)
    for (i=0; i<6; i++)
      if (j<3)
        r2[j] += m[j][i]*(i<3?r[i]:v[i-3]);
      else
        v1[j-3] += m[j][i]*(i<3?r[i]:v[i-3]);

  for (i=0; i<3; i++)
    r1[i] = r2[i]+v1[i]*ARCSEC*t1;

  dotp = 0.0;
  for (i=0; i<3; i++)
    {
    a1[i] = a[i]+ap[i]*ARCSEC*t1;
    dotp += a1[i]*(r1[i]+a1[i]);
    }
  dotp = 2.0/(sqrt(1+4.0*dotp)+1.0);
  rmod = 0.0;
  for (i=0; i<3; i++)
    {
    ro[i] = dotp*(r1[i]+a1[i]);
    rmod += ro[i]*ro[i];
    }
  rmod = sqrt(rmod);
  delta = asin(ro[2]/rmod);
  alpha = acos(ro[0]/cos(delta)/rmod);
  if (ro[1]<0)
    alpha = 2.0*PI - alpha;
  *alphaout = alpha/DEG;
  *deltaout = delta/DEG;

  return;
  }


/******************************** degtosexal *********************************/
/*
Convert degrees to hh mm ss.xx alpha coordinates.
*/
char    *degtosexal(double alpha)

  {
   static char	str[12];
   int		hh, mm;
   double	ss;

  if (alpha>=0.0 && alpha <360.0)
    {
    hh = (int)(alpha/15.0);
    mm = (int)(60.0*(alpha/15.0 - hh));
    ss = 60.0*(60.0*(alpha/15.0 - hh) - mm);
    }
  else
    hh = mm = ss = 0.0;
  sprintf(str,"%2d:%02d:%05.2f", hh, mm, ss);

  return str;
  }


/******************************** degtosexde *********************************/
/*
Convert degrees to dd dm ds.x delta coordinates.
*/
char    *degtosexde(double delta)

  {
   static char	str[12];
   char		sign;
   double	ds;
   int		dd, dm;

  sign = delta<0.0?'-':'+';
  delta = fabs(delta);
  if (delta>=-90.0 && delta <=90.0)
    {
    dd = (int)delta;
    dm = (int)(60.0*(delta - dd));
    ds = 60.0*fabs(60.0*(delta - dd) - dm);
    }
  else
    dd = dm = ds = 0.0;
  sprintf(str,"%c%02d:%02d:%04.1f", sign, dd, dm, ds);
  return str;
  }


/******************************** sextodegal *********************************/
/*
Convert hh mm ss.xxx alpha coordinates to degrees.
*/
double  sextodegal(char *hms)
  {
   double	val;

  val = atof(strtok(hms, ": \t"))*15.0;		/* Hours */
  val += atof(strtok(NULL, ": \t"))/4.0;	/* Minutes */
  val += atof(strtok(NULL, ": \t"))/240.0;	/* Seconds */

  return val;
  }


/******************************** sextodegde *********************************/
/*
Convert dd dm ds.xxx alpha coordinates to degrees.
*/
double  sextodegde(char *dms)
  {
   double	val, sgn;

  sgn = (val = atof(strtok(dms, ": \t")))<0.0? -1.0:1.0;/* Degrees */
  val += atof(strtok(NULL, ": \t"))*sgn/60.0;		/* Minutes */
  val += atof(strtok(NULL, ": \t"))*sgn/3600.0;		/* Seconds */

  return val;
  }