/*
 *  pipeline.c
 *
 *  Project 1640 Data Processing Pipeline
 *  version 2.0
 *  Neil Zimmerman and Douglas S. Brenner
 *  
 *
 */

#include <pcxp.h>

#define STRSIZE 500
#define BIGSTRSIZE STRSIZE*5

int check_dir_exist(char *);
int check_file_exist(char *);
int check_file_type(char *);
int interpret_flags(int, char **, int, float *, float *, int *, int *);
void replaceblanks(char *, char *);
int getfinalcubename(char *, char *, char *, char *, char *, char *, int, int, char *);
void fixdate(char *, char *);
void abbrev_fitsname(char *, char *);
void chop_off_leading_path(char *, char *);
void copy_array_2d_double(double *, double *, int, int, int, int, int, int, int, int, int, int);
void copy_array_2d_float(float *, float *, int, int, int, int, int, int, int, int, int, int);
void copy_array_2d_float_to_double(float *, double *, int, int, int, int, int, int, int, int, int, int);
int multiply_images(float *, float *, float *, int);
int subtract_images(float *, float *, float *, int);
int apply_dewarflat(float *, float *);
int median_filter(float *, int, int);
int weightedsum_cubextract(float *, float *, float ***, float *, float *, float *, float *, float *, float *, FrameTok);
int write_cube(float *, char *);
int write_cube_header(FrameTok, char *, char *, char *, char **, char *, int);
int write_collapsed_headers(FrameTok, char *, char **, char *, int, int, char *);
int copy_fits_header(char *, char *, FrameTok, char **, char *, int);
int write_resbgmap(float *, FrameTok *, char *);
int subtract_residual_DC(float *, float *);
int write_collapsed(float *, FrameTok, int, float, float);
int make_speccal_cube(float *, float *, float *);
int smooth_badlenslets(float *, int, char, float [CUBEWIDTH][CUBEWIDTH][5]);
void load_focplsol(float ***, char *);
void make_psf(double *, double *, int);
void get_int_date(char *, int *, int *, int *);
int get_peak_region(float *, int, int *, int *);
int get_image_coords(int, int, int *, int *);
int get_cube_coords(int, int, int *, int *);
int get_peak_corr(float **, float **, int, int, int, int, int, int *, int *);
int get_crude_offset(float *, int *, int *, int, int, int, int, int, char *);
int get_fine_vert_offset(float *, float ***, float *, float *, FrameTok, double *, double *, double *, int, int *, int *);
int get_fine_horiz_offset(float *, float *, double *, int *, int, int, int, int);
int specsim(int, int, double *, double *, double **, void *);
void FT_convolve_2dreal(double *, fftw_complex *, double *, fftw_complex *, double *, fftw_complex *, int, int);
void count_badpix(int *, char *);
void load_badpix_list(int **, int, char *);
int clean_badpix(float *, int **, int);
void load_psf_params(double *, char *);

int owswitch = 0;
int fine_vert_reg_switch = 1;
int applydewarflat = 1;
int applyspeccal = 0;
int regoverride = 0;
int doallswitch = 0;
int extract_switch = 1;
int procd_switch = 0;
int darksubswitch = 1;
int localbgsubswitch = 1;
int alignanchor = 0;

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

time_t starttime = time(NULL);

char *name, *datadir, *librarydir, *pipelinedir, *irafscriptsdir, *intermsubdir, *rootdatadir, *intermdir, *procddir, *date, *cmd, *objectdir, *datedir, *auxdir, *frameseriesdir, *cubedir, *collapsedir, *aligndir, *laserdir, *psfdir, *maskdir, *darkdir, *longmaskname, *tranmaskname, *default_pedmaskname, *pedmaskname, *shiftedmaskname, *flatdir, *focplmoddir, *lensletflatname, *base_dewarflatname, *dewarflatname, *darkname; 
char *fitsname, *combinelistname, *combinename, *combinelog, *cubename, *calcubename, *fitsold, *keyword, *keyval, *comment, *record, *oldfilename, *newfilename, *shortname, *alignrefname, *lasertemplatename, *collapsedname, *dummystr, *meanfocplname, *tempstr;
char *testdir;
char *framelistname, *darklistname;
char *sepstr;
char *datestr, *framestr, *typestr, *flagstr;
//FILE *speccalfile;
FrameTok *TokArr = NULL;
DIR *dirptr;
struct dirent *entry;

fitsfile *fptr;
int status = 0;

float *focplbuf;
float *cube = NULL;
//float *speccalcube = NULL;

int i, j, p, r, length;
int imtotal, proctotal;
char *nullstr = NULL;
int dirarg, flagarg;

float shiftx, shifty;
int anchorx, anchory;

char *medianfilterfocpl = (char *)calloc(STRSIZE, sizeof(char));

long axissize_focpl[2];
axissize_focpl[0] = FOCPLWIDTH;
axissize_focpl[1] = FOCPLWIDTH;
long axissize_cube[3];
axissize_cube[0] = CUBEWIDTH;
axissize_cube[1] = CUBEWIDTH;
axissize_cube[2] = NCHAN;
int ndatas, nd;              /* number of elements returned  */
long loaded_naxis, loaded_image_axislength[2];

/*variables for new cube extraction algorithm*/
int Nbadpix; int **badpix_list = NULL;
float ***focplsol = (float ***)calloc(CUBEWIDTH, sizeof(float **));
for(i = 0; i < CUBEWIDTH; i++)
{
	focplsol[i] = (float **)calloc(CUBEWIDTH, sizeof(float *));
	for(j = 0; j < CUBEWIDTH; j++)
	{
		focplsol[i][j] = (float *)calloc(5, sizeof(float));
	}
}

float *crude_template = NULL;
float *fine_template = NULL;
float *lensletflat = NULL;
float *heightmap = NULL;
float *tiltmap = NULL;
double *pix_resp = NULL ;
double *vpix_resp = NULL;
int psf_width = 9;
double jlaser_psf_params[9];
double hlaser_psf_params[9];
double *jlaser_psf = (double *)calloc(psf_width*psf_width*VBINDENS*VBINDENS, sizeof(double));
double *hlaser_psf = (double *)calloc(psf_width*psf_width*VBINDENS*VBINDENS, sizeof(double));
float *weights_J = NULL, *weights_H = NULL;
float *master_dark = NULL;
float *master_bulbflat = NULL;

//allocate variables for cube extraction
//float spec_response[NCHAN];
float *resbgframe;
laserdir = (char *) calloc(STRSIZE, sizeof(char));
psfdir = (char *) calloc(STRSIZE, sizeof(char));
aligndir = (char *) calloc(STRSIZE, sizeof(char));
char *speccalname = (char *) calloc(STRSIZE, sizeof(char));
char *speccal_tok = (char *) calloc(STRSIZE, sizeof(char));
char *guess_spectrum_tok = (char *) calloc(STRSIZE, sizeof(char));
char *guess_spectrum_name = (char *) calloc(STRSIZE, sizeof(char));

//string variable allocation for file processing and organization 

char *calib_library_str = (char *)calloc(STRSIZE, sizeof(char));
char *focplsol_fname = (char *)calloc(STRSIZE, sizeof(char));
char *jlaser_psf_fname = (char *)calloc(STRSIZE, sizeof(char));
char *hlaser_psf_fname = (char *)calloc(STRSIZE, sizeof(char));
char *jlaser_xtractweights_fname = (char *)calloc(STRSIZE, sizeof(char));
char *hlaser_xtractweights_fname = (char *)calloc(STRSIZE, sizeof(char));
char *badpix_list_fname = (char *)calloc(STRSIZE, sizeof(char));
char *master_dark_fname = (char *)calloc(STRSIZE, sizeof(char));
char *master_bulbflat_fname = (char *)calloc(STRSIZE, sizeof(char));
char *lenslet_flat_fname = (char *)calloc(STRSIZE, sizeof(char));
char *height_map_fname = (char *)calloc(STRSIZE, sizeof(char));
char *tilt_map_fname = (char *)calloc(STRSIZE, sizeof(char));
char *crudetemplate_fname = (char *) calloc(STRSIZE, sizeof(char));
char *finetemplate_fname = (char *) calloc(STRSIZE, sizeof(char));
char *pix_resp_fname = (char *) calloc(STRSIZE, sizeof(char));
char *vpix_resp_fname = (char *) calloc(STRSIZE, sizeof(char));

testdir = (char *) calloc(STRSIZE, sizeof(char));
datadir = (char *) calloc(STRSIZE, sizeof(char));
intermsubdir = (char *) calloc(STRSIZE, sizeof(char));
objectdir = (char *) calloc(STRSIZE, sizeof(char));
auxdir = (char *) calloc(STRSIZE, sizeof(char));
datedir = (char *) calloc(STRSIZE, sizeof(char));
frameseriesdir = (char *) calloc(STRSIZE, sizeof(char));
cubedir = (char *) calloc(STRSIZE, sizeof(char));
collapsedir = (char *) calloc(STRSIZE, sizeof(char));
flatdir = (char *) calloc(STRSIZE, sizeof(char));
focplmoddir = (char *) calloc(STRSIZE, sizeof(char));
maskdir = (char *) calloc(STRSIZE, sizeof(char));
darkdir = (char *) calloc(STRSIZE, sizeof(char));
tranmaskname = (char *) calloc(STRSIZE, sizeof(char));
longmaskname = (char *) calloc(STRSIZE, sizeof(char));
pedmaskname = (char *) calloc(STRSIZE, sizeof(char));
default_pedmaskname = (char *) calloc(STRSIZE, sizeof(char));
shiftedmaskname = (char *) calloc(STRSIZE, sizeof(char));
cubename = (char *) calloc(STRSIZE, sizeof(char));
calcubename = (char *) calloc(STRSIZE, sizeof(char));
name = (char *) calloc(STRSIZE, sizeof(char));
date = (char *) calloc(STRSIZE, sizeof(char));
cmd = (char *) calloc(BIGSTRSIZE, sizeof(char));
fitsname = (char *) calloc(STRSIZE, sizeof(char));
combinelistname = (char *) calloc(STRSIZE, sizeof(char));
combinelog = (char *) calloc(STRSIZE, sizeof(char));
framelistname = (char *) calloc(STRSIZE, sizeof(char));
darklistname = (char *) calloc(STRSIZE, sizeof(char));
combinename = (char *) calloc(STRSIZE, sizeof(char));
fitsold = (char *) calloc(STRSIZE, sizeof(char));
keyword = (char *) calloc(STRSIZE, sizeof(char));
comment = (char *) calloc(STRSIZE, sizeof(char));
lensletflatname = (char *) calloc(STRSIZE, sizeof(char));
dewarflatname = (char *) calloc(STRSIZE, sizeof(char));
base_dewarflatname = (char *) calloc(STRSIZE, sizeof(char));
darkname = (char *) calloc(STRSIZE, sizeof(char));
keyval = (char *) calloc(STRSIZE, sizeof(char));
record = (char *) calloc(STRSIZE, sizeof(char));
oldfilename = (char *) calloc(STRSIZE, sizeof(char));
newfilename = (char *) calloc(STRSIZE, sizeof(char));
shortname = (char *) calloc(STRSIZE, sizeof(char));
alignrefname = (char *) calloc(STRSIZE, sizeof(char));
lasertemplatename = (char *) calloc(STRSIZE, sizeof(char));
meanfocplname = (char *) calloc(STRSIZE, sizeof(char));
rootdatadir = (char *) calloc(STRSIZE, sizeof(char));
intermdir = (char *) calloc(STRSIZE, sizeof(char));
procddir = (char *) calloc(STRSIZE, sizeof(char));
librarydir = (char *) calloc(STRSIZE, sizeof(char));
pipelinedir = (char *) calloc(STRSIZE, sizeof(char));
irafscriptsdir = (char *) calloc(STRSIZE, sizeof(char));
datestr = (char *) calloc(STRSIZE, sizeof(char));
framestr = (char *) calloc(STRSIZE, sizeof(char));
typestr = (char *) calloc(STRSIZE, sizeof(char));
flagstr = (char *) calloc(STRSIZE, sizeof(char));
sepstr = (char *) calloc(STRSIZE, sizeof(char));
collapsedname = (char *) calloc(STRSIZE, sizeof(char));
dummystr = (char *) calloc(STRSIZE, sizeof(char));
tempstr = (char *) calloc(STRSIZE, sizeof(char));

sprintf(rootdatadir, "%s", STRINGIFY(ROOTDATADIR));
if(check_dir_exist(rootdatadir) == 0)
{
	fprintf(stderr, "Exit on error: root data directory %s (specified at compilation time) does not exist.\n", rootdatadir);
	exit(-1);
}
sprintf(pipelinedir, "%s", STRINGIFY(PIPELINEDIR));
if(check_dir_exist(pipelinedir) == 0)
{
	fprintf(stderr, "Exit on error: pipeline directory %s (specified at compilation time) does not exist.\n", pipelinedir);
	exit(-1);
}
sprintf(irafscriptsdir, "%s", STRINGIFY(IRAFSCRIPTSDIR));
if(check_dir_exist(pipelinedir) == 0)
{
	fprintf(stderr, "Exit on error: pipeline iraf script directory %s (specified at compilation time) does not exist.\n", irafscriptsdir);
	exit(-1);
}
sprintf(librarydir, "%s", STRINGIFY(LIBRARYDIR));
if(check_dir_exist(librarydir) == 0)
{
	fprintf(stderr, "Exit on error: pipeline library directory %s (specified at compilation time) does not exist.\n", librarydir);
	exit(-1);
}
sprintf(intermdir, "%s/INTERMEDIATE", rootdatadir);
if(check_dir_exist(intermdir) == 0)
{
	snprintf(cmd, BIGSTRSIZE, "mkdir %s", intermdir);
	system(cmd);
}
sprintf(procddir, "%s/PROCESSED", rootdatadir);
if(check_dir_exist(procddir) == 0)
{
	snprintf(cmd, BIGSTRSIZE, "mkdir %s", procddir);
	system(cmd);
}

printf("Root data directory: %s\n", rootdatadir);
printf("Pipeline path: %s\n", pipelinedir);
printf("Pipeline iraf script path: %s\n", irafscriptsdir);
printf("Pipeline library path: %s\n", librarydir);
fflush(stdout);
flagarg = 2;
if(argc == 1 || (argc == 2 && strncmp(argv[1], "-", 1) == 0))
{
	printf("***************************************************************************************\n");
	printf("Will run a test (no raw data directory was specified at the command line)...\n");
	printf("***************************************************************************************\n");
	if(argc == 2 && strncmp(argv[1], "-", 1) == 0)
	{
		flagarg = 1;
		interpret_flags(argc, argv, flagarg, &shiftx, &shifty, &anchorx, &anchory);
	}
	sprintf(datadir, "%s/TEST/", librarydir);
}
else if(argc == 2)
{
	strcpy(datadir, argv[1]);
	dirptr = opendir(datadir);
	if(dirptr)
	{
		closedir(dirptr);
	}
	else
	{ 	
		fprintf(stderr, "Exit on error: raw data directory %s does not exist.\n", datadir);
		exit(-1);
	}
}
else if(argc >= 3)
{
	if(strncmp(argv[1], "-", 1) == 0)
	{
		dirarg = 2;
		flagarg = 1;
	}
	else if(strncmp(argv[2], "-",1) == 0)
	{
		dirarg = 1;
		flagarg = 2;
	}
	else
	{
		fprintf(stderr, "Exit on error: too many arguments.\n");
		exit(-1);
	}
	strcpy(datadir, argv[dirarg]);
	dirptr = opendir(datadir);
	if(dirptr)
	{
		closedir(dirptr);
	}
	else
	{ 	
		fprintf(stderr, "Exit on error: raw data directory %s does not exist.\n", datadir);
		exit(-1);
	}
//check user's command line flags: o = overwrite mode, f = moon flat-fielding, d = use dat files, s = focal plane shift override 
	interpret_flags(argc, argv, flagarg, &shiftx, &shifty, &anchorx, &anchory);	
}

//remove trailing forward slash of the directory given by the user
length = strlen(datadir);
if(datadir[length-1] == '/')
{
	datadir[length-1] = 0;
}

printf("Processing the data in %s\n", datadir);
chop_off_leading_path(datadir, name);
sprintf(intermsubdir, "%s/%s", intermdir, name);

if(check_dir_exist(intermsubdir) == 0 || owswitch)
{
	if(check_dir_exist(intermsubdir) != 0)
	{
		printf("Deleting files in existing intermediate directory (OVERWRITE mode is on)\n");
		sprintf(cmd, "rm %s/*.fits", intermsubdir);
		system(cmd);
	}
	else
	{
		printf("Making a new directory for intermediate files: %s\n", intermsubdir);
		sprintf(cmd, "mkdir %s", intermsubdir);
		system(cmd);
	}
}

imtotal = 0;//tally of raw images examined
proctotal = 0;//tally of raw images actually processed
dirptr = opendir(datadir);

