/********************************/
/*	setcathdr		*/
/*				*/
/*	Set default catalog	*/
/********************************/
#include "focas1.h"


long	oldsz;		/* Old catalog size */

/*Set initial header	*/
setcathdr () {
    int     i, j;

    *sp.ptfl = 0;
    *sp.arfl = 0;
    for (i = 0; i < 3; i++)
	for (j = 0; j < 3; j++)
	    sp.xfrm[i][j] = 0.;
    sp.xfrm[0][0] = sp.xfrm[1][1] = sp.xfrm[2][2] = 1.;
    sp.comment = 0;
    sp.maglim = 100.;
    sp.magfst = 30.;
    sp.satr = 30000.;
    sp.sa =.1;
    sp.sb =.1;
    sp.nupdte = 3;
    sp.skwdth = 3;
    sp.skybw =.1;
    sp.skyhw = 0.;
    sp.sptbw =.2;
    sp.minudt = 500;
    sp.thrln0 = 2.5;
    sp.thrln1 = 10.;
    sp.buf = 5;
    sp.grow = 2.0;
    sp.minarea = 6;
    sp.sig = -100;
    sp.rfca = 5;
    sp.fx = sp.fy = 1;
    dfltr[0][0] = 1.;
    sp.nx = sp.ny = sp.nden = 0;
    sp.pflags = LNR;
    sp.nrules = 5;
    crules[0].mag1 = -100.;
    crules[0].mag2 = 100.;
    crules[0].scl1 =.1;
    crules[0].scl2 = 0.6;
    crules[0].frc1 = 0.;
    crules[0].frc2 = 1.;
    strcpy (crules[0].class, "n");
    crules[1].mag1 = -100.;
    crules[1].mag2 = 100.;
    crules[1].scl1 =.61;
    crules[1].scl2 = 1.2;
    crules[1].frc1 = 0.;
    crules[1].frc2 = 1.;
    strcpy (crules[1].class, "s");
    crules[2].mag1 = -100.;
    crules[2].mag2 = 100.;
    crules[2].scl1 = 1.21;
    crules[2].scl2 = 10.0;
    crules[2].frc1 = 0.;
    crules[2].frc2 = 0.24;
    strcpy (crules[2].class, "sf");
    crules[3].mag1 = -100.;
    crules[3].mag2 = 100.;
    crules[3].scl1 = 1.21;
    crules[3].scl2 = 10.0;
    crules[3].frc1 = 0.25;
    crules[3].frc2 = 1.;
    strcpy (crules[3].class, "g");
    crules[4].mag1 = -100.;
    crules[4].mag2 = 100.;
    crules[4].scl1 = 10.01;
    crules[4].scl2 = 100.;
    crules[4].frc1 = 0.;
    crules[4].frc2 = 1.;
    strcpy (crules[4].class, "d");
    ncmmnts = COLNS;
}



