/*
 *  p1640.h
 *  
 *
 *  Created by Douglas S. Brenner on 10/6/08.
 *  Copyright 2008 __MyCompanyName__. All rights reserved.
 *
 */



#include "/cfitsio/include/fitsio.h"
#include "/cfitsio/include/fitsio2.h"
#include "/cfitsio/include/longnam.h"
#include <fcntl.h>
#include <unistd.h>
#include <math.h>
#include <sys/types.h>
#include <dirent.h>
#include <string.h>
#include <stdlib.h>
#include <math.h>
#include <sys/dir.h>
#include <sys/param.h>
#include <pthread.h>
#include "./memwatch.h"


#define fpsize 2048
#define lenslets_size  250
#define imagesize (lenslets_size*lenslets_size)
#define centersize  407
#define ctrpixel 125
#define ctrradius 10
#define low_res 25
#define high_res 67
#define magnifyfpfactor 3
#define pixelsperlenset 10.568  
 //lenslet rotation 
#define sinr .322 
#define cosr .947

#define reduced_size  175
#define red_ctr  82
#define reduced_image (reduced_size*reduced_size)
#define edge  (lenslets_size-reduced_size)/2

#define mag 10.568
#define waves 23

#define no 0
#define yes 1
//#define psfdir	"/extract/psfs/"
#define psfdir "/proj1640_pipeline-debug5/library/laser/psf/"
#define psf24dir "/data/psfs24/"

#define psf30dir "/data/psfs30/"
#define psfalldir "/proj1640_pipeline-debug5/library/laser/psfall/"
#define mapdir  "/data/maps25/"

#define data "/DATA/"
#define filtersize  15
#define fwhm1 5
#define fwhm2 16
#define processed_data "/data/processed/"

#define sizeofcube lenslets_size*lenslets_size*waves

#define MAXAXES   3     
#define Quad_Threads 4



float * opencube (char *av);
float * makefilter();
long naxis;
long naxes[MAXAXES];
int ndatas;
int row, col,fpcol, fprow;

struct input{
        float  *varimage;
        float  *fixedimage;
        int r1;
        int r2;
        int c1;
        int c2;
        int quart;
        int magnif;
        int xs;
        int ys;
        float sd;
        } ;

void ludcmp( float **a, int n, int *indx, float *d);
void lubksb(float **a, int n, int *indx, float b[]);
void inverse(float  **a, float **b, int N);
float ** MatrixMult(float **a, float **b, int N);
float ABAT(float **a, float *b, int n);

int findlaserspots(char * infile, char * outfile);
int makelaserpsf(char * spotfile, char *infile, char * outfile);
int centroidpsfs (char * spotfile, char * psffile);
void makefilters(int* fpcols, int* fprows,float ** spots, int col, int row);
//char *laserdirectory, *cmd1, *codedirectory,*psfdirectory,*darkdirectory, *infile1, *outfile1, *spotfile1, *spotsdirectory ;
//char  *cmd2,  *infile2, *outfile2, *spotfile2 ;



//void spline(float x[], float y[], int n, float yp1, float ypn, float y2[]);
//void splint(float xa[], float ya[], float y2a[], int n, float x, float *y);
//void splie2 (float x1a[], float x2a[], float **ya, int m, int n, float **y2a);
//void splin2(float x1a[], float x2a[], float **ya, float **y2a, int m, int n, float x1, float x2, float *y);


float * load_simple_fits_float_data(char *fn, long *naxis, long *naxes, int *ndatas);
int write_simple_fits_float_data(char *fn, long naxis, long *naxes, float *dbuf);
static void printerror(const char *msg, int status);