while((entry = readdir(dirptr)) != 0)
{
	sprintf(fitsname, "%s/%s", datadir, entry->d_name);
	if(check_file_type(fitsname))
	{
		strcpy(tempstr, fitsname);
		abbrev_fitsname(tempstr, shortname);
		getfinalcubename(procddir, fitsname, newfilename, calcubename, collapsedir, dummystr, FITS, 0, shortname);
		if((check_file_exist(newfilename) == 0) || owswitch)
		{
			proctotal++;
		}
		imtotal++;
	}
}
closedir(dirptr);

if(proctotal == 0)
{
	printf("No new images to process in the given directory. Use the -o switch to overwrite existing cubes.\n");
	exit(0);
}

TokArr = (FrameTok *)malloc(imtotal*sizeof(FrameTok));
dirptr = opendir(datadir);
i = 0;
printf("\n********************************* Spectrograph Focal Plane Image Cleaning Loop ********************************\n\n");
while((entry = readdir(dirptr)) != 0)
{//make list of detector images to extract cubes from
	sprintf(fitsname, "%s/%s", datadir, entry->d_name);
	if(check_file_type(fitsname))
	{
		TokArr[i].FocplInitName = (char *) calloc(STRSIZE, sizeof(char));
		TokArr[i].FocplRootName = (char *) calloc(STRSIZE, sizeof(char));
		TokArr[i].FocplFinalName = (char *) calloc(STRSIZE, sizeof(char));
		TokArr[i].FocplIntermName = (char *) calloc(STRSIZE, sizeof(char));
		TokArr[i].FocplShortName = (char *) calloc(STRSIZE, sizeof(char));
		TokArr[i].CubeName = (char *) calloc(STRSIZE, sizeof(char));
		TokArr[i].CollapsedCubeName = (char *) calloc(STRSIZE, sizeof(char));
		TokArr[i].Object = (char *) calloc(STRSIZE, sizeof(char));
		TokArr[i].Type = (char *) calloc(STRSIZE, sizeof(char));
		TokArr[i].Date = (char *) calloc(STRSIZE, sizeof(char));
		TokArr[i].DetectorImageDone = 0;
		strcpy(TokArr[i].FocplInitName, fitsname);
		strcpy(TokArr[i].FocplFinalName, fitsname);
		abbrev_fitsname(TokArr[i].FocplInitName, TokArr[i].FocplShortName);
		printf("********************************** %s (%d of %d) *********************************\n", TokArr[i].FocplShortName, i+1, imtotal);
		sprintf(TokArr[i].FocplRootName, "%s/%s", intermsubdir, TokArr[i].FocplShortName);
		fits_open_file(&fptr, TokArr[i].FocplInitName, READONLY, &status);
		if(status)
		{
			fits_report_error(stdout, status);
		}
		//fits_read_key(fptr, TINT, "NUM_READ", &TokArr[i].NReads, comment, &status);
		fits_read_keyword(fptr, "IMG_TYPE", TokArr[i].Type, comment, &status);
		strcpy(keyword, "OBJECT");
		fits_read_keyword(fptr, keyword, keyval, comment, &status);               
		if(status)
		{
			fits_report_error(stdout, status);
			fprintf(stdout, "WARNING: No OBJECT keyword in header of %s\n", fitsname); 
			abbrev_fitsname(fitsname, TokArr[i].Object);
		}
		else
		{
			replaceblanks(keyval, TokArr[i].Object);
		}
		strcpy(keyword, "UT_DATE");
		fits_read_keyword(fptr, keyword, keyval, comment, &status);               
		if(status)
		{
			fits_report_error(stdout, status);
			fprintf(stdout, "WARNING: No UT_DATE keyword in header of %s\n", fitsname); 
			strcpy(TokArr[i].Date, "0000-00-00");
		}
		else
		{
			fixdate(keyval, TokArr[i].Date);
		}
		get_int_date(TokArr[i].Date, &TokArr[i].Year, &TokArr[i].Month, &TokArr[i].Day);
		fits_close_file(fptr, &status);
		if(i == 0)
		{
			sprintf(testdir, "%s/%.4d-%.2d", librarydir, TokArr[0].Year, TokArr[0].Month);
			if(check_dir_exist(testdir) != 0)
			{
				sprintf(calib_library_str, "%.4d-%.2d", TokArr[0].Year, TokArr[0].Month);
			}
			else
			{//default set of calibration data and focal plane model
				sprintf(calib_library_str, "2011-09");
			}
			sprintf(focplsol_fname, "%s/%s/focplsol_%s.txt", librarydir, calib_library_str, calib_library_str);
			sprintf(jlaser_psf_fname, "%s/%s/binned_dpgfit_jlaser_%s.fits", librarydir, calib_library_str, calib_library_str);
			sprintf(hlaser_psf_fname, "%s/%s/binned_dpgfit_hlaser_%s.fits", librarydir, calib_library_str, calib_library_str);
			sprintf(badpix_list_fname, "%s/%s/badpix_list_%s.txt", librarydir, calib_library_str, calib_library_str);
			sprintf(master_dark_fname, "%s/%s/master_dark_%s.fits", librarydir, calib_library_str, calib_library_str);
			sprintf(master_bulbflat_fname, "%s/%s/master_bulbflat_%s.fits", librarydir, calib_library_str, calib_library_str);
			sprintf(crudetemplate_fname, "%s/%s/detector_focpl_%s.fits", librarydir, calib_library_str, calib_library_str);
			sprintf(finetemplate_fname, "%s/%s/hires_focpl_%s.fits", librarydir, calib_library_str, calib_library_str);
			sprintf(pix_resp_fname, "%s/COMMON/pix_resp_bindens7.fits", librarydir);
			sprintf(vpix_resp_fname, "%s/COMMON/pix_resp_bindens5.fits", librarydir);
			sprintf(lenslet_flat_fname, "%s/%s/med_amp_map_%s.fits", librarydir, calib_library_str, calib_library_str);
			sprintf(height_map_fname, "%s/%s/med_height_map_%s.fits", librarydir, calib_library_str, calib_library_str);
			sprintf(tilt_map_fname, "%s/%s/med_tilt_map_%s.fits", librarydir, calib_library_str, calib_library_str);
			
			assert(check_file_exist(focplsol_fname));
			assert(check_file_exist(jlaser_psf_fname));
			assert(check_file_exist(hlaser_psf_fname));
			assert(check_file_exist(badpix_list_fname));
			assert(check_file_exist(master_dark_fname));
			assert(check_file_exist(master_bulbflat_fname));
			assert(check_file_exist(crudetemplate_fname));
			assert(check_file_exist(finetemplate_fname));
			assert(check_file_exist(pix_resp_fname));
			assert(check_file_exist(vpix_resp_fname));
			assert(check_file_exist(lenslet_flat_fname));
			assert(check_file_exist(height_map_fname));
			assert(check_file_exist(tilt_map_fname));

			Nbadpix = 0;
			count_badpix(&Nbadpix, badpix_list_fname);	
			badpix_list = (int **)calloc(Nbadpix, sizeof(int *));
			for(p = 0; p < Nbadpix; p++)
			{
				badpix_list[p] = (int *)calloc(2, sizeof(int));					
			}
			load_badpix_list(badpix_list, Nbadpix, badpix_list_fname);
			master_dark = load_simple_fits_float_data(master_dark_fname, &loaded_naxis, loaded_image_axislength, &ndatas);			
			master_bulbflat = load_simple_fits_float_data(master_bulbflat_fname, &loaded_naxis, loaded_image_axislength, &ndatas);			
		}
		//clean bad pixels, dark sub, cosmic ray clean, apply bulb flat
		sprintf(TokArr[i].FocplIntermName, "%s-interm.fits", TokArr[i].FocplRootName);
		sprintf(TokArr[i].FocplFinalName, "%s-final.fits", TokArr[i].FocplRootName);
		if(check_file_exist(TokArr[i].FocplFinalName))
		{
			if(owswitch)
			{
				sprintf(cmd, "rm %s", TokArr[i].FocplFinalName);		
				system(cmd);
			}
			else
			{
				printf("%s already processed; skipping.\n", TokArr[i].FocplShortName);
				TokArr[i].DetectorImageDone = 1;
			}
		}	
		if(!TokArr[i].DetectorImageDone)
		{
			focplbuf = load_simple_fits_float_data(fitsname, &loaded_naxis, loaded_image_axislength, &ndatas);
			printf("Cleaning bad pixels...\n");
			clean_badpix(focplbuf, badpix_list, Nbadpix);
			printf("Dark subtraction...\n");
			subtract_images(focplbuf, master_dark, focplbuf, FOCPLWIDTH);
			TokArr[i].DarkSubFlag = 1;
			printf("Applying dewar bulb flat...\n");
			apply_dewarflat(focplbuf, master_bulbflat);
			write_simple_fits_float_data(TokArr[i].FocplFinalName, 2, axissize_focpl, focplbuf);
			free(focplbuf);
		}
		/*
		printf("Cleaning cosmic rays...\n");
		snprintf(cmd, BIGSTRSIZE, "%s/cosmicrays.cl input = %s output = %s", irafscriptsdir, TokArr[i].FocplIntermName, TokArr[i].FocplFinalName);
		printf("cmd: %s\n", cmd);
		system(cmd);
		*/
		i++;
	}
}
free(master_dark);
free(master_bulbflat);
printf("\n********************************* Spectrograph Focal Plane Image Registration Loop ********************************\n\n");
load_focplsol(focplsol, focplsol_fname);
lensletflat = load_simple_fits_float_data(lenslet_flat_fname, &loaded_naxis, loaded_image_axislength, &nd);
heightmap = load_simple_fits_float_data(height_map_fname, &loaded_naxis, loaded_image_axislength, &nd);
tiltmap = load_simple_fits_float_data(tilt_map_fname, &loaded_naxis, loaded_image_axislength, &nd);
if(!regoverride)
{
	pix_resp = load_simple_fits_double_data(pix_resp_fname, &loaded_naxis, loaded_image_axislength, &ndatas);
	vpix_resp = load_simple_fits_double_data(vpix_resp_fname, &loaded_naxis, loaded_image_axislength, &ndatas);
	load_psf_params(jlaser_psf_params, jlaser_psf_fname);
	load_psf_params(hlaser_psf_params, hlaser_psf_fname);
	make_psf(jlaser_psf, jlaser_psf_params, VBINDENS);
	make_psf(hlaser_psf, hlaser_psf_params, VBINDENS);
	printf("Loading crude focal plane template %s\n", crudetemplate_fname);
	crude_template = load_simple_fits_float_data(crudetemplate_fname, &loaded_naxis, loaded_image_axislength, &ndatas);
	printf("Loading high resolution focal plane template %s\n", finetemplate_fname);
	fine_template = load_simple_fits_float_data(finetemplate_fname, &loaded_naxis, loaded_image_axislength, &ndatas);
}
for(i = 0; i < imtotal; i++)
{	//Registration loop
	printf("********************************** %s (%d of %d) *********************************\n", TokArr[i].FocplShortName, i+1, imtotal);

	if(TokArr[i].DetectorImageDone)
	{
		printf("%s already processed; skipping.\n", TokArr[i].FocplShortName);
	}
	else
	{
		focplbuf = load_simple_fits_float_data(TokArr[i].FocplFinalName, &loaded_naxis, loaded_image_axislength, &ndatas);

		if(!regoverride)
		{
			if(alignanchor)
			{
				TokArr[i].AnchorX = anchorx;
				TokArr[i].AnchorY = anchory;
				get_image_coords(TokArr[i].AnchorX, TokArr[i].AnchorY, &TokArr[i].PeakCol, &TokArr[i].PeakRow);
			}
			else
			{
				get_peak_region(focplbuf, 16, &TokArr[i].PeakCol, &TokArr[i].PeakRow);
				printf("Peak region: %d, %d\n", TokArr[i].PeakCol, TokArr[i].PeakRow);
				get_cube_coords(TokArr[i].PeakCol, TokArr[i].PeakRow, &TokArr[i].AnchorX, &TokArr[i].AnchorY);
			}
			if(i > 0)
			{
				if(strcmp(TokArr[i].Object,TokArr[i-1].Object) == 0 && strcmp(TokArr[i].Date,TokArr[i-1].Date) == 0)
				{
					if(abs(TokArr[i].PeakCol - TokArr[i-1].PeakCol) < 300 && abs(TokArr[i].PeakRow - TokArr[i-1].PeakRow) < 300)
					{
						printf("Using previously determined offsets as a guess for the alignment...\n");
						get_crude_offset(focplbuf, &TokArr[i].crudeXoffset, &TokArr[i].crudeYoffset, TokArr[i-1].crudeXoffset, TokArr[i-1].crudeXoffset, 1, TokArr[i].PeakCol, TokArr[i].PeakRow, crudetemplate_fname); 
					}
					else
					{
						get_crude_offset(focplbuf, &TokArr[i].crudeXoffset, &TokArr[i].crudeYoffset, 0, 0, 0, TokArr[i].PeakCol, TokArr[i].PeakRow, crudetemplate_fname);
					}
				}
				else
				{
					get_crude_offset(focplbuf, &TokArr[i].crudeXoffset, &TokArr[i].crudeYoffset, 0, 0, 0, TokArr[i].PeakCol, TokArr[i].PeakRow, crudetemplate_fname);
				}
			}
			else
			{
				get_crude_offset(focplbuf, &TokArr[i].crudeXoffset, &TokArr[i].crudeYoffset, 0, 0, 0, TokArr[i].PeakCol, TokArr[i].PeakRow, crudetemplate_fname);
			}
			printf("Crude registration offset solution: %d columns, %d rows\n", TokArr[i].crudeXoffset, TokArr[i].crudeYoffset);
			/*GET FINE HORIZONTAL OFFSET*/
			get_fine_horiz_offset(focplbuf, fine_template, pix_resp, &TokArr[i].fineXoffset, TokArr[p].crudeXoffset, TokArr[p].crudeYoffset, TokArr[p].PeakCol, TokArr[p].PeakRow);
			/*GET FINE VERTICAL OFFSET*/
			if(fine_vert_reg_switch)
			{
				get_fine_vert_offset(focplbuf, focplsol, heightmap, tiltmap, TokArr[i], vpix_resp, jlaser_psf, hlaser_psf, psf_width, &TokArr[i].crudeYoffset, &TokArr[i].fineYoffset);
			}
			else
			{
				TokArr[i].fineYoffset = 0;
			}
		}
		else
		{
			TokArr[i].crudeXoffset = (int)lroundf(shiftx);
			TokArr[i].crudeYoffset = (int)lroundf(shifty);
			TokArr[i].fineXoffset = (int)lroundf((shiftx - TokArr[i].crudeXoffset)*BINDENS);
			TokArr[i].fineYoffset = (int)lroundf((shifty - TokArr[i].crudeYoffset)*BINDENS);
			printf("Assuming user-specified registration offset: %0.2f columns, %0.2f rows (shift override ON)\n", shiftx, shifty);  
			printf("Coarse offset: %d, %d; fine offset: %d, %d\n", TokArr[i].crudeXoffset, TokArr[i].crudeYoffset, TokArr[i].fineXoffset, TokArr[i].fineYoffset);
		}
		TokArr[i].xshift = (float)TokArr[i].crudeXoffset + TokArr[i].fineXoffset/(float)BINDENS;
		TokArr[i].yshift = (float)TokArr[i].crudeYoffset + TokArr[i].fineYoffset/(float)BINDENS;
		printf("Binned deltaX, deltaY = %0.2f, %0.2f\n", TokArr[i].xshift, TokArr[i].yshift);
		free(focplbuf);
		TokArr[i].DetectorImageDone = 1;
		copy_fits_header(TokArr[i].FocplInitName, TokArr[i].FocplFinalName, TokArr[i], argv, datadir, flagarg);
	}
}
if(!regoverride)
{
	free(pix_resp);
	free(vpix_resp);
	free(jlaser_psf);
	free(hlaser_psf);
	free(crude_template);
	free(fine_template);
}
closedir(dirptr);

