/* header file for pipeline.c
*/

#define _STRINGIFY(x) #x
#define STRINGIFY(x) _STRINGIFY(x) 

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <dirent.h>
#include <fcntl.h>
#include <time.h>
#include <unistd.h>
#include <assert.h>
#include <math.h>
#include <fitsio.h>
#include <fitsio2.h>
#include <longnam.h>
#include <sys/types.h>
#include <complex.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_statistics.h>
#include <gsl/gsl_interp.h>
#include <gsl/gsl_sort.h>
#include <gsl/gsl_fit.h>
#include <gsl/gsl_multifit.h>
#include <fftw3.h>
#include <mpfit.h>

#include <memwatch.h>

#define PI 3.14159265358979

#define BINDENS 7
#define VBINDENS 5
#define FOCPLWIDTH 2048
#define QUADWIDTH 1024
#define BORDER 4
#define CUBEWIDTH 250
#define NCHAN 28
#define LAMBDA_MIN 980
#define LAMBDA_MAX 1790
#define LAMBDA_REF 1370
#define LONGLEG 10.10
#define SHORTLEG 3.25
#define PSFWIDTH 9
#define FITWIDTH 5
#define FITHEIGHT 31
#define SIMWIDTH 9
#define SIMHEIGHT 31
#define REFROW 16

#define SIZEOFFOCPL FOCPLWIDTH*FOCPLWIDTH 
#define SIZEOFCUBE CUBEWIDTH*CUBEWIDTH*NCHAN
#define HRFOCPLWIDTH FOCPLWIDTH*BINDENS

#define TRUE 1
#define FALSE 0
#define UNKNOWN -1
#define IN_DEADQUAD 3
#define IN_CORNERMASK 2
#define AT_EDGE 1
#define CLEAR 0

#define low_res 25
#define high_res 75
#define magnifyfpfactor 3
#define pixelsperlenset 10.568  
 //lenslet rotation 
#define sinr .322 
#define cosr .947

#define mag 10.568
#define sinr .322
#define cosr .947
#define no 0
#define yes 1

#define EXTRACTBOXWIDTH 9
#define filtersize  15
#define fwhm1 5
#define fwhm2 16
#define JBEGIN 1130
#define JEND 1340
#define HBEGIN 1460
#define HEND 1730
#define FITS 0
#define READDIFF 1
#define DAT 2
#define FITSMEAN 3
#define DATMEAN 4
#define HIFI 5
#define PROCD 6

#define MAXAXES        7

#define MISSINGVAL 100000
#define UNDEFINED -99999
#define e 2.71828
#define E .0024788				//smallest S/N ratio bin
#define K  8					//number of signal (S) cuts
#define M 64					//Max number of reads
#define FULL_WELL 65000
#define percent_of_full_well 0.975
#define too_slopy 2.0
#define saturated -999
#define MAXLINE 100

#define XSIZE 1024  
#define YSIZE 1024  

typedef struct FrameTok{
	char *FocplInitName;
	char *FocplRootName;
	char *FocplFinalName;
	char *FocplIntermName;
	char *FocplShortName;
	char *CubeName;
	char *CollapsedCubeName;
	char *Object;
	char *Type;
	char *Date;
	int Year;
	int Month;
	int Day;
	int crudeXoffset, crudeYoffset; 
	int fineXoffset, fineYoffset; 
	int NReads;
	int DarkSubFlag;
	int DetectorImageDone;
	float xshift;
	float yshift;
	int PeakCol;
	int PeakRow;
	int LastFlag;
	float EdgeCnt;	
	int AnchorX;
	int AnchorY;
} FrameTok;

