/***************************************************************\
*  This subroutine makes calls functions that make the models   *
*  and convolve the models with the PSF.                        *
*                                                               *
*  The model making happens in two steps.  First, models will   *
*  be created centered on the values specified in fpar->a[1]    *
*  and fpar->a[2] for all objects.  This model will not be      *
*  used in convolution because doing so would broaden out the   *
*  objects at the center.  While this is ok for galaxies it's   *
*  rather bad for fitting point sources.  Therefore, a second   *
*  model is created which is centered precisely on the pixel    *
*  center.                                                      *
*                                                               *
*  Note that regions inside "df->imgsect" will be convolved by  *
*  the PSF, but regions outside it will not.  First, to do      *
*  this correctly, one has to worry about how border padding    *
*  will influence the convolution within "df->imgsect".         *
*                                                               *
*  Specifically, since we're using FFT to do convolution, we    *
*  need to extend the image so that each side is 2^N in         *
*  length and width.  If we pad the extension with 0, then      *
*  around the edges of the model, convolution will effectively  *
*  convolve the model flux values with 0, hence corrupting the  *
*  model there.  But, instead of padding with 0, we pad a       *
*  square annulus around that region with real model data,      *
*  then only the region inside the annulus will be corrupted,   *
*  because the thickness of the annulus will be exactly half    *
*  the width and length of the convolution PSF.  The corrupted  *
*  region (annulus) will then be thrown away.                   *
*                                                               *
*  Therefore we devise the following scheme:                    *
*    1) The model inside the region "df->imgsect" will          *
*       be created in a separate array.                         *
*    2) That array will be extended on all sides by 1/2 the     *
*       length/width of the PSF, and that region will be        *
*       filled with model values.                               *
*    3) The rest of the image will be extended so that each     *
*       side has 2^N pixels in length and width.  This region   *
*       can be filled with 0 (or whatever, actually) padding.   *
*    4) Convolution takes place in this 2-D array.              *
*    5) Once convolution is done, copy this array back into     *
*       the final model image, being careful to leave out the   *
*       PSF extended annulus and the 0 padding.                 *
*                                                               *
\***************************************************************/

#include <math.h>
#include <string.h>
#include "mymath.h"
#include "structs.h"
#include "nrutil.h"
#include "debug.h"

void realign_trunc_mirrors (struct fitpars *fpar);
void realign_compnum (int compnum, struct fitpars *fpar);
struct fitpars *dup_fitpars_struct (struct fitpars *fpar);
struct derivs *create_work_arrays (struct fitpars *fpar, struct image *model,
    int convolve, struct convpars *cpar);
void point_fpar_at_df (struct fitpars *fpar, struct derivs *df);
void trunc_cent (struct fitpars *mod, struct fitpars *orig, float *refpos,
    int os);
int istrunc (char *objtype);
void oversamp (struct fitpars *fpar, int os);
void strong_links (struct fitpars *fpar, struct derivs *df, 
    struct cons *constr);
void assemble_trunc_df (struct fitpars *fpar, struct derivs *df);
void delete_work_arrays (struct fitpars *fpar, struct derivs *df);
void del_fitpars_struct (struct fitpars *fpar);
void objselect (float a[], int ia[], int normtype, struct fitpars *high, 
    char objtype[], struct image *model, struct derivs *df);
double **get_fftpsf (struct fitpars *fpar, struct image model, struct image psf,
    long complex_naxes[]);
void do_conv (struct fitpars *fpar, struct image model, double **psfft,
    long complex_naxes[], struct derivs *df, struct derivs *cdf,
    struct convpars *cpar, int output);
void conv_highpars (struct fitpars *fpar, struct image model, double **psfft,
    long complex_naxes[], struct derivs *df, struct derivs *cdf,
    struct convpars *cpar, int output);
void clearmodel (struct image *model, struct derivs *df, struct fitpars *fpar);
void delete_derivs (struct derivs *df);
int istrunc (char *objtype);