if(extract_switch)
{
	printf("\n$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ Cube Extraction Loop $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$\n\n");

	cube = (float *) calloc(SIZEOFCUBE, sizeof(float));
	/*
	if(applyspeccal)
	{
		speccalcube = (float *) calloc(SIZEOFCUBE, sizeof(float));
	}
	*/
	resbgframe = (float *)calloc(CUBEWIDTH*CUBEWIDTH, sizeof(float));

	long Nax, axlength[3];
	int nd;
	sprintf(jlaser_xtractweights_fname, "%s/%s/xtract_weights_jlaser_%s.fits", librarydir, calib_library_str, calib_library_str);
	sprintf(hlaser_xtractweights_fname, "%s/%s/xtract_weights_hlaser_%s.fits", librarydir, calib_library_str, calib_library_str);
	weights_J = load_simple_fits_float_data(jlaser_xtractweights_fname, &Nax, axlength, &nd);	
	weights_H = load_simple_fits_float_data(hlaser_xtractweights_fname, &Nax, axlength, &nd);	
	free(jlaser_xtractweights_fname);
	free(hlaser_xtractweights_fname);

	/* to be fixed
	//read in spectral response
	sprintf(speccalname, "%s/SPECCAL/speccal_phaseI.txt", librarydir);
	speccalfile = fopen(speccalname, "r");
	fgets(speccal_tok, STRSIZE, speccalfile);
	speccal_tok[strlen(speccal_tok)-1] = 0;//get rid of newline character
	printf("Spectral calibration reference source: %s\n", speccal_tok);
	for(i=0; i < NCHAN; i++)
	{
		fscanf(speccalfile, "%f\n", &spec_response[i]); 	
	//	printf("channel %d: %.3f\n", i+1, spec_response[i]);
	}
	fclose(speccalfile);
	*/

	for(r=0; r<imtotal; r++)
	{	//make cubes out of the processed data focal planes
		printf("********************************** %s (%d of %d) *********************************\n", TokArr[r].FocplShortName, r+1, imtotal);
		getfinalcubename(procddir, TokArr[r].FocplInitName, TokArr[r].CubeName, calcubename, collapsedir, TokArr[r].CollapsedCubeName, FITS, 0, TokArr[r].FocplShortName);
		if(((check_file_exist(TokArr[r].CubeName) == 0) || owswitch) && TokArr[r].DetectorImageDone)
		{	//load processed detector image
			if ((focplbuf = load_simple_fits_float_data(TokArr[r].FocplFinalName, &loaded_naxis, loaded_image_axislength, &ndatas)) == 0)
			{
				(void) stdout, fprintf(stdout, "main: load_simple_fits_float_data %s failed\n", TokArr[r].FocplFinalName);
				return -1;
			}
			//clear cube array, extract new cube
			memset(cube, 0, SIZEOFCUBE*sizeof(float));
			if(weightedsum_cubextract(focplbuf, cube, focplsol, weights_J, weights_H, lensletflat, heightmap, tiltmap, resbgframe, TokArr[r]) != 0)
			{
				fprintf(stdout, "Exit on error, main: extract_cube for %s failed\n", TokArr[r].FocplFinalName);
				return -1;
			}
			free(focplbuf);
			/*
			smooth_badlenslets(cube, TokArr[r].NReads, TokArr[r].Type[0], default_lookuptable);//replace missing lenslet pixels with average of neighbors
			*/

			//prepare directories and set cube file name
			sprintf(objectdir, "%s/%s", procddir, TokArr[r].Object);
			if(check_dir_exist(objectdir) == 0)
			{
				sprintf(cmd, "mkdir %s", objectdir);
				system(cmd);
			}
			sprintf(datedir, "%s/%s", objectdir, TokArr[r].Date);
			if(check_dir_exist(datedir) == 0)
			{
				sprintf(cmd, "mkdir %s", datedir);
				system(cmd);
			}
			sprintf(cubedir, "%s/CUBES", datedir);
			if(check_dir_exist(cubedir) == 0)
			{
				sprintf(cmd, "mkdir %s", cubedir);
				system(cmd);
			}
			sprintf(auxdir, "%s", intermsubdir);
			if(check_dir_exist(auxdir) == 0)
			{
				sprintf(cmd, "mkdir %s", auxdir);
				system(cmd);
			}
			//write cube to .fits file; apply spectral calibration if requested
			if(write_cube(cube, TokArr[r].CubeName) == 0)
			{
				/*
				if(applyspeccal)
				{//make cube with spectral calibration and write to .fits file
					printf("Applying crude spectral calibration...\n");
					make_speccal_cube(cube, speccalcube, spec_response);
					write_cube(speccalcube, calcubename);
				}
				*/
				nullstr = NULL;
				write_cube_header(TokArr[r], TokArr[r].CubeName, calib_library_str, nullstr, argv, datadir, flagarg);
				if(applyspeccal)
				{
					write_cube_header(TokArr[r], calcubename, calib_library_str, speccal_tok, argv, datadir, flagarg);
				}
				write_resbgmap(resbgframe, &TokArr[r], auxdir);
			}
			else
			{
				fprintf(stderr, "Problem writing cube file %s\n", TokArr[r].CubeName);
			}
			fflush(stdout);
		} //end of "if" block checking for processed detector image
		else
		{
			printf("%s already exists, skipping. Use the -o command line switch to overwrite.\n", TokArr[r].CubeName);
		}	
	} //end of cube extraction loop (r)
	free(cube);
	free(weights_J);
	free(weights_H);
	free(resbgframe);
	/*	
	if(applyspeccal)
	{
		free(speccalcube);
	}
	*/
}

time_t endtime = time(NULL);

double exectime = difftime(endtime, starttime);
double avgexectime = 0;
int hrs, mins;
double secs;
hrs = exectime/3600;
if(hrs > 0)
{
	mins = (int)exectime%60;
}
else
{
	mins = (int)exectime/60;
}
if(mins > 0)
{
	secs = exectime - 60*mins; 
}
else
{
	secs = exectime;
}
if(proctotal > 0)
{
	avgexectime = exectime/proctotal;
}

printf("***************************************************************************************\n");
printf("The Project 1640 Pipeline has finished processing the data in %s.\n", datadir);
printf("***************************************************************************************\n");
if(proctotal == 0)
{
	printf("No unprocessed data in %s. To overwrite existing cubes, use the -o command line switch.\n", datadir);	
}
else if(extract_switch)
{ //list completed data cubes
	printf("Created the following data cubes:\n");
	for(i=0; i<proctotal; i++)
	{
		printf("%s\n", TokArr[i].CubeName);
	}
}
printf("\n");
if(proctotal > 0)
{
	printf("Processed %d detector images in %dh %dm %1.1fs (an average of %1.1f seconds per image).\n", proctotal, hrs, mins, secs, avgexectime);
}
else
{
	printf("Processed %d focal planes in %1.1f seconds.\n", proctotal, exectime);
}
printf("Freeing allocated memory...\n");

if(lensletflat != NULL)
{
	free(lensletflat);
}
if(heightmap != NULL)
{
	free(heightmap);
}
if(tiltmap != NULL)
{
	free(tiltmap);
}

for(p=0; p < imtotal; p++)
{
	free(TokArr[p].FocplInitName);
	free(TokArr[p].FocplRootName);
	free(TokArr[p].FocplFinalName);
	free(TokArr[p].FocplIntermName);
	free(TokArr[p].CubeName);
	free(TokArr[p].CollapsedCubeName);
	free(TokArr[p].Object);
	free(TokArr[p].Type);
	free(TokArr[p].Date);
}
free(TokArr);

free(crudetemplate_fname);
free(finetemplate_fname);
free(lenslet_flat_fname);
free(height_map_fname);
free(tilt_map_fname);
free(testdir);
free(medianfilterfocpl);
free(laserdir);
free(psfdir);
free(aligndir);
free(speccalname);
free(speccal_tok);
free(datadir);
free(intermsubdir);
free(objectdir);
free(auxdir);
free(datedir);
free(frameseriesdir);
free(cubedir);
free(collapsedir);
free(flatdir);
free(focplmoddir);
free(maskdir);
free(tranmaskname);
free(longmaskname);
free(pix_resp_fname);
free(vpix_resp_fname);
free(pedmaskname);
free(default_pedmaskname);
free(shiftedmaskname);
free(cubename);
free(calcubename);
free(name);
free(date);
free(cmd);
free(fitsname);
free(combinelistname);
free(combinelog);
free(framelistname);
free(darklistname);
free(combinename);
free(fitsold);
free(keyword);
free(comment);
free(dewarflatname);
free(base_dewarflatname);
free(lensletflatname);
free(keyval);
free(record);
free(oldfilename);
free(newfilename);
free(lasertemplatename);
free(meanfocplname);
free(shortname);
free(alignrefname);
free(rootdatadir);
free(intermdir);
free(procddir);
free(librarydir);
free(pipelinedir);
free(irafscriptsdir);
free(datestr);
free(framestr);
free(typestr);
free(flagstr);
free(sepstr);
free(collapsedname);
free(dummystr);
free(tempstr);
free(darkname);
free(darkdir);
free(guess_spectrum_tok);
free(guess_spectrum_name);

return 0;

}

int subtract_residual_DC(float *cube, float *resbgmap)
{
	float quadI_med, quadII_med, quadIII_med, quadIV_med, avg;
	int boxwidth = 7;
	int boxsize = boxwidth*boxwidth;

	float *sample = (float *)calloc(boxsize, sizeof(float)); 

	int row, col, index;

	int quadI_xbegin = 182;
	int quadI_ybegin = 238;	

	int quadII_xbegin = 9;
	int quadII_ybegin = 181;
	
	int quadIII_xbegin = 64;
	int quadIII_ybegin = 8;

	int quadIV_xbegin = 237;
	int quadIV_ybegin = 66;

//get quad I estimate	
	index = 0;
	for(row = quadI_ybegin; row < quadI_ybegin + boxwidth; row++)
	{
		for(col = quadI_xbegin; col < quadI_xbegin + boxwidth; col++)
		{
			sample[index] = resbgmap[row*CUBEWIDTH + col];	 
			index++;
		}
	}
	qsort(sample, boxsize, sizeof(float), compare_floats);
	quadI_med = sample[boxsize/2]; 
//get quad II estimate	
	index = 0;
	for(row = quadII_ybegin; row < quadII_ybegin + boxwidth; row++)
	{
		for(col = quadII_xbegin; col < quadII_xbegin + boxwidth; col++)
		{
			sample[index] = resbgmap[row*CUBEWIDTH + col];	 
			index++;
		}
	}
	qsort(sample, boxsize, sizeof(float), compare_floats);
	quadII_med = sample[boxsize/2]; 
//get quad III estimate	
	index = 0;
	for(row = quadIII_ybegin; row < quadIII_ybegin + boxwidth; row++)
	{
		for(col = quadIII_xbegin; col < quadIII_xbegin + boxwidth; col++)
		{
			sample[index] = resbgmap[row*CUBEWIDTH + col];	 
			index++;
		}
	}
	qsort(sample, boxsize, sizeof(float), compare_floats);
	quadIII_med = sample[boxsize/2]; 
//get quad IV estimate	
	index = 0;
	for(row = quadIV_ybegin; row < quadIV_ybegin + boxwidth; row++)
	{
		for(col = quadIV_xbegin; col < quadIV_xbegin + boxwidth; col++)
		{
			sample[index] = resbgmap[row*CUBEWIDTH + col];	 
			index++;
		}
	}
	qsort(sample, boxsize, sizeof(float), compare_floats);
	quadIV_med = sample[boxsize/2]; 

	avg = (quadI_med + quadII_med + quadIII_med + quadIV_med)/4;

//	printf("DC estimates:\n");
//	printf("quad I: %.3f, quad II: %.3f, quad III: %.3f, quad IV: %.3f\n", quadI_med, quadII_med, quadIII_med, quadIV_med);	
//	printf("avg of four DC estimates = %.3f\n", avg);

	for(index = 0; index < SIZEOFCUBE; index++)
	{
		cube[index] -= avg;
	} 

	free(sample);
	return 0;
}

int interpret_flags(int argc, char **argv, int flagarg, float *shiftx, float *shifty, int *anchorx, int *anchory)
{
	int i, length;
	char *flagstr = (char *)calloc(STRSIZE, sizeof(char));
	strcpy(flagstr, argv[flagarg]);

	length = strlen(flagstr);

	for(i=0; i<length; i++)
	{
		if(flagstr[i] == 'l')
		{
			printf("Local background subtraction is OFF (on by default).\n");
			localbgsubswitch = 0;
			break;
		}
	}
	for(i=0; i<length; i++)
	{
		if(flagstr[i] == 'p')
		{
			printf("Already Processed switch ON, will treat focal plane images as ready for extraction.\n");
			procd_switch = 1;
			break;
		}
	}
	for(i=0; i<length; i++)
	{
		if(flagstr[i] == 'y')
		{
			printf("Bias/dark-subtraction OFF (on by default).\n");
			darksubswitch = 0;
			break;
		}
	}
	for(i=0; i<length; i++)
	{
		if(flagstr[i] == 'o')
		{
			printf("OVERWRITE switch is ON.\n");
			owswitch = 1;
			break;
		}
	}
	for(i=0; i<length; i++)
	{
		if(flagstr[i] == 'q')
		{
			printf("Dewar flat-fielding is OFF (on by default).\n");
			applydewarflat = 0;
			break;
		}
	}
	for(i=0; i<length; i++)
	{
		if(flagstr[i] == 'w')
		{
			printf("Spectral calibration is ON, will produce cubes with and without crude spectral calibration (off by default).\n");
			applyspeccal = 1;
			break;
		}		
	}
	for(i=0; i<length; i++)
	{
		if(flagstr[i] == 'a')
		{
			printf("Process all image type switch ON.\n");
			doallswitch = 1;
			break;
		}
	}
	for(i=0; i<length; i++)
	{
		if(flagstr[i] == 's')
		{
			if(argc >= 5)
			{
				sscanf(argv[3], "%f", shiftx);	
				sscanf(argv[4], "%f", shifty);
				printf("Spectrograph focal plane image registration override ON. User-specified offset from calibration image:\n\t%0.2f columns, %0.2f rows\n", *shiftx, *shifty);
			}
			else
			{
				*shiftx = 0;
				*shifty = 0;
				printf("Spectrograph focal plane image registration override ON. No offset specified; assuming (0,0).\n\t(Example of command line syntax to specify an override vector: './pipeline ../DATA/goodstuff -os 2.2 -4.6')\n");
			}
			regoverride = 1;
			break;
		}		
	}
	for(i=0; i<length; i++)
	{
		if(flagstr[i] == 'r')
		{
			if(argc >= 5)
			{
				sscanf(argv[3], "%d", anchorx);
				sscanf(argv[4], "%d", anchory);
				printf("Alignment anchor ON. User-specified alignment anchor (cube coordinates):\n\t%d, %d\n", *anchorx, *anchory);
				alignanchor = 1;
			}
			break;
		}		
	}
	for(i=0; i<length; i++)
	{
		if(flagstr[i] == 'v')
		{
			fine_vert_reg_switch = 0;
			printf("Fine vertical registration OFF.\n"); 
		}
	}
	for(i=0; i<length; i++)
	{
		if(flagstr[i] == 'x')
		{
			extract_switch = 0;
			printf("Will not extract cubes from the detector images.\n"); 
		}
	}

	free(flagstr);
	return 0;	
}

