#define VERSION "3.0 -- June 27, 2009"

/***********************************************************************\
* Note that all functions whose prototypes are declared before the      *
* start of the first function code will be found in an external file    *
* (usually because more than one function calls them).  Functions whose *
* prototypes are declared inside another function will be found in the  *
* same file as the calling function (usually because they are used only *
* by one calling function)						*
\***********************************************************************/

#include <string.h>
#include <stdio.h>
#include <stdlib.h>
#include <curses.h>
#include "nrutil.h"
#include "structs.h"
#include "debug.h"

#define NINTERP 100

float Y2[NINTERP+1];
int ninterp = NINTERP;

void help ();
void read_input (struct inpars *input, struct fitpars *fpar);
int assoc_trunc_comps (struct fitpars *fpar);
void link_trunc_fpar (struct fitpars *fpar);
void realign_trunc_mirrors (struct fitpars *fpar);
int assign_pars (struct inpars *input, struct image *, struct image *, 
    struct image *, struct image *, struct fitpars *, struct convpars *,
    struct cons *constr);
int readfits (struct image *);
float sigcheck (struct image sigma);
void psfcheck (struct image *psf);
struct image *getkernel (char *kernel_file);
float psf_fwhm (struct image psf);
int badpixlist (struct image *badpix);
int read_cons (struct inpars input, struct cons *constr, struct fitpars *fpars);
void errorcheck (struct inpars *input, struct image *data, 
    struct image *badpix, struct image *psf, struct image *sigma);
void getmu (float *, int);
void search_min (struct image *, struct image *, struct image *, 
    struct image *, struct fitpars *, struct cons *constr, struct convpars *);
void initcurses (char device[]);
void writefits (char *, struct image *, char *object, int);
void detach_trunc (struct fitpars *fpar);
void del_fitpars_struct (struct fitpars *fpar);
void del_constr_struct (struct cons *constr);

struct inpars input = {"", "", "", "none", "", "", "none", 1, "none", "", 
                     {0., 0., 0.}, {0,0,0}, 0., {0.,1., 1.}, "both", 0};
struct sampars sample = {0, 0};
struct runflags flags = {0, 0., -1., 0, 0};
float fwhmpsf, sigval;

int (*pfunc)(const char *, ...);
struct image psf, *kernel;
char device[10];