void mkmodel (struct image *model, struct derivs *df, struct fitpars *fpar, 
    struct convpars *cpar, struct cons *constr, int output) 
{

    void set_ia_zero (struct fitpars *fpar);
    void copy_model (struct fitpars *fpar, struct image *model, 
        struct derivs *df);
    extern struct image psf;
    extern struct inpars input;

    float xorig, yorig;
    unsigned long ndata;
    struct image tmodel;
    struct derivs *cdf, *tdfptr, *dfptr, *dfhigh, *cdfhigh;
    struct fitpars *fptr, *cfpar, *cfptr, *high, *hfptr;
    int i, j, xmin, xmax, ymin, ymax, os;

    long complex_naxes[3];
    double **psfft;

    /*----------------------------------------------------------------\
    |  After each iteration, if there's a truncation function, its    |
    |  parameters in the base structure will be misaligned with the   |
    |  mirror parameters in the fpar->high structures.  Therefore,    |
    |  we need to realign them.                                       |
    \----------------------------------------------------------------*/

    realign_trunc_mirrors (fpar);

    /*----------------------------------------------------------------\
    |  "cfpar" is a clone of fpar, created because the convolution    |
    |  process has to modify the position, and sometimes the sizes    |
    |  of objects if the PSF is oversampled compared to the data.     |
    |  The clone also critical for the truncation functions, because  |
    |  the centroid of a specific truncation function may change      |
    |  from object to object.  So "cfpar" is modified in lieu of the  |
    |  original fpar structure.                                       |
    \----------------------------------------------------------------*/

    cfpar = dup_fitpars_struct (fpar);
    point_fpar_at_df (cfpar, df); 
    clearmodel (model, df, fpar);

    /*--------------------------------------------------------------\
    |  This sets all the ia[]'s to 0 so that the derivative images  |
    |  are not produced.  Saves time when we just want to create    |
    |  an image block.                                              |
    \--------------------------------------------------------------*/

    if (output) 
	set_ia_zero (cfpar);

    /*----------------------------------------------------------------\
    |  Create models for the entire fitting region.  The models       |
    |  created here are centered on the parameter centroids in        |
    |  fpar->a[1] and fpar->a[2], as opposed to exactly on the pixel  |
    |  center (below, in the next block of model generation).  This   |
    |  two step process is necessary so we don't broaden out the      |
    |  model profiles unnecessarily at the center-most pixel.         |
    \----------------------------------------------------------------*/

    fptr = fpar;
    dfptr = df;
    cfptr = cfpar;

    os = input.sampfac;

    while (cfptr != NULL) {   

	/*-----------------------------------------------------------\
        |  Note that fptr is used in trunc_cent because trunc_cent   |
	|  figures out from the original fpar what position to       |
        |  assign the truncation function and stores it into cfpar   |
        |  so as not to mess with the original fpar pointer.  The    |
        |  benefit of doing so is more obvious below when            |
        |  convolution is involved.                                  |
	\-----------------------------------------------------------*/

	trunc_cent (cfptr, fptr, fptr->a, 1);

	if (output == 0 || (output == 1 && cfptr->outtype == 0)) {

	    if (strncmp(cfptr->objtype, "psf", 3)!=0)
	        input.sampfac = 1;

            objselect (cfptr->a, cfptr->ia, cfptr->normtype, cfptr->high, 
		cfptr->objtype, model, dfptr);
	    input.sampfac = os;
	};

	cfptr = cfptr->next;
	dfptr = dfptr->next;
	fptr = fptr->next;
    };


    /*----------------------------------------------------------------\
    |  Now create a temporary model image of *just* the convolution   |
    |  region to be used in convolution.  The temporary model is      |
    |  created exactly on pixel center, as opposed to the models      |
    |  that are not used in convolution (above), which are centered   |
    |  on actual object centers.  The fractional pixel shift in the   |
    |  convolved model comes about by convolving the temporary        |
    |  model with a shifted PSF.  The convolution region in the       |
    |  temporary model is then copied back to the unconvolved model.  |
    |								      |
    |  Notice the model used that's being convolved is created to     |
    |  be larger than the convolution box, by the size of the PSF,    |
    |  or one pixel larger if the size is an odd number.	      |
    \----------------------------------------------------------------*/