int smooth_badlenslets(float *cube, int NReads, char Type, float lookuptable[CUBEWIDTH][CUBEWIDTH][5])
{//replace zero-valued lenslet pixels with an average of neighboring pixels 
	//also, replace anomalously high or low pixels 
	int row, col, w;
	float here, medneighb;
	float surround[8];
	int surround_count = 0;
	int high_count = 0;
	int low_count = 0;
	/*	thresh = 10*sigma where blank sky sigma = 0.11 cts/sec for 10-read exposure	*/
	float thresh;
	float lfr = 2/3.0;
	float hfr = 3.0;

	if(NReads <= 2)
	{
		thresh = 25;
	}
	else
	{
		thresh = 1.1;
	}

	for(w=0; w<NCHAN; w++)
	{
		for(row=1; row<CUBEWIDTH-1; row++)
		{
			for(col=1; col<CUBEWIDTH-1; col++)
			{
				if(lookuptable[col][row][2] > 0)
				{
					//if current lenslet is defined, form array of surrounding ones and check for high pixel
					here = cube[w*CUBEWIDTH*CUBEWIDTH + row*CUBEWIDTH + col];
					surround_count = 0;

					if(lookuptable[col][row+1][2] > 0)
					{//above
						surround[surround_count] = cube[w*CUBEWIDTH*CUBEWIDTH + (row+1)*CUBEWIDTH + col];
						surround_count++;
					}
					if(lookuptable[col][row-1][2] > 0)
					{//below
						surround[surround_count] = cube[w*CUBEWIDTH*CUBEWIDTH + (row-1)*CUBEWIDTH + col];
						surround_count++;
					}
					if(lookuptable[col-1][row][2] > 0)
					{//left
						surround[surround_count] = cube[w*CUBEWIDTH*CUBEWIDTH + row*CUBEWIDTH + col-1];
						surround_count++;
					}
					if(lookuptable[col+1][row][2] > 0)
					{//right
						surround[surround_count] = cube[w*CUBEWIDTH*CUBEWIDTH + row*CUBEWIDTH + col+1];
						surround_count++;
					}
					if(lookuptable[col-1][row+1][2] > 0)
					{//above left
						surround[surround_count] = cube[w*CUBEWIDTH*CUBEWIDTH + (row+1)*CUBEWIDTH + col-1];
						surround_count++;
					}
					if(lookuptable[col+1][row+1][2] > 0)
					{//above right
						surround[surround_count] = cube[w*CUBEWIDTH*CUBEWIDTH + (row+1)*CUBEWIDTH + col+1];
						surround_count++;
					}
					if(lookuptable[col-1][row-1][2] > 0)
					{//below left
						surround[surround_count] = cube[w*CUBEWIDTH*CUBEWIDTH + (row-1)*CUBEWIDTH + col-1];
						surround_count++;
					}
					if(lookuptable[col+1][row-1][2] > 0)
					{//below right
						surround[surround_count] = cube[w*CUBEWIDTH*CUBEWIDTH + (row-1)*CUBEWIDTH + col+1];
						surround_count++;
					}
					
					if(surround_count > 2)
					{	
						qsort(surround, surround_count, sizeof(float), compare_floats);
						if(surround_count % 2 == 0)
						{
							medneighb = (surround[surround_count/2 - 1] + surround[surround_count/2])/2;
						}
						else
						{
							medneighb = surround[surround_count/2]/2;
						}
						/* original condition
						if((here > medneighb + thresh) && ((medneighb > 0 && here > hfr*medneighb)))
						*/
						if(here > 0 && here > (medneighb + thresh) && (here - medneighb) > hfr*medneighb)
						{//replace high pixel
							cube[w*CUBEWIDTH*CUBEWIDTH + row*CUBEWIDTH + col] = medneighb;
							high_count++;
						}
					}
				}
			}
		}
	}
	for(w=0; w<NCHAN; w++)
	{
		for(row=1; row<CUBEWIDTH-1; row++)
		{
			for(col=1; col<CUBEWIDTH-1; col++)
			{
				if(lookuptable[col][row][2] > 0)
				{
					//if current lenslet is defined, form array of surrounding ones and check for low pixel
					here = cube[w*CUBEWIDTH*CUBEWIDTH + row*CUBEWIDTH + col];
					surround_count = 0;

					if(lookuptable[col][row+1][2] > 0)
					{//above
						surround[surround_count] = cube[w*CUBEWIDTH*CUBEWIDTH + (row+1)*CUBEWIDTH + col];
						surround_count++;
					}
					if(lookuptable[col][row-1][2] > 0)
					{//below
						surround[surround_count] = cube[w*CUBEWIDTH*CUBEWIDTH + (row-1)*CUBEWIDTH + col];
						surround_count++;
					}
					if(lookuptable[col-1][row][2] > 0)
					{//left
						surround[surround_count] = cube[w*CUBEWIDTH*CUBEWIDTH + row*CUBEWIDTH + col-1];
						surround_count++;
					}
					if(lookuptable[col+1][row][2] > 0)
					{//right
						surround[surround_count] = cube[w*CUBEWIDTH*CUBEWIDTH + row*CUBEWIDTH + col+1];
						surround_count++;
					}
					if(lookuptable[col-1][row+1][2] > 0)
					{//above left
						surround[surround_count] = cube[w*CUBEWIDTH*CUBEWIDTH + (row+1)*CUBEWIDTH + col-1];
						surround_count++;
					}
					if(lookuptable[col+1][row+1][2] > 0)
					{//above right
						surround[surround_count] = cube[w*CUBEWIDTH*CUBEWIDTH + (row+1)*CUBEWIDTH + col+1];
						surround_count++;
					}
					if(lookuptable[col-1][row-1][2] > 0)
					{//below left
						surround[surround_count] = cube[w*CUBEWIDTH*CUBEWIDTH + (row-1)*CUBEWIDTH + col-1];
						surround_count++;
					}
					if(lookuptable[col+1][row-1][2] > 0)
					{//below right
						surround[surround_count] = cube[w*CUBEWIDTH*CUBEWIDTH + (row-1)*CUBEWIDTH + col+1];
						surround_count++;
					}
					
					if(surround_count > 2)
					{	
						qsort(surround, surround_count, sizeof(float), compare_floats);
						if(surround_count % 2 == 0)
						{
							medneighb = (surround[surround_count/2 - 1] + surround[surround_count/2])/2;
						}
						else
						{
							medneighb = surround[surround_count/2]/2;
						}
						/*original
						if(here < medneighb - thresh && (here < 0 || medneighb > lfr*here))
						*/
						if(here < medneighb - thresh && (here < 0 || (medneighb - here) > lfr*medneighb))
						{//replace low pixel
							cube[w*CUBEWIDTH*CUBEWIDTH + row*CUBEWIDTH + col] = medneighb;
							low_count++;
						}
					}
				}
			}
		}
	}

	printf("Replaced %d high and %d low pixels\n", high_count, low_count);
	return 0;
}

int make_speccal_cube(float *cube, float *speccalcube, float *speccal_array)
{
	int w, i, j;

	for(w = 0; w < NCHAN; w++)
	{
		for(i = 0; i < CUBEWIDTH; i++)
		{
			for(j = 0; j < CUBEWIDTH; j++)
			{
				speccalcube[w*CUBEWIDTH*CUBEWIDTH + j*CUBEWIDTH + i] = cube[w*CUBEWIDTH*CUBEWIDTH + j*CUBEWIDTH + i]/speccal_array[w];   
			}	
		}
	}		

	return 0;
} 

int write_resbgmap(float *resbgframe, FrameTok *Tok, char *auxdir)
{
	char *resbgname, *cmd;
	resbgname = (char *) calloc(STRSIZE, sizeof(char));
	cmd = (char *) calloc(STRSIZE, sizeof(char));
	int i,j;

	long naxes[2];
	naxes[0] = CUBEWIDTH;
	naxes[1] = CUBEWIDTH;

	for(i=0; i<CUBEWIDTH; i++)
	{//set values of rescount image with UNDEFINED values to zero
		for(j=0; j<CUBEWIDTH; j++)
		{
			if(resbgframe[j*CUBEWIDTH + i] == UNDEFINED)
			{
				resbgframe[j*CUBEWIDTH + i] = 0;
			}	
		}
	}

	sprintf(resbgname, "%s/%s-res.fits", auxdir, Tok->FocplShortName);
//	printf("resbgname = %s\n", resbgname);
	if(check_file_exist(resbgname) != 0)
	{
		sprintf(cmd, "rm %s", resbgname);
		system(cmd);
	}
	write_simple_fits_float_data(resbgname, 2, naxes, resbgframe);

	free(resbgname);
	free(cmd);
	return 0;
}

int copy_fits_header(char *oldfitsname, char *newfitsname, FrameTok Tok, char **argv, char *datadir, int flagarg)
{
	fitsfile *fptr, *old;
	int *nullptr = NULL;
	int status = 0;
	int i, nkeys;
	char *record;
	record = (char *) calloc(STRSIZE, sizeof(char));

	//add header from original focal plane image	
	fits_open_file(&fptr, newfitsname, READWRITE, &status);
	if(status)
	{
		fits_report_error(stdout, status);
		return(status);
	}
	fits_open_file(&old, oldfitsname, READONLY, &status);
	if(status)
	{
		fits_report_error(stdout, status);
		return(status);
	}
	fits_get_hdrspace(old, &nkeys, nullptr, &status); 
	if(status)
	{
		fits_report_error(stdout, status);
		return(status);
	}
	for(i=0; i<=nkeys; i++)
	{
		fits_read_record(old, i, record, &status);	
		fits_write_record(fptr, record, &status);
		if(status)
		{
			fits_report_error(stdout, status);
			return(status);
		}
	}
	if(status)
	{
		fits_report_error(stdout, status);
		return(status);
	}
	//add detector image processing keywords
	fits_write_key(fptr, TSTRING, "EXECFILE", argv[0], "executable file", &status);
	fits_write_key(fptr, TSTRING, "DATADIR", datadir, "raw data directory", &status);
	fits_write_key(fptr, TSTRING, "EXECARGS", argv[flagarg], "command line arguments", &status);
	fits_write_key(fptr, TINT, "DARKSUB", &Tok.DarkSubFlag, "= 1 if did dark-subtraction", &status);
	fits_write_key(fptr, TFLOAT, "XSHIFT", &Tok.xshift, "Computed column offset from laser reference", &status);
	fits_write_key(fptr, TFLOAT, "YSHIFT", &Tok.yshift, "Computed row offset from laser reference", &status);
	fits_close_file(fptr, &status);
	fits_close_file(old, &status);

	free(record);

	return 0;
}

int write_cube_header(FrameTok Tok, char *cubename, char *calib_library_str, char *speccal_tok, char **argv, char *datadir, int flagarg)
{
	fitsfile *fptr, *old;
	int *nullptr = NULL;
	int status = 0;
	int i, nkeys, nslices;
	char *record, *keyval;
	record = (char *) calloc(STRSIZE, sizeof(char));
	keyval = (char *) calloc(STRSIZE, sizeof(char));

	//add header from focal plan image
	fits_open_file(&fptr, cubename, READWRITE, &status);
	if(status)
	{
		fits_report_error(stdout, status);
		return(status);
	}
	fits_open_file(&old, Tok.FocplInitName, READONLY, &status);
	if(status)
	{
		fits_report_error(stdout, status);
		return(status);
	}
	fits_get_hdrspace(old, &nkeys, nullptr, &status); 
	if(status)
	{
		fits_report_error(stdout, status);
		return(status);
	}
	for(i=6; i<=nkeys; i++)
	{
		fits_read_record(old, i, record, &status);	
		fits_write_record(fptr, record, &status);
		if(status)
		{
			fits_report_error(stdout, status);
			return(status);
		}
	}
	if(status)
	{
		fits_report_error(stdout, status);
		return(status);
	}
	//add detector image processing keywords
	fits_write_key(fptr, TSTRING, "EXECFILE", argv[0], "executable file", &status);
	fits_write_key(fptr, TSTRING, "DATADIR", datadir, "raw data directory", &status);
	fits_write_key(fptr, TSTRING, "EXECARGS", argv[flagarg], "command line arguments", &status);
	fits_write_key(fptr, TINT, "DARKSUB", &Tok.DarkSubFlag, "= 1 if did dark-subtraction", &status);
	fits_write_key(fptr, TFLOAT, "XSHIFT", &Tok.xshift, "Computed column offset from laser reference", &status);
	fits_write_key(fptr, TFLOAT, "YSHIFT", &Tok.yshift, "Computed row offset from laser reference", &status);
	if(status)
	{
		fits_report_error(stdout, status);
		return(status);
	}
	fits_write_key(fptr, TINT, "XANCHOR", &Tok.AnchorX, "Cube anchor x position for offset calc (approx)", &status);
	fits_write_key(fptr, TINT, "YANCHOR", &Tok.AnchorY, "Cube anchor y position for offset calc (approx)", &status);
	
	//add cube extraction keywords
	fits_write_key(fptr, TSTRING, "CUBEMODE", "Weighted sum", "Cube extraction mode", &status);
	fits_write_key(fptr, TSTRING, "FOCPLSOL", calib_library_str, "Epoch of focal plane solution", &status);
	nslices = NCHAN;
	fits_write_key(fptr, TINT, "NCHAN", &nslices, "Number of wavelength channels", &status);
	if(speccal_tok != NULL)
	{
		fits_write_key(fptr, TSTRING, "SPECCAL", speccal_tok, "spectral calibration token", &status);
	}
	else
	{
		fits_write_key(fptr, TSTRING, "SPECCAL", "NONE", "spectral calibration token", &status);
	}
	if(status)
	{
		fits_report_error(stdout, status);
		return(status);
	}
	fits_close_file(fptr, &status);
	if(status)
	{
		fits_report_error(stdout, status);
		return(status);
	}
	fits_close_file(old, &status);
	if(status)
	{
		fits_report_error(stdout, status);
		return(status);
	}

	free(record);
	free(keyval);

	return 0;
}

int write_cube(float *cube, char *cubename)
{
	char *cmd;
	int owflag;
	cmd = (char *) calloc(BIGSTRSIZE, sizeof(char));
	long naxes_cube[3];
	naxes_cube[0] = CUBEWIDTH;
	naxes_cube[1] = CUBEWIDTH;
	naxes_cube[2] = NCHAN;

	owflag = 0;
	if(check_file_exist(cubename) == 0 || owswitch)
	{
		if(check_file_exist(cubename) != 0)
		{
			sprintf(cmd, "rm %s", cubename);
			system(cmd);
			owflag = 1;
		}
		write_simple_fits_float_data(cubename, 3, naxes_cube, cube);
		if(owflag)
		{
			printf("Overwriting existing data cube %s\n", cubename);
		}
		else
		{
			printf("Writing data cube %s\n", cubename);
		}
		free(cmd);
		return 0;
	}
	else
	{ //not writing cube
		free(cmd);
		return -1;
	}
}

void load_focplsol(float ***focplsol, char *focplsol_fname)
{
	int lensx, lensy;
	char *line = (char *)calloc(STRSIZE, sizeof(char));
	float X, Y, height, tilt, amp, eucliddist, deltaX, deltaY, DClevel;
	
	FILE *fptr = fopen(focplsol_fname, "r");
	int linecount = 0;

	while(fgets(line, STRSIZE, fptr) != NULL)
	{
		sscanf(line, "%d, %d, %f, %f, %f, %f, %f, %f, %f, %f, %f\n", &lensx, &lensy, &X, &Y, &deltaX, &deltaY, &height, &tilt, &amp, &DClevel, &eucliddist);
		focplsol[lensx][lensy][0] = X;
		focplsol[lensx][lensy][1] = Y;
		focplsol[lensx][lensy][2] = height;
		focplsol[lensx][lensy][3] = tilt;
		focplsol[lensx][lensy][4] = amp;
		linecount++;
	}
	fclose(fptr);
	printf("Loaded %d lines from focal plane solution %s\n", linecount, focplsol_fname);
	free(line);
}

void load_psf_params(double *psf_params, char *psf_model_fname)
{
	fitsfile *psf_fptr;
	char *comment = NULL;
	int status = 0;
	
	fits_open_file(&psf_fptr, psf_model_fname, READONLY, &status);
	if(status)
	{
		fits_report_error(stdout, status);
	}
	fits_read_key(psf_fptr, TDOUBLE, "A", &psf_params[0], comment, &status);
	if(status)
	{
		fits_report_error(stdout, status);
	}
	fits_read_key(psf_fptr, TDOUBLE, "B", &psf_params[1], comment, &status);
	fits_read_key(psf_fptr, TDOUBLE, "SIGX_A", &psf_params[2], comment, &status);
	if(status)
	{
		fits_report_error(stdout, status);
	}
	fits_read_key(psf_fptr, TDOUBLE, "SIGY_A+", &psf_params[3], comment, &status);
	fits_read_key(psf_fptr, TDOUBLE, "SIGY_A-", &psf_params[4], comment, &status);
	fits_read_key(psf_fptr, TDOUBLE, "SIGX_B", &psf_params[5], comment, &status);
	fits_read_key(psf_fptr, TDOUBLE, "SIGY_B+", &psf_params[6], comment, &status);
	fits_read_key(psf_fptr, TDOUBLE, "SIGY_B-", &psf_params[7], comment, &status);
	fits_read_key(psf_fptr, TDOUBLE, "THETA", &psf_params[8], comment, &status);
	fits_close_file(psf_fptr, &status);
}

void count_badpix(int *Nbadpix, char *badpix_list_fname)
{
	char *line = (char *)calloc(STRSIZE, sizeof(char));

	FILE *fptr = fopen(badpix_list_fname, "r");
	while(fgets(line, STRSIZE, fptr) != NULL)
	{
		(*Nbadpix) += 1;	
	}
	fclose(fptr);

	free(line);
}

void load_badpix_list(int **badpix_list, int Nbadpix, char *badpix_list_fname)
{
	int i;
	char *line = (char *)calloc(STRSIZE, sizeof(char));

	FILE *fptr = fopen(badpix_list_fname, "r");
	i = 0;
	while(fgets(line, STRSIZE, fptr) != NULL)
	{
		sscanf(line, "%d, %d\n", &badpix_list[i][0], &badpix_list[i][1]);	
		i++;
	}
	fclose(fptr);
	//printf("last bad pixel: %d, %d\n", badpix_list[Nbadpix-1][0], badpix_list[Nbadpix-1][1]);
	free(line);
}

void make_psf(double *psf, double *P, int bindens)
{
	int c, r, i, j, m, n;
	int width = 9*bindens;
	int xcent = 4;
	int ycent = 4;

	double *Xfine = (double *)calloc(width*width, sizeof(double));
	double *Yfine = (double *)calloc(width*width, sizeof(double));

	double xpix, ypix;

	for(c = 0; c < width; c++)
	{
		for(r = 0; r < width; r++)
		{
			Xfine[r*width + c] = (c*cos(-P[8]*PI/180) - r*sin(-P[8]*PI/180))/bindens;
			Yfine[r*width + c] = (c*sin(-P[8]*PI/180) + r*cos(-P[8]*PI/180))/bindens;
		}
	}

	for(i = 0; i < 9; i++)
	{
		for(j = 0; j < 9; j++)
		{
			for(m = -bindens/2; m <= bindens/2; m++)
			{
				for(n = -bindens/2; n <= bindens/2; n++)
				{
					xpix = Xfine[(j*bindens + bindens/2 + n)*width + i*bindens + bindens/2 + m] - Xfine[(ycent*bindens + bindens/2)*width + xcent*bindens + bindens/2];	
					ypix = Yfine[(j*bindens + bindens/2 + n)*width + i*bindens + bindens/2 + m] - Yfine[(ycent*bindens + bindens/2)*width + xcent*bindens + bindens/2];	


					if(ypix > 0)
					{
						psf[(j*bindens + bindens/2 + n)*width + i*bindens + bindens/2 + m] = P[0]*exp(-(gsl_pow_2(xpix/P[2]) + gsl_pow_2(ypix/P[3]))/2) + P[1]*exp(-(gsl_pow_2(xpix/P[5]) + gsl_pow_2(ypix/P[6]))/2);
					}
					else
					{
						psf[(j*bindens + bindens/2 + n)*width + i*bindens + bindens/2 + m] = P[0]*exp(-(gsl_pow_2(xpix/P[2]) + gsl_pow_2(ypix/P[4]))/2) + P[1]*exp(-(gsl_pow_2(xpix/P[5]) + gsl_pow_2(ypix/P[7]))/2);
					}
/*
					if(i == 6 && m == 0 && j == 3 && n == 0)
					{
						printf("xpix, ypix = %0.2f, %0.2f\n", xpix, ypix);
						printf("* %0.2f\n", gsl_pow_2(xpix/P[2]));
						printf("** %0.2f\n", P[0]*exp(-1/2*(gsl_pow_2(xpix/P[2]) + gsl_pow_2(ypix/P[4]))));
						printf("*** %0.2f\n", gsl_pow_2(xpix/P[2]));
						printf("**** %f\n", (float)exp((double)(-gsl_pow_2(xpix/P[2]))/2.0));
						printf("%f\n", psf[(j*bindens + bindens/2 + n)*width + i*bindens + bindens/2 + m]);
					}
*/
				}
			}
		} 
	}

	free(Xfine);
	free(Yfine);
}

