fitswcs.c

/*-- 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 = world[lng];
  else
    {
    lc = 0.0;   /* to avoid gcc -Wall warnings */
    lng = -1;
    }

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


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

  radmax = 0.0;
  worldc = wcs->wcsscalepos;

  for (j=npoints; j--;)
    {
    raw_to_wcs(wcs, raw, world);
/*-- Compute maximum distance to center */
    if ((rad=wcs_dist(wcs, world, worldc)) > radmax)
      radmax = rad;
    for (i=0; i<naxis; i++)
      {
/*---- Handle longitudes around 0 */
      if (i==lng)
        {
        world[i] -= lc;
        if (world[i]>180.0)
          world[i] -= 360.0;
        else if (world[i] <= -180.0)
          world[i] += 360.0;
        }
      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];
        }
      }
    }

  wcs->wcsmaxradius = radmax;

  if (lng!=lat)
    {
    worldmin[lng] = fmod_0_p360(worldmin[lng]+lc);
    worldmax[lng] = fmod_0_p360(worldmax[lng]+lc);
    if (worldmax[lat]<-90.0)
      worldmax[lat] = -90.0;
    if (worldmax[lat]>90.0)
      worldmax[lat] = 90.0;
    }

  return;
  }


/******* frame_wcs ***********************************************************
PROTO	int 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	1 if frames overlap, 0 otherwise.
NOTES	-.
AUTHOR	E. Bertin (IAP)
VERSION	26/03/2012
 ***/
int	frame_wcs(wcsstruct *wcsin, wcsstruct *wcsout)

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

  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);
    linecount[i] = 0;
    }

/* Check if lng and lat are swapped between in and out wcs (vicious idea!) */
  swapflag = (((wcsin->lng != wcsout->lng) || (wcsin->lat != wcsout->lat))
	&& (wcsin->lng != wcsin->lat) && (wcsout->lng != wcsout->lat));

  for (j=npoints; j--;)
    {
    if (!raw_to_wcs(wcsin, rawin, world))
      {
      if (swapflag)
        {
        worldc = world[wcsout->lat];
        world[wcsout->lat] = world[wcsin->lat];
        world[wcsin->lat] = worldc;
        }
      if (!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;
        }
      }
    }

/* Add a little margin, just in case... */
  for (i=0; i<naxis; i++)
    {
    if (min[i]>-2147483647)
      min[i] -= 2;
    if (max[i]>2147483647)
      max[i] += 2;
    }

/* Check overlap */
  overlapflag = 1;
  for (i=0; i<naxis; i++)
    if (min[i]>wcsout->naxisn[i] || max[i]<0)
      {
      overlapflag = 0;
      break;
      }

  return overlapflag;
  }


/******* reaxe_wcs ***********************************************************
PROTO	int reaxe_wcs(wcsstruct *wcs, int lng, int lat)
PURPOSE	Reformulate a wcs structure to match lng and lat axis indices
INPUT	WCS structure,
	longitude index,
	latitude index.
OUTPUT	RETURN_OK if successful, RETURN_ERROR otherwise.
NOTES	.
AUTHOR	E. Bertin (IAP)
VERSION	20/11/2003
 ***/
int	reaxe_wcs(wcsstruct *wcs, int lng, int lat)

  {
   char		strlng[80], strlat[80];
   double	dlng,dlat,dlng2,dlat2;
   int		l, ilng,ilat,olng,olat, naxis;

  olng = wcs->lng;
  olat = wcs->lat;
  if (lng<0 || lat<0 || olng<0 || olat<0)
    return RETURN_ERROR;

  ilng = wcs->naxisn[olng];
  ilat = wcs->naxisn[olat];
  wcs->naxisn[lng] = ilng;
  wcs->naxisn[lat] = ilat;
  strcpy(strlng, wcs->ctype[olng]);
  strcpy(strlat, wcs->ctype[olat]);
  strcpy(wcs->ctype[lng], strlng);
  strcpy(wcs->ctype[lat], strlat);
  dlng = wcs->crval[olng];
  dlat = wcs->crval[olat];
  wcs->crval[lng] = dlng;
  wcs->crval[lat] = dlat;
  naxis = wcs->naxis;
  dlng =  wcs->cd[olng+olng*naxis];
  dlng2 = wcs->cd[olng+olat*naxis];
  dlat =  wcs->cd[olat+olat*naxis];
  dlat2 = wcs->cd[olat+olng*naxis];
  wcs->cd[lng+lng*naxis] = dlng2;
  wcs->cd[lng+lat*naxis] = dlng;
  wcs->cd[lat+lat*naxis] = dlat2;
  wcs->cd[lat+lng*naxis] = dlat;
  for (l=0; l<100; l++)
    {
    dlng = wcs->projp[l+olng*100];
    dlat = wcs->projp[l+olat*100];
    wcs->projp[l+lng*100] = dlng;
    wcs->projp[l+lat*100] = dlat;
    }

/*-- Reinitialize wcs */
    wcsset(wcs->naxis,(const char(*)[9])wcs->ctype, wcs->wcsprm);

/*-- Initialize other WCS structures */
    init_wcs(wcs);
/*-- Find the range of coordinates */
    range_wcs(wcs);

  return RETURN_OK;
  }