modcathdr () {
    int     i;
    char    in[80];
    float   magfst;

    printf ("\n  Field image file (%s): ", sp.ptfl); fflush (stdout);
    if (strlen (gets (in))) {
	strncpy (sp.ptfl, in, SZ_FNAME);
	imopen (sp.ptfl, &pthdr, 0);
	strncpy (sp.pltnm, pthdr.object, 32);
	strncpy (sp.pltepoch, pthdr.date, 32);
	strncpy (sp.observer, pthdr.observer, 24);
	strncpy (sp.origin, pthdr.origin, 24);
	sp.xs = pthdr.crpix1;
	sp.ys = pthdr.crpix2;
	sp.kpts = pthdr.naxis1;
	sp.kscan = pthdr.naxis2;
	sprintf (sp.arfl, "%s.ar", sp.ptfl);
    }
    if (feof (stdin))
	return;

    printf ("  Field name (%s): ", sp.pltnm); fflush (stdout);
    if (strlen (gets (in)))
	strncpy (sp.pltnm, in, 32);
    if (feof (stdin))
	return;
    printf ("  Field epoch (%s): ", sp.pltepoch); fflush (stdout);
    if (strlen (gets (in)))
	strncpy (sp.pltepoch, in, 32);
    if (feof (stdin))
	return;
    printf ("  Field passband (%s): ", sp.band); fflush (stdout);
    if (strlen (gets (in)))
	strncpy (sp.band, in, 2);
    if (feof (stdin))
	return;
    printf ("  Field coordinates in decimal (%g %g): ", sp.pltra, sp.pltdec);
    fflush (stdout);
    if (strlen (gets (in)))
	sscanf (in, "%hf%hf", &sp.pltra, &sp.pltdec);
    if (feof (stdin))
	return;
    printf ("  Exposure or integration (%d): ", sp.pltexp); fflush (stdout);
    if (strlen (gets (in)))
	sscanf (in, "%hd", &sp.pltexp);
    if (feof (stdin))
	return;
    printf ("  Observer (%s): ", sp.observer); fflush (stdout);
    if (strlen (gets (in)))
	strncpy (sp.observer, in, 24);
    if (feof (stdin))
	return;
    printf ("  Origin (%s): ", sp.origin); fflush (stdout);
    if (strlen (gets (in)))
	strncpy (sp.origin, in, 24);
    if (feof (stdin))
	return;
    printf ("  Saturation value (%d): ", sp.satr); fflush (stdout);
    if (strlen (gets (in)))
	sscanf (in, "%d", &sp.satr);
    if (feof (stdin))
	return;
    printf ("  Magnitude zero point (%g): ", sp.magfst); fflush (stdout);
    if (strlen (gets (in))) {
	sscanf (in, "%hf", &magfst);
	fseek (cfd, oldsz, 0);
	chmag (magfst);
	sp.magfst = magfst;
    }
    if (feof (stdin))
	return;
    printf ("  Catalog magnitude limit (%g): ", sp.maglim); fflush (stdout);
    if (strlen (gets (in)))
	sscanf (in, "%hf", &sp.maglim);
    if (feof (stdin))
	return;
    printf ("  Radius of fixed circular aperature (%g): ", sp.rfca);
    fflush (stdout);
    if (strlen (gets (in)))
	sscanf (in, "%hf", &sp.rfca);
    if (feof (stdin))
	return;
    printf ("  Sigma density above sky for detection (%g): ", sp.thrln0);
    fflush (stdout);
    if (strlen (gets (in)))
	sscanf (in, "%hf", &sp.thrln0);
    if (feof (stdin))
	return;
    printf ("  Sigma density below sky for detection (%g): ", sp.thrln1);
    fflush (stdout);
    if (strlen (gets (in)))
	sscanf (in, "%hf", &sp.thrln1);
    if (feof (stdin))
	return;
    printf ("  Sigma of sky (0 = automatic determination) (%g): ", sp.skyhw);
    fflush (stdout);
    if (strlen (gets (in)))
	sscanf (in, "%hf", &sp.skyhw);
    if (feof (stdin))
	return;
    printf ("  Minimum area for detection (%d): ", sp.minarea); fflush (stdout);
    if (strlen (gets (in)))
	sscanf (in, "%d", &sp.minarea);
    if (feof (stdin))
	return;
    printf ("  Significance level for evaluation and splitting (%g): ", sp.sig);
    fflush (stdout);
    if (strlen (gets (in)))
	sscanf (in, "%hf", &sp.sig);
    if (feof (stdin))
	return;
    printf ("  Area description filename (%s): ", sp.arfl); fflush (stdout);
    if (strlen (gets (in)))
	strncpy (sp.arfl, in, 11);
    if (access (sp.arfl, 0) == 0) {
	printf ("  Warning: File %s already exists.  Repeat to accept this filename.\n", sp.arfl);
	printf ("  Area description filename (%s): ", sp.arfl); fflush (stdout);
	if (strlen (gets (in)))
	    strncpy (sp.arfl, in, 11);
    }
    if (feof (stdin))
	return;
    printf ("  Sky updating constants (%.6g %.6g): ", sp.sa, sp.sb);
    fflush (stdout);
    if (strlen (gets (in)))
	sscanf (in, "%hf%hf", &sp.sa, &sp.sb);
    if (feof (stdin))
	return;
    printf ("  Comments:\n");
    for (i = 0; i < sp.comment % ncmmnts; i++)
	printf ("  %s\n", comments[i]);
    printf ("  Additional comments (Exit with < cr > or EOF):\n");
    fflush (stdout);
    for (;; sp.comment++)
	if (strlen (gets (comments[sp.comment % ncmmnts])) == 0)
	    break;
}