typedef struct specsim_datapack{
	int psf_width;
	int simwidth;
	int simheight;
	int fitwidth;
	int fitheight;
	int psfwidth;
	double *focpl_cutout;
	double *pix_resp;
	double *spec_shape;
	double *J_psf;
	double *H_psf;
	fftw_complex *FT_J_psf;
	fftw_complex *FT_H_psf;
	double *J_spec_trace;
	double *H_spec_trace;
	fftw_complex *FT_J_spec_trace;
	fftw_complex *FT_H_spec_trace;
	fftw_complex *FT_J_spec_image;
	fftw_complex *FT_H_spec_image;
	double *J_spec_image;
	double *H_spec_image;
	double *binned_spec_image;
	double longleg;
	double shortleg;
} specsim_datapack;

typedef struct Coord{
	int Col;
	int Row;
} Coord;

typedef struct FocplBox{
	Coord TR;
	Coord TL;
	Coord BL;
	Coord BR;
} FocplBox;

static void printerror(char *msg, int status) {

    if (status != 0) {
        fprintf(stderr, "%s ", msg);
        fits_report_error(stderr, status);
    }
}

int compare_unsignedints (const void *a, const void *b)
{
	const unsigned int *da = (const unsigned int *) a;
	const unsigned int *db = (const unsigned int *) b;
     
	return (*da > *db) - (*da < *db);
}

int compare_floats (const void *a, const void *b)
{
	const float *da = (const float *) a;
	const float *db = (const float *) b;
     
	return (*da > *db) - (*da < *db);
}

int compare_float_ptrs (const void *a, const void *b)
{
	const float **fa = (const float **) a;
	const float **fb = (const float **) b;

	return (**fa > **fb) - (**fa < **fb);
}

double * load_simple_fits_double_data(char *fn, long *naxis, long *naxes, int *ndatas)
{

    fitsfile   *fptr;        /* fits ptr    */
    int        status=0;
    long       nelem;

    char       cmnt[128];

    float      nullval = 0.0;    /* fits arcanum */
    int        anynull = 0;    /* fits arcanum */


    char        *axisname[7] = { "NAXIS1",  "NAXIS2",  "NAXIS3",  
                    "NAXIS4",  "NAXIS5",  "NAXIS6",  "NAXIS7" };
    int        i;
    double      *ldbuf;        /* local space */


    /* Open fits file... */
    if (fits_open_file(&fptr, fn, READONLY, &status)) {
        printerror("load_simple_fits_double_data: fits_open_file: ", status);
        return NULL;
    }

    /* get a few keywords */

    /* find image dimensionality */
    if (fits_read_key(fptr, TLONG, "NAXIS", naxis, cmnt, &status)) {
        printerror("load_simple_fits_double_data: fits_read_key: NAXIS", status);
        return NULL;
    }
    /*(void) fprintf(stderr, "    NAXIS = %d\n", (int) *naxis);*/


    /* 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_double_data: fits_read_key:",
                            status);
        return NULL;
        }
        /*(void) fprintf(stderr, "    %s = %d\n", axisname[i],  (int) naxes[i]);*/
        nelem *= naxes[i];
    }
    if (*naxis > MAXAXES) {
        /*(void) fprintf(stderr, "    %d axes > allowable dimensionality (%d)\n", (int) *naxis, MAXAXES);*/
        return NULL;
        }



    /* (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 = (double *) calloc(nelem, sizeof(double))) == NULL) {
        (void) fprintf(stderr, "load_simple_fits_float_data: Out of mem for ldbuf\n");
        return NULL;
    }
     if (fits_read_img(fptr, TDOUBLE, (long)1, nelem, &nullval,
                          ldbuf, &anynull, &status)) {
        printerror("load_simple_fits_double_data: fits_read_image:",
                            status);
        return NULL;
    }


    /* close the table */
    if (fits_close_file(fptr, &status)) {
        printerror("load_simple_fits_double_data: fits_close_file: ", status);
        return NULL;
    }


	*ndatas = (int)nelem;

    return ldbuf;
}

