astrom.c

/*
*				astrom.c
*
* High level astrometric computations.
*
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
*
*	This file part of:	SExtractor
*
*	Copyright:		(C) 1993-2010 Emmanuel Bertin -- IAP/CNRS/UPMC
*
*	License:		GNU General Public License
*
*	SExtractor 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.
*	SExtractor 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 SExtractor. If not, see <http://www.gnu.org/licenses/>.
*
*	Last modified:		11/10/2010
*
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/

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

#include	<math.h>
#include	<stdlib.h>
#include	<string.h>

#include 	"wcs/wcs.h"
#include	"define.h"
#include	"globals.h"
#include	"prefs.h"
#include	"astrom.h"
#include	"fitswcs.h"
#include	"wcs/tnx.h"

static obj2struct	*obj2 = &outobj2;

/****************************** initastrom **********************************/
/*
Initialize astrometrical structures.
*/
void	initastrom(picstruct *field)

  {
   wcsstruct	*wcs;

  wcs = field->wcs;

/* Test if the WCS is in use */
  if (wcs->lng != wcs->lat)
    {
    if (FLAG(obj2.dtheta2000) || FLAG(obj2.dtheta1950))
      {
      if (fabs(wcs->equinox-2000.0)>0.003)
        precess(wcs->equinox, 0.0, 90.0, 2000.0, &wcs->ap2000, &wcs->dp2000);
      else
        {
        wcs->ap2000 = 0.0;
        wcs->dp2000 = 90.0;
        }

      if (FLAG(obj2.theta1950) || FLAG(obj2.poserr_theta1950))
        j2b(wcs->equinox, wcs->ap2000, wcs->dp2000, &wcs->ap1950, &wcs->dp1950);
      }
    }

/* Override astrometric definitions only if user supplies a pixel-scale */
  if (prefs.pixel_scale == 0.0)
    field->pixscale = wcs->pixscale*3600.0;	/* in arcsec */
  else
    field->pixscale = prefs.pixel_scale;

  return;
  }


/***************************** astrom_pos **********************************/
/*
Compute real FOCAL and WORLD coordinates according to FITS info.
*/
void	astrom_pos(picstruct *field, objstruct *obj)

  {
   wcsstruct	*wcs;
   double	rawpos[NAXIS], wcspos[NAXIS],
		da,dd;
   int		lng,lat;

  wcs = field->wcs;
  lng = wcs->lng;
  lat = wcs->lat;

/* If working with WCS, compute FOCAL coordinates and local matrix */
  if (FLAG(obj2.mxf))
    {
    rawpos[0] = obj2->posx;
    rawpos[1] = obj2->posy;
    raw_to_red(wcs, rawpos, wcspos);
    obj2->mxf = wcspos[0];
    obj2->myf = wcspos[1];
    }

/* If working with WCS, compute WORLD coordinates and local matrix */
  if (FLAG(obj2.mxw))
    {
    rawpos[0] = obj2->posx;
    rawpos[1] = obj2->posy;
    raw_to_wcs(wcs, rawpos, wcspos);
    obj2->mxw = wcspos[0];
    obj2->myw = wcspos[1];
    if (lng != lat)
      {
      obj2->alphas = lng<lat? obj2->mxw : obj2->myw;
      obj2->deltas = lng<lat? obj2->myw : obj2->mxw;
      if (FLAG(obj2.alpha2000))
        {
        if (fabs(wcs->equinox-2000.0)>0.003)
          precess(wcs->equinox, wcspos[lng<lat?0:1], wcspos[lng<lat?1:0],
		2000.0, &obj2->alpha2000, &obj2->delta2000);
        else
          {
          obj2->alpha2000 = lng<lat? obj2->mxw : obj2->myw;
          obj2->delta2000 = lng<lat? obj2->myw : obj2->mxw;
          }
        if (FLAG(obj2.dtheta2000))
          {
          da = wcs->ap2000 - obj2->alpha2000;
          dd = (sin(wcs->dp2000*DEG)
		-sin(obj2->delta2000*DEG)*sin(obj2->deltas*DEG))
		/(cos(obj2->delta2000*DEG)*cos(obj2->deltas*DEG));
          dd = dd<1.0? (dd>-1.0?acos(dd)/DEG:180.0) : 0.0;
          obj2->dtheta2000 = (((da>0.0 && da<180.0) || da<-180.0)?-dd:dd);
          }
        if (FLAG(obj2.alpha1950))
          {
          j2b(wcs->equinox, obj2->alpha2000, obj2->delta2000,
		&obj2->alpha1950, &obj2->delta1950);
          if (FLAG(obj2.dtheta1950))
            {
            da = wcs->ap1950 - obj2->alpha1950;
            dd = (sin(wcs->dp1950*DEG)
		-sin(obj2->delta1950*DEG)*sin(obj2->deltas*DEG))
		/(cos(obj2->delta1950*DEG)*cos(obj2->deltas*DEG));
            dd = dd<1.0? (dd>-1.0?acos(dd)/DEG:180.0) : 0.0;
            obj2->dtheta1950 = (((da>0.0 && da<180.0) || da<-180.0)?-dd:dd);
           }
          }
        }
      }
    }

/* Custom coordinate system for the MAMA machine */
  if (FLAG(obj2.mamaposx))
    {
    rawpos[0] = obj2->posx - 0.5;
    rawpos[1] = obj2->posy - 0.5;
    raw_to_wcs(wcs, rawpos, wcspos);
    obj2->mamaposx = wcspos[1]*(MAMA_CORFLEX+1.0);
    obj2->mamaposy = wcspos[0]*(MAMA_CORFLEX+1.0);
    }

  return;
  }


/***************************** astrom_peakpos *******************************/
/*
Compute real FOCAL and WORLD peak coordinates according to FITS info.
*/
void	astrom_peakpos(picstruct *field, objstruct *obj)

  {
   wcsstruct	*wcs;
   double	rawpos[NAXIS], wcspos[NAXIS];
   int		lng,lat;

  wcs = field->wcs;
  lng = wcs->lng;
  lat = wcs->lat;

  if (FLAG(obj2.peakxf))
    {
    rawpos[0] = obj->peakx;
    rawpos[1] = obj->peaky;
    raw_to_red(wcs, rawpos, wcspos);
    obj2->peakxf = wcspos[0];
    obj2->peakyf = wcspos[1];
    }
  if (FLAG(obj2.peakxw))
    {
    rawpos[0] = obj->peakx;