intensity () {
    int     i, j, k;
    float   temp;
    FILE * fd;
    char    str[80];

    printf ("  Linear relation? "); fflush (stdout);
    gets (str);
    if (str[0] == 'y') {
	sp.pflags &= ~INTP;
	sp.pflags &= ~DLUT;
	sp.pflags &= ~FUNC;
	sp.pflags |= LNR;
	return;
    }
    printf ("  Use calibration spots? "); fflush (stdout);
    gets (str);
    if (str[0] == 'y') {
	sp.pflags &= ~LNR;
	sp.pflags &= ~INTP;
	sp.pflags &= ~DLUT;
	sp.pflags &= ~FUNC;
	spotcalib ();
	printf ("\n");
	return;
    }
    printf ("  Input file (if < cr > input from the terminal): ");
    fflush (stdout);
    if (strlen (gets (str)) == 0)
	fd = getstdin;
    else if ((fd = fopen (str, "r")) == NULL) {
	printf ("  Can't read file %s\n", str);
	fd = getstdin;
    }
    if (fd == stdin) {
	sp.pflags &= ~LNR;
	sp.pflags &= ~INTP;
	sp.pflags &= ~FUNC;
	sp.pflags |= DLUT;
	for (sp.nden = 0; sp.nden < MXDEN; sp.nden++) {
	    printf ("  Pixel value: "); fflush (stdout);
	    if (strlen (gets (str)) == 0)
		break;
	    sscanf (str, "%d", &flmcv[sp.nden].d);
	    printf ("  Instrumental intensity: "); fflush (stdout);
	    if (strlen (gets (str)) == 0)
		continue;
	    sscanf (str, "%hf", &flmcv[sp.nden].i);
	}
    } else {
	sp.pflags &= ~LNR;
	sp.pflags &= ~INTP;
	sp.pflags &= ~FUNC;
	sp.pflags |= DLUT;
	for (sp.nden = 0; sp.nden < MXDEN; sp.nden++) {
	    fscanf (fd, "%d", &i);
	    if (feof (fd))
		break;
	    flmcv[i].d = i;
	    fscanf (fd, "%hf", &flmcv[i].i);
	}
	fclose (fd);
	return;
    }

    for (i = 0; i < sp.nden - 1; i++)
	for (j = i + 1; j < sp.nden; j++) {
	    if (flmcv[i].d > flmcv[j].d) {
		k = flmcv[i].d;
		flmcv[i].d = flmcv[j].d;
		flmcv[j].d = k;
		temp = flmcv[i].i;
		flmcv[i].i = flmcv[j].i;
		flmcv[j].i = temp;
	    }
	}
    for (i = 0; i < sp.nden - 1; i++)
	flmcv[i].slope = (flmcv[i + 1].i - flmcv[i].i) / (flmcv[i + 1].d - flmcv[i].d);
    printf ("\n");
}

/*Define the detection filter*/
dflt () {
    int     i, j, fx, fy;
    FILE * fd;
    char    str[80];

    fx = sp.fx;
    fy = sp.fy;
    printf ("  Use builtin filter? "); fflush (stdout);
    gets (str);
    if (str[0] == 'y') {
	sp.pflags |= FLTR;
	sp.fx = sp.fy = 5;
	dfltr[0][0] = 0.;
	dfltr[0][1] = 1.;
	dfltr[0][2] = 2.;
	dfltr[0][3] = 1.;
	dfltr[0][4] = 0.;
	dfltr[1][0] = 1.;
	dfltr[1][1] = 2.;
	dfltr[1][2] = 3.;
	dfltr[1][3] = 2.;
	dfltr[1][4] = 1.;
	dfltr[2][0] = 2.;
	dfltr[2][1] = 3.;
	dfltr[2][2] = 4.;
	dfltr[2][3] = 3.;
	dfltr[2][4] = 2.;
	dfltr[3][0] = 1.;
	dfltr[3][1] = 2.;
	dfltr[3][2] = 3.;
	dfltr[3][3] = 2.;
	dfltr[3][4] = 1.;
	dfltr[4][0] = 0.;
	dfltr[4][1] = 1.;
	dfltr[4][2] = 2.;
	dfltr[4][3] = 1.;
	dfltr[4][4] = 0.;
	return;
    }
    sp.pflags &= ~FLTR;
    printf ("  Input file (if < cr > input from the terminal): ");
    fflush (stdout);
    if (strlen (gets (str)) == 0)
	fd = getstdin;
    else
	if ((fd = fopen (str, "r")) == NULL) {
	    printf ("  Can't read file %s\n", str);
	    fd = getstdin;
	}
    if (fd == stdin) {
	printf ("  Size of filter (%d %d): ", fx, fy); fflush (stdout);
	if (strlen (gets (str)))
	    sscanf (str, "%d%d", &fx, &fy);
	if ((fx < 1) || (fx > NXFLT))
	    goto DFLT_ERR;
	if ((fy < 1) || (fy > NYFLT))
	    goto DFLT_ERR;
	for (i = 0; i < fy; i++) {
	    printf ("  y = %2d: ", i + 1); fflush (stdout);
	    for (j = 0; j < fx; j++)
		scanf ("%hf", &dfltr[i][j]);
	}
	getchar ();
    }
    else {
	fscanf (fd, "%d%d", &fx, &fy);
	if ((fx < 1) || (fx > NXFLT))
	    goto DFLT_ERR;
	if ((fy < 1) || (fy > NYFLT))
	    goto DFLT_ERR;
	for (i = 0; i < fy; i++)
	    for (j = 0; j < fx; j++)
		fscanf (fd, "%hf", &dfltr[i][j]);
	fclose (fd);
    }
    sp.fx = fx;
    sp.fy = fy;
    printf ("\n");
    return;
DFLT_ERR: 
    fprintf (stderr, "  Error in detection filter dimensions %d %d\n", fx, fy);
    printf ("\n");
    return;
}