    if ( strncmp (psf.name, "none", 4) != 0) {

        /*-------------------------------------------------------------\
        |  tmodel is just a pointer to df->dp[i] matrices, used for    |
        |  convolution.                                                |
        \-------------------------------------------------------------*/

        tmodel.dp[1] = model->dp[1];
        tmodel.dp[2] = model->dp[2];
        tmodel.muzpt = model->muzpt;
        tmodel.magzpt = model->magzpt;
        tmodel.naxes[1] = 0;
        tmodel.naxes[2] = 0;

        cdf = create_work_arrays (fpar, model, 1, cpar);
	point_fpar_at_df (cfpar, cdf);
        clearmodel (&tmodel, cdf, fpar);

        /*--------------------------------------------------\
        |  Cycle through all the objects and create models  |
        \--------------------------------------------------*/

        fptr = fpar;
	dfptr = df;
        tdfptr = cdf;

	if (os > 1)       /*  Scale sizes by oversample factor  */
            oversamp (cfpar, os);

        cfptr = cfpar;

        while (fptr != NULL) {
            if (strncmp(fptr->objtype, "sky", 3)!=0 && 
		strncmp(fptr->objtype, "psf", 3)!=0 &&
		!istrunc(fptr->objtype)) {

                sscanf (tdfptr->imgsect, "[%d:%d,%d:%d]", &xmin, &xmax, &ymin,
		    &ymax); 

		/*--------------------------------------------------\
		|  Bypass convolution if img size < 1 pixels which  |
		|  means the convolution box is outside the image.  |
		\--------------------------------------------------*/

		if (tdfptr->naxes[1] > 1 && tdfptr->naxes[2] > 1) {
	            strcpy (tmodel.imgsect, tdfptr->imgsect);

                    tmodel.naxes[1] = tdfptr->naxes[1];
                    tmodel.naxes[2] = tdfptr->naxes[2];

		    /*----------------------------------------------------\
		    |  Offset the object center and size when oversampled |
		    \----------------------------------------------------*/

		    xorig = fptr->a[1] - xmin + 1;
		    yorig = fptr->a[2] - ymin + 1;

	            cfptr->a[1] = NINT(xorig*os - 0.5*(os-1)) + NINT (cpar->psfsz[1]/2.);
                    cfptr->a[2] = NINT(yorig*os - 0.5*(os-1)) + NINT (cpar->psfsz[2]/2.);

		    /*-------------------------------------------------------\
		    |  Now, how to handle the truncation function centroid:  |
		    |  The truncation function can also have a centroid.     |
		    |  BUT, the peak position of the light profile model     |
		    |  that's being truncated essentially does NOT change    |
		    |  when truncated.  At least, if the centroid does       |
		    |  change when truncated, it is virtually impossible to  |
		    |  tell what then the offset should be, a priori.        |
		    |  Even then, the offset of the truncation must follow   |
		    |  that of the light profile centroid in a relative      |
		    |  sense, and the convolution should be done in exactly  |
		    |  the same way as the light profile component.          |
		    \-------------------------------------------------------*/

		    trunc_cent (cfptr, fptr, cfptr->a, os);

		    /*-------------------------------------------------------\
                    |  Now create the model, unless it's a truncation.       |
		    |  Truncation parameters are stored in high order        |
		    |  structures instead of the main object structure, and  |
		    |  then added together below in 'trunc_link_df'.  Note   |
		    |  that the derivative images of the truncation 'mirror' |
                    |  parameters are created at the same time as the        |
                    |  regular parameters in objselect.                      |
		    \-------------------------------------------------------*/

                    if (output && fptr->outtype == 0&& tdfptr->naxes[1] != 0 && 
		        tdfptr->naxes[2] != 0)

                        objselect (cfptr->a, cfptr->ia, cfptr->normtype, cfptr->high, 
		            cfptr->objtype, &tmodel, tdfptr);

		    else if (!output && tdfptr->naxes[1] != 0 && 
		        tdfptr->naxes[2] != 0)

                        objselect (cfptr->a, fptr->ia, cfptr->normtype, cfptr->high, 
		            cfptr->objtype, &tmodel, tdfptr);

                    /*-------------------------------------------------------\
                    |  Now do convolution, cycle through all the parameters  |
                    \-------------------------------------------------------*/

		    psfft = get_fftpsf (fptr, tmodel, psf, complex_naxes);
	            do_conv (fptr, tmodel, psfft, complex_naxes, dfptr, 
		        tdfptr, cpar, output);

		    if (fptr->high != NULL)
		        conv_highpars (fptr->high, tmodel, psfft, 
		            complex_naxes, dfptr->high, tdfptr->high, cpar, 
			    output);

                    if (tmodel.naxes[1] > 0 && tmodel.naxes[2] > 0) 
                        free_dmatrix (psfft, 0, complex_naxes[2]-1, 0,
                            complex_naxes[1]*2-1);
		};
	    };

	    fptr = fptr->next;
	    cfptr = cfptr->next; 
	    dfptr = dfptr->next;
	    tdfptr = tdfptr->next;
	};
    };

