/* try to correlate two lists of stellar positions, finding the offset
	   which produces the most matches between stars in the two lists */

#include <stdio.h>
#include <math.h>
#include <ctype.h>

#define LINELEN      512
#define MAX_STAR    1000		/* max number of stars in each list */
#define MAX_LIST      10		/* max number of lists we can handle at once */

#define BIG       1.0e6			/* a number >> max number of rows, cols */
#define RADIUS      5.0			/* distance (in pixels) within which two */
                      			/*   stars count as 'matched' */
#define RAD2  (RADIUS*RADIUS)	/* square of the RADIUS distance */

struct s_star {
	int index;
	double row, col, peak;
};

static int num_star[MAX_LIST];	/* actual number of stars in each list */
static struct s_star star[MAX_LIST][MAX_STAR];
static struct s_star temp[MAX_STAR];


main(argc, argv)
int argc;
char *argv[];
{
	int i, j, num_lists;
	
	if (argc < 3) {
		printf("usage: match list1 list2 ... listN\n");
		exit(0);
	}
	for (i = 1; i < argc; i++) {
		if (slurp(argv[i], i - 1, &(num_star[i - 1])) != 0) {
			fprintf(stderr, "slurp returns non-zero\n");
			exit(1);
		}
	}
	num_lists = argc - 1;

	for (i = 0; i < num_lists - 1; i++) {	
		for (j = i + 1; j < num_lists; j++) {
			if (find_best(i, j) != 0) {
				fprintf(stderr, "find best fails for lists %d and %d ?!\n", i, j); 
			}
		}
	}

	exit(0);
}


	/* read index, row, col values from a STARS or PHOT output
	   file into the given row of the array.  
	   place the number of items read into the passed parameter 'num_read'.

	   return 0 if all is well, 1 if there's a
	   problem. */

int
slurp(filename, list, num_read)
char *filename;
int list, *num_read;
{
	int i;
	char line[LINELEN];
	FILE *fp;

	if ((fp = fopen(filename, "r")) == NULL) {
		fprintf(stderr, "can't open file %s\n", filename);
		return(1);
	}
	i = 0;
	skip_header(fp);
	while (fgets(line, LINELEN, fp) != NULL) {
		if (sscanf(line, "%d %lf %lf", &(star[list][i].index),
				&(star[list][i].row), &(star[list][i].col)) != 3) {
			fprintf(stderr, "error reading the following line from file %s\n",
					filename);
			fprintf(stderr, "%s\n", line);
			return(1);
		}
		i++;
	}
	*num_read = i;

	return(0);
}

	/* skip over any pseudo-FITS file header there might be */

int 
skip_header(fp)
FILE *fp;
{
	int done;
	long pos;
	char line[LINELEN];

	done = 0;
	pos = 0;
	while (done == 0) {
		pos = ftell(fp);
		if (fgets(line, LINELEN, fp) == NULL) {
			done = 1;
		}
		else if ((line[0] == '#') || (isalpha(line[0]))) {
			if (strncmp(line, "END", 3) == 0) {
				pos = ftell(fp);
				done = 1;
			}
		}
		else
			done = 1;
	}
	fseek(fp, pos, 0);
}

	

	/* find the best offset between the two given star lists.
	   when found, print it out. 

	   return 0 if all OK, or 1 if there's a problem. */

int
find_best(list1, list2)
int list1, list2;
{
	int i, j, max, n, min;
	double dr, dc, bestdr, bestdc;

	max = 0;
	min = (num_star[list1] < num_star[list2] ? num_star[list1] : num_star[list2]);
	for (i = 0; i < num_star[list1]; i++) {
		for (j = 0; j < num_star[list2]; j++) {
printf("i is %5d, j is %5d\n", i, j);
			n = try_align(list1, list2, &(star[list1][i]), &(star[list2][j]),
							&dr, &dc);
			if (n > max) {
				max = n;
				bestdr = dr;
				bestdc = dc;
			}
		}
	}
	if (n == 0)
		return(1);
	else {
		printf("lists %2d %2d  dr = %6.2lf dc = %6.2lf  %4d of %4d\n",
				list1, list2, bestdr, bestdc, max, min);
		return(0);
	}
}
	


	/* assume that the two given stars should be superposed to create
	   matching images.  using the shifts in row and col necessary to 
	   match them, compare all the other stars in both lists to find
	   the number of other matches.  return the number of 'good' matches
	   (stars which fall within one critical distance of each other),
	   as well as the average 'dr' and 'dc' for all the good matches.

	   returning 0 means that there were no matches, which shouldn't
	   occur, since we explicitly match the two given.  */

int
try_align(list1, list2, s1, s2, dr, dc)
int list1, list2;
struct s_star *s1, *s2;
double *dr, *dc;
{
	int i, j, k, num_match, bestj;
	double x, y, z, sumdr, sumdc, tdr, tdc, bestdr, bestdc;

	x = s1->row - s2->row;
	y = s1->col - s2->col;
	for (i = 0; i < num_star[list2]; i++) {
		temp[i].index = star[list2][i].index;
		temp[i].row = star[list2][i].row + x;
		temp[i].col = star[list2][i].col + y;
		temp[i].peak = star[list2][i].peak;
	}
		
	sumdr = 0.0;
	sumdc = 0.0;
	num_match = 0;
	for (i = 0; i < num_star[list1]; i++) {
		bestdr = BIG;
		bestdc = BIG;
		bestj = 0;
		for (j = 0; j < num_star[list2]; j++) {
			if (is_a_match(&(star[list1][i]), &(temp[j]), &tdr, &tdc)) {
				if (tdr*tdr + tdc*tdc < bestdr*bestdr + bestdc*bestdc) {
					bestdr = tdr;
					bestdc = tdc;
					bestj = j;
				}
			}
		}
		if (bestdr < BIG) {
			sumdr += bestdr;
			sumdc += bestdc;
			num_match++;
		}
	}

	if (num_match > 0) {
		*dr = x + sumdr/(double)num_match;
		*dc = y + sumdc/(double)num_match;
	}
	return(num_match);
}

	/* given two stars, find the difference in their row and col positions,
	   placing the results in the passed parameters 'dr' and 'dc'.  Then
	   compare the distance between the two stars to the critical distance for
	   a match; if the distance is less (and we have a match), return 1.
	   otherwise, return 0. */

int
is_a_match(s1, s2, dr, dc)
struct s_star *s1, *s2;
double *dr, *dc;
{
	*dr = (s1->row - s2->row);
	*dc = (s1->col - s2->col);
	if (((*dr)*(*dr) + (*dc)*(*dc)) <= RAD2) 
		return(1);
	else
		return(0);
}