#include "stdio.h"
#include "string.h"
#include "stdlib.h"
#include "dos.h"
#include "alloc.h"
#include "float.h"
#include "math.h"

char in_file[80];
char out_file[80];

#define MAXCOL 2200
#define MAXROW 2200
char buffer[MAXCOL*2];
unsigned int data_point[MAXROW];
int FITSrows, FITScols; /* rows and columns in FITS data file */

#define TRUE 1
#define FALSE 0

int c_break(void)  {
	printf("Control-break, ending program.\n");
	return(0);
	}

/* make a new file name with old+extension */
void strmfe(char *new,char *old,char *ext) {

	/* make new string from old+extension if necessary */
	while(*old != 0 && *old != '.') *new++=*old++;
	if (*old == 0) { /* need default extension */
		*new++='.';
		while(*ext) *new++=*ext++;
		*new=0;
	} else while(*old != 0) *new++=*old++;
}

/* input file: FITS Mark IV with 2043 or 2037 columns of 16-bit integers */
/* preceeded by FITS header of 80 columns and LF */
/* version 0.5: read lines of FITS header of input data file, find rows and columns */
/*   someday return eventualy row and column value */
/* for now process FITS header to get to data */
/* 28 oct 99: FITS header has several data lines, then END, then */
/* a number of blank (ASCII 20h) lines. Strategy: after END, look for no */
/* blanks. */
/*  version 0.6, 03 Nov 99: seems to produce reasonable values, */
/*           fixed bins for full range by floating increment */
/*        findings on Mark IV data disk 5:
			   col 0 modest standard deviation
			   col 1 very small deviation
			   col 2 small deviation
			   col 3 very small deviation
			   col 4 medium deviation
			   col 5 very small deviation
			   col 6-11,   modest standard deviation
			   col 2037-2040, modest st. dev
			   col 2041 very large st. dev for h3 images
			   col 2042 modest st. dev, in h4 images IDENTICAL to col 2041 */

/* ver 0.7 NOv 4 1999 varient "rows" to read and process rows. */
/* version 0.8 13 Jan 00: added 3-sigma filter for new average, st.dev */
/*            changes results file to list avr, std dev, and new ones. */

/* first, process FITS header */

void read_FITS(FILE *source) {
	int i;
	char *endptr;
	char *s_value; /*FITS value */
	long l_value;
    int endheader;

	FITSrows = 0; FITScols = 0; endheader = FALSE; i = 0;
    while ( !(feof(source)) || (endheader == TRUE) ) {
		fgets(buffer,81,source); /* read a line, terminated by newline */
		s_value = buffer + 10; /* point just past equals sign */
		i++;
		if (strncmp(buffer, "END     ",8) == 0) {
			endheader = TRUE;
            break;
        }
		else if (strncmp(buffer, "NAXIS1  ",8) == 0) {
			l_value = atol(s_value);
			FITScols = l_value;
		}
		else if (strncmp(buffer, "NAXIS2  ",8) == 0) {
			l_value = atol(s_value);
			FITSrows = l_value;
		}
	}
	if (!endheader) {
        printf("Did not find FITS record END.\n");
        exit (0);
    }
	printf("%i FITS records read to END, rows= %i, cols = %i.\n",
		   i,FITSrows,FITScols);
    /* now run through blank FITS header lines */
	while  (!(feof(source))) {
		fgets(buffer,81,source); /* read a line, terminated by newline */
		if (strncmp(buffer, "        ",8) != 0) {
			printf("Found nonblank line.\n");
			fread(buffer+80,(FITScols*2)-80,1,source); /* read rest of line */
			return;
		}
        /* found end of FITS header, nonblank data */
		/* so pass buffer to next routine */
    }
    printf("Did not find end of FITS header.\n");
    exit (0);
}