/******* celsys_to_eq *********************************************************
PROTO	int celsys_to_eq(wcsstruct *wcs, double *wcspos)
PURPOSE	Convert arbitrary celestial coordinates to equatorial.
INPUT	WCS structure,
	Coordinate vector.
OUTPUT	RETURN_OK if mapping successful, RETURN_ERROR otherwise.
NOTES	-.
AUTHOR	E. Bertin (IAP)
VERSION	08/02/2007
 ***/
int	celsys_to_eq(wcsstruct *wcs, double *wcspos)

  {
   double	*mat,
		a2,d2,sd2,cd2cp,sd,x,y;
   int		lng, lat;

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


/******* eq_to_celsys *********************************************************
PROTO	int eq_to_celsys(wcsstruct *wcs, double *wcspos)
PURPOSE	Convert equatorial to arbitrary celestial coordinates.
INPUT	WCS structure,
	Coordinate vector.
OUTPUT	RETURN_OK if mapping successful, RETURN_ERROR otherwise.
NOTES	-.
AUTHOR	E. Bertin (IAP)
VERSION	08/02/2007
 ***/
int	eq_to_celsys(wcsstruct *wcs, double *wcspos)

  {
   double	*mat,
		a,d,sd2,cdcp,sd,x,y;
   int		lng, lat;

  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;

  return RETURN_OK;
  }


/******* raw_to_wcs ***********************************************************
PROTO	int 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	08/02/2007
 ***/
int	raw_to_wcs(wcsstruct *wcs, double *pixpos, double *wcspos)

  {
   double	imgcrd[NAXIS],
		phi,theta;
   int		i;

  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)
    celsys_to_eq(wcs, wcspos);

  return RETURN_OK;
  }


/******* wcs_to_raw ***********************************************************
PROTO	int 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	08/02/2007
 ***/
