#include <stdio.h>
#include <string.h>
#include "nrutil.h"
#include "structs.h"

#define NEW 0
#define ADD 1

int istrunc (char *objtype);
void writefits(char *phead, struct image *img, char *object, int append);
void writeheader (struct image *img, struct fitpars *fpar, double chisq,
    int ndof, int nfree, int nfixed, int hdu, int offx, int offy);
void subimg (struct image *data, struct image *model, struct image *sub);
void mkmodel (struct image *model, struct derivs *df, struct fitpars *fpar,
    struct convpars *cpar, struct cons *constr, int output);
double chi2calc (int nfree, float **data, float **sigma, float **model, 
    long naxes[], int *ndof);

void outmodel (struct fitpars *fpar, struct convpars *cpar, 
    struct image *d, struct image *model, struct image *sigma, struct derivs *df, 
    double *chisq, int *ndof, int nfree, int nfixed, int offx, int offy)
{
    struct image sub;
    extern struct inpars input;
    char obj[25], infile[80];
    extern unsigned long *pos, *cpos;
    struct fitpars *tfpar, *current;
    int hdu;

    /****************************************************\
    *  Output original image into the output image block *
    \****************************************************/

    /* Output original image into the output image block */

    if (input.create != 1 && input.create != 3) {
        strcpy (infile, d->name);
        strcpy (d->name, input.output);
        writefits (infile, d, "original image", NEW);
    };

    /**********************************************************\
    *  Create and output final models into 2nd image in block  *
    \**********************************************************/

    sprintf (obj, "model");
    strcpy (model->name, input.output);
    mkmodel (model, df, fpar, cpar, (struct cons *)NULL, 1);

    if (input.create != 1 && input.create != 3) {
	if (input.create == 2)
            *chisq = chi2calc (nfree, d->z, sigma->z, model->z, model->naxes, 
		    ndof);
        writefits (infile, model, obj, ADD);
	writeheader (model, fpar, *chisq, *ndof, nfree, nfixed, 3, offx, offy);
    };

    /****************************************************\
    *  Output residual image into third image in block   *
    \****************************************************/

    if (input.create != 1 && input.create != 3) {
        sub = *d;
	sub.z = matrix (1, d->naxes[2], 1, d->naxes[1]);
        subimg (d, model, &sub);
        sprintf (obj, "residual map");
        writefits (infile, &sub, obj, ADD);
        free_matrix(sub.z, 1, d->naxes[2], 1, d->naxes[1]);
    } else if (input.create == 1) {
        writefits (infile, model, obj, NEW);
	writeheader (model, fpar, *chisq, *ndof, nfree, nfixed, 1, offx, offy); 
    };

    if (input.create == 3) {

        strcpy (infile, "subcomps.fits");
        strcpy (model->name, "subcomps.fits");
        strcpy (d->name, "subcomps.fits");
        writefits (infile, d, "original image", NEW);

	hdu = 3;
	current = fpar;
	while (current != NULL) {
	    if (!istrunc (current->objtype)) {

	        tfpar = fpar;
	        while (tfpar != NULL) {
	            if (tfpar->compnum == current->compnum)
		        tfpar->outtype = 0;
		    else
		        tfpar->outtype = 1;

	            tfpar = tfpar->next;
	        };

                sprintf (obj, current->objtype);
                mkmodel (model, df, fpar, cpar, (struct cons *)NULL, 1);
                writefits (infile, model, obj, ADD);
                writeheader (model, fpar, *chisq, *ndof, nfree, nfixed, hdu, 
		    offx, offy); 
	        hdu++;
	    };

	    current = current->next;
	};


	/*--------------------------------------------\
        |  Now that we're done, set all outtype to 0  |
	\--------------------------------------------*/

	current = fpar;
	while (current != NULL) {
	    current->outtype = 0;
	    current = current->next;
	};
    };

}