# include <stdio.h>
# include <stdlib.h>		/* calloc */
# include <math.h>		/* fabs, sqrt */

# include <c_iraf.h>
# include <xtables.h>
# include "../stis.h"
# include "calstis7.h"
# include "../stiserr.h"

typedef struct {
	IRAFPointer tp;			/* pointer to table descriptor */
	IRAFPointer cp_opt_elem;	/* column descriptors */

	IRAFPointer cp_xsize;
	IRAFPointer cp_ysize;
	IRAFPointer cp_xref;
	IRAFPointer cp_yref;

	IRAFPointer cp_pedigree;
	IRAFPointer cp_descrip;
	int nrows;			/* number of rows in table */
} TblInfo;

typedef struct {
	char opt_elem[STIS_CBUF+1];	/* optical element name */
} TblRow;

static int OpenDistTab (char *, TblInfo *);
static int ReadDistTab (TblInfo *, int, TblRow *);
static int ReadCoefficients (TblInfo *, int, CoordInfo **, DistInfo *);
static int CloseDistTab (TblInfo *);

/* This routine reads the distortion information from the imaging
   distortion table IDCTAB.  This is only used for obstype=IMAGING.

   The distortion information table should contain the following:
	header parameters:
		NORDER:  polynomial order (int)
	columns:
		OPT_ELEM:  grating (or mirror) name (string)
		XSIZE, YSIZE:  size for output image (int)
				(also saved in dist->npix[0], npix[1])
		XREF, YREF:  zero points in X and Y (double)
		C00, C10, etc:  arrays of coefficients (double)

   The table is read to find the first row for which the value of
   OPT_ELEM is the same as in the input image header.

   The output image size is read from XSIZE & YSIZE, and this is the
   only (useful) information saved in the CoordInfo structure.  It's
   copied there because Do2Dx uses it when creating the output image.

   Note:
	Memory is allocated for the coords list; it should be freed
	by calling FreeCoord.
	coords should have been initialized to NULL.
	Memory is allocated for the distortion coefficients, dist->xcoeff
	and dist->ycoeff.  The memory should be freed by calling freeDist.

   Phil Hodge, 2000 Jan 13:
	Add one to aperture and opt_elem buffer sizes.

   Phil Hodge, 2000 Oct 9:
	Remove references to aperture; select row only on opt_elem.
	Get columns C00, C10, ... D00, D10, ... .
*/

int GetIDC (StisInfo7 *sts, CoordInfo **coords, DistInfo *dist) {

/* arguments:
StisInfo7 *sts      i: calibration switches and info
CoordInfo **coords  o: coordinate parameters for output image
DistInfo *dist      o: distortion coefficients
*/

	int status;

	TblInfo tabinfo;	/* pointer to table descriptor, etc */
	TblRow tabrow;		/* values read from a table row */

	int row;		/* loop index */
	int foundit = 0;	/* true if parameters found in table */
	int RowPedigree (RefTab *, int, IRAFPointer, IRAFPointer, IRAFPointer);
	int SameString (char *, char *);

	/* Open the distortion information table. */
	if (status = OpenDistTab (sts->distntab.name, &tabinfo))
	    return (status);

	/* Check for a match with opt_elem. */

	for (row = 1;  row <= tabinfo.nrows;  row++) {

	    if (status = ReadDistTab (&tabinfo, row, &tabrow))
		return (status);

	    if (SameString (tabrow.opt_elem, sts->opt_elem)) {

		/* This is the row; read the info. */
		foundit = 1;

		/* Get pedigree & descrip from the row. */
		if (status = RowPedigree (&sts->distntab, row,
			tabinfo.tp, tabinfo.cp_pedigree, tabinfo.cp_descrip))
		    return (status);
		if (sts->distntab.goodPedigree == DUMMY_PEDIGREE) {
		    printf ("Warning  DUMMY pedigree in row %d of %s.\n",
			row, sts->distntab.name);
		    sts->x2dcorr_o = DUMMY;
		    CloseDistTab (&tabinfo);
		    return (0);
		}

		if (status = ReadCoefficients (&tabinfo, row, coords, dist))
		    return (status);

		break;
	    }
	}

	if (!foundit) {
	    printf ("Warning  Matching row not found in IDCTAB %s; \\\n",
				sts->distntab.name);

	    printf ("Warning  OPT_ELEM %s.\n", sts->opt_elem);
	    sts->x2dcorr_o = OMIT;
	}

	if (status = CloseDistTab (&tabinfo))
	    return (status);

	return (0);
}