float * load_simple_fits_float_data(char *fn, long *naxis, long *naxes, int *ndatas)
{
	int DEBUG = 0;
    fitsfile   *fptr;        /* fits ptr    */
    int        status=0;
    long       nelem;
	
    char       cmnt[128];
	
    float      nullval = 0.0;    /* fits arcanum */
    int        anynull = 0;    /* fits arcanum */
	
	
    const char        *axisname[7] = { "NAXIS1",  "NAXIS2",  "NAXIS3",  
		"NAXIS4",  "NAXIS5",  "NAXIS6",  "NAXIS7" };
    int        i;
    float      *ldbuf;        /* local space */
	if(DEBUG)
		printf("%s\n", "0");	

	
    /* Open fits file... */
    if (fits_open_file(&fptr, fn, READONLY, &status)) {
        printerror("load_simple_fits_float_data: fits_open_file: ", status);
        return NULL;
    }
	if(DEBUG)
		printf("%s\n", "A");	
    /* get a few keywords */
	
    /* find image dimensionality */
    if (fits_read_key(fptr, TLONG, "NAXIS", naxis, cmnt, &status)) {
        printerror("load_simple_fits_float_data: fits_read_key: NAXIS", status);
        return NULL;
    }
    /*(void) fprintf(stderr, "    NAXIS = %d\n", (int) *naxis);*/
	if(DEBUG)
		printf("%s\n", "B");
	
    /* find image size for as many axes as exist */
    for (i=0, nelem=1; i<*naxis; i++) {
        if (fits_read_key(fptr, TLONG, axisname[i], naxes+i, cmnt, &status)) {
			printerror("load_simple_fits_float_data: fits_read_key:",
					   status);
			return NULL;
        }
        /*(void) fprintf(stderr, "    %s = %d\n", axisname[i],  (int) naxes[i]);*/
        nelem *= naxes[i];
    }
    
	if(DEBUG)
		printf("%s\n", "b");
	if (*naxis > MAXAXES) {
        /*(void) fprintf(stderr, "    %d axes > allowable dimensionality (%d)\n", (int) *naxis, MAXAXES);*/
        return NULL;
	}
	
	if(DEBUG)
		printf("%s\n", "C");
	
    /* (void) fprintf(stderr, "    There are %d pixels in %s, %d * %d is %d\n", (int)nelem, fn,
	 (int)naxes[0], (int)naxes[1],  (int)(naxes[0] * naxes[1])); */
	
	
    /* allocate zeroed space for data in the buffer ptr ldbuf */
    if ((ldbuf = (float *) calloc(nelem, sizeof(float))) == NULL) {
        (void) fprintf(stderr, "load_simple_fits_float_data: Out of mem for ldbuf\n");
        return NULL;
    }
	if(DEBUG)
		printf("%s\n","D");
	if (fits_read_img(fptr, TFLOAT, (long)1, nelem, &nullval,
					  ldbuf, &anynull, &status)) {
        printerror("load_simple_fits_float_data: fits_read_image:",
				   status);
        return NULL;
    }
	if(DEBUG)
		printf("%s\n", "E");
	
    /* close the table */
    if (fits_close_file(fptr, &status)) {
        printerror("load_simple_fits_float_data: fits_close_file: ", status);
        return NULL;
    }
	
	*ndatas = (int)nelem;
	if(DEBUG)
		printf("%s\n", "F\n");
	
    return ldbuf;
}

int write_simple_fits_float_data(char *fn, long naxis, long *naxes, float *dbuf)
{
	
    fitsfile   *fptr;        /* fits ptr    */
    int        status=0;
    long       nelem;
    long       fpixel=1;
	
    char       cmnt[128];
	
	float       meaning = 42.0;
    int        i;
	
   // (void) printf("    NAXIS = %d\n", (int) naxis);
	nelem = 1;
	for (i=0; i<naxis; i++) {
		nelem *= naxes[i];
    	(void) printf("%d.", (int) naxes[i]);
	}
    (void) printf(" = %ld pixels\n", (long int) nelem);
	
	
	/* open up a fits file handle.. */
	if (fits_create_file(&fptr, fn, &status)) {
		printerror("write_simple_fits_float_data: fits_create_file: ", status);
		return -1; 
	}
	
	
    /* Create fits file... */
    if (fits_create_img(fptr, FLOAT_IMG, naxis, naxes, &status)) {
		printerror("write_simple_fits_float_data: fits_create_file: ", status);
		return -1; }
	
    /* set a keyword value for fun */
	strcpy(cmnt, "Life, the Universe and everything");
    if (fits_update_key(fptr, TFLOAT, "LUNVNVRT", &meaning, cmnt, &status)) {
        printerror("write_simple_fits_float_data: fits_update_key: NAXIS", status);
        return -1;
    }
	
	
    /* dump data to fits file */
    if (fits_write_img(fptr, TFLOAT, fpixel, nelem, dbuf, &status)) {
        printerror("write_simple_fits_float_data: fits_write_file: ", status);
        return -1;
    }
	
	
    /* close the fits file */
    if (fits_close_file(fptr, &status)) {
        printerror("write_simple_fits_float_data: fits_close_file: ", status);
        return -1;
    }
	return nelem;
	
}


static void printerror(const char *msg, int status) {
	
    if (status != 0) {
        fprintf(stderr, "%s ", msg);
        fits_report_error(stderr, status);
    }
}



//float      *dbuf1;        /* local space */