int	wcs_to_raw(wcsstruct *wcs, double *wcspos, double *pixpos)

  {
   double	imgcrd[NAXIS],
		phi,theta;
   int		i;

/* If needed, convert to a coordinate system different from equatorial */
  if (wcs->celsysconvflag)
    eq_to_celsys(wcs, wcspos);

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


/******* red_to_raw **********************************************************
PROTO	int red_to_raw(wcsstruct *, double *, double *)
PURPOSE	Convert reduced (World Coordinate System) coords to raw (pixel)
	coords.
INPUT	WCS structure,
	Pointer to the array of input (reduced) coordinates,
	Pointer to the array of output (pixel) coordinates.
OUTPUT	RETURN_OK if mapping successful, RETURN_ERROR otherwise.
NOTES	-.
AUTHOR	E. Bertin (IAP)
VERSION	23/10/2003
 ***/
int	red_to_raw(wcsstruct *wcs, double *redpos, double *pixpos)

  {
   struct wcsprm	*wcsprm;
   double		offset;

  wcsprm = wcs->wcsprm;
/* Initialize if required */
  if (wcsprm && wcsprm->flag != WCSSET)
    {
    if (wcsset(wcs->naxis, (const char(*)[9])wcs->ctype, wcsprm))
      return RETURN_ERROR;
    }

  if (wcsprm && wcsprm->flag != 999 && wcsprm->cubeface != -1)
    {
/*-- Separation between faces */
    offset = (wcs->prj->r0 == 0.0 ? 90.0 : wcs->prj->r0*PI/2.0);
/*-- Stack faces in a cube */
    if (redpos[wcs->lat] < -0.5*offset)
      {
      redpos[wcs->lat] += offset;
      redpos[wcsprm->cubeface] = 5.0;
      }
    else if (redpos[wcs->lat] > 0.5*offset)
      {
      redpos[wcs->lat] -= offset;
      redpos[wcsprm->cubeface] = 0.0;
      }
    else if (redpos[wcs->lng] > 2.5*offset)
      {
      redpos[wcs->lng] -= 3.0*offset;
      redpos[wcsprm->cubeface] = 4.0;
      }
    else if (redpos[wcs->lng] > 1.5*offset)
      {
      redpos[wcs->lng] -= 2.0*offset;
      redpos[wcsprm->cubeface] = 3.0;
      }
    else if (redpos[wcs->lng] > 0.5*offset)
      {
      redpos[wcs->lng] -= offset;
      redpos[wcsprm->cubeface] = 2.0;
      }
    else
      redpos[wcsprm->cubeface] = 1.0;
    }

/* Apply forward linear transformation */
  if (linfwd(redpos, wcs->lin, pixpos))
      return RETURN_ERROR;

  return RETURN_OK;
  }


/******* raw_to_red **********************************************************
PROTO	int raw_to_red(wcsstruct *, double *, double *)
PURPOSE	Convert raw (pixel) coordinates to reduced WCS coordinates.
INPUT	WCS structure,
	Pointer to the array of input (pixel) coordinates,
	Pointer to the array of output (reduced) coordinates.
OUTPUT	RETURN_OK if mapping successful, RETURN_ERROR otherwise.
NOTES	-.
AUTHOR	E. Bertin (IAP)
VERSION	23/10/2003
 ***/
int	raw_to_red(wcsstruct *wcs, double *pixpos, double *redpos)

  {
   struct wcsprm	*wcsprm;
   double		offset;
   int			face;

  wcsprm = wcs->wcsprm;
/* Initialize if required */
  if (wcsprm && wcsprm->flag != WCSSET)
    {
    if (wcsset(wcs->naxis, (const char(*)[9])wcs->ctype, wcsprm))
      return RETURN_ERROR;
    }
   
/* Apply reverse linear transformation */
   if (linrev(pixpos, wcs->lin, redpos))
     return RETURN_ERROR;

  if (wcsprm && wcsprm->flag != 999 && wcsprm->cubeface != -1)
    {
/*-- Do we have a CUBEFACE axis? */
    face = (int)(redpos[wcsprm->cubeface] + 0.5);
    if (fabs(redpos[wcsprm->cubeface]-face) > 1e-10)
      return RETURN_ERROR;

/*-- Separation between faces. */
    offset = (wcs->prj->r0 == 0.0 ? 90.0 : wcs->prj->r0*PI/2.0);
/*-- Lay out faces in a plane. */
    switch (face)
      {
      case 0:
        redpos[wcs->lat] += offset;
        break;
      case 1:
        break;
      case 2:
        redpos[wcs->lng] += offset;
        break;
      case 3:
        redpos[wcs->lng] += offset*2;
        break;
      case 4:
        redpos[wcs->lng] += offset*3;
        break;
      case 5:
        redpos[wcs->lat] -= offset;
        break;
      default:
        return RETURN_ERROR;
      }
    }

  return RETURN_OK;
  }


/******* wcs_dist ***********************************************************
PROTO	double wcs_dist(wcsstruct *wcs, double *wcspos1, double *wcspos2)
PURPOSE	Compute the angular distance between 2 points on the sky.
INPUT	WCS structure,
	Pointer to the first array of world coordinates,
	Pointer to the second array of world coordinates.
OUTPUT	Angular distance (in degrees) between points.
NOTES	-.
AUTHOR	E. Bertin (IAP)
VERSION	24/07/2002
 ***/
double	wcs_dist(wcsstruct *wcs, double *wcspos1, double *wcspos2)

  {
  double	d, dp;
  int		i, lng, lat;

  lng = wcs->lng;
  lat = wcs->lat;
  if (lat!=lng)
    {
/*-- We are operating in angular coordinates */
    d = sin(wcspos1[lat]*DEG)*sin(wcspos2[lat]*DEG)
	+ cos(wcspos1[lat]*DEG)*cos(wcspos2[lat]*DEG)
		*cos((wcspos1[lng]-wcspos2[lng])*DEG);
    return d>-1.0? (d<1.0 ? acos(d)/DEG : 0.0) : 180.0;
    }
  else
    {
    d = 0.0;
    for (i=0; i<wcs->naxis; i++)
      {
      dp = wcspos1[i] - wcspos2[i];
      d += dp*dp;
      }
    return sqrt(d);
    }
  }


/******* 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	03/01/2008
 ***/
double	wcs_scale(wcsstruct *wcs, double *pixpos)

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

  if (raw_to_wcs(wcs, pixpos, wcspos))
    return 0.0;

  lng = wcs->lng;
  lat = wcs->lat;
  if (lng == lat)
    {
    lng = 0;
    lat = 1;
    }

/* Compute pixel scale */
  pixpos2[lng] = pixpos[lng] + 1.0;
  pixpos2[lat] = pixpos[lat];
  if (raw_to_wcs(wcs, pixpos2, wcspos1))
    return 0.0;
  pixpos2[lng] -= 1.0;
  pixpos2[lat] += 1.0;
  if (raw_to_wcs(wcs, pixpos2, wcspos2))
    return 0.0;
  dpos1 = wcspos1[lng]-wcspos[lng];
  dpos2 = wcspos2[lng]-wcspos[lng];
  if (wcs->lng!=wcs->lat)
    {
    if (dpos1>180.0)
      dpos1 -= 360.0;
    else if (dpos1<-180.0)
      dpos1 += 360.0;
    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));
    }
  else
    return fabs((dpos1*(wcspos2[lat]-wcspos[lat])
		-(wcspos1[lat]-wcspos[lat])*dpos2));
  }


