#include "p1640.h"

int main (int argc, char * const argv[]) {
    
//open spot images
	int  i, n = 0, m = 0;
	FILE *spots[waves];
	char *fitsfile;
	float *spotbuf[waves];
	float *sumspots;
	struct dirent  *laserfile;
	int row, col, fprow, fpcol;
	float *spect;
	int *indx;
	float *dbuf = NULL;
	char *testfile = "/data/intermediate/leo/hr8799_O_90-ds.fits";
	float *d,  **matrixa;
	
	matrixa = (float **)calloc(1, sizeof(float*));
	for (n=0; n < (waves + 4); n ++)
		matrixa[n] = (float * )calloc((waves+ 4), sizeof(float));
	indx = (int *)calloc((waves + 4), sizeof(int));
	sumspots = (float *)calloc(3*waves,sizeof(float));
	d  = (float *)calloc(1,sizeof(float));
	spect = (float*)calloc(waves + 4, sizeof(float));
	DIR *ld;
	ld = (DIR *) calloc (1, sizeof(DIR));
	ld = opendir(psfdir);
	naxes[0]=naxes[1] = fpsize;
	

	// load test file
	if ( (dbuf = load_simple_fits_float_data(testfile, &naxis, naxes, &ndatas) ) == NULL) 
	{
		(void) stderr, fprintf(stderr, "main: load_simple_fits_float_data %s failed\n", testfile);
		return 0; 
	}
	
	//	spots = (FILE **)calloc(1,sizeof(FILE*));
	
	fitsfile = (char *) calloc(200, sizeof(char*));
	while ( (laserfile = readdir(ld)) != NULL)
	{	
		if ( strstr(laserfile->d_name, "fits") != NULL)
		{
			strcpy(fitsfile, psfdir);
			strcat(fitsfile, laserfile->d_name);
		//	printf("loaded %s\n ", fitsfile);
			if ( (spotbuf[n++] = load_simple_fits_float_data(fitsfile, &naxis, naxes, &ndatas) ) == NULL) 
			{
				(void) stderr, fprintf(stderr, "main: load_simple_fits_float_data %s failed\n", fitsfile);
				return 0; 
			}
			else 
			{
				printf("loaded %s, size %d\n ", fitsfile, ndatas);
		//		for (i = 0; i < 100; i ++)
		//			printf( "%f \n", spotbuf[n][2110466+i]);
				
			}
		
		}
		else if ( strstr(laserfile->d_name, "txt") != NULL)
		{
			strcpy(fitsfile, psfdir);
			strcat(fitsfile, laserfile->d_name);
			if ( (spots[m++] = fopen(fitsfile, "r"))== NULL) 
			{
				printf("failed to open txt file %s\n", fitsfile);
				return 0; 
			}
			else 
			{
				printf("opened %s\n ", fitsfile);
				
			}
		}

		
		
	}	// end loading files.
	
		
	// make 1 - d spectrum and sumspots
	
	printf("here\n");
	for (n=0; n<waves; n++)
		{
		if(n==0)
			{
			fscanf(spots[n],"%d %d, %d %d\n", &col, &row, &fpcol, &fprow);
			for(i=-1; i<2; i++)
				{
				spect[2] += dbuf[(fprow)*fpsize + fpcol - 1+ i];
				spect[3] += dbuf[(fprow-1)*fpsize + fpcol - 1 + i];	
				printf("%d\n",(fprow)*fpsize + fpcol - 1+ i);
				sumspots[3*n] += spotbuf[n][(fprow)*fpsize + fpcol - 1+ i];
				sumspots[3*n + 1] += 4;//spotbuf[n][(fprow - 1)*fpsize + fpcol - 1+ i];
				sumspots[3*n + 2] += 4;//spotbuf[n][(fprow -2)*fpsize + fpcol - 1 + i];
				}
			}
		/*	else if ((n>0) && n<(waves-1))
			{
				fscanf(spots[n],"%d %d, %d %d\n", &col, &row, &fpcol, &fprow);
				for(i=-1; i<2; i++)
					{
					spect[n+3] += dbuf[(fprow-1)*fpsize + fpcol - 1 + i];
					sumspots[3*n] += spotbuf[n][(fprow - 2)*fpsize + fpcol - 1+ i];
					sumspots[3*n + 1] += spotbuf[n][(fprow - 1)*fpsize + fpcol - 1+ i];
					sumspots[3*n + 2] += spotbuf[n][(fprow)*fpsize + fpcol - 1 + i];
					}

			}
			else if (n== (waves - 1))
			{
				fscanf(spots[n],"%d %d, %d %d\n", &col, &row, &fpcol, &fprow);
				for(i=-1; i<2; i++)
					{
					spect[n+3] += dbuf[(fprow-1)*fpsize + fpcol - 1 + i];
					spect[n+4] += dbuf[(fprow-2)*fpsize + fpcol - 1 + i];
					sumspots[3*n] += spotbuf[n][(fprow - 2)*fpsize + col - 1+ i];
					sumspots[3*n + 1] += spotbuf[n][(fprow - 1)*fpsize + col - 1+ i];
					sumspots[3*n + 2] += spotbuf[n][(fprow)*fpsize + col - 1 + i];
				}
			}
		//printf("%d %d, %d %d\n", col, row, fpcol, fprow);										  
*/		}
	printf("here\n");	
/*	for (n=0; n<(waves+2); n++)
			{
			matrixa[3*n+2][n + 2] = sumspots[3*n];
			matrixa[3*n+3][n + 2] = sumspots[3*n+1];
			matrixa[3*n+4][n + 2] = sumspots[3*n+2];
			}	
				
	ludcmp(matrixa, waves + 4, indx, d);
	lubksb(matrixa, waves + 4, indx, spect);
	
*/	
					
						
								
	for (n=0; n<waves + 4; n++)
		printf("%f\n", spect[n]);
	free(sumspots);
	for (n = 0; n < waves; n++)
		free(spotbuf[n]);
	free(matrixa);
	for (n = 0; n < waves + 4; n++)
		free(matrixa[n]);
	free(fitsfile);
	free(ld);
	free(d);
	free(indx);
printf("done\n ");	
//open spot map files

	
//open -ds data file
//for pixel put data into arrays
//do least square fits

	
	return 0;
}