    /*----------------------------------------------------------------\
    |  Now add all the subcomponents together into one net model      |
    |  image, merging together all the convolution pieces, and        |
    |  taking care of strong constraints.  The derivative images are  |
    |  taken care of by do_conv above.  			      |
    \----------------------------------------------------------------*/

    copy_model (fpar, model, df);
    assemble_trunc_df (fpar, df);
    strong_links (fpar, df, constr);

    if ( strncmp (psf.name, "none", 4) != 0) {
        tmodel.naxes[1] = 0;
        tmodel.naxes[2] = 0;

	/*-----------------------------------------------------------\
	|  Note that in the following, it is important that fpar is  |
        |  used instead of cfpar, because the truncation function    |
	|  in cfpar may have ia[1] and ia[2] modified to 2 (i.e.     |
	|  they are not free parameters.  This would cause galfit    |
	|  to try and delete work arrays that don't exist.           |
	\*----------------------------------------------------------*/

        delete_work_arrays (fpar, cdf);
    };

    del_fitpars_struct (cfpar);

}


/*
    sprintf (tmodel.name, "tmodel.fits");
    writefits ("tmodel.fits", &tmodel, "tmodel", 0);
*/

/****************************************************************************/

void copy_model (struct fitpars *fpar, struct image *model, struct derivs *df)
{
    struct derivs *dptr;
    struct fitpars *fptr;
    int ix, iy, xmin, xmax, ymin, ymax;

    fptr = fpar;
    dptr = df;

    while (fptr != NULL) {
	for (iy = 1; iy <= model->naxes[2]; iy++)
	    for (ix = 1; ix <= model->naxes[1]; ix++)
		    model->z[iy][ix] += dptr->dpm[0][iy][ix];

        fptr = fptr->next;
	dptr = dptr->next;
    };

}

/****************************************************************************\
*  Set all fitting parameters toggle to 0, which bypasses the convolution    *
*  for derivative images, which makes for a faster output.                   *
\****************************************************************************/

void set_ia_zero (struct fitpars *fpar)
{
    int i;
    static int high=0;


    /*  Go to the end of the current list, whether high order or main  */

    if (fpar->next != NULL)
	set_ia_zero (fpar->next);

    for (i=1; i <= fpar->npars; i++)
	if (fpar->ia[i] == 1)
	    fpar->ia[i] = 0;

    /*  If there are higher order terms, go into that list */

    if (fpar->high != NULL && high == 0) {
	high += 1;
	set_ia_zero (fpar->high);
	high -= 1;
    };
}