/*
  This file is part of the Astrometry.net suite.
  Copyright 2006, 2007 Keir Mierle, David W. Hogg, Sam Roweis and Dustin Lang.

  The Astrometry.net suite 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, version 2.

  The Astrometry.net suite 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 the Astrometry.net suite ; if not, write to the Free Software
  Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301 USA
*/

#ifndef ANWCS_H
#define ANWCS_H

#include <stdio.h>
#include "an-bool.h"
#include "keywords.h"

#define SIP_MAXORDER 10

// WCS TAN header.
typedef struct tan_s {

	// World coordinate of the tangent point, in ra,dec.
	double crval[2];

	// Tangent point location in pixel (CCD) coordinates
	// This may not be in the image; consider a telescope with an array of
	// CCD's, where the tangent point is in one or none of the CCD's.
	double crpix[2];

	// Matrix for the linear transformation of relative pixel coordinates
	// (u,v) onto "intermediate world coordinates", which are in degrees
	// (x,y). The x,y coordinates are on the tangent plane. If the SIP
	// terms are all zero, then the equation to get from pixel coordinates
	// to intermediate world coordinates is:
	//
	//   u = pixel_x - crpix0
	//   v = pixel_y - crpix1
	// 
	//   x  = [cd00 cd01] * u
	//   y    [cd10 cd11]   v
	// 
	// where x,y are in intermediate world coordinates (i.e. x points
	// along negative ra and y points to positive dec) and u,v are in pixel
	// coordinates.
	double cd[2][2];

    // size of the image in pixels.  Not strictly part of the WCS, but useful!
    double imagew;
    double imageh;

} tan_t;

// Flat structure for minimal SIP wcs. This structure should contain enough
// information to effectively represent any image, provided it is possible to
// convert that image's projecton to a TAN projection and the distortion to SIP
// distortion.
typedef struct sip_s {

	// A basic TAN header.
	tan_t wcstan;

	// Forward SIP coefficients
	// The transformation from relative pixel coordinates to intermediate
	// world coordinates[1] is:
	// 
	//   x  = [cd00 cd01] * (u + f(u,v))          x,y are intermediate coordinates on the sky (in deg)
	//   y    [cd10 cd11]   (v + g(u,v))            which are just a "translation" away from final WCS
	//                                            u,v are original (unwarped) pixel coordinates
	// where
	//                                       p    q
	//   U = u + f(u,v) = u + SUM a[p][q] * u  * v  ,  p+q <= a_order
  	//                        p,q
   //
	//                                       p    q
	//   V = v + g(u,v) = v + SUM b[p][q] * u  * v  ,  p+q <= b_order
	//                        p,q
	// 
	// [1] The SIP convention for representing distortion in FITS image
	// headers. D. L. Shupe, M.Moshir, J. Li, D. Makovoz, R. Narron, R. N.
	// Hook. Astronomical Data Analysis Software and Systems XIV.
	// http://ssc.spitzer.caltech.edu/postbcd/doc/shupeADASS.pdf
	//
	// Note: These matricies are larger than they strictly need to be
	// because aij = 0 if i+j > a_order and similarily for b.
	// 
	// Note: The convention for indicating that no SIP polynomial is
	// present is to simply set [ab]_order to zero.
	int a_order, b_order;
	double a[SIP_MAXORDER][SIP_MAXORDER];
	double b[SIP_MAXORDER][SIP_MAXORDER];

	// Inverse SIP coefficients
	// To convert from world coordinates back into image coordinates, the
	// inverse transformation may be stored. To convert from intermediate
	// world coordinates, first we calculate the linear pixel coordinates:
	// 
	//                   -1
	//    U = [cd00 cd01]    * x
	//    V   [cd10 cd11]      y
	//
	// Then, the original pixel coordinates are computed as:
	// 
	//                            p    q
	//   u  = U + SUM ap[p][q] * U  * V  ,  p+q <= ap_order
	//            p,q
        //
	//                            p    q
	//   v  = V + SUM bp[p][q] * U  * V  ,  p+q <= ap_order
	//            p,q
	// 
	// Note: ap_order does not necessarily equal a_order, because the
	// inverse of a nth-order polynomial may be of higer order than n.
	// 
	// Note: The convention for indicating that no inverse SIP polynomial
	// is present is to simply set [ab]p_order to zero.
	int ap_order, bp_order;
	double ap[SIP_MAXORDER][SIP_MAXORDER];
	double bp[SIP_MAXORDER][SIP_MAXORDER];
} sip_t;