/****** wcs jacobian *********************************************************
PROTO	double wcs_jacobian(wcsstruct *wcs, double *pixpos, double *jacob)
PURPOSE	Compute the local Jacobian matrix of the astrometric deprojection.
INPUT	WCS structure,
	Pointer to the array of local raw coordinates,
	Pointer to the jacobian array (output).
OUTPUT	Determinant over spatial coordinates (=pixel area), or -1.0 if mapping
	was unsuccesful.
NOTES   Memory must have been allocated (naxis*naxis*sizeof(double)) for the
        Jacobian array.
AUTHOR	E. Bertin (IAP)
VERSION	11/10/2007
 ***/
double	wcs_jacobian(wcsstruct *wcs, double *pixpos, double *jacob)
  {
   double	pixpos0[NAXIS], wcspos0[NAXIS], wcspos[NAXIS],
		dpos;
   int		i,j, lng,lat,naxis;

  lng = wcs->lng;
  lat = wcs->lat;
  naxis = wcs->naxis;
  for (i=0; i<naxis; i++)
    pixpos0[i] = pixpos[i];
  if (raw_to_wcs(wcs, pixpos0, wcspos0) == RETURN_ERROR)
    return -1.0;
  for (i=0; i<naxis; i++)
    {
    pixpos0[i] += 1.0;
    if (raw_to_wcs(wcs, pixpos0, wcspos) == RETURN_ERROR)
      return -1.0;
    pixpos0[i] -= 1.0;
    for (j=0; j<naxis; j++)
      {
      dpos = wcspos[j]-wcspos0[j];
      if (lng!=lat && j==lng)
        {
        if (dpos>180.0)
          dpos -= 360.0;
        else if (dpos<-180.0)
          dpos += 360.0;
        dpos *= cos(wcspos0[lat]*DEG);
        }
      jacob[j*naxis+i] = dpos;
      }
    }

  if (lng==lat)
    {
    lng = 0;
    lat = 1;
    }

  return fabs(jacob[lng+naxis*lng]*jacob[lat+naxis*lat]
		- jacob[lat+naxis*lng]*jacob[lng+naxis*lat]);
  }