int read_data(FILE *source, FILE *dest, int row_val) {
/* read to desired row lines of input data file, find int's on line,
   assign row  values into data_point array,
   returns array and number of columns read */

	int i, cols, rows;
	int index;
	char *buf_ptr;
	unsigned int datum;
	unsigned char twochar[2];

    fprintf(dest,"%i ",row_val); /* write values to results file */

	rows = 0;
	while (!(feof(source))) { /* on entry, buffer already read */
		if (rows == row_val) {
			for (cols = 0; cols < FITScols; cols++){
				index = cols * 2;  /* index to col-th int, from zero */
				twochar[0] = buffer[index]; twochar[1] = buffer[index+1];
				datum = (unsigned int)twochar[1] +  256*(unsigned int)twochar[0];
				data_point[cols] = datum;
				/* construct a 16 bit int */
				/* assumes first byte is high, second low */
			}
			break;
		}
		if (rows == FITSrows) {
			printf("Read to last row for some reason.\n");
			exit(0);
		}
		if (fread(buffer,(FITScols*2),1,source)== NULL) {; /* read next line, terminated by newline */
			printf ("bad read of source.\n");
			exit(0);
		} else rows++;
	}
	printf("%i data points read from row %i\n",cols, rows);
	return(cols);

}

#define BINNUMB 20    /* number of histogram bins */
#define DRIFTNUMB 50   /* size of each drift average group */
#define MAXK 100  /* max number of drift groups */
void process_data(FILE *dest, int cols) {

	int i,j,k,dist[BINNUMB],binned;
	double sum, curbin;
	unsigned int mindat, maxdat, dat, minbin, bin[BINNUMB];
	float incr;

	double avr;
	double dif, dev, stdev;
	double davr, ddif[MAXK], meddat, dsum;
    double newavr, newdev, threesig, threemin, threemax, newstdev;

    /* read row data into array, find min and max */
	mindat = maxdat = (int) data_point[0];
	sum = 0.0;
	printf("Will read first %i cols of data.\n",FITScols);
	for (i=0; i<cols; i++) {
		dat = data_point[i];
		/* putw(dat, dest); */ /* use if you want to write out datapoints */
		sum += dat*1.0;
		if (mindat>dat) mindat=dat;
		if (maxdat<dat) maxdat=dat;
	}

    /* compute average, median */
	avr = (sum*1.0)/(cols*1.0); /* force float conversion */
	meddat = ((maxdat*1.0)+(mindat*1.0))/2.0;
    printf("mean of %i values is %6.2f\n",cols,avr);
    printf("min is %u, max is %u, median is %6.2f\n",mindat, maxdat,meddat);

	/*count up for deviation */
	dev = 0.0;
    for (i=0; i<cols; i++) {
		dat = data_point[i];
		dif = ((float)dat)-avr;
		dev += dif*dif; /* sum of squared differences */
	}

    j=BINNUMB; stdev = (dev/cols); stdev = sqrt(stdev);
	printf("the standard deviation is %6.2f\n", stdev);

	fprintf(dest,"%6.2f %6.2f ",avr, stdev); /* write values to results file */

	/* ver 0.7: compute histogram by std. dev. */
	/* set up bins for distribution histogram */
	/* ver 0.6: incr must be float or roundoff error accumulates */
	incr = stdev; minbin = (avr-(stdev*(BINNUMB/2))) + .5;
	if (minbin < 0) {
		printf("WARNING: 10 std dev below avr is below zero, setting min to zero.\n");
		minbin = 0;
	}
	for (i=0; i<BINNUMB; i++) {
		curbin = minbin + (int)((incr*i)+ .5);
		dist[i] = 0;
		if (curbin > 65535) bin[i] = 65535;
		else bin[i] = (unsigned int) curbin;
		if (avr - (bin[i]*1.0) < 1.0) bin[i]++;
	}
    for (i=0; i<cols; i++) {
		dat = data_point[i]; binned = FALSE;
		for (j=0; j<BINNUMB; j++) {
			if (bin[j] >=dat) { /* if datum in the bin */
				dist[j]++; binned = TRUE;     /* count it */
				break;
			}
		}
		if (binned==FALSE)	dist[BINNUMB-1]++; /* count overflow here */
	}

	printf("the distribution over %i less than or equal bins\n",BINNUMB);
	printf("from %u to %u by std. dev. %6.2f is below.\n",bin[0], bin[BINNUMB-1],incr);
	printf("note: first and last bins accumlate under and overflow.\n");
	printf("the bin with a preceeding * includes the average and below.\n");
	for (i=0; i<BINNUMB; i++) {
		printf("%i ",dist[i]);
		if (((float)bin[i+1] > avr) && ((float)bin[i] < avr)) printf("*");
		if ((float)bin[i+1] == avr) printf("*");
	}
	printf("\n");

    /* compute variation along row based on local average */
	k=0;
	/* * for i groups of size DRIFTNUMB */
    /* for (i=(DRIFTNUMB-1); i<cols; i+=DRIFTNUMB) { */
    for (i=0; i<(cols-DRIFTNUMB); i+=DRIFTNUMB) {
		dsum = 0;
		for (j=0; j<(DRIFTNUMB); j++) {
            dsum +=data_point[j+i];
        }
        davr = dsum/DRIFTNUMB;
		/* if (avr>davr) ddif[k++] = avr - davr; else ddif[k++] = davr - avr; */
		ddif[k++] = davr - avr;
		if (k>MAXK) break;
    }
    k--;
	printf("running differences between average and local average,\n");
    printf("computed every %i cols as localavr - avr is:\n",DRIFTNUMB);
	dsum = 0;
	for (i=0;i<k; i++) {
		printf("%6.2f ",ddif[i]);
		dsum += ddif[i];
	}
	davr = dsum/k;
	printf("\naverage difference is %6.2f\n",davr);

    /* recompute deviation by removing three sigma outliers */
	newavr = 0.0; newdev = 0.0; threesig = stdev * 3.0;
    threemin = avr-threesig; threemax = avr+threesig;
    j = 0;
    for (i=0; i<cols; i++) {
		dat = data_point[i];
		if ((dat >= threemin) && (dat <= threemax)) {
			j++;
            newavr+=dat;
			dif = ((float)dat)-avr;
			newdev += dif*dif; /* sum of squared differences */
		}
	}

    newavr = newavr/j;
	newstdev = (newdev/j); newstdev = sqrt(newstdev);
    printf("If %i 3-sigma outliers are removed, that leaves %i datapoints.\n",cols-j,j);
	printf("New average is %6.2f, std dev is %6.2f\n", newavr,newstdev);

    fprintf(dest,"%i %6.2f %6.2f\n",cols-j,newavr, newstdev); /* write values to results file */

}

void main(int argc,char *argv[]) {

	FILE *source, *dest;
	int cols, rownum;

	ctrlbrk(c_break);	/* allow control break to end loops */
    puts("Mkivrows.c ver 0.8 Jan 15 2000");

	if (argc >1) {
		/* make file names */
		strmfe(in_file,argv[1],"DAT"); /* input file  */
		strmfe(out_file,argv[2],"RES");  /* output file for image */
		if (argc>3) rownum = atol(argv[3]); else rownum = 0;

		/* open files */
        if((source=fopen(in_file,"rb")) != NULL) {
			/* create dest file */
			if((dest=fopen(out_file,"ab")) != NULL) {
				printf("Input=%s, Output=%s, Row=%i\n",in_file,out_file,rownum);
				fprintf(dest,"%s ",in_file); /* output filename */
				read_FITS(source);
				if (rownum > FITSrows) {
					printf("row value exceeds row size of image %i\n",FITSrows);
					exit(0);
				}
                cols = read_data(source, dest, rownum); /* returns row read */
				process_data(dest, cols);
				fclose(dest);
			} else printf("Error creating %s\n",out_file);
			fclose(source);
		} else printf("Error opening %s\n",in_file);
	} else puts("I need input filename, output filename, row number");
}