/********************************/
/*	xfrm	valdes 1 12 83	*/
/*				*/
/*	Coordinate transform	*/
/*	routines		*/
/********************************/
#include	"focas1.h"

/* Definitions for hfti solution */

#define	PTS	1250		/* Maximum number of data points */

int     mda = PTS, mdb = PTS, m, n = 3, nb = 2, krank, ip[3];
float   tau, rnorm[2], h[3], g[3];

/* Flat field parameters */
static int  flatflag = 0;
static float    flatparam[4];

/*Determine transform from object reference points*/
refpts (type)
int     type;
{
    float   a[3][PTS], b[2][PTS], xy[2], xyf[2];
    struct objrec   ob;

 /* Collect reference points */

    if (type == 0) {		/* Use catalog reference points */
	fseek (cfd, szcathdr (), 0);
	for (m = 0; m < PTS;) {
	    if (rdcatob (cfd, 0L, &ob))
		break;
	    if (ob.eflgs & REFP) {
		flat (&ob.xc, xyf);
		a[0][m] = xyf[0];
		a[1][m] = xyf[1];
		a[2][m] = 1.;
		b[0][m] = ob.ra;
		b[1][m] = ob.dec;
		m++;
	    }
	}
    }
    else
	if (type == 1) {	/* Enter reference points */
	    for (m = 0; m < PTS; m++) {
		printf ("    X - Y scan coordinate of point  %1d: ", m + 1);
		scanf ("%hf%hf", &xy[0], &xy[1]);
		if (feof (stdin))
		    break;
		printf ("    Reference coordinate of point %1d: ", m + 1);
		scanf ("%hf%hf", &b[0][m], &b[1][m]);
		if (feof (stdin))
		    break;
		flat (xy, xyf);
		a[0][m] = xyf[0];
		a[1][m] = xyf[1];
		a[2][m] = 1.;
	    }
	}
    if (m < 3) {
	printf ("coords: Not enough reference points for solution\n");
	return (1);
    }

 /* Compute transform */

    tau = 0.;
    hfti (a, &mda, &m, &n, b, &mdb, &nb, &tau, &krank, rnorm, h, g, ip);

    sp.xfrm[0][0] = b[0][0];
    sp.xfrm[0][1] = b[0][1];
    sp.xfrm[0][2] = b[0][2];
    sp.xfrm[1][0] = b[1][0];
    sp.xfrm[1][1] = b[1][1];
    sp.xfrm[1][2] = b[1][2];
    sp.xfrm[2][0] = 0;
    sp.xfrm[2][1] = 0;
    sp.xfrm[2][2] = 0;
    return (0);
}

/*Transform scan coordinates to reference coordinates*/

coords (a, b, flag)
float   a[2];
float   b[2];
int     flag;
{
    float   c[2];

    flat (a, c);
    b[0] = c[0] * sp.xfrm[0][0] + c[1] * sp.xfrm[0][1] + sp.xfrm[0][2];
    b[1] = c[0] * sp.xfrm[1][0] + c[1] * sp.xfrm[1][1] + sp.xfrm[1][2];
}

/* Set flat field terms */
setflat (flag, xc, yc, b3, b5)
int     flag;
float   xc, yc, b3, b5;
{
    flatflag = flag;
    if (flatflag) {
	flatparam[0] = xc;
	flatparam[1] = yc;
	flatparam[2] = b3;
	flatparam[3] = b5;
    }
}

/* 4-meter Distortion corrections */
flat (a, b)
float   a[2];
float   b[2];
{
    float   x, y, r2, s;

    if (flatflag) {
	x = a[0] - flatparam[0];
	y = a[1] - flatparam[1];
	r2 = x * x + y * y;
	s = 1.+ flatparam[2] * r2 + flatparam[3] * r2 * r2;
	b[0] = x * s + flatparam[0];
	b[1] = y * s + flatparam[1];
    }
    else {
	b[0] = a[0];
	b[1] = a[1];
    }
}


/* Logical to physical coordinate transformations */

float
lpcoords (x, axis, type)
float	x;
int	axis, type;
{
    switch (type) {
    case 1:	/* logical to physical */
	switch (axis) {
	case 1:
	    return (x + sp.xs);
	case 2:
	    return (x + sp.ys);
	}
    case 2:	/* physical to logical */
	switch (axis) {
	case 1:
	    return (x - sp.xs);
	case 2:
	    return (x - sp.ys);
	}
    }
}