/* This routine opens the 2-D distortion table, finds the columns
   that we need, and gets the total number of rows in the table.
*/

static int OpenDistTab (char *tname, TblInfo *tabinfo) {

	tabinfo->tp = c_tbtopn (tname, IRAF_READ_ONLY, 0);
	if (c_iraferr()) {
	    printf ("ERROR    IDCTAB `%s' not found.\n", tname);
	    return (OPEN_FAILED);
	}

	tabinfo->nrows = c_tbpsta (tabinfo->tp, TBL_NROWS);

	/* Find the columns. */

	c_tbcfnd1 (tabinfo->tp, "OPT_ELEM", &tabinfo->cp_opt_elem);

	c_tbcfnd1 (tabinfo->tp, "XSIZE", &tabinfo->cp_xsize);
	c_tbcfnd1 (tabinfo->tp, "YSIZE", &tabinfo->cp_ysize);
	c_tbcfnd1 (tabinfo->tp, "XREF", &tabinfo->cp_xref);
	c_tbcfnd1 (tabinfo->tp, "YREF", &tabinfo->cp_yref);

	if (tabinfo->cp_opt_elem == 0 ||
	    tabinfo->cp_xsize == 0  || tabinfo->cp_ysize == 0 ||
	    tabinfo->cp_xref == 0  || tabinfo->cp_yref == 0) {

	    c_tbtclo (tabinfo->tp);
	    printf ("ERROR    Column not found in IDCTAB.\n");
	    return (COLUMN_NOT_FOUND);
	}

	/* Pedigree and descrip are optional columns. */
	c_tbcfnd1 (tabinfo->tp, "PEDIGREE", &tabinfo->cp_pedigree);
	c_tbcfnd1 (tabinfo->tp, "DESCRIP", &tabinfo->cp_descrip);

	return (0);
}

/* This routine reads the column (OPT_ELEM) used to select
   the correct row.
*/

static int ReadDistTab (TblInfo *tabinfo, int row, TblRow *tabrow) {

	c_tbegtt (tabinfo->tp, tabinfo->cp_opt_elem, row,
			tabrow->opt_elem, STIS_CBUF);
	if (c_iraferr())
	    return (TABLE_ERROR);

	return (0);
}

/* This routine reads the data from one row into the coords and dist
   structures.  The size and coordinate parameters are gotten into coords.
   The number of elements in the arrays of distortion coefficients is
   gotten into dist, memory is allocated, and the array values are read.
   The order of the polynomial coefficients is gotten from the NORDER
   header keyword.
*/