sip_t* sip_create(void);
void   sip_free(sip_t* sip);

void sip_copy(sip_t* dest, const sip_t* src);

// Set the given SIP wcs to the given TAN wcs.
void sip_wrap_tan(const tan_t* tan, sip_t* sip);

double sip_imagew(sip_t* sip);
double sip_imageh(sip_t* sip);

// Pixels to RA,Dec in degrees.
void   sip_pixelxy2radec(const sip_t* sip, double px, double py, double *a, double *d);

// Pixels to XYZ unit vector.
void   sip_pixelxy2xyzarr(const sip_t* sip, double px, double py, double *xyz);

// RA,Dec in degrees to Pixels.
// Returns FALSE if the point is on the opposite side of the sphere (and hence the point
// does not project onto the tangent plane)
WarnUnusedResult
bool sip_radec2pixelxy(const sip_t* sip, double a, double d, double *px, double *py);

WarnUnusedResult
bool sip_radec2pixelxy_check(const sip_t* sip, double ra, double dec, double *px, double *py);

WarnUnusedResult
bool sip_xyzarr2pixelxy(const sip_t* sip, const double* xyz, double *px, double *py);

WarnUnusedResult
bool sip_xyz2pixelxy(const sip_t* sip, double x, double y, double z, double *px, double *py);

// Pixels to Intermediate World Coordinates in degrees.
void sip_pixelxy2iwc(const sip_t* sip, double px, double py,
					 double *iwcx, double* iwcy);


double tan_det_cd(const tan_t* tan);
double sip_det_cd(const sip_t* sip);
// returns pixel scale in arcseconds/pixel (NOT arcsec^2)
double tan_pixel_scale(const tan_t* tn);
double sip_pixel_scale(const sip_t* sip);

// these take *relative* pixel coords (WRT crpix)
void   sip_calc_inv_distortion(const sip_t* sip, double U, double V, double* u, double *v);
void   sip_calc_distortion(const sip_t* sip, double u, double v, double* U, double *V);
      
// applies forward SIP distortion to pixel coords.
void sip_pixel_distortion(const sip_t* sip, double x, double y, double* X, double *Y);

// Pixels to XYZ unit vector.
void   tan_pixelxy2xyzarr(const tan_t* tan, double px, double py, double *xyz);

// Pixels to RA,Dec in degrees.
void   tan_pixelxy2radec(const tan_t* wcs_tan, double px, double py, double *ra, double *dec);

// RA,Dec in degrees to Pixels.
// Returns FALSE if the point is on the opposite side of the sphere.
WarnUnusedResult
bool   tan_radec2pixelxy(const tan_t* wcs_tan, double ra, double dec, double *px, double *py);

// xyz unit vector to Pixels.
WarnUnusedResult
bool   tan_xyzarr2pixelxy(const tan_t* wcs_tan, const double* xyz, double *px, double *py);

void tan_iwc2pixelxy(const tan_t* tan, double iwcx, double iwcy,
					 double *px, double* py);
void tan_iwc2xyzarr(const tan_t* tan, double x, double y, double *xyz);

/**
 Subtracts off CRPIX, multiplies by CD matrix.
 Results are in degrees.
 */
void tan_pixelxy2iwc(const tan_t* tan, double px, double py, double *iwcx, double* iwcy);

void   sip_print(const sip_t*);
void   sip_print_to(const sip_t*, FILE* fid);

void tan_print(const tan_t* tan);
void tan_print_to(const tan_t* tan, FILE* f);

// for python
void sip_get_crval(const sip_t* sip, double* ra, double* dec);

double sip_get_orientation(const sip_t* sip);

#endif