int main (int argc, char *argv[])
{
    void blank_img (struct image *data, struct inpars input);

    char phead[FLEN_FILENAME], gname[FLEN_FILENAME];
    struct image data, sigma, badpix;
    struct fitpars fpar = {0, "", 0, NULL, NULL, NULL, NULL, NULL, 0, 0, 0, NULL, NULL, NULL};
    struct convpars cpar = {{0., 0., 0.}, {0, 0, 0}, "", 0, 0};
    struct cons constr = {NULL, "", NULL, 0, NULL, 0, NULL};
    int returnval, err, i, output=0;

    printf ("\nGALFIT Version %s\n\n", VERSION);

    for (i= 1; i < argc; i++) {
	if (strncmp(argv[i], "-noskyest", 9) == 0)
	    flags.noestsky = 1;
	else if (strncmp(argv[i], "-skyped", 7) == 0) {
	    flags.skyped = atof(argv[++i]);
	    flags.noestsky = 2;
	} else if (strncmp(argv[i], "-skyrms", 7) == 0) {
	    flags.skyrms = atof(argv[++i]);
	} else if (strncmp(argv[i], "-outsig", 6) == 0)
	    flags.outsig = 1;
	else if (strncmp(argv[i], "-i", 2) == 0)
	    flags.interact = 1;
	else if (strncmp(argv[i], "-help", 5) == 0) {
	    help ();
	    return(0);
	} else if (strncmp(argv[i], "-o1", 3) == 0)
	    output = 1;
	else if (strncmp(argv[i], "-o2", 3) == 0)
	    output = 2;
	else if (strncmp(argv[i], "-o3", 3) == 0)
	    output = 3;
	else
            strcpy (input.initparfile, argv[i]);
    };


    /*************************************************************************\
    *  Read in data, psf, and sigma images as well as initial fit parameters  *
    \*************************************************************************/

    read_input (&input, &fpar); 

    /******************************************************\
    *  Set up output to either be in regular window mode,  *
    *  curses mode, or combination of both.                *
    \******************************************************/

    strcpy (device, input.device);
    if (strncmp (device, "regular", 7)!= 0 && output == 0) {
        initcurses(device);  /* Initialize the curses window */
        if (strncmp (device, "curses", 6)== 0)
            pfunc = printw;
        else
            pfunc = printf;
    } else
        pfunc = printf;

    pfunc ("\n\n");

    /*************************************************************************\
    *  Figure out which profiles are truncated by which truncation component  *
    \*************************************************************************/

    err = assoc_trunc_comps (&fpar);
    if (err != 0) {
	pfunc ("-- ERROR: Component #%d is being truncated but it is not pointing \n", err);
	pfunc ("          to a truncation function.  Must quit now....\n\n");
	exit (1);
    }

    link_trunc_fpar (&fpar); 
    realign_trunc_mirrors (&fpar); 

    if (output != 0) input.create = output;

    /*******************************************************\
    *  Look for a file that contains parameter constraints  *
    \*******************************************************/

    err = read_cons (input, &constr, &fpar); 

    if (err == 1) 
	pfunc ("-- No constraint file found.\n");
    else if (err == 2) {
	pfunc ("\n");
	pfunc ("-- WARNING: there's a problem in the constraint file.\n\n");
    } else if (err == 3) {
	pfunc ("\n");
	pfunc ("-- QUIT REASON: trying to constrain a component which doesn't exist.\n\n");
	exit (1);
    }

    err = assign_pars (&input, &data, &sigma, &psf, &badpix, &fpar, &cpar, &constr);

    if (err == 4)
        pfunc ("-- WARNING: trying to constrain a parameter that is being held fixed.\n");

    strcpy (phead, data.name);     /* Output images will have the same  */
				   /* header as the primary input image */
  
    /*********************************************************************\
    *  Check to see that the input images exist, else deal with failures  *
    \*********************************************************************/

    if ((returnval = readfits (&data)) == 1) {
        pfunc ("-- WARNING: No input data image found.  Generating a model based\n ");
        pfunc ("  on input parameters, with exposure time of 1 second.\n\n");
	input.create = 1;
        sprintf (input.data, "none");
        blank_img (&data, input);
    } else if (returnval == 2) {
	pfunc ("-- WARNING: The exposure time header is missing.  Default to 1 second.\n\n");
    } else if (returnval == 3) {
	pfunc ("-- WARNING: The exposure time is zero seconds.  Default to 1 second.\n\n");
    };

    if (readfits (&psf) == 1) {
	pfunc ("-- NO CONVOLUTION DONE: PSF image not found. \n"); 
        sprintf (psf.name, "none");
    } else {
	/*  Make sure the PSF naxes are even and flux normalized to unity.  */
	psfcheck (&psf);
	kernel = getkernel (input.kernel);

	if (kernel ==  NULL) {
            sprintf (input.kernel, "");
	    pfunc ("-- No CCD charge diffusion kernel found or applied.\n");
        };
	
        /* Keep track of PSF size to figure out minimal padding later  */
        cpar.psfsz[1] = psf.naxes[1];
        cpar.psfsz[2] = psf.naxes[2];
	fwhmpsf = psf_fwhm (psf);
    };

    /*************************************************************\
    *  Look for a badpix mask image.  Even if none found, create  *
    *  an array of 0s (by badpixlist) anyway for later use as     *
    *  a lookup table to prevent the same pixel from being        *  
    *  used twice.                                                *
    \*************************************************************/
   
    if (readfits (&badpix) == 1 || badpix.naxes[1] != data.naxes[1] ||
				   badpix.naxes[2] != data.naxes[2]) {
	if (badpix.err == 0) {
	    pfunc ("-- WARNING:  The pixel mask is not the same size as the data.\n\n");
	} else {
            badpix.naxes[1] = data.naxes[1];
            badpix.naxes[2] = data.naxes[2];
	    if (badpixlist(&badpix) == 1) {
	        pfunc ("-- No bad pixel mask image used.\n");
	        sprintf (badpix.name, "none");
	    };
	};
    };


    /******************************************************\
    *  Catch problems with the user input, specifically    *
    *  with image region sizes, that might cause the       *
    *  program to crash or cause unanticipated problems.   *
    \******************************************************/

    errorcheck (&input, &data, &badpix, &psf, &sigma);

    /*************************\
    *  Check for sigma image  *
    \*************************/

    if (readfits (&sigma) == 1) {
	pfunc ("-- No sigma image.  Creating one using: GAIN=%.2f, RDNOISE=%.2f, NCOMBINE=%.1f.\n", data.gain, data.rdnoise, data.ncombine);
    };

    /***************************************************************\
    *  For Sersic profile only -- interpolates to get the relation  *
    *                           between the powerlaw index and mu   *
    \***************************************************************/

    getmu (Y2, NINTERP);

    /*****************************\
    *  Now the real fun begins!!  *
    \*****************************/

    search_min (&data, &sigma, &psf, &badpix, &fpar, &constr, &cpar); 

    /************************************************\
    *  DONE!  Now free up all the input matrices.    *
    \************************************************/

    if (sigma.err == 0)
        free_matrix (sigma.z, 1, sigma.naxes[2], 1, sigma.naxes[1]);
    if ( strncmp (psf.name, "none", 4) != 0)
        free_matrix (psf.z, 1, psf.naxes[2], 1, psf.naxes[1]);
    if (kernel != NULL) 
	free_matrix (kernel->z, 1, kernel->naxes[2], 1, kernel->naxes[1]);

    free_matrix (badpix.z, 1, badpix.naxes[2], 1, badpix.naxes[1]);

    /* Free fpar */

    if (fpar.next != NULL)
        del_fitpars_struct (fpar.next);

    if (fpar.high != NULL) 
        del_fitpars_struct (fpar.high);
    free ((float *)fpar.a);
    free ((float *)fpar.sig);
    free ((int *)fpar.ia);
    if (fpar.outer != NULL)
	free ((int *)fpar.outer);
    if (fpar.inner != NULL)
	free ((int *)fpar.inner);
    if (constr.next != NULL) 
        del_constr_struct (constr.next);
    if (constr.comp != NULL) 
        free ((int *) constr.comp);
    if (constr.par != NULL) 
        free ((int *) constr.par);
    if (constr.cval != NULL) 
        free ((int *) constr.cval);

    /* Beep when done */
    printf ("%c", 7);

    if (strncmp (device, "regular", 7)!= 0)
        endwin ();

}


void blank_img (struct image *data, struct inpars input)
{
    int xmin, xmax, ymin, ymax;

    sscanf (input.imgsect, "[%d:%d,%d:%d]", &xmin, &xmax, &ymin, &ymax);
    sprintf (data->imgsect, "[%d:%d,%d:%d]", &xmin, &xmax, &ymin, &ymax);
    data->naxes[1] = xmax;
    data->naxes[2] = ymax;
    data->z = matrix (1, data->naxes[2], 1, data->naxes[1]);
}