int check_cutout_clearance(int Xwater, int Ywater)
{
	/*	assumes simulated region with 3 pixel movement in each direction */
	if(Xwater < BORDER + 3 + SIMWIDTH/2 || Xwater > FOCPLWIDTH - BORDER - 3 - SIMWIDTH/2 || Ywater < BORDER + 3 + REFROW  || Ywater > FOCPLWIDTH - BORDER - 3 - (SIMHEIGHT - (REFROW + 1)))
	{
		return AT_EDGE;
	}
	else
	{
		return CLEAR;
	}	
}

int get_fine_vert_offset(float *imagebuf, float ***lookuptable, float *heightmap, float *tiltmap, FrameTok Token, double *pix_resp, double *J_psf, double *H_psf, int psfwidth, int *crudeYoffset, int *fineYoffset)
{
	int samplebox_width = 11;
	//int samplebox_width = 3;
	int Nfit = samplebox_width*samplebox_width;
	int focpl_x, focpl_y, lensx, lensy, i, j;
	int *lensxarray = (int *)calloc(Nfit, sizeof(int));
	int *lensyarray = (int *)calloc(Nfit, sizeof(int));
	double *Yoffarray = (double *)calloc(Nfit, sizeof(double));
	double med_Yoff, stddev_Yoff;
	int tries;

	/* debug variables */
	/*
	int p;
	char *cmd = (char *)calloc(STRSIZE, sizeof(char));
	char *name = (char *)calloc(STRSIZE, sizeof(char));
	long axis_length[2];
	*/

	specsim_datapack specpack;
	specpack.simwidth = SIMWIDTH;
	specpack.simheight = SIMHEIGHT;
	specpack.fitwidth = 3;
	specpack.fitheight = 11;
	specpack.psfwidth = PSFWIDTH;
	int hires_psfwidth = specpack.psfwidth*VBINDENS; 
	int hires_paddedheight = (specpack.simheight + specpack.psfwidth)*VBINDENS;
	int hires_paddedwidth = (specpack.simwidth + specpack.psfwidth)*VBINDENS;
	double *deviates = (double *)calloc(specpack.fitwidth*specpack.fitheight, sizeof(double));	
	double **derivs = (double **)calloc(specpack.fitwidth*specpack.fitheight, sizeof(double));
	specpack.focpl_cutout = (double *)calloc(specpack.fitwidth*specpack.fitheight, sizeof(double));
	specpack.pix_resp = pix_resp;
	specpack.J_psf = (double *)calloc(hires_paddedwidth*hires_paddedheight, sizeof(double));
	specpack.H_psf = (double *)calloc(hires_paddedwidth*hires_paddedheight, sizeof(double));
	specpack.FT_J_psf = (fftw_complex *)fftw_malloc(hires_paddedwidth*hires_paddedheight*sizeof(fftw_complex));
	specpack.FT_H_psf = (fftw_complex *)fftw_malloc(hires_paddedwidth*hires_paddedheight*sizeof(fftw_complex));
	specpack.J_spec_trace = (double *)calloc(hires_paddedwidth*hires_paddedheight, sizeof(double));
	specpack.H_spec_trace = (double *)calloc(hires_paddedwidth*hires_paddedheight, sizeof(double));
	specpack.FT_J_spec_trace = (fftw_complex *)fftw_malloc(hires_paddedwidth*hires_paddedheight*sizeof(fftw_complex));
	specpack.FT_H_spec_trace = (fftw_complex *)fftw_malloc(hires_paddedwidth*hires_paddedheight*sizeof(fftw_complex));
	specpack.FT_J_spec_image = (fftw_complex *)fftw_malloc(hires_paddedwidth*hires_paddedheight*sizeof(fftw_complex));
	specpack.FT_H_spec_image = (fftw_complex *)fftw_malloc(hires_paddedwidth*hires_paddedheight*sizeof(fftw_complex));
	specpack.J_spec_image = (double *)calloc(hires_paddedwidth*hires_paddedheight, sizeof(double));
	specpack.H_spec_image = (double *)calloc(hires_paddedwidth*hires_paddedheight, sizeof(double));
	specpack.binned_spec_image = (double *)calloc(specpack.fitwidth*specpack.fitheight, sizeof(double));

	specpack.shortleg = SHORTLEG;		
	specpack.longleg = LONGLEG;
	
	printf("Determining fine vertical offset...\n");

	copy_array_2d_double(J_psf, specpack.J_psf, hires_psfwidth, hires_psfwidth, hires_paddedwidth, hires_paddedheight, 0, 0, hires_paddedwidth/2 - hires_psfwidth/2, hires_paddedheight/2 - hires_psfwidth/2, hires_psfwidth, hires_psfwidth);
	copy_array_2d_double(H_psf, specpack.H_psf, hires_psfwidth, hires_psfwidth, hires_paddedwidth, hires_paddedheight, 0, 0, hires_paddedwidth/2 - hires_psfwidth/2, hires_paddedheight/2 - hires_psfwidth/2, hires_psfwidth, hires_psfwidth);

	size_t *sorted_indices = (size_t *)calloc(FOCPLWIDTH*FOCPLWIDTH, sizeof(long));
	gsl_sort_float_index(sorted_indices, imagebuf, 1, FOCPLWIDTH*FOCPLWIDTH);
	long ninetyninth_cut_index = lroundf(FOCPLWIDTH*FOCPLWIDTH*0.99) - 1;
	float ninetyninth_pct = imagebuf[sorted_indices[ninetyninth_cut_index]];
	/*printf("99th percentile pixel value = %.2f\n", ninetyninth_pct);*/
	float Xwater, Ywater;

	int Nparams = 7; /*color, height (fixed), tilt (fixed), amplitude, deltaX (fixed), deltaY, DC offset*/
	int status;
	int Npoints = specpack.fitwidth*specpack.fitheight;
	double max_vertoff = 1.5;
	double *P = (double *)calloc(Nparams, sizeof(double));
	mp_par *P_constraints = (mp_par *)calloc(Nparams, sizeof(mp_par));
	memset(P_constraints, 0, sizeof(P_constraints));
	mp_result result;
	/*	structure of parameter (P) vector
	 *	0		color (J/H ratio)
     *	1		height
     *	2		tilt
     *	3		amplitude
     *	4		deltaX
	 *	5		deltaY
	 *	6		DC offset
     */
	/*color	*/
	P_constraints[0].fixed = 0;	
	P_constraints[0].limited[0] = 1;
	P_constraints[0].limited[1] = 0;
	P_constraints[0].limits[0] = 0;
	P_constraints[0].relstep = 0.1;
	/*height, tilt: */
	P_constraints[1].fixed = 1;
	P_constraints[2].fixed = 1;
	/*amplitude:	*/
	P_constraints[3].limited[0] = 1;
	P_constraints[3].limited[1] = 0;
	P_constraints[3].limits[0] = 0;
	P_constraints[3].relstep = 0.1;
	/*deltaX:		*/
	P_constraints[4].fixed = 1;
	/*deltaY:		*/
	P_constraints[5].fixed = 0;
	P_constraints[5].limited[0] = 1;
	P_constraints[5].limited[1] = 1;
	P_constraints[5].limits[0] = -max_vertoff;
	P_constraints[5].limits[1] = max_vertoff;
	P_constraints[5].step = 0.2;
	/*DC offset:	*/
	P_constraints[6].limited[0] = 1;
	P_constraints[6].limited[1] = 0;
	P_constraints[6].limits[0] = 0;
	P_constraints[6].relstep = 0.1;

	/*
	int patch_width = FOCPLWIDTH/DISTMAPWIDTH*2;
	*/
	int patch_width = 600;
	int fit_gap = patch_width/(samplebox_width + 1);
	/**/
	int focpl_xstart = Token.PeakCol - patch_width/2;
	int focpl_ystart = Token.PeakRow - patch_width/2;
	/**/
	/*
	int focpl_xstart = 1010 - patch_width/2;
	int focpl_ystart = 1075 - patch_width/2;
	*/
	
	*crudeYoffset = Token.crudeYoffset;
	for(tries = 0; tries < 4; tries++)
	{
		int fit_attempt_count = 0, fit_success_count = 0;
		for(i = 0; i < samplebox_width; i++)
		{
			for(j = 0; j < samplebox_width; j++)
			{
				focpl_x = focpl_xstart + fit_gap*(1 + i);	
				focpl_y = focpl_ystart + fit_gap*(1 + j);	
				get_cube_coords(focpl_x, focpl_y, &lensx, &lensy);
				/* debug begin */
				/*
				lensx = 146;
				lensy = 123;
				*/
				/* debug end */

				Xwater = lookuptable[lensx][lensy][0] + Token.crudeXoffset + Token.fineXoffset/(float)BINDENS;
				Ywater = lookuptable[lensx][lensy][1] + *crudeYoffset;
							
				/*debug begin*/
				/*
				printf("%d, %d: %d, %d <-> %d, %d\n", i, j, focpl_x, focpl_y, lensx, lensy);
				printf("Xwater, Ywater = %0.1f, %0.1f => %ld, %ld\n", Xwater, Ywater, lroundf(Xwater), lroundf(Ywater));
				*/
				/*debug end*/
		
				if(lookuptable[lensx][lensy][2] > 0 && check_cutout_clearance(Xwater, Ywater) == CLEAR)
				{
					P[0] = 1.0; /* color */
					P[1] = heightmap[lensy*CUBEWIDTH + lensx]; /* height */
					P[2] = tiltmap[lensy*CUBEWIDTH + lensx]; /* tilt */
					P[3] = ninetyninth_pct/10; /* amplitude */
					P[4] = Xwater - lroundf(Xwater); /* delta X */
					P[5] = Ywater - lroundf(Ywater); /*delta Y */
					P[6] = 0.1; /* DC */
					copy_array_2d_float_to_double(imagebuf, specpack.focpl_cutout, FOCPLWIDTH, FOCPLWIDTH, specpack.fitwidth, specpack.fitheight, (int)lroundf(Xwater) - specpack.fitwidth/2, (int)lroundf(Ywater) - 7, 0, 0, specpack.fitwidth, specpack.fitheight);
					memset(&result, 0, sizeof(result));
					/**/
					status = mpfit(specsim, Npoints, Nparams, P, P_constraints, 0, (void *)&specpack, &result);
					/**/
					/*
					specsim(Npoints, Nparams, P, deviates, derivs, (void *)&specpack);
					*/
					/*debug begin*/
					/*
					printf("\tP = ");
					for(p = 0; p < Nparams; p++)
					{
						printf("%.2f, ", P[p]);
					}
					printf("\n");
					*/
					/*debug end*/
	
					//if(abs(P[5]) < max_vertoff && P[0] < 2.0 && P[0] > 0.5)
					if(P[0] < 2.0 && P[0] > 0.5)
					{ /*exclude extreme fit cases*/
						Yoffarray[fit_success_count] = (P[5] + lroundf(Ywater)) - Ywater;
						fit_success_count++;
					}
					fit_attempt_count++;

					/*debug begin*/
					/*
					axis_length[0] = specpack.fitwidth;
					axis_length[1] = specpack.fitheight;
					sprintf(name, "final_model_%d_%d.fits", i, j);
					if(check_file_exist(name))
					{
						sprintf(cmd, "rm %s", name);
						system(cmd);
					}
					write_simple_fits_double_data(name, 2, axis_length, specpack.binned_spec_image);
					sprintf(name, "data_cutout_%d_%d.fits", i, j);
					if(check_file_exist(name))
					{
						sprintf(cmd, "rm %s", name);
						system(cmd);
					}
					write_simple_fits_double_data(name, 2, axis_length, specpack.focpl_cutout);
					*/
					/*debug end*/
				}
			}
		}

		gsl_sort(Yoffarray, 1, fit_success_count);
		med_Yoff = gsl_stats_median_from_sorted_data(Yoffarray, 1, fit_success_count);
		stddev_Yoff = gsl_stats_sd_m(Yoffarray, 1, Nfit, med_Yoff);
		*crudeYoffset += (int)lround(med_Yoff);
		*fineYoffset = (int)lround((med_Yoff - (int)lround(med_Yoff))*BINDENS);

		printf("Vertical registration: %d successful fits. Median, std. dev. of Y offset array = %0.2f, %0.2f\n", fit_success_count, med_Yoff, stddev_Yoff);
	
		if(fabs(med_Yoff) < 0.5)
		{
			break;
		}
		else
		{
			printf("Repeating offset fit with crude Y offset set to %d.\n", *crudeYoffset);
		}
	}

	printf("Final vertical offsets (crude, fine): %d, %d\n", *crudeYoffset, *fineYoffset);
	
	free(P);
	free(P_constraints);
	free(sorted_indices);
	free(specpack.focpl_cutout);
	free(specpack.J_psf);
	free(specpack.H_psf);
	free(specpack.J_spec_trace);
	free(specpack.H_spec_trace);
	fftw_free(specpack.FT_J_psf);
	fftw_free(specpack.FT_H_psf);
	fftw_free(specpack.FT_J_spec_trace);
	fftw_free(specpack.FT_H_spec_trace);
	fftw_free(specpack.FT_J_spec_image);
	fftw_free(specpack.FT_H_spec_image);
	free(specpack.J_spec_image);
	free(specpack.H_spec_image);
	free(specpack.binned_spec_image);
	free(deviates);
	free(derivs);
	free(Yoffarray);
	free(lensxarray);
	free(lensyarray);

	return 0;
}