/****** wcs rawtoraw *********************************************************
PROTO	double wcs_rawtoraw(wcsstruct *wcsin, wcsstruct *wcsout,
		double *pixpos, double *jacob)
PURPOSE	Convert raw coordinates from input wcs structure to raw coordinates
	from output wcs structure, and the local Jacobian matrix of the
	reprojection.
INPUT	WCS input structure,
	WCS output structure,
	pointer to the array of local input raw coordinates,
	pointer to the array of local output raw coordinates (output),
	pointer to the jacobian array (output).
OUTPUT	Determinant over spatial coordinates (ratio of pixel areas),
	0.0 if jacob is NULL (Jacobian not computed in that case),
	or -1.0 if mapping was unsuccesful.
NOTES   Memory must have been allocated (naxis*naxis*sizeof(double)) for the
        Jacobian array.
AUTHOR	E. Bertin (IAP)
VERSION	12/06/2012
 ***/
double	wcs_rawtoraw(wcsstruct *wcsin, wcsstruct *wcsout,
		double *pixposin, double *pixposout, double *jacob)
  {
   double	pixpos0[NAXIS], pixpos2[NAXIS], wcspos[NAXIS];
   int		i,j, lng,lat,naxis;

  naxis = wcsin->naxis;
  for (i=0; i<naxis; i++)
    pixpos0[i] = pixposin[i];

/* Coordinate transformation */
  if (raw_to_wcs(wcsin, pixposin, wcspos) == RETURN_ERROR)
    return -1.0;
  if (wcs_to_raw(wcsout, wcspos, pixposout) == RETURN_ERROR)
    return -1.0;
 
/* Jacobian */
  if (!jacob)
    return 0.0;

  lng = wcsin->lng;
  lat = wcsin->lat;
  for (i=0; i<naxis; i++)
    {
    pixpos0[i] += 1.0;
    if (raw_to_wcs(wcsin, pixpos0, wcspos) == RETURN_ERROR)
      return -1.0;
    if (wcs_to_raw(wcsout, wcspos, pixpos2) == RETURN_ERROR)
      return -1.0;
    pixpos0[i] -= 1.0;
    for (j=0; j<naxis; j++)
      jacob[j*naxis+i] = pixpos2[j] - pixposout[j];
    }

  if (lng==lat)
    {
    lng = 0;
    lat = 1;
    }

  return fabs(jacob[lng+naxis*lng]*jacob[lat+naxis*lat]
		- jacob[lat+naxis*lng]*jacob[lng+naxis*lat]);
  }


/******* wcs_chirality *******************************************************
PROTO	int wcs_chirality(wcsstruct *wcs)
PURPOSE	Compute the chirality of a WCS projection.
INPUT	WCS structure.
OUTPUT	+1 if determinant of matrix is positive, -1 if negative, 0 if null.
NOTES	-.
AUTHOR	E. Bertin (IAP)
VERSION	26/09/2006
 ***/
int	wcs_chirality(wcsstruct *wcs)

  {
   double	a;
   int		lng,lat, naxis;

  lng = wcs->lng;
  lat = wcs->lat;
  naxis = wcs->naxis;
  if (lng==lat && naxis>=2)
    {
    lng = 0;
    lat = 1;
    }

  a = wcs->cd[lng*naxis+lng]*wcs->cd[lat*naxis+lat]
	- wcs->cd[lng*naxis+lat]*wcs->cd[lat*naxis+lng];
  return a>TINY? 1 : (a<-TINY? -1 : 0);
  }


/****** precess_wcs **********************************************************
PROTO	void precess_wcs(wcsstruct *wcs, double yearin, double yearout)
PURPOSE	Precess the content of a WCS structure according to the equinox.
INPUT	WCS structure,
	Input year,
	Output year.
OUTPUT	-.
NOTES	Epoch for coordinates should be J2000 (FK5 system).
AUTHOR	E. Bertin (IAP)
VERSION	04/01/2008
 ***/