float * load_simple_fits_float_data(char *fn, long *naxis, long *naxes, int *ndatas)
{

    fitsfile   *fptr;        /* fits ptr    */
    int        status=0;
    long       nelem;

    char       cmnt[128];

    float      nullval = 0.0;    /* fits arcanum */
    int        anynull = 0;    /* fits arcanum */


    char        *axisname[7] = { "NAXIS1",  "NAXIS2",  "NAXIS3",  
                    "NAXIS4",  "NAXIS5",  "NAXIS6",  "NAXIS7" };
    int        i;
    float      *ldbuf;        /* local space */


    /* Open fits file... */
    if (fits_open_file(&fptr, fn, READONLY, &status)) {
        printerror("load_simple_fits_float_data: fits_open_file: ", status);
        return NULL;
    }

    /* 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);*/


    /* 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 (*naxis > MAXAXES) {
        /*(void) fprintf(stderr, "    %d axes > allowable dimensionality (%d)\n", (int) *naxis, MAXAXES);*/
        return NULL;
        }



    /* (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 (fits_read_img(fptr, TFLOAT, (long)1, nelem, &nullval,
                          ldbuf, &anynull, &status)) {
        printerror("load_simple_fits_float_data: fits_read_image:",
                            status);
        return NULL;
    }


    /* close the table */
    if (fits_close_file(fptr, &status)) {
        printerror("load_simple_fits_float_data: fits_close_file: ", status);
        return NULL;
    }


	*ndatas = (int)nelem;

    return ldbuf;
}