/*Define the point spread function*/
psf (psffile)
char	*psffile;
{
    int     i, j, nx, ny;
    FILE * fd;
    char    str[80];

    nx = sp.nx;
    ny = sp.ny;
    if (strlen (psffile) == 0)
	fd = 0;
    else if (strcmp (psffile, "STDIN") == 0)
	fd = getstdin;
    else {
	if ((fd = fopen (psffile, "r")) == NULL) {
	    printf ("  Can't read file %s\n", psffile);
	    fd = 0;
	}
    }
    if (fd == 0) {
	printf ("  Size of template (%d %d): ", nx, ny); fflush (stdout);
	if (strlen (gets (str)))
	    sscanf (str, "%d%d", &nx, &ny);
	if ((nx < 1) || (nx > NXMAX))
	    goto TMP_ERR;
	if ((ny < 1) || (ny > NYMAX))
	    goto TMP_ERR;
	for (i = 0; i < ny; i++) {
	    printf ("  y = %2d: ", i + 1); fflush (stdout);
	    for (j = 0; j < nx; j++)
		scanf ("%hf", &template[i][j]);
	}
	getchar ();
	printf ("\n");
    } else {
	fscanf (fd, "%d%d", &nx, &ny);
	if ((nx < 1) || (nx > NXMAX))
	    goto TMP_ERR;
	if ((ny < 1) || (ny > NYMAX))
	    goto TMP_ERR;
	for (i = 0; i < ny; i++)
	    for (j = 0; j < nx; j++)
		fscanf (fd, "%hf", &template[i][j]);
	fclose (fd);
    }
    sp.nx = nx;
    sp.ny = ny;
    return;
TMP_ERR: 
    fprintf (stderr, "setcat: Error in template dimensions %d %d\n", nx, ny);
    printf ("\n");
    return;
}

/*Define the classification rules*/
clssrls () {
    char    str[80];

    printf ("  Exit with 'quit':\n"); fflush (stdout);
    str[0] = 'q';
    for (sp.nrules = 0; sp.nrules < NRULES; sp.nrules++) {
	printf ("  Minimum and maximum magnitudes (%f %f): ",
		crules[sp.nrules].mag1 = 0., crules[sp.nrules].mag2 = 100.);
	fflush (stdout);
	if (strlen (gets (str)) == 0);
	else
	    if (str[0] == 'q')
		break;
	    else
		sscanf (str, "%hf%hf", &crules[sp.nrules].mag1, &crules[sp.nrules].mag2);
	printf ("  Minimum and maximum scales (%f %f): ",
		crules[sp.nrules].scl1 = 0.1, crules[sp.nrules].scl2 = 100.);
	fflush (stdout);
	if (strlen (gets (str)) == 0);
	else
	    if (str[0] == 'q')
		break;
	    else
		sscanf (str, "%hf%hf", &crules[sp.nrules].scl1, &crules[sp.nrules].scl2);
	printf ("  Minimum and maximum fraction (%f %f): ",
		crules[sp.nrules].frc1 = 0., crules[sp.nrules].frc2 = 1.);
	fflush (stdout);
	if (strlen (gets (str)) == 0);
	else
	    if (str[0] == 'q')
		break;
	    else
		sscanf (str, "%hf%hf", &crules[sp.nrules].frc1, &crules[sp.nrules].frc2);
	strcpy (crules[sp.nrules].class, "u");
	printf ("  Class (%s): ", crules[sp.nrules].class); fflush (stdout);
	if (strlen (gets (str)) == 0);
	else
	    if (str[0] == 'q')
		break;
	    else
		sscanf (str, "%5s", crules[sp.nrules].class);
    }
    printf ("\n");
    fseek (cfd, oldsz, 0);
    chclass ();
}