int specsim(int Npoints, int Nparams, double *P, double *deviates, double **derivs, void *private)
{
	specsim_datapack *datapack = (specsim_datapack *)private;
	int x, y, i, j, m, n, w;
	int hires_paddedheight = (datapack->simheight + datapack->psfwidth)*VBINDENS;
	int hires_paddedwidth = (datapack->simwidth + datapack->psfwidth)*VBINDENS;

	memset(datapack->binned_spec_image, 0, Npoints*sizeof(double));
	memset(datapack->J_spec_trace, 0, hires_paddedwidth*hires_paddedheight*sizeof(double));
	memset(datapack->H_spec_trace, 0, hires_paddedwidth*hires_paddedheight*sizeof(double));

	/*	structure of parameter (P) vector
	 *	0		color (J/H ratio)
     *	1		height
     *	2		tilt
     *	3		amplitude
     *	4		deltaX
	 *	5		deltaY
	 *	6		DC offset
     */
	
	double color = P[0];
	double height = P[1];
	double tilt = P[2];
	double amp = P[3];
	double Xoffset = P[4];
	double Yoffset = P[5];
	double DC = P[6];

	double watershape[21] = {0.95, 0.89, 0.84, 0.68, 0.70, 0.38, 0.04, 0.06, 0.12, 0.15, 0.18, 0.20, 0.34, 0.39, 0.46, 0.63, 0.69, 0.48, 0.62, 0.83, 0.91};
	int waterindex;

	double xdist, xdist_lowerleft, xdist_upperright;
	double vertspecscale = (1760 - 1100)/(height*VBINDENS); 
	int w_jedge = (1280 - LAMBDA_MIN)/30;
	int w_hedge = (1520 - LAMBDA_MIN)/30;

	double refx = (datapack->simwidth/2)*VBINDENS + VBINDENS/2 + Xoffset*VBINDENS;
	double refy = REFROW*VBINDENS + VBINDENS/2 + Yoffset*VBINDENS;
	double bluey = refy + (LAMBDA_REF - LAMBDA_MIN)/vertspecscale;

	double refx_lowerleft = refx - datapack->shortleg*VBINDENS;
	double refy_lowerleft = refy - datapack->longleg*VBINDENS;
	double bluey_lowerleft = bluey - datapack->longleg*VBINDENS;

	double refx_upperright = refx + datapack->shortleg*VBINDENS;
	double refy_upperright = refy + datapack->longleg*VBINDENS;
	double bluey_upperright = bluey + datapack->longleg*VBINDENS; 

	for(x = 0; x < datapack->simwidth*VBINDENS; x++)
	{
		for(y = 0; y < datapack->simheight*VBINDENS; y++)
		{
			xdist = ((y - refy)*tilt + refx) - x;
			xdist_lowerleft = ((y - refy_lowerleft)*tilt + refx_lowerleft) - x;
			xdist_upperright = ((y - refy_upperright)*tilt + refx_upperright) - x;
			if(xdist > -0.5 && xdist <= 0.5)
			{
				w = (int)lround((bluey - y)*vertspecscale/30);
				if(w >= 0 && w <= w_jedge)
				{
					if(w == 0)
					{
						datapack->J_spec_trace[y*hires_paddedwidth + x] = 0.2*color*amp;
					}
					else if(w == 1)
					{
						datapack->J_spec_trace[y*hires_paddedwidth + x] = 0.5*color*amp;
					}
					else
					{
						datapack->J_spec_trace[y*hires_paddedwidth + x] = color*amp;
					}
				}
				else if(w >= w_hedge && w < NCHAN)
				{
					if(w == NCHAN-1)
					{
						datapack->H_spec_trace[y*hires_paddedwidth + x] = 0.2*amp;
					}
					else if(w == NCHAN-2)
					{
						datapack->H_spec_trace[y*hires_paddedwidth + x] = 0.5*amp;
					}
					else
					{
						datapack->H_spec_trace[y*hires_paddedwidth + x] = amp;
					}
				}
				else if(w > w_jedge && w < w_hedge)
				{
					waterindex = (int)lround(((bluey - y)*vertspecscale + LAMBDA_MIN - 1300)/10);
					if(waterindex >= 0 && waterindex <= 9)
					{
						datapack->J_spec_trace[y*hires_paddedwidth + x] = color*amp*watershape[waterindex];
					}
					else if(waterindex > 9 && waterindex <= 20)
					{
						datapack->H_spec_trace[y*hires_paddedwidth + x] = amp*watershape[waterindex];
					}
				}
			}
			if(xdist_lowerleft > -0.5 && xdist_lowerleft <= 0.5)
			{
				w = (int)lround((bluey_lowerleft - y)*vertspecscale/30);
				if(w >= 0 && w <= w_jedge)
				{
					if(w == 0)
					{
						datapack->J_spec_trace[y*hires_paddedwidth + x] = 0.2*color*amp;
					}
					else if(w == 1)
					{
						datapack->J_spec_trace[y*hires_paddedwidth + x] = 0.5*color*amp;
					}
					else
					{
						datapack->J_spec_trace[y*hires_paddedwidth + x] = color*amp;
					}
				}
				else if(w >= w_hedge && w < NCHAN)
				{
					if(w == NCHAN-1)
					{
						datapack->H_spec_trace[y*hires_paddedwidth + x] = 0.2*amp;
					}
					else if(w == NCHAN-2)
					{
						datapack->H_spec_trace[y*hires_paddedwidth + x] = 0.5*amp;
					}
					else
					{
						datapack->H_spec_trace[y*hires_paddedwidth + x] = amp;
					}
				}
				else if(w > w_jedge && w < w_hedge)
				{
					waterindex = (int)lround(((bluey_lowerleft - y)*vertspecscale + LAMBDA_MIN - 1300)/10);
					if(waterindex >= 0 && waterindex <= 9)
					{
						datapack->J_spec_trace[y*hires_paddedwidth + x] = color*amp*watershape[waterindex];
					}
					else if(waterindex > 9 && waterindex <= 20)
					{
						datapack->H_spec_trace[y*hires_paddedwidth + x] = amp*watershape[waterindex];
					}
				}
			}
			if(xdist_upperright > -0.5 && xdist_upperright <= 0.5)
			{
				w = (int)lround((bluey_upperright - y)*vertspecscale/30);
				if(w >= 0 && w <= w_jedge)
				{
					if(w == 0)
					{
						datapack->J_spec_trace[y*hires_paddedwidth + x] = 0.2*color*amp;
					}
					else if(w == 1)
					{
						datapack->J_spec_trace[y*hires_paddedwidth + x] = 0.5*color*amp;
					}
					else
					{
						datapack->J_spec_trace[y*hires_paddedwidth + x] = color*amp;
					}
				}
				else if(w >= w_hedge && w < NCHAN)
				{
					if(w == NCHAN-1)
					{
						datapack->H_spec_trace[y*hires_paddedwidth + x] = 0.2*amp;
					}
					else if(w == NCHAN-2)
					{
						datapack->H_spec_trace[y*hires_paddedwidth + x] = 0.5*amp;
					}
					else
					{
						datapack->H_spec_trace[y*hires_paddedwidth + x] = amp;
					}
				}
				else if(w > w_jedge && w < w_hedge)
				{
					waterindex = (int)lround(((bluey_upperright - y)*vertspecscale + LAMBDA_MIN - 1300)/10);
					if(waterindex >= 0 && waterindex <= 9)
					{
						datapack->J_spec_trace[y*hires_paddedwidth + x] = color*amp*watershape[waterindex];
					}
					else if(waterindex > 9 && waterindex <= 20)
					{
						datapack->H_spec_trace[y*hires_paddedwidth + x] = amp*watershape[waterindex];
					}
				}
			}
		}
	}

	/* debug */
	/*
	printf("* %0.3f\n", datapack->H_spec_image[10*hires_paddedwidth + 4*VBINDENS]);
	*/

	FT_convolve_2dreal(datapack->J_spec_trace, datapack->FT_J_spec_trace, datapack->J_psf, datapack->FT_J_psf, datapack->J_spec_image, datapack->FT_J_spec_image, hires_paddedwidth, hires_paddedheight);
	FT_convolve_2dreal(datapack->H_spec_trace, datapack->FT_H_spec_trace, datapack->H_psf, datapack->FT_H_psf, datapack->H_spec_image, datapack->FT_H_spec_image, hires_paddedwidth, hires_paddedheight);

	/* debug */
	/*
	long nax[2];
	nax[0] = hires_paddedwidth;
	nax[1] = hires_paddedheight;
	write_simple_fits_double_data("J_spec_trace.fits", 2, nax, datapack->J_spec_trace);
	write_simple_fits_double_data("J_psf.fits", 2, nax, datapack->J_psf);
	write_simple_fits_double_data("J_spec_image.fits", 2, nax, datapack->J_spec_image);
	*/
	/* end debug */

	int widthmargin = (datapack->simwidth - datapack->fitwidth)/2;
	//int ybegin = (datapack->simheight - datapack->fitheight)/2;
	int ybegin = REFROW - 7;
	int yend = ybegin + datapack->fitheight;
	for(i = widthmargin;  i < datapack->simwidth - widthmargin; i++)
	{
		for(j = ybegin; j < yend; j++)
		{
			for(m = 0; m < VBINDENS; m++)
			{
				for(n = 0; n < VBINDENS; n++)
				{
					datapack->binned_spec_image[(j - ybegin)*datapack->fitwidth + i - widthmargin] += (datapack->J_spec_image[(j*VBINDENS + n)*hires_paddedwidth + i*VBINDENS + m] + datapack->H_spec_image[(j*VBINDENS + n)*hires_paddedwidth + i*VBINDENS + m])*datapack->pix_resp[n*VBINDENS + m]/(VBINDENS*VBINDENS);
				}
			}
		}
	}
	/* debug */
	/*
	printf("** %0.2f\n", datapack->binned_spec_image[10*fitwidth + 2]);
	*/
	/* end debug */

	for(i = 0; i < datapack->fitwidth; i++)
	{
		for(j = 0; j < datapack->fitheight; j++)
		{
			datapack->binned_spec_image[j*datapack->fitwidth + i] += DC;
			deviates[j*datapack->fitwidth + i] = datapack->focpl_cutout[j*datapack->fitwidth + i] - datapack->binned_spec_image[j*datapack->fitwidth + i];
		}
	}
	/* weights */
	/*
	for(j = 2; j <= 8; j++)
	{
		deviates[j*3 + 1] *= 5;
	}
	*/

	return 0;
}

void FT_convolve_2dreal(double *A, fftw_complex *FT_A, double *B, fftw_complex *FT_B, double *C, fftw_complex *FT_C, int width, int height)
{
	int i, j;
	float temp;
	int Npoints = width*height;

	fftw_plan A_forward, B_forward, C_backward;
	A_forward = fftw_plan_dft_r2c_2d(height, width, A, FT_A, FFTW_ESTIMATE); 
	B_forward = fftw_plan_dft_r2c_2d(height, width, B, FT_B, FFTW_ESTIMATE); 
	C_backward = fftw_plan_dft_c2r_2d(height, width, FT_C, C, FFTW_ESTIMATE); 

	fftw_execute(A_forward);	
	fftw_execute(B_forward);	

	for(i = 0; i < width; i++)
	{
		for(j = 0; j < height; j++)
		{
			FT_C[j*width + i] = FT_A[j*width + i]*FT_B[j*width + i]; 	
		}
	}

	fftw_execute(C_backward);
	
	fftw_destroy_plan(A_forward);
	fftw_destroy_plan(B_forward);
	fftw_destroy_plan(C_backward);

	if(width % 2 != 0 || height % 2 != 0)
	{
		printf("FT_convolve_2dreal: WARNING, one of the array dimensions (%d x %d) is odd\n", width, height);
	}

	/* diagoonally swap quadrants (equivalent to shift by width/2 in x and height/2 in y) and rescale */
	for(i = 0; i < width/2; i++)
	{
		for(j = 0; j < height/2; j++)
		{/* swap lower left with upper right */
			temp = C[(j + height/2)*width + i + width/2]; 
			C[(j + height/2)*width + i + width/2] = C[j*width + i]/Npoints; 
			C[j*width + i] = temp/Npoints;
		}
		for(j = height/2; j < height; j++)
		{/* swap upper left with lower right */
			temp = C[(j - height/2)*width + i + width/2];
			C[(j - height/2)*width + i + width/2] = C[j*width + i]/Npoints;
			C[j*width + i] = temp/Npoints;
		}
	}
}

int get_fine_horiz_offset(float *imagebuf, float *hires_template, double *pix_resp, int *fine_col_offset, int crude_col_offset, int crude_row_offset, int anchor_col, int anchor_row)
{
	int i, j, m, n, r, s, col, row;
	int patch_width = FOCPLWIDTH/8;
	int box_width = 200;

	float *sampled_offset_template = (float *)calloc(box_width*box_width, sizeof(float));

	/* debug variables */
	/*
	int start_col, start_row;
	float *data_cutout = (float *)calloc(box_width*box_width, sizeof(float));
	char *cmd = (char *)calloc(STRSIZE, sizeof(char));
	char *name = (char *)calloc(STRSIZE, sizeof(char));
	*/
	/* */

	double dotproduct[BINDENS];
	long int HR_x, HR_y;
	int max_dotproduct_index;

	for(r = 0; r < BINDENS; r++)
	{	/*loop through fractional offsets */
		for(s = BINDENS/2; s <= BINDENS/2; s++)
		{	
			for(i = 0; i < box_width; i++)
			{
				for(j = 0; j < box_width; j++)
				{
					for(m = 0; m < BINDENS; m++)
					{
						for(n = 0; n < BINDENS; n++)
						{	/*form detector-sampled template based on current combined crude + fine offset */
							HR_x = ((anchor_col/patch_width)*patch_width + patch_width/2 - box_width/2 + i - crude_col_offset)*BINDENS - (r - BINDENS/2) + m;
							HR_y = ((anchor_row/patch_width)*patch_width + patch_width/2 - box_width/2 + j - crude_row_offset)*BINDENS - (s - BINDENS/2) + n;
							sampled_offset_template[j*box_width + i] += hires_template[HR_y*HRFOCPLWIDTH + HR_x]*pix_resp[n*BINDENS + m]/(BINDENS*BINDENS);
							/*debug begin*/
							/*
							if(i == 0 && j == 0 && m == 0 && n == 0 && r == BINDENS/2 && s == BINDENS/2)
							{
								printf("data cutout start x,y = %d, %d\n", (anchor_col/patch_width)*patch_width + patch_width/2 - box_width/2, (anchor_row/patch_width)*patch_width + patch_width/2 - box_width/2);
								printf("HR_x, HR_y = %ld, %ld, ", HR_x, HR_y);
								printf("hires_template[%ld, %ld] = %.3e\n", HR_x, HR_y, hires_template[HR_y*HRFOCPLWIDTH + HR_x]);
							}
							*/
							/*
							if(i == 117 && j == 177 && r == BINDENS/2)
							{
								printf("hires_template[%ld, %ld] = %.3e\n", HR_x, HR_y, hires_template[HR_y*HRFOCPLWIDTH + HR_x]);
							}
							*/
							/*debug end*/
						}
					}
					/* debug begin */
					/*
					if(i == 117 && j == 177 && r == 3)
					{
						printf("sampled_offset_template[%d, %d] = %.3e\n", i, j, sampled_offset_template[j*box_width + i]);
					}
					*/
					/* debug end */
				}
			}
			/* debug begin */
			/*
			long axislength[2];
			axislength[0] = box_width;
			axislength[1] = box_width;
			if(1)
			{
				start_col = (anchor_col/patch_width)*patch_width + patch_width/2 - box_width/2;
				start_row = (anchor_row/patch_width)*patch_width + patch_width/2 - box_width/2;
				if(check_file_exist("template_cutout.fits"))
				{
					sprintf(cmd, "rm template_cutout.fits"); 
					system(cmd);
				}
				sprintf(name, "template_cutout_%d_%d.fits", r, s);
				if(check_file_exist(name))
				{
					sprintf(cmd, "rm %s", name);
					system(cmd);
				}
				write_simple_fits_float_data(name, 2, axislength, sampled_offset_template);
				copy_array_2d_float(imagebuf, data_cutout, FOCPLWIDTH, FOCPLWIDTH, box_width, box_width, start_col, start_row, 0, 0, box_width, box_width);
				if(check_file_exist("data_cutout.fits"))
				{
					sprintf(cmd, "rm data_cutout.fits");
					system(cmd);
				}
				write_simple_fits_float_data("data_cutout.fits", 2, axislength, data_cutout);
			}
			*/
			/* debug end */
			/*determine sum of products between the template and image cutouts for the current offset*/
			dotproduct[r] = 0;
			for(i = 0; i < box_width; i++)
			{
				for(j = 0; j < box_width; j++)
				{
					col = (anchor_col/patch_width)*patch_width + patch_width/2 - box_width/2 + i;
					row = (anchor_row/patch_width)*patch_width + patch_width/2 - box_width/2 + j;
					dotproduct[r] += (double)(sampled_offset_template[j*box_width + i]*imagebuf[row*FOCPLWIDTH + col]);
				}
			}
			/*debug begin*/
			/*
			printf("dotproduct for r = %d, s = %d: %0.2f\n", r, s, dotproduct[s*BINDENS + r]); 

			if(r == 3 && s == 4)
			{
				printf("start_col, start_row = %d, %d\n", start_col, start_row);
				multiply_images(data_cutout, sampled_offset_template, product_array, box_width);
				axislength[0] = box_width;
				axislength[1] = box_width;
				write_simple_fits_float_data("data_cutout.fits", 2, axislength, data_cutout);
				char *product_name = (char *)calloc(STRSIZE, sizeof(char));
				sprintf(product_name, "product_%d_%d.fits", r, s);
				if(check_file_exist(product_name))
				{
					sprintf(cmd, "rm %s", product_name);
					system(cmd);
				}
				write_simple_fits_float_data(product_name, 2, axislength, product_array);
				free(product_name);
				float sum = 0;
				for(i = 0; i < box_width*box_width; i++)
				{
					sum += product_array[i];
				}
				printf("sum of product array for r = %d, s = %d: %0.2f\n", r, s, sum); 
				memset(product_array, 0, box_width*box_width*sizeof(float));
				memset(data_cutout, 0, box_width*box_width*sizeof(float));
			}
			*/
			/*debug end */
			memset(sampled_offset_template, 0, box_width*box_width*sizeof(float));
		}
	}

	/* debug begin */
	/*
	printf("dot product values: \n");
	for(i = 0; i < BINDENS; i++)
	{
		printf("%3e ", dotproduct[i]);
	}
	printf("\n");
	*/
	/* debug end */

	max_dotproduct_index = gsl_stats_max_index(dotproduct, 1, BINDENS);

	*fine_col_offset = max_dotproduct_index - BINDENS/2; 	

	printf("Fine x offset: %d bins = %0.2f detector pixels\n", *fine_col_offset, (float)(*fine_col_offset)/BINDENS);

	/* debug begin */
	/*
	long axislength[2];
	axislength[0] = BINDENS;
	axislength[1] = BINDENS;
	if(check_file_exist("dotproduct.fits"))
	{
		sprintf(cmd, "rm dotproduct.fits"); 
		system(cmd);
	}
	write_simple_fits_double_data("dotproduct.fits", 2, axislength, dotproduct);	
	*/
	/*
	free(data_cutout);
	*/
	/* debug end */

	free(sampled_offset_template);

	/* debug begin */
	/*
	free(cmd);
	free(name);
	*/
	/*
	free(product_array);
	*/
	/* debug end */

	return 0;
}

