#include	"focas1.h"
#include	"match1.h"


/* MATCH -- Match objects in catalogs */

/* This array is not local because the SUN3/V3.x compiler can't handle it */
/* struct objrec   moba[PMAX][MAXMATCH]; */

match (mcatalog, catalog, ncatalogs, delta)
char	*mcatalog;
char	*catalog[];
int	*ncatalogs;
float	*delta;
{
    FILE    *fin, *fout, *fmatch, *fnomatch;
    float   mindec[PMAX], maxdec[PMAX];
    int	    index[PMAX];
    struct mentry   mentrya[PMAX];
    struct objrec   moba[PMAX][MAXMATCH];

    char    syscmd[80];
    int     i, j, k, l, nbuf1, nbuf2;

    for (k = 0; k < *ncatalogs; k++) {
	i = 'd';
	catsort (catalog[k], "match1.tmp", &i);
	fin = catopen ("match1.tmp", 0);

	if (access (mcatalog, 0) != 0) {/* create new match file */
	    printf ("match: Creating new match file %s from %s\n",
		mcatalog, catalog[k]);
	    fflush (stdout);
	    ncats = 1;
	    fout = mcatopen (mcatalog, 1);
	    cmmnt ("match");
	    wtcathdr (fout, 1);
	    mentry.entries = 1;
	    mentry.wt = 0;
	    for (;;) {
		if (rdcatob (fin, 0L, mob))
		    break;
		if (ffilter (mob))
		    wtmcat (fout, &mentry, mob);
	    }
	    fclose (fout);
	    fclose (fin);
	    continue;
	} else {
	    fmatch = mcatopen (mcatalog, 0);
	    if (ncats >= MAXMATCH)
		focaserr (1, "match:  Matched catalog full", "");
	    printf ("match: adding file %s to match file %s\n",
		catalog[k], mcatalog);
	    fflush (stdout);
	    ncats++;
	    fout = mcatopen ("match2.tmp", 1);
	    for (j = 1; j < ncats; j++) {
		rdcathdr (fmatch, 1);
		wtcathdr (fout, 1);
	    }
	    rdcathdr (fin, 0);
	    cmmnt ("match");
	    wtcathdr (fout, 1);
	    fnomatch = mcatopen ("match3.tmp", 1);
	    near (fin, fout, fmatch, fnomatch,
		mentrya, moba, mindec, maxdec, *delta);
	}
	fclose (fmatch);
	fclose (fin);
	fclose (fout);
	fclose (fnomatch);
	c_delete ("match1.tmp");
	fin = mcatopen ("match2.tmp", 0);
	fnomatch = mcatopen ("match3.tmp", 0);
	c_delete (mcatalog);
	fout = mcatopen (mcatalog, 1);
	for (j = 0; j < ncats; j++) {
	    rdcathdr (fin, 1);
	    wtcathdr (fout, 1);
	}

	nbuf1 = 0;
	for (i = 0; i < PMAX; i++) {
	    if (rdmcatob (fin, &mentrya[i], moba[i]))
		break;
	    declimits (moba[i], ncats, &mindec[i], &maxdec[i]);
	    for (j = i; j > 0; j--) {
		if (maxdec[i] >= maxdec[index[j-1]])
		    break;
		index[j] = index[j-1];
	    }
	    index[j] = i;
	    nbuf1++;
	}

	nbuf2 = 0;
	if (!rdmcatob (fnomatch, &mentry, mob))
	    nbuf2++;

	l = 0;
	while (nbuf1 || nbuf2) {
	    i = index[l];
	    if (nbuf1 && nbuf2) {
		if (maxdec[i] < mob[ncats-1].dec) {
		    wtmcat (fout, &mentrya[i], moba[i]);
		    if (rdmcatob (fin, &mentrya[i], moba[i])) {
			l++;
			nbuf1--;
		    } else {
			declimits (moba[i], ncats, &mindec[i], &maxdec[i]);
			for (j = 0; j < nbuf1-1; j++) {
			    if (maxdec[i] < maxdec[index[j+1]])
				break;
			    index[j] = index[j+1];
			}
			index[j] = i;
		    }
		} else {
		    wtmcat (fout, &mentry, mob);
		    if (rdmcatob (fnomatch, &mentry, mob))
			nbuf2--;
		}
	    } else if (nbuf1) {
		wtmcat (fout, &mentrya[i], moba[i]);
		if (rdmcatob (fin, &mentrya[i], moba[i])) {
		    l++;
		    nbuf1--;
		} else {
		    declimits (moba[i], ncats, &mindec[i], &maxdec[i]);
		    for (j = 0; j < nbuf1-1; j++) {
			if (maxdec[i] < maxdec[index[j+1]])
			    break;
			index[j] = index[j+1];
		    }
		    index[j] = i;
		}
	    } else if (nbuf2) {
		wtmcat (fout, &mentry, mob);
		if (rdmcatob (fnomatch, &mentry, mob))
		    nbuf2--;
	    }
	}
	fclose (fin);
	fclose (fout);
	fclose (fnomatch);
	c_delete ("match2.tmp");
	c_delete ("match3.tmp");
    }
}