trnsfrm () {
    int     i, j;
    char    str[80];
    float   xc, yc, s, b3, b5, a[2], b[2];

    printf ("  Apply 4-meter flat field corrections? "); fflush (stdout);
    gets (str);
    if (str[0] == 'y') {
	xc = lpcoords (sp.kpts / 2., 1, 1);
	yc = lpcoords (sp.kscan / 2., 2, 1);
	s = 0.02;
	b3 = -3.65e-6;
	b5 = -2.5e-11;
	printf ("  Field center scan coordinates (%g %g): ", xc, yc);
	fflush (stdout);
	if (strlen (gets (str)))
	    sscanf (str, "%hf%hf", &xc, &yc);
	printf ("  Milimeters per pixel (%g): ", s); fflush (stdout);
	if (strlen (gets (str)))
	    sscanf (str, "%hf", &s);
	printf ("  Distortion constants in mm units (%g %g): ", b3, b5);
	fflush (stdout);
	if (strlen (gets (str)))
	    sscanf (str, "%hf%hf", &b3, &b5);
	b3 *= s * s;
	b5 *= s * s * s * s;
	setflat (1, xc, yc, b3, b5);
    }
    else
	setflat (0, xc, yc, b3, b5);
    printf ("  Use catalog object reference points? ");
    fflush (stdout);
    gets (str);
    if (str[0] == 'y') {
	fseek (cfd, oldsz, 0);
	if (refpts (0))
	    return;
	fseek (cfd, oldsz, 0);
	chcoords ();
	return;
    }
    printf ("  Enter reference points manually (if yes exit with EOF)? ");
    fflush (stdout);
    gets (str);
    if (str[0] == 'y') {
	if (refpts (1))
	    return;
	fseek (cfd, oldsz, 0);
	chcoords ();
	return;
    }
    printf ("  Enter transform matrix? "); fflush (stdout);
    gets (str);
    if (str[0] == 'y') {
	for (i = 0; i < 3; i++)
	    for (j = 0; j < 3; j++)
		scanf ("%lf", &sp.xfrm[i][j]);
	getchar ();
	fseek (cfd, oldsz, 0);
	chcoords ();
	return;
    }
    fseek (cfd, oldsz, 0);
    chcoords ();
}

chcoords () {
    long    pos, err;
    struct objrec   ob;

    if ((sp.xfrm[0][0] == 0.) && (sp.xfrm[0][1] == 0.) &&
	    (sp.xfrm[0][2] == 0.) && (sp.xfrm[1][0] == 0.) &&
	    (sp.xfrm[1][1] == 0.) && (sp.xfrm[1][2] == 0.) &&
	    (sp.xfrm[2][0] == 0.) && (sp.xfrm[2][1] == 0.) &&
	    (sp.xfrm[2][2] == 0.)) {
	fprintf (stderr, "Transform not defined\n");
	return;
    }
    for (;;) {
	pos = ftell (cfd);
	if (rdcatob (cfd, 0L, &ob))
	    break;
	if (!(ob.eflgs & REFP))
	    coords (&ob.xc, &ob.ra);
	fseek (cfd, pos, 0);
	err = ftell (cfd);
	wtcat (&ob, cfd);

	/* Debugging */
	fseek (cfd, pos, 0);
	err = ftell (cfd);
	rdcatob (cfd, 0L, &ob);
	if (!(ob.eflgs & REFP))
	    if (ob.ra == 0. && ob.dec == 0.)
		printf ("%d.%-d: Coordinate error\n", ob.entnum, ob.subent);
    }
}

chclass () {
    long    pos;
    struct objrec   ob;

    for (;;) {
	pos = ftell (cfd);
	if (rdcatob (cfd, 0L, &ob))
	    break;
	classes (&ob);
	fseek (cfd, pos, 0);
	wtcat (&ob, cfd);
    }
}

chmag (magfst)
float   magfst;
{
    float   dmag;
    long    pos;
    struct objrec   ob;

    dmag = magfst - sp.magfst;
    if ((cattype == 2) || (dmag == 0.))
	return;
    for (;;) {
	pos = ftell (cfd);
	if (rdcatob (cfd, 0L, &ob))
	    break;
	if (ob.Li > 0.) {
	    ob.mag += dmag;
	    fseek (cfd, pos, 0);
	    wtcat (&ob, cfd);
	}
    }
}