int get_crude_offset(float *imagebuf, int *col_offset, int *row_offset, int init_col_offset, int init_row_offset, int init_flag, int anchor_col, int anchor_row, char *laser_template_name)
{
	float *laser_template;
	long naxes, axis_length[2];
	int nd, u, v, r, s, box_width, patch_width, max_x_lag, max_y_lag;
	box_width = 200;
	patch_width = FOCPLWIDTH/8;
	if(init_flag)
	{
		max_x_lag = 2;
		max_y_lag = 2;
	}
	else
	{
		init_col_offset = 0;
		init_row_offset = 0;
		max_x_lag = 4;
		max_y_lag = 4;
	}
//allocate arrays
	float **template_box = (float **)calloc(box_width, sizeof(float *));
	float **data_box = (float **)calloc(box_width, sizeof(float *));
	for(u = 0; u < box_width; u++)
	{
		template_box[u] = (float *)calloc(box_width, sizeof(float));
		data_box[u] = (float *)calloc(box_width, sizeof(float));
	}
//copy data to arrays
	if ((laser_template = load_simple_fits_float_data(laser_template_name, &naxes, axis_length, &nd)) == 0)
	{
		(void) stdout, fprintf(stdout, "get_crude_offset: load_simple_fits_float_data %s failed\n", laser_template_name);
		return -1;
	}
	r = (anchor_col/patch_width)*patch_width + patch_width/2 - box_width/2;
	//r = anchor_col - box_width/2;
	for(u = 0; u < box_width; u++)
	{
		s = (anchor_row/patch_width)*patch_width + patch_width/2 - box_width/2;
		/* debug begin */
		/*
		if(u == 0)
		{
			printf("crude align cutouts start at x,y = %d, %d\n", r, s);
		}
		*/
		/* debug end */
		//s = anchor_row - box_width/2;
		for(v = 0; v < box_width; v++)
		{
			template_box[u][v] = laser_template[s*FOCPLWIDTH + r];
			data_box[u][v] = imagebuf[s*FOCPLWIDTH + r];
			s++;
		}
		r++;
	}
	free(laser_template);	
	get_peak_corr(template_box, data_box, box_width, max_x_lag, max_y_lag, init_col_offset, init_row_offset, col_offset, row_offset);

//free allocated arrays
	for(u = 0; u < box_width; u++)
	{
		free(template_box[u]);
		free(data_box[u]);
	}
	free(template_box);
	free(data_box);

	return 0;
}

int get_peak_region(float *orig, int grid_size, int *peak_u, int *peak_v)
{
	int u, v, r, s, grid_box_width, max_index;
	grid_box_width = FOCPLWIDTH/grid_size;
	double *gridsum = (double *)calloc(grid_size*grid_size, sizeof(double));
	double sum;

	for(u = 0; u < grid_size; u++)
	{
		for(v = 0; v < grid_size; v++)
		{
			sum = 0;
			for(r = 0; r < grid_box_width; r++)
			{
				for(s=0 ; s < grid_box_width; s++)
				{
					sum += orig[(v*grid_box_width + s)*FOCPLWIDTH + u*grid_box_width + r];
				}
			}
			gridsum[v*grid_size + u] = sum;
		}
	}

	max_index = gsl_stats_max_index(gridsum, 1, grid_size*grid_size);	

	*peak_v = (max_index/grid_size)*grid_box_width + grid_box_width/2;
	*peak_u = (max_index%grid_size)*grid_box_width + grid_box_width/2;

	free(gridsum);

	return 0;
}

int get_cube_coords(int xfocpl, int yfocpl, int *xlens, int *ylens)
{
	float cubeimagewidth = 193.0;
	float cubetilt = 18.5*PI/180;
	int x_cube_origin = 65;
	int y_cube_origin = 3;

	*xlens = (int)round((xfocpl*cos(cubetilt) - yfocpl*sin(cubetilt))*cubeimagewidth/FOCPLWIDTH) + x_cube_origin;
	*ylens = (int)round((xfocpl*sin(cubetilt) + yfocpl*cos(cubetilt))*cubeimagewidth/FOCPLWIDTH) + y_cube_origin;

	return 0;
}

int get_image_coords(int xlens, int ylens, int *xfocpl, int *yfocpl)  
{	
	float x_cube, y_cube;
	float cubeimagewidth = 193.0;
	float cubetilt = 18.5*PI/180;
	int x_cube_origin = 65;
	int y_cube_origin = 3;

	x_cube = (float)(xlens + 1 - x_cube_origin);
	y_cube = (float)(ylens + 1 - y_cube_origin);
	
	*xfocpl = (int)round((x_cube*cos(-cubetilt) - y_cube*sin(-cubetilt))*FOCPLWIDTH/cubeimagewidth);
	*yfocpl = (int)round((x_cube*sin(-cubetilt) + y_cube*cos(-cubetilt))*FOCPLWIDTH/cubeimagewidth);
	
	if(*xfocpl < 0)
	{
		*xfocpl = 0;
	}
	else if(*xfocpl > FOCPLWIDTH-1)
	{
		*xfocpl = FOCPLWIDTH-1;
	}
	if(*yfocpl < 0)
	{
		*yfocpl = 0;
	}
	else if(*yfocpl > FOCPLWIDTH-1)
	{
		*yfocpl = FOCPLWIDTH-1;
	}

	return 0;
}

int get_peak_corr(float **a, float **b, int width, int x_max_lag, int y_max_lag, int init_x_lag, int init_y_lag, int *x_offset, int *y_offset)
{
	int x_lag, y_lag, x_lag_index, y_lag_index, x_a, y_a, x_lag_range, y_lag_range;
	int max_corr_index;
	x_lag_range = x_max_lag*2 + 1;
	y_lag_range = y_max_lag*2 + 1;
	double corr[x_lag_range*y_lag_range];
	memset(corr, 0, x_lag_range*y_lag_range*sizeof(double));

	for(x_lag_index = 0; x_lag_index < x_lag_range; x_lag_index++)
	{
		for(y_lag_index = 0; y_lag_index < y_lag_range; y_lag_index++)
		{
			x_lag = x_lag_index - x_max_lag + init_x_lag;
			y_lag = y_lag_index - y_max_lag + init_y_lag;
			for(x_a = x_max_lag + abs(init_x_lag); x_a < width - x_max_lag - abs(init_x_lag); x_a++)
			{
				for(y_a = y_max_lag + abs(init_y_lag); y_a < width - y_max_lag - abs(init_y_lag); y_a++)
				{
					corr[y_lag_index*x_lag_range + x_lag_index] += (double)a[x_a][y_a]*b[x_a + x_lag][y_a + y_lag];
				}
			}
		}
	}

	max_corr_index = gsl_stats_max_index(corr, 1, x_lag_range*y_lag_range);
	//printf("max corr index = %d\n", max_corr_index);

	*x_offset = max_corr_index%x_lag_range - x_max_lag + init_x_lag;
	*y_offset = max_corr_index/x_lag_range - y_max_lag + init_y_lag;

	return 0;
}

int weightedsum_cubextract(float *dbuf, float *cube, float ***lookuptable, float *weights_J, float *weights_H, float *lensletflat, float *heightmap, float *tiltmap, float *resbgframe, FrameTok Token)
{
	int lensx, lensy, crudeXcntr, crudeYcntr, i, j, w;
	int fineXoffIndex, fineYoffIndex, weightPageIndex;
	float Xref, Yref, height, tilt, vertscale;
	/*float normamp;*/
	float Xcntr, Ycntr, extractsum;
	int reflambda = LAMBDA_REF;
	int weightboxwidth = 3;
	float lambda, lambdainc = (LAMBDA_MAX - LAMBDA_MIN)/(NCHAN - 1);
	int lenstally = 0;
	float flattally = 0, cubetally = 0;
	float meanflatval;

	/*variables for background subtraction*/
	int lambda_left = LAMBDA_MIN + lambdainc*8;
	int lambda_right = LAMBDA_MIN + lambdainc*20;
	float *resbg;
	resbg = (float *)calloc(8, sizeof(float));
	int resbgcount;
	int upperleft_darkcol, upperleft_darkrow, lowerright_darkcol, lowerright_darkrow, upperleft_col_begin, upperleft_row_begin, lowerright_col_begin, lowerright_row_begin;
	float image_edge_medvals[4];
	float edge_resbg_estimate;
	float local_resval;
				
	for(lensx = 0; lensx < CUBEWIDTH; lensx++)
	{//form residual count array from dark pixels in between spectra 
		for(lensy = 0; lensy < CUBEWIDTH; lensy++)
		{
			if(lookuptable[lensx][lensy][2] > 0)
			{//make sure the lenslet coordinates have an entry in the lookup table
				resbgcount = 0;
				//set coordinates of the corners of the dark boxes
				Xref = lookuptable[lensx][lensy][0] + Token.xshift;
				Yref = lookuptable[lensx][lensy][1] + Token.yshift;
				height = heightmap[lensy*CUBEWIDTH + lensx];
				tilt = tiltmap[lensy*CUBEWIDTH + lensx];
				//tilt = 0.03;
				vertscale = (1760 - 1100)/height;
				/*fix for FPM shadow in NRM data*/
				
				upperleft_col_begin = (int)lroundf(Xref + (reflambda - lambda_left)/vertscale*tilt) - 4;
				upperleft_row_begin = (int)lroundf(Yref + (reflambda - lambda_left)/vertscale);
				lowerright_col_begin = (int)lroundf(Xref + (reflambda - lambda_right)/vertscale*tilt) + 3;
				lowerright_row_begin = (int)lroundf(Yref + (reflambda - lambda_right)/vertscale);

				/*debug begin*/
				/*
				if(lensx == 110 && lensy == 120)
				{
					printf("upperleft_row_begin = %d, lowerright_row_begin = %d\n", upperleft_row_begin, lowerright_row_begin);
				}
				*/
				/*debug end*/
	
				for(i = 0; i < 2; i++)
				{
					for(j = 0; j < 2; j++)
					{//grab dark pixels for this lenslet
						upperleft_darkcol = upperleft_col_begin + i;
						upperleft_darkrow = upperleft_row_begin + j;
						if(upperleft_darkcol > BORDER && upperleft_darkcol < FOCPLWIDTH - BORDER - 1  && upperleft_darkrow > BORDER && upperleft_darkrow < FOCPLWIDTH - BORDER - 1)
						{//if the dark pixel is in range, grab it
							resbg[resbgcount] = dbuf[upperleft_darkrow*FOCPLWIDTH + upperleft_darkcol];
							resbgcount++;
						}
						lowerright_darkcol = lowerright_col_begin + i;
						lowerright_darkrow = lowerright_row_begin + j;
						if(lowerright_darkcol > BORDER && lowerright_darkcol < FOCPLWIDTH - BORDER - 1 && lowerright_darkrow > BORDER && lowerright_darkrow < FOCPLWIDTH - BORDER - 1)
						{//if the dark pixel is in range, grab it
							resbg[resbgcount] = dbuf[lowerright_darkrow*FOCPLWIDTH + lowerright_darkcol];
							resbgcount++;
						}
					}
				}//end of dark pixel grabbing loop
				if(resbgcount >= 4)
				{
					qsort(resbg, resbgcount, sizeof(float), compare_floats);
					/*resbgframe[lensy*CUBEWIDTH + lensx] = resbg[resbgcount/2];*/
					resbgframe[lensy*CUBEWIDTH + lensx] = resbg[lroundf(0.2*resbgcount)];
				}
				else
				{
					/*	
					printf("WARNING only %d residual bias pixels found for lenslet %d, %d\n", resbgcount, lensx, lensy);			
					*/
					resbgframe[lensy*CUBEWIDTH + lensx] = UNDEFINED;
				}
			}
			else
			{
				resbgframe[lensy*CUBEWIDTH + lensx] = UNDEFINED;
			}
		}
	}
	//median filter the resbg map
	median_filter(resbgframe, CUBEWIDTH, 3);

	//form edge resbg estimate from corners
	image_edge_medvals[0] = resbgframe[238*CUBEWIDTH + 183];
	image_edge_medvals[1] = resbgframe[181*CUBEWIDTH + 5];
	image_edge_medvals[2] = resbgframe[7*CUBEWIDTH + 63];
	image_edge_medvals[3] = resbgframe[63*CUBEWIDTH + 239];
	//sort all quadrants
	qsort(image_edge_medvals, 4, sizeof(float), compare_floats);	
	edge_resbg_estimate = (image_edge_medvals[1] + image_edge_medvals[2])/2;   
	printf("extract_cube: residual count/sec estimate at edge of image = %0.2f\n", edge_resbg_estimate);

	for(lensx = 0; lensx < CUBEWIDTH; lensx++)
	{
		for(lensy = 0; lensy < CUBEWIDTH; lensy++)
		{
			if(lookuptable[lensx][lensy][2] > 0)
			{
				Xref = lookuptable[lensx][lensy][0] + Token.xshift;
				Yref = lookuptable[lensx][lensy][1] + Token.yshift;
				height = heightmap[lensy*CUBEWIDTH + lensx];
				tilt = tiltmap[lensy*CUBEWIDTH + lensx];
				vertscale = (1760 - 1100)/height;
				for(w = 0; w < NCHAN; w++)
				{
					lambda = LAMBDA_MIN + lambdainc*w;
					crudeXcntr = (int)lroundf(Xref + (LAMBDA_REF - lambda)/vertscale*tilt);
					crudeYcntr = (int)lroundf(Yref + (LAMBDA_REF - lambda)/vertscale);
					if(crudeXcntr > BORDER && crudeXcntr < FOCPLWIDTH - BORDER - 1 && crudeYcntr > BORDER && crudeYcntr < FOCPLWIDTH - BORDER - 1)
					{
						extractsum = 0;
						/* debug begin */
						/*
						if(lensx == 248 && lensy == 82)
						{
							printf("col 2 = %0.2f, height = %0.2f, tilt = %0.2f\n", lookuptable[lensx][lensy][2], height, tilt);
						}
						*/
						/* debug end */
						Xcntr = Xref + (LAMBDA_REF - lambda)/vertscale*tilt;
						Ycntr = Yref + (LAMBDA_REF - lambda)/vertscale;

						/* debug begin */
						/*
						if(lensx == 82 && lensy == 152 && (w == 0 || w == 22))
						{
							printf("w = %d: Extraction center X, Y = %.1f, %.1f\n", w, Xcntr, Ycntr);
	 					}
						*/	
						/* debug end */

						fineXoffIndex = (int)lroundf((Xcntr - crudeXcntr)*BINDENS) + BINDENS/2;
						fineYoffIndex = (int)lroundf((Ycntr - crudeYcntr)*BINDENS) + BINDENS/2;
						if(fineXoffIndex > BINDENS - 1)
						{
							fineXoffIndex = BINDENS - 1;
						}
						else if(fineXoffIndex < 0)
						{
							fineXoffIndex = 0;
						}
						if(fineYoffIndex > BINDENS - 1)
						{
							fineYoffIndex = BINDENS - 1;
						}
						else if(fineYoffIndex < 0)
						{
							fineYoffIndex = 0;
						}
						weightPageIndex = fineYoffIndex*BINDENS + fineXoffIndex;

						if(localbgsubswitch)
						{
							if(resbgframe[lensy*CUBEWIDTH + lensx] != UNDEFINED)
							{
								local_resval = resbgframe[lensy*CUBEWIDTH + lensx];
							}
							else
							{
								printf("extract_cube: WARNING undefined value in resbg array at %d,%d\n", lensx, lensy);
								local_resval = edge_resbg_estimate;
							}
							for(i = 0; i <= 2; i++)
							{
								for(j = 0; j<= 2; j++)
								{
									if(lambda <= LAMBDA_REF)
									{
										extractsum += weights_J[weightPageIndex*weightboxwidth*weightboxwidth + weightboxwidth*j + i]*(dbuf[(crudeYcntr + j - 1)*FOCPLWIDTH + crudeXcntr + i - 1] - local_resval);
									}
									else
									{
										extractsum += weights_H[weightPageIndex*weightboxwidth*weightboxwidth + weightboxwidth*j + i]*(dbuf[(crudeYcntr + j - 1)*FOCPLWIDTH + crudeXcntr + i - 1] - local_resval);	
									}
								}
							}
						}/*end of block for computing extraction sum with local b.g. subtraction*/
						else
						{//not using local background subtraction
							for(i = 0; i <= 2; i++)
							{
								for(j = 0; j<= 2; j++)
								{
									if(lambda <= LAMBDA_REF)
									{
										extractsum += weights_J[weightPageIndex*weightboxwidth*weightboxwidth + weightboxwidth*j + i]*dbuf[(crudeYcntr + j - 1)*FOCPLWIDTH + crudeXcntr + i - 1];
									}
									else
									{
										extractsum += weights_H[weightPageIndex*weightboxwidth*weightboxwidth + weightboxwidth*j + i]*dbuf[(crudeYcntr + j - 1)*FOCPLWIDTH + crudeXcntr + i - 1];	
									}
								}
							}
						}
						assert(lensletflat[lensy*CUBEWIDTH + lensx] > 0);
						cube[w*CUBEWIDTH*CUBEWIDTH + lensy*CUBEWIDTH + lensx] = extractsum/lensletflat[lensy*CUBEWIDTH + lensx];
						lenstally++;
						flattally += lensletflat[lensy*CUBEWIDTH + lensx];
					}/*end of block testing for valid extraction location*/
				}
			}
		}
	}

	meanflatval = flattally/lenstally;
	for(lensx = 0; lensx < CUBEWIDTH; lensx++)
	{//rescale cube based on mean-normalized flat field, and tally the full signal
		for(lensy = 0; lensy < CUBEWIDTH; lensy++)
		{
			for(w = 0; w < NCHAN; w++)
			{
				if(cube[w*CUBEWIDTH*CUBEWIDTH + lensy*CUBEWIDTH + lensx] > 0)
				{
					cube[w*CUBEWIDTH*CUBEWIDTH + lensy*CUBEWIDTH + lensx] *= meanflatval; 
					cubetally += cube[w*CUBEWIDTH*CUBEWIDTH + lensy*CUBEWIDTH + lensx];
				}
			}
		}
	}
	printf("Cube tally: %.2e\n", cubetally);

	free(resbg);
	return 0;
}