void	precess_wcs(wcsstruct *wcs, double yearin, double yearout)

  {
   double	crval[NAXIS],a[NAXIS*NAXIS],b[NAXIS*NAXIS],
		*c,*at,
		val, cas, sas, angle, dalpha;
   int		i,j,k, lng,lat, naxis;

  lng = wcs->lng;
  lat = wcs->lat;
  if (lat==lng || yearin==yearout)
    return;
  naxis = wcs->naxis;
/* Precess to year out */
  precess(yearin, wcs->crval[lng], wcs->crval[lat], yearout,
	&crval[lng], &crval[lat]);

  dalpha = (crval[lng] - wcs->crval[lng])*DEG;

/* Compute difference angle with the north axis between start and end */
  angle = (dalpha!=0.0 && (crval[lat] - wcs->crval[lat])*DEG != 0.0) ?
	180.0 - (atan2(sin(dalpha),
	cos(crval[lat]*DEG)*tan(wcs->crval[lat]*DEG)
	- sin(crval[lat]*DEG)*cos(dalpha))
	+ atan2(sin(dalpha),
	cos(wcs->crval[lat]*DEG)*tan(crval[lat]*DEG)
	- sin(wcs->crval[lat]*DEG)*cos(dalpha)))/DEG
	: 0.0;

/* A = C*B */
  c = wcs->cd;
/* The B matrix is made of 2 numbers */

  cas = cos(angle*DEG);
  sas = sin(-angle*DEG);
  for (i=0; i<naxis; i++)
    b[i+i*naxis] = 1.0;
  b[lng+lng*naxis] = cas;
  b[lat+lng*naxis] = -sas;
  b[lng+lat*naxis] = sas;
  b[lat+lat*naxis] = cas;
  at = a;
  for (j=0; j<naxis; j++)
    for (i=0; i<naxis; i++)
      {
      val = 0.0;
      for (k=0; k<naxis; k++)
        val += c[k+j*naxis]*b[i+k*naxis];
      *(at++) = val;
      }

  at = a;

  for (i=0; i<naxis*naxis; i++)
    *(c++) = *(at++);

  wcs->crval[lng] = crval[lng];
  wcs->crval[lat] = crval[lat];
  wcs->equinox = yearout;

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

  return;
  }


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


/********************************* b2j ***********************************/
/*
conver equatorial coordinates from equinox and epoch B1950 to equinox and
epoch J2000 for extragalactic sources (from Aoki et al. 1983).
*/
void	b2j(double yearobs, double alphain, double deltain,
		double *alphaout, double *deltaout)
  {
   int		i,j;
   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.9999256782,     -0.0111820611,     -0.0048579477,
    0.00000242395018, -0.00000002710663, -0.00000001177656},
  { 0.0111820610,      0.9999374784,     -0.0000271765,
    0.00000002710663,  0.00000242397878, -0.00000000006587},
  { 0.0048579479,     -0.0000271474,      0.9999881997,
    0.00000001177656, -0.00000000006582,  0.00000242410173},
  {-0.000551,        -0.238565,           0.435739,
    0.99994704,	     -0.01118251,        -0.00485767},
  { 0.238514,        -0.002662,          -0.008541,
    0.01118251,	      0.99995883,        -0.00002718},
  {-0.435623,         0.012254,           0.002117,
    0.00485767,      -0.00002714,         1.00000956}},
 		a1[3], r[3], ro[3], r1[3], r2[3], v1[3], v[3];
   double		cai, sai, cdi, sdi, dotp, rmod, alpha, delta,
			t1 = (yearobs - 1950.0)/100.0;

  alphain *= PI/180.0;
  deltain *= PI/180.0;
  cai = cos(alphain);
  sai = sin(alphain);
  cdi = cos(deltain);
  sdi = sin(deltain);
  ro[0] = cdi*cai;
  ro[1] = cdi*sai;
  ro[2] = sdi;
  dotp = 0.0;
  for (i=0; i<3; i++)
    {
    a1[i] = a[i]+ap[i]*ARCSEC*t1;
    dotp += a1[i]*ro[i];
    }
  for (i=0; i<3; i++)
    {
    r1[i] = ro[i] - a1[i] + dotp*ro[i];
    r[i] = v[i] = v1[i] = 0.0;
    }
  for (j=0; j<6; j++)
    for (i=0; i<6; i++)
      {
      if (j<3)
        r[j] += m[j][i]*(i<3?r1[i]:v1[i-3]);
      else
         v[j-3] += m[j][i]*(i<3?r1[i]:v1[i-3]);
      }
  rmod = 0.0;
  for (i=0; i<3; i++)
    {
    r2[i] = r[i]+v[i]*ARCSEC*(t1-0.5);
    rmod += r2[i]*r2[i];
    }
  rmod = sqrt(rmod);
  delta = asin(r2[2]/rmod);
  alpha = acos(r2[0]/cos(delta)/rmod);
  if (r2[1]<0)
    alpha = 2*PI - alpha;
  *alphaout = alpha*180.0/PI;
  *deltaout = delta*180.0/PI;

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