static int ReadCoefficients (TblInfo *tabinfo, int row,
		CoordInfo **coords, DistInfo *dist) {

	char colname[SZ_COLNAME];	/* a column name */
	IRAFPointer cp;		/* for getting coefficients */
	int i, j, k;
	int kx, ky;		/* indexes for coefficients of x and y */
	int npix[2];		/* the size of the output image */
	double temp;		/* for computing the image scale */
	int status;

	CoordInfo *newrec;
	int NewCoord (CoordInfo **, CoordInfo *);

	if ((newrec = malloc (sizeof (CoordInfo))) == NULL) {
	    printf ("ERROR    Can't allocate memory.\n");
	    return (OUT_OF_MEMORY);
	}
	newrec->sporder = 1;	/* use nominal values for imaging mode */
	newrec->a2center = 1.;
	newrec->next = NULL;

	/* Get size info and coordinate parameters. */
	c_tbegti (tabinfo->tp, tabinfo->cp_xsize, row, &npix[0]);
	c_tbegti (tabinfo->tp, tabinfo->cp_ysize, row, &npix[1]);
	if (c_iraferr())
	    return (TABLE_ERROR);

	newrec->npix[0] = npix[0];
	newrec->npix[1] = npix[1];

	/* dummy values (won't be used) */
	newrec->crpix[0] = 0.;
	newrec->crpix[1] = 0.;
	newrec->cdelt[0] = 1.;
	newrec->cdelt[1] = 1.;
	newrec->crval[0] = 0.;
	newrec->crval[1] = 0.;

	/* Insert in the coords list. */
	if (status = NewCoord (coords, newrec))
	    return (status);

	free (newrec);

	/* Allocate memory for the coefficients. */
	dist->xcoeff = calloc (MAX_NCOEFF, sizeof(double));
	dist->ycoeff = calloc (MAX_NCOEFF, sizeof(double));
	if (dist->xcoeff == NULL || dist->ycoeff == NULL) {
	    c_tbtclo (tabinfo->tp);
	    return (OUT_OF_MEMORY);
	}
	dist->allocated = 1;			/* set flag */

	/* Initialize.  It is not required that all columns be present;
	   missing columns result in default values for coefficients.
	*/

	for (k = 0;  k < MAX_NCOEFF;  k++) {
	    dist->xcoeff[k] = 0.;
	    dist->ycoeff[k] = 0.;
	}
	kx = WHICH_COEFF (1, 1);	/* D11, coefficient of x */
	ky = WHICH_COEFF (1, 0);	/* C10, coefficient of y */
	dist->xcoeff[kx] = 1.;		/* note:  just default values */
	dist->ycoeff[ky] = 1.;

	/* Copy image size. */
	dist->npix[0] = npix[0];
	dist->npix[1] = npix[1];

	/* The polynomial order is a header keyword. */
	dist->norder = c_tbhgti (tabinfo->tp, "NORDER");
	if (c_iraferr()) {
	    printf ("Warning  Can't read NORDER from IDCTAB table header.\n");
	    dist->norder = MAX_ORDER;
	    clear_cvoserr();
	}
	if (dist->norder > MAX_ORDER) {
	    printf (
	"ERROR    Polynomial order = %d in IDCTAB exceeds maximum of %d\n",
		dist->norder, MAX_ORDER);
	    return (GENERIC_ERROR_CODE);
	}

	/* Get the zero-point offsets, and convert to zero indexing. */
	c_tbegtd (tabinfo->tp, tabinfo->cp_xref, row, &dist->xref);
	if (c_iraferr())
	    return (TABLE_ERROR);
	dist->xref--;
	c_tbegtd (tabinfo->tp, tabinfo->cp_yref, row, &dist->yref);
	if (c_iraferr())
	    return (TABLE_ERROR);
	dist->yref--;

	/* Get the coefficients. */

	for (i = 0;  i <= dist->norder;  i++) {

	    for (j = 0;  j <= i;  j++) {

		k = WHICH_COEFF (i, j);

		/* C coefficient for x coordinate. */

		sprintf (colname, "C%d%d", i, j);

		c_tbcfnd1 (tabinfo->tp, colname, &cp);
		if (cp != 0) {
		    c_tbegtd (tabinfo->tp, cp, row, &dist->xcoeff[k]);
		    if (c_iraferr())
			return (TABLE_ERROR);
		}

		/* D coefficient for y coordinate. */

		sprintf (colname, "D%d%d", i, j);

		c_tbcfnd1 (tabinfo->tp, colname, &cp);
		if (cp != 0) {
		    c_tbegtd (tabinfo->tp, cp, row, &dist->ycoeff[k]);
		    if (c_iraferr())
			return (TABLE_ERROR);
		}
	    }
	}

	/* Compute a scale that preserves area at the reference pixel.
	   (kx and ky were assigned during the initialization step for
	   the coefficients.)
	*/
	temp = dist->xcoeff[kx] * dist->ycoeff[ky];
	dist->scale = 1. / sqrt (fabs (temp));

	return (0);
}

/* This routine closes the IDCTAB table. */

static int CloseDistTab (TblInfo *tabinfo) {

	c_tbtclo (tabinfo->tp);
	if (c_iraferr())
	    return (TABLE_ERROR);

	return (0);
}