int median_filter(float *image, int imwidth, int winwidth)
{
	int col, row, i, j, p;
	float *window = (float *)calloc(winwidth*winwidth, sizeof(float)); 
	float *medianimage = (float *)calloc(imwidth*imwidth, sizeof(float));

	if(winwidth % 2 == 0)
	{//winwidth is even
		printf("median_filter: WARNING changing window width (%d) to an odd number (%d)\n", winwidth, winwidth-1); 
		winwidth = winwidth - 1;
	}
	for(col = 0; col < imwidth; col++)
	{
		for(row = 0; row < imwidth; row++)
		{
			p = 0;
			for(i = -winwidth/2; i <= winwidth/2; i++)
			{//loop through peripheral pixels in the window around the current position
				for(j = -winwidth/2; j <= winwidth/2; j++)
				{
					if(row+j>0 && row+j<imwidth && col+i>0 && col+i<imwidth)
					{//this peripheral pixel is inside image bounds 
						if(image[(row + j)*imwidth + (col + i)] != UNDEFINED)
						{//this peripheral pixel has a valid value
							window[p] = image[(row + j)*imwidth + (col + i)];
							p++;
						}
					}
				}
			}
			if(p >= 3)
			{
				qsort(window, p, sizeof(float), compare_floats);
				medianimage[row*imwidth + col] = window[p/2];
			}
			else
			{
				medianimage[row*imwidth + col] = UNDEFINED;
			}
		}
	}
	//copy over median-filtered result to image
	for(col = 0; col < imwidth; col++)
	{
		for(row = 0; row < imwidth; row++)
		{
			image[row*imwidth + col] = medianimage[row*imwidth + col];
		}
	}

	free(medianimage);
	free(window);

	return 0;
}

int clean_badpix(float *image, int **badpix_list, int Nbadpix)
{
	int i, neighb_count, row, col, m, n;
	float *neighb_arr = (float *)calloc(8, sizeof(float));
	for(i = 0; i < Nbadpix; i++)
	{
		col = badpix_list[i][0];
		row = badpix_list[i][1];

		neighb_count = 0;

		for(m = col-1; m <= col+1; m += 2)
		{
			for(n = row-1; n <= row+1; n += 2)
			{
				if(m >= BORDER && m < FOCPLWIDTH - BORDER && n >= BORDER && n < FOCPLWIDTH - BORDER)
				{
					neighb_arr[neighb_count] = image[n*FOCPLWIDTH + m];
					neighb_count++;
				}
			}
		}
		qsort(neighb_arr, neighb_count, sizeof(float), compare_floats);
		if(neighb_count % 2 == 0)
		{
			image[row*FOCPLWIDTH + col] = (neighb_arr[neighb_count/2 - 1] + neighb_arr[neighb_count/2])/2;
		}
		else
		{
			image[row*FOCPLWIDTH + col] = neighb_arr[neighb_count/2];
		}
	}

	free(neighb_arr);

	return 0;
}

int apply_dewarflat(float *image, float *flat)
{
	int col, row;

	for(col = BORDER; col < FOCPLWIDTH - BORDER; col++)
	{
		for(row = BORDER; row < FOCPLWIDTH - BORDER; row++)
		{
			if(flat[row*FOCPLWIDTH + col] > 0.05)//threshold to exclude dead pixels
			{
				image[row*FOCPLWIDTH + col] = image[row*FOCPLWIDTH + col]/flat[row*FOCPLWIDTH + col];
			}
		}
	}

	return 0;
}

void copy_array_2d_double(double *src, double *dest, int width_src, int height_src, int width_dest, int height_dest, int xsrcbeg, int ysrcbeg, int xdestbeg, int ydestbeg, int width_copy, int height_copy)
{
	int col, row;

	for(col = xdestbeg; col < xdestbeg + width_copy; col++)
	{
		for(row = ydestbeg; row < ydestbeg + height_copy; row++)
		{
			if(col < width_dest && row < height_dest)
			{/* boundary check */
				dest[width_dest*row + col] = src[width_src*(row - ydestbeg + ysrcbeg) + col - xdestbeg + xsrcbeg];
			}
		}
	}
}

void copy_array_2d_float(float *src, float *dest, int width_src, int height_src, int width_dest, int height_dest, int xsrcbeg, int ysrcbeg, int xdestbeg, int ydestbeg, int width_copy, int height_copy)
{
	int col, row;

	for(col = xdestbeg; col < xdestbeg + width_copy; col++)
	{
		for(row = ydestbeg; row < ydestbeg + height_copy; row++)
		{
			if(col < width_dest && row < height_dest)
			{/* boundary check */
				dest[width_dest*row + col] = src[width_src*(row - ydestbeg + ysrcbeg) + col - xdestbeg + xsrcbeg];
			}
		}
	}
}

void copy_array_2d_float_to_double(float *src, double *dest, int width_src, int height_src, int width_dest, int height_dest, int xsrcbeg, int ysrcbeg, int xdestbeg, int ydestbeg, int width_copy, int height_copy)
{
	int col, row;

	for(col = xdestbeg; col < xdestbeg + width_copy; col++)
	{
		for(row = ydestbeg; row < ydestbeg + height_copy; row++)
		{
			if(col < width_dest && row < height_dest)
			{/* boundary check */
				dest[width_dest*row + col] = (double)(src[width_src*(row - ydestbeg + ysrcbeg) + col - xdestbeg + xsrcbeg]);
			}
		}
	}
}

int subtract_images(float *imageone, float *imagetwo, float *result, int width)
{
	int i;

	for(i = 0; i < width*width; i++)
	{
		result[i] = imageone[i] - imagetwo[i];
	}

	return 0;
}

int multiply_images(float *imageone, float *imagetwo, float *result, int width)
{
	int i;

	for(i = 0; i < width*width; i++)
	{
		result[i] = imageone[i]*imagetwo[i];
	}

	return 0;
}

void get_int_date(char *datestr, int *year, int *month, int *day)
{
    char *yearstr, *monthstr, *daystr;
    char *sepstr = (char *) calloc(STRSIZE, sizeof(char));
    char *datestr_copy = (char *) calloc(STRSIZE, sizeof(char));
	strcpy(datestr_copy, datestr);

/*	extract year, month, and day */
    strcpy(sepstr, "-");
    yearstr = strtok(datestr_copy, sepstr);
    monthstr = strtok(0, sepstr);
    daystr = strtok(0, sepstr);
    *year = atoi(yearstr);
    *month = atoi(monthstr);
    *day = atoi(daystr);
/*
	printf("year = %.2d, month = %.2d, day = %.2d\n", *year, *month, *day);
*/
    free(sepstr);
    free(datestr_copy);
}

void fixdate(char *olddate, char *newdate)
{
	int i, length, year, month, day;
	char *sepstr, *yearstr, *monthstr, *daystr; 
	sepstr = (char *) calloc(STRSIZE, sizeof(char));

	length = strlen(olddate);
	
	for(i=1; i < length - 1; i++)
	{
		if(olddate[i] == '/')
		{
			olddate[i] = '_';
		}
		if(olddate[i] == ' ')
		{
			break;
		}
		newdate[i-1] = olddate[i];
	}		
	newdate[i-1] = 0;

//extract year, month, and day
	strcpy(sepstr, "_");
	yearstr = strtok(newdate, sepstr); 
	monthstr = strtok(0, sepstr);
	daystr = strtok(0, sepstr);
	year = atoi(yearstr);
	month = atoi(monthstr);
	day = atoi(daystr);

//	printf("year = %.2d, month = %.2d, day = %.d\n", year, month, day);
	
	sprintf(newdate, "%.4d-%.2d-%.2d", year, month, day);
	free(sepstr);
}

int getfinalcubename(char *processeddatadir, char *focalfits, char *cubefits, char *calcubefits, char *collapsedir, char *collapsedname, int cubetype, int frame, char *shortname)
{
	int status, i, startindex, length, exp, blankheaderflag;
	fitsfile *old;
	char *keyword, *keyval, *comment, *datestr, *typestr, *expstr, *objectstr, *cubedir;

	status = 0;
	blankheaderflag = 0;
	startindex = 0;
	keyword = (char *) calloc(STRSIZE, sizeof(char));
	keyval = (char *) calloc(STRSIZE, sizeof(char));
	datestr = (char *) calloc(STRSIZE, sizeof(char));
	expstr = (char *) calloc(STRSIZE, sizeof(char));
	typestr = (char *) calloc(STRSIZE, sizeof(char));
	objectstr = (char *) calloc(STRSIZE, sizeof(char));
	cubedir = (char *) calloc(STRSIZE, sizeof(char));
	comment = (char *) calloc(STRSIZE, sizeof(char));
	
//create the cube filename based on the object, date, type, and frame # of the original focal plane image
	fits_open_file(&old, focalfits, READONLY, &status);
//	fits_get_hdrspace(old, &nkeys, nullptr, &status); 
	if(status)
	{
		fits_report_error(stdout, status);
		return(status);
	}
	strcpy(keyword, "UT_DATE");
	fits_read_keyword(old, keyword, keyval, comment, &status);               
	if(status)
	{
		fits_report_error(stdout, status);
		fprintf(stdout, "getfinalcubename WARNING: No UT_DATE keyword in header of %s\n", focalfits); 
		sprintf(datestr, "0000-00-00");
		blankheaderflag = 1;
	}
	else
	{
		fixdate(keyval, datestr);
	}

	length = strlen(focalfits);
	if(blankheaderflag)
	{
		abbrev_fitsname(focalfits, expstr);
	}
	else
	{
		fits_read_key(old, TINT, "EXP_NUM", &exp, comment, &status);
	}
	if(!blankheaderflag)	
	{
		strcpy(keyword, "IMG_TYPE");
		fits_read_keyword(old, keyword, keyval, comment, &status);               
		if(status)
		{
			fits_report_error(stdout, status);
			fprintf(stdout, "getfinalcubename WARNING: No IMG_TYPE keyword in header of %s\n", focalfits); 
			sprintf(typestr, "NOTYPE");
			blankheaderflag = 1;
		}
		else
		{
			length = strlen(keyval);
			for(i = 1; i < length; i++)
			{
				if(keyval[i] == ' ')
				{
					break;
				}	
				else
				{
					typestr[i-1] = keyval[i];
				}
			}
			typestr[i-1] = 0;
		}
	}

	strcpy(keyword, "OBJECT");
	fits_read_keyword(old, keyword, keyval, comment, &status);
	if(status)
	{
		fits_report_error(stdout, status);
		fprintf(stdout, "getfinalcubename WARNING: No OBJECT keyword in header of %s\n", focalfits); 
		strcpy(objectstr, expstr);
		blankheaderflag = 1;
	}
	else
	{
		replaceblanks(keyval, objectstr);
	}

//Create the .fits data cube name
	if(!blankheaderflag)
	{
		if(cubetype == FITS)
		{
			sprintf(cubedir, "%s/%s/%s/CUBES", processeddatadir, objectstr, datestr);
			sprintf(cubefits, "%s/%s_%s_%s_%.3d.fits", cubedir, objectstr, typestr, datestr, exp);
			sprintf(calcubefits, "%s/%s_%s_%s_%.3d_speccal.fits", cubedir, objectstr, typestr, datestr, exp);
			sprintf(collapsedname, "%s/%s_%s_%s_%.3d", collapsedir, objectstr, typestr, datestr, exp);
		}
		else
		{
			fprintf(stdout, "getfinalcubename: unrecognized cubetype %d\n", cubetype);
			return -1;
		}
	}
	else
	{
		if(cubetype == FITS)
		{
			sprintf(cubedir, "%s/%s/%s/CUBES", processeddatadir, objectstr, datestr);
			sprintf(cubefits, "%s/%s_cube.fits", cubedir, expstr);
			sprintf(calcubefits, "%s/%s_cube_speccal.fits", cubedir, expstr);
			sprintf(collapsedname, "%s/%s", collapsedir, expstr);
		}
		else
		{
			fprintf(stdout, "getfinalcubename: unrecognized cubetype %d\n", cubetype);
			return -1;
		}
	}

	fits_close_file(old, &status);

	free(keyword); free(keyval); free(datestr); free(typestr); free(expstr); free(objectstr);
	free(cubedir); free(comment);

	return 0;
}

void abbrev_fitsname(char *srcstr, char *deststr)
{
	int i, j, startindex, length;
	length = strlen(srcstr);
	startindex = 0;

	for(i=length-1; i>0; i--)
	{
		if(srcstr[i] == '/')
		{
			startindex = i+1;
			break;
		}
	}
	if(startindex == 0)
	{
		printf("abbrev_fitsname: WARNING no '/' character found in given filename, copying entire name to abbreviated version.\n");
	}

	for(i=0, j=startindex; j<length-5; i++, j++)
	{
		deststr[i] = srcstr[j];	
	} 	
	deststr[i+1] = 0;
}

void chop_off_leading_path(char *srcstr, char *deststr)
{
	int i, j, startindex, length;
	length = strlen(srcstr);
	startindex = 0;

	for(i=length-1; i>0; i--)
	{
		if(srcstr[i] == '/')
		{
			startindex = i+1;
			break;
		}
	}
	if(startindex == 0)
	{
		printf("chop_off_leading_path: WARNING no '/' character found in given filename, copying entire name to abbreviated version.\n");
	}

	for(i=0, j=startindex; j<length; i++, j++)
	{
		deststr[i] = srcstr[j];	
	} 	
	deststr[i+1] = 0;
}

void replaceblanks(char *srcstr, char *deststr)
{
	int j, k, enddist, offset;

	offset = 1;	
	char blankstr[STRSIZE];
	char endsubstr[STRSIZE];

	for(j=1; j < (strlen(srcstr)-1); j++)
	{
		if(srcstr[j] == ' ')
		{
			enddist = (strlen(srcstr)-2)-j;
			if(enddist > 0)
			{
				for(k = 0; k < enddist; k++)
				{
					blankstr[k] = ' ';										
					endsubstr[k] = srcstr[j+k];
				}	
				blankstr[k] = '\0';
				endsubstr[k] = '\0';
				if(strcmp(blankstr, endsubstr) == 0) //only blanks left?
				{
		//			printf("found trailing blank. j = %d\n", j);
					break;
				}
				else
				{
					offset++;
				}
			}
			else 
			{
				break;
			}
		}
		else //copy over normally
		{
			deststr[j-offset] = srcstr[j];
		}
	}
	deststr[j-1] = '\0';
}

int check_file_type(char *filename)
{
	fitsfile *fptr;	
	int status = 0;
	int naxis, axisone_width, axistwo_width;
	char *nullstr = NULL;
	if(strstr(filename, ".fits") != 0 && ((strstr(filename, "_P_") == 0 && strstr(filename, "_F_") == 0 && strstr(filename, "_D_") == 0) || doallswitch))
	{
		//name is ok, now check image size
		fits_open_file(&fptr, filename, READONLY, &status);
		if(status)
		{
			fits_report_error(stdout, status);
			return 0;
		}
		fits_read_key(fptr, TINT, "NAXIS", &naxis, nullstr, &status);
		fits_read_key(fptr, TINT, "NAXIS1", &axisone_width, nullstr, &status);
		fits_read_key(fptr, TINT, "NAXIS2", &axistwo_width, nullstr, &status);
		fits_close_file(fptr, &status);	
		if(status)
		{
			fits_report_error(stdout, status);
			return 0;
		}
		else
		{
			if(naxis == 2 && axisone_width == FOCPLWIDTH && axistwo_width == FOCPLWIDTH)	
			{
				return 1;
			}
			else
			{
				printf("check_file_type: %s has invalid image size %d x %d\n", filename, axisone_width, axistwo_width);
				return 0;
			}
		}
	}
	else
	{
		return 0;
	}
}

int check_dir_exist(char *dirname)
{
	DIR *dirptr = opendir(dirname);
	if(dirptr)
	{
		closedir(dirptr);
		return 1;
	}
	else
	{
		return 0;
	}
}

int check_file_exist(char *filename)
{
	FILE *fileptr = fopen(filename, "r");
	if(fileptr)
	{
		fclose(fileptr);
		return 1;
	}
	else
	{
		return 0;
	}
}