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