int write_simple_fits_double_data(char *fn, long naxis, long *naxes, double *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(" = %d pixels\n", (int) nelem);


	/* open up a fits file handle.. */
	if (fits_create_file(&fptr, fn, &status)) {
       printerror("write_simple_fits_double_data: fits_create_file: ", status);
       return -1; 
	}


    /* Create fits file... */
    if (fits_create_img(fptr, DOUBLE_IMG, naxis, naxes, &status)) {
       printerror("write_simple_fits_double_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, TDOUBLE, "LUNVNVRT", &meaning, cmnt, &status)) {
        printerror("write_simple_fits_double_data: fits_update_key: NAXIS", status);
        return -1;
    }


    /* dump data to fits file */
    if (fits_write_img(fptr, TDOUBLE, fpixel, nelem, dbuf, &status)) {
        printerror("write_simple_fits_double_data: fits_write_file: ", status);
        return -1;
    }


    /* close the fits file */
    if (fits_close_file(fptr, &status)) {
        printerror("write_simple_fits_double_data: fits_close_file: ", status);
        return -1;
    }
	return nelem;
}

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(" = %d pixels\n", (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;
}

float *clean_bad_pixels(float *image, int num_images, char *maskpath)
{
	int j,k,r,c;
	int i;
	float  *dbuf1     = NULL;    /* ptr to physical space in mem with data */
	long   naxis     = -1;      /* number of axes in image: should be 2 */	
	long   naxes[2];      /* each axis length in # elements */
	int    ndatas;              /* number of elements returned  */
	float sum;
	char badpixmapfilename[100];

	sprintf(badpixmapfilename, "%s/badpixelmap.fits", maskpath);

	if ((dbuf1 = load_simple_fits_float_data(badpixmapfilename, &naxis, naxes, &ndatas)) == NULL)
			{
			(void) stderr, fprintf(stderr, "main: load_simple_fits_float_data %s failed\n", badpixmapfilename);
			return 0;
			}
	
	for( r = 0; r < FOCPLWIDTH; r++)	
		for (c = 0; c < FOCPLWIDTH; c++)
			if (dbuf1[r*FOCPLWIDTH + c] == 1)
				{
				for (i = 0; i < num_images; i ++)
					{
					sum = 0;
					for (j = -2; j < 1; j++)
						for (k = -2; k < 1; k++)
							if ( !((j == 0) && (k == 0)) )
								{
								sum += image[i*FOCPLWIDTH*FOCPLWIDTH + (r+j)*FOCPLWIDTH + c + k];
							//	if( (c == 1156) && (r == 1475) )
							//		printf("%d %d %d %f %f\n", i, j , k,image[i][(r+j)*FOCPLWIDTH + c + k], sum);
								}
					image[i*FOCPLWIDTH*FOCPLWIDTH + r*FOCPLWIDTH + c]   = sum/8;
					}
				}

free(dbuf1);

return image;
}

float ** clean_bad(float **image, int num_images, char *maskpath)
{
	int j,k,r,c;
	int i;
	float  *dbuf1     = NULL;    /* ptr to physical space in mem with data */
	long   naxis     = -1;      /* number of axes in image: should be 2 */	
	long   naxes[2];      /* each axis length in # elements */
	int    ndatas;              /* number of elements returned  */
	float sum;
	char badpixmapfilename[100];

	sprintf(badpixmapfilename, "%s/badpixelmap.fits", maskpath);

	if ((dbuf1 = load_simple_fits_float_data(badpixmapfilename, &naxis, naxes, &ndatas)) == NULL)
			{
			(void) stderr, fprintf(stderr, "main: load_simple_fits_float_data %s failed\n", badpixmapfilename);
			return 0;
			}
	
	for( r = 0; r < FOCPLWIDTH; r++)	
		for (c = 0; c < FOCPLWIDTH; c++)
			if (dbuf1[r*FOCPLWIDTH + c] == 1)
				{
				for (i = 0; i < num_images; i ++)
					{
					sum = 0;
					for (j = -2; j < 1; j++)
						for (k = -2; k < 1; k++)
							if ( !((j == 0) && (k == 0)) )
								{
								sum += image[i][(r+j)*FOCPLWIDTH + c + k];
							//	if( (c == 1156) && (r == 1475) )
							//		printf("%d %d %d %f %f\n", i, j , k,image[i][(r+j)*FOCPLWIDTH + c + k], sum);
								}
					image[i][r*FOCPLWIDTH + c]   = sum/8;
					}
				}

free(dbuf1);

return image;
}

float *clean_firstread_sats(float *image, int num_reads, int sat_level)
{
	int i, j, k,r,c;
	float sum = 0;
// take out saturated or almost saturated pixels.
	for( r = 0; r < FOCPLWIDTH; r++)	
		for (c = 0; c < FOCPLWIDTH; c++)
			if ( image[0*FOCPLWIDTH*FOCPLWIDTH + r*FOCPLWIDTH + c] > percent_of_full_well*sat_level)
				{ 
				sum = 0;
				for ( i = 0; i < num_reads; i ++)
					{
					for ( k = -1; k <= 1; k++)
						for (j=-1; j <= 1; j++)
							if ( !((i == 0 ) && (j == 0)))
							sum += image[i*FOCPLWIDTH*FOCPLWIDTH + (r+k)*FOCPLWIDTH + c + j];
					image[i*FOCPLWIDTH*FOCPLWIDTH + r*FOCPLWIDTH + c] = sum/8;
					}
				}			
	return image;
}

// take out saturated or almost saturated pixels. 	
float ** clean_sats(float **image, int num_reads)
{
	int i, j, k,r,c;
	float sum = 0;
// take out saturated or almost saturated pixels. 	
	for( r = 0; r < FOCPLWIDTH; r++)	
		for (c = 0; c < FOCPLWIDTH; c++)
			if ( image[0][r*FOCPLWIDTH + c] > percent_of_full_well*FULL_WELL)
				{ 
				sum = 0;
				for ( i = 0; i < num_reads; i ++)
					{
					for ( k = -1; k <= 1; k++)
						for (j=-1; j <= 1; j++)
							if ( !((i == 0 ) && (j == 0)))
							sum += image[i][(r+k)*FOCPLWIDTH + c + j];
					image[i][r*FOCPLWIDTH + c] = sum/8;
					}
				}			
return image;
}

float ** clean_crs(float **image, int num_images)
{ 
	int i,r,c,m, k;
	float sum, /*sum2, sd,*/ p1, p2;
	float slope;
	
	for( r = 0; r < FOCPLWIDTH; r++)	
		for (c = 0; c < FOCPLWIDTH; c++)
			{ 
			sum = 0;
	//		sum2 = 0;
		
			m = 0;
			if( (p1 = image[0][r*FOCPLWIDTH + c] ) != saturated )
				{ 
				sum = p1;
				k = 1;
				for (i = 1; i < num_images; i++, k++)
					if ((p2 = image[i][r*FOCPLWIDTH + c]) == saturated)
							{
							slope = sum/i;
							i = num_images;
							}
							else 
							sum += p2;
				slope = sum/num_images;	
				
			for (m = 0; m < num_images; m++)
					{
					p2 = image[m][r*FOCPLWIDTH + c];
						if (p2 > too_slopy * slope)
							{
							sum -= p2;
							num_images--;
							}
					}

			slope = sum/num_images;	
			for (m = 0; m < num_images; m++)
					{
					p2 = image[m][r*FOCPLWIDTH + c];
					if (p2 > too_slopy * slope)
						image[m][r*FOCPLWIDTH + c] = slope;
					}

				
				}
			}
															
	return(image);
}

float *load_dot_dat(char *filename, int num_reads)
{
	float *image;
	unsigned short int *data;
	int fd, num_floats, numfloats;
	int i, c, r, counter = 12;
	
	image = (float *) calloc(num_reads*FOCPLWIDTH*FOCPLWIDTH, sizeof(float));
	
	if ((fd = open(filename, O_RDONLY)) == -1){
		(void) printf("file %s not opened\n", filename);
		(void) exit(20);
		}

	num_floats = (int) lseek(fd, 0L, SEEK_END)/sizeof(unsigned short int);
	(void) lseek(fd, 0L, SEEK_SET);
	
//	printf("%d\n",num_floats);
	data = (unsigned short int *) calloc(num_floats, sizeof(unsigned short int));

	if ( (numfloats=read(fd, (char *)data,  num_floats*sizeof(unsigned short int) )) != num_floats*sizeof(unsigned short int))
		{
		(void) fprintf(stderr, "oops... %d bytes read\n", numfloats);
		(void) perror("data not read");
		(void) exit(25);
		}
//		printf("num chars read %d\n",numfloats);

	(void) close(fd);
	
	for(i = 0; i < num_reads; i ++)
		{
	// Skip first 8 pixels of each frame 
	// The skip 8 part is needed for NDRs but not for CDSs.

//		printf("getting frame %d\n ", i + 1);
			counter += 8;  
			for (c = 0; c < XSIZE - 1;c++)									// c is the column index for a quadrant, r is the row index FOR 2 WHILE
				{
				for (r = 0; r < YSIZE; counter += 4, r++)					// for all the columns except the last one put the data values into first_frame
					{
					image[i*FOCPLWIDTH*FOCPLWIDTH + (r + 1024)*FOCPLWIDTH + c] = (float) data[counter];				//quad 1 lower left
					image[i*FOCPLWIDTH*FOCPLWIDTH + (2047 - c)*FOCPLWIDTH + r + 1024] = (float)data[counter + 1];	//quad 2 upper upper
					image[i*FOCPLWIDTH*FOCPLWIDTH + (1023 - r) *FOCPLWIDTH + 2047 - c] = (float)data[counter + 2];	//quad 3 upper right
					image[i*FOCPLWIDTH*FOCPLWIDTH + c*FOCPLWIDTH + (1023 - r)] =(float) data[counter + 3];			//quad 4
					}
				}
				
    // now deal with the last column
			c = XSIZE - 1;
			for (r = 0; r < YSIZE - 2; counter += 4, r++)	
				{   
				image[i*FOCPLWIDTH*FOCPLWIDTH + (r + 1024)*FOCPLWIDTH + c] = (float)data[counter];				//quad 1
				image[i*FOCPLWIDTH*FOCPLWIDTH + (2047 - c)*FOCPLWIDTH + r + 1024] =(float) data[counter + 1];	//quad 2
				image[i*FOCPLWIDTH*FOCPLWIDTH + (1023 - r) *FOCPLWIDTH + 2047 - c] = (float)data[counter+ 2];	//quad 3
				image[i*FOCPLWIDTH*FOCPLWIDTH + c*FOCPLWIDTH + (1023 - r)] = (float)data[counter +  3];			//quad 4
				}
			// end for i == 0
			counter += 12;												// jump 12 values to get to the start of the next frame

		} //end for image
	//***********************************************************************************************************************************

	free(data);
	
	return image;
}

float ** process_dot_dat_file(char *filename, int num_reads)
	   // trial c conversion of ian's idl .dat file conversion idl file
{
	float **image;
	unsigned short int *data;
	int fd, num_floats, numfloats;
	int i, c, r, counter = 12;
	
	image = (float **) calloc(num_reads , sizeof(float *));
	
	if ((fd = open(filename, O_RDONLY)) == -1){
		(void) printf("file %s not opened\n", filename);
		(void) exit(20);
		}

	num_floats = (int) lseek(fd, 0L, SEEK_END)/sizeof(unsigned short int);
	(void) lseek(fd, 0L, SEEK_SET);
	
//	printf("%d\n",num_floats);
	data = (unsigned short int *) calloc(num_floats, sizeof(unsigned short int));

	if ( (numfloats=read(fd, (char *)data,  num_floats*sizeof(unsigned short int) )) != num_floats*sizeof(unsigned short int))
		{
		(void) fprintf(stderr, "oops... %d bytes read\n", numfloats);
		(void) perror("data not read");
		(void) exit(25);
		}
//		printf("num chars read %d\n",numfloats);

	(void) close(fd);
	
	for(i = 0; i < num_reads; i ++)
		{
		image[i] = (float *) calloc (FOCPLWIDTH*FOCPLWIDTH, sizeof(float));
	// Skip first 8 pixels of each frame 
	// The skip 8 part is needed for NDRs but not for CDSs.

//		printf("getting frame %d\n ", i + 1);
			counter += 8;  
			for (c = 0; c < XSIZE - 1;c++)									// c is the column index for a quadrant, r is the row index FOR 2 WHILE
				{
				for (r = 0; r < YSIZE; counter += 4, r++)					// for all the columns except the last one put the data values into first_frame
					{
					image[i][(r + 1024)*FOCPLWIDTH + c] = (float) data[counter];				//quad 1 lower left
					image[i][(2047 - c)*FOCPLWIDTH + r + 1024] = (float)data[counter + 1];	//quad 2 upper upper
					image[i][(1023 - r) *FOCPLWIDTH + 2047 - c] = (float)data[counter + 2];	//quad 3 upper right
					image[i][c*FOCPLWIDTH + (1023 - r)] =(float) data[counter + 3];			//quad 4
					}
				}
				
    // now deal with the last column
			c = XSIZE - 1;
			for (r = 0; r < YSIZE - 2; counter += 4, r++)	
				{   
				image[i][(r + 1024)*FOCPLWIDTH + c] = (float)data[counter];				//quad 1
				image[i][(2047 - c)*FOCPLWIDTH + r + 1024] =(float) data[counter + 1];	//quad 2
				image[i][(1023 - r) *FOCPLWIDTH + 2047 - c] = (float)data[counter+ 2];	//quad 3
				image[i][c*FOCPLWIDTH + (1023 - r)] = (float)data[counter +  3];			//quad 4
				}
			// end for i == 0
			counter += 12;												// jump 12 values to get to the start of the next frame

		} //end for image
	//***********************************************************************************************************************************

	free(data);
	
	return image;
}