/**************************************************************************
 *                                                                        *
 * Copyright (c) 1996 Michael Richmond and Richard Treffers               *
 *                                                                        *
 *                    This software may be copied and distributed for     *
 *                    educational, research and not for profit services   *
 *                    provided that this copyright and statement are      *
 *                    included in all such copies.                        *
 *                                                                        *
 **************************************************************************/



/*  
	creates a radial profile of the object whose centroid is
    given on the command line from the given file 

		PROFILE file [axr=] [axc=] [radius=] [outfile=] [eef] [help] 

    for all pixels within the given radius from the centroid,
    lists the distance of the pixel and its value in a two-column
    output.

	9/22/92 - usage, exit(0) added -rrt

	3/19/92 - write eef (encircled energy function stuff ) -rrt
	4/22/94 - correct structure error.
	1/8/96  - add 'help' option.  MWR
*/

#include <stdio.h>
#include <math.h>
#include <string.h>
#include <unistd.h>
#include <stdlib.h>
#include "pcvista.h"
#include "fits.h"

#ifdef PROTO
int main(int, char **);
static void do_profile(int, int, int, int, int);
static void do_eef(int, int, int, int, int);
#else
static void do_profile();
static void do_eef();
#endif


#define RADIUS  5		/* default radius of profile in pixels */

struct d_type{
	double dist; 
	int val;
}; 
static double radius;
static int rad;
static double axr, axc;
static FITS_HANDLE fh;
static FILE *outfp;

static char *usage = 
     "usage: profile file [axr=] [axc=] [radius=] [outfile=] [eef] [help]";

int
main(argc,argv)
int argc;
char *argv[];
{
	static int i;
	static int sc, sr, nr, nc, nrow, ncol;
	char *cp, *gotit;
	int eef_flag=0;

	if (argc < 2 ) {
		error(-1, usage);
	}
	if (strcmp(argv[1], "help") == 0) {
		printf("%s\n", usage);
		exit(0);
	}

	radius = (double) RADIUS;
	axr = -1.0;
	axc = -1.0;
	outfp = stdout;

	fh = fits_open(argv[1], "r", &nrow, &ncol);
	nr = nrow;
	nc = ncol;

	for (i = 2; i < argc; i++) {
		if (keyword("help", argv[i])) {
			printf("%s\n", usage);
			exit(0);
		}
		if ((gotit = find("axr", argv[i])) != NULL) {
			axr = evaluate(gotit);
			continue;
		}
		if ((gotit = find("axc", argv[i])) != NULL) {
			axc = evaluate(gotit);
			continue;
		}
		if ((gotit = find("radius", argv[i])) != NULL) {
			radius = evaluate(gotit);
			continue;
		}
		if ((gotit = find("outfile", argv[i])) != NULL) {
			if ((outfp = fopen(gotit, "w")) == NULL)
				error(-1, "can't open output file ");
			continue;
		}
		if(keyword("eef",argv[i])){
			eef_flag=1;
			continue;
		}
		fprintf(stderr,"unknown option %s\n",argv[i]);
		exit(1);
	}
	if (axr == -1.0)
		if ((cp = getsym("axr")) != NULL)
			sscanf(cp, "%lf", &axr);
	if (axc == -1.0)
		if ((cp = getsym("axc")) != NULL)
			sscanf(cp, "%lf", &axc);

	if ((axr < 0.0) || (axr >= (double) nrow) || 
	    (axc < 0.0) || (axc >= (double) ncol))
		error(-1, "centroid is negative or outside image");
	if ((radius < 0.0) || (2*radius > (double)nrow) || (2*radius > (double)ncol))
	    error(-1, "radius is negative or larger than image");

	rad = (int) (radius + 0.5);

	/* determine the starting points */
	if ((int)(axr - radius) < 0)
		sr = 0;
	else if ((int)(axr + radius) > nrow)
		sr = nrow - rad;
	else
		sr = (int)(axr - radius);
	if ((int)(axc - radius) < 0)
		sc = 0;
	else if ((int)(axc + rad) > ncol)
		sc = ncol - rad;
	else
		sc = (int)(axc - rad);

	if(eef_flag)
		do_eef(sr, sc, nr, nc, rad);
	else
		do_profile(sr, sc, nr, nc, rad);

	exit(0);
}

	/* do the actual computations on the data  */

static void
do_profile(sr, sc, nr, nc, rad)
int sr, sc, nr, nc, rad;
{
	static int row, col, n;
	static int16 data[NMAX];
	static double drow, dcol, dist;

	n = 2*rad;

	for (row = sr; (row < sr + n) && (row < nr); row++) {
		if (fits_get_data(fh, row, 0, data, nc) == FITS_FAIL)
			error(-1, "error reading data from file");
		for (col = sc; (col < sc + n) && (col < nc); col++) {
			drow = axr - (double) row;
			dcol = axc - (double) col;
			dist = sqrt(drow*drow + dcol*dcol);
			if (dist <= radius)
				fprintf(outfp, "%6.2f  %5d\n", dist, data[col]);
		}
	}
}

static int compar();

static void
do_eef(sr, sc, nr, nc, rad)
int sr, sc, nr, nc, rad;
{
	int i,step;
	static int row, col, n,count,little_count=0;
	static int16 data[NMAX];
	static double drow, dcol, dist;
	struct d_type *d;
	double little_sum=0.,sum=0.;

	n = 2*rad;

	d=(struct d_type *)malloc(n*n*sizeof(struct d_type));
	if( d == NULL)
		error(1,"malloc error");

	for (row = sr; (row < sr + n) && (row < nr); row++) {
		if (fits_get_data(fh, row, 0, data, nc) == FITS_FAIL)
			error(-1, "error reading data from file");
		for (col = sc; (col < sc + n) && (col < nc); col++) {
			drow = axr - (double) row;
			dcol = axc - (double) col;
			dist = sqrt(drow*drow + dcol*dcol);
			if (dist <= radius){
				d[count].dist=dist;
				d[count].val=data[col];
				count++;
			}
		}
	}
	qsort(d,count,sizeof(struct d_type),compar);
	for(i=0,step=0; i<count; i++){
		sum += d[i].val;
		little_sum += d[i].val;
		little_count++;
		if(d[i].dist > (double)step){
			fprintf(outfp, "%d %6.1f  %.0f\n", 
				step, little_sum/little_count,sum);
			little_count=0;
			little_sum=0.;
			step++;
		}
	}
}

/* compar on distance alone */

static int compar(a,b)
struct d_type *a,*b;
{
	if(a->dist > b->dist)
		return 1;
	return -1;
}