near (fin, fout, fmatch, fnomatch, mentr, m, miny, maxy, delta)
FILE    *fin, *fout, *fmatch, *fnomatch;
struct mentry   mentr[PMAX];
struct objrec   m[PMAX][MAXMATCH];
float   miny[PMAX], maxy[PMAX];
float	delta;
{
    struct objrec   mtemp;
    int     i, j, n, first, last, nent, imatch, nbuf, eof = 0;
    unsigned    wt[PMAX], wt1;
    float   x, y, ymin, ymax, d, dmin;

    nent = ncats - 1;
    mentry.entries = 1 << nent;
    mentry.wt = 0;

    for (nbuf = 0; nbuf < PMAX; nbuf++) {/* Fill wraparound buffer */
	wt[nbuf] = 0;
	if (rdmcatob (fmatch, &mentr[nbuf], m[nbuf]))
	    break;
	declimits (m[nbuf], nent, &miny[nbuf], &maxy[nbuf]);
    }

    first = 0;
    last = nbuf - 1;

 /* Read through input catalog */
    for (;;) {
	if (rdcatob (fin, 0L, &mob[nent]))
	    break;
	if (ffilter (&mob[nent]) == 0)
	    continue;

A: 
	dmin = 1000000000.0;
	ymin = mob[nent].dec - delta - 1;/* Band of possible matches */
	ymax = mob[nent].dec + delta + 1;

	for (i = first; nbuf != 0; i = (++i) % PMAX) {
				/* Test for matches in buffer */

	    if ((maxy[i] < ymin) && (i == first)) {
		wtmcat (fout, &mentr[i], m[i]);
		wt[i] = 0;
		first = (i + 1) % PMAX;
		if (rdmcatob (fmatch, &mentr[i], m[i])) {
		    nbuf--;
		    eof = 1;
		} else {
		    declimits (m[i], nent, &miny[i], &maxy[i]);
		    last = i;
		}
		continue;
	    }
	    if (miny[i] > ymax)
		break;

	    for (j = 0; j < nent; j++) {/* Test overlapping objects */
		if (m[i][j].entnum < 0)
		    continue;
		x = mob[nent].ra - m[i][j].ra;
		y = mob[nent].dec - m[i][j].dec;
		d = sqrt (x * x + y * y);
		if ((d <= delta) && (d < dmin)) {
		    wt1 = 254 * (1 - (d / delta)) + 1;
		    if (wt1 > wt[i]) {
			dmin = d;
			imatch = i;
		    }
		}
	    }
	    if (i == last)
		break;
	}

	if ((i == last) && (eof == 0))
	    printf ("WARNING:  insufficient storage for full matching\n");


	if (dmin <= delta) {	/* Match found */

	/* Calculate weight */
	    for (j = 0, d = 0., n = 0; j < nent; j++) {
		if (m[imatch][j].entnum < 0)
		    continue;
		x = mob[nent].ra - m[imatch][j].ra;
		y = mob[nent].dec - m[imatch][j].dec;
		d += x * x + y * y;
		n++;
	    }
	    mentr[imatch].wt = 254 * (1 - (sqrt (d / n) / delta)) + 1;
	    wt[imatch] = 254 * (1 - (dmin / delta)) + 1;

	    if (mentr[imatch].entries & mentry.entries) {
		mtemp = m[imatch][nent];
		m[imatch][nent] = mob[nent];
		mob[nent] = mtemp;
		goto A;
	    } else {		/* Write match */
		mentr[imatch].entries |= mentry.entries;
		m[imatch][nent] = mob[nent];
	    }
	} else
	    wtmcat (fnomatch, &mentry, mob);/* No match */
    }

    for (i = first; nbuf != 0; i = (++i) % PMAX) {
				/* Write remaining unmatched objects */
	wtmcat (fout, &mentr[i], m[i]);
	if (rdmcatob (fmatch, &mentr[i], m[i]))
	    nbuf--;
    }
}


declimits (mob, nent, mindec, maxdec)
struct objrec  *mob;
int	nent;
float	*mindec, *maxdec;
{
	int	i;

	for (i = 0; mob[i].entnum < 0; i++);
	*mindec = *maxdec = mob[i].dec;
	for (i++; i < nent; i++) {
	    if (mob[i].entnum >= 0) {
		if (mob[i].dec < *mindec)
		    *mindec = mob[i].dec;
		if (mob[i].dec > *maxdec)
		    *maxdec = mob[i].dec;
	    }
	}
}