/****************************************************************************\
*  Create a structured list of derivative images, one image for each         *
*  parameter that needs to be fitted.                                        *
*                                                                            *
*  Note that each of the subcomponents involved in truncation has structs    *
*  in fpar->high which mirror all the main truncation parameters.  However   *
*  *each* mirrored parameter, though it only refers back to the truncation   *
*  and thus not an independent parameter, will nevertheless come with a      *
*  derivative image, stored in the df->high struct of the light model.       *
*  These individual derivative images will be summed together to create the  *
*  overal image corresponding to each truncation parameter.  This is done    *
*  in mkmodel, by the subroutine "trunc_links."                              *
*                                                                            *
*  BUT, because the truncation profile can either have its own centroid or   *
*  have the centroid tied to the objects, implicitly, we have to remember    *
*  the centroid position of the parent light profile if there's an implicit  *
*  centroid linking, when creating an array for convolution....  If not      *
*  convolving, then it doesn't matter.                                       *
*                                                                            *
*  convolve = create arrays for convolution region.  The size will depend    *
*             on where an object center is relative to the entire fitted     *
*             region.  If part of the convolution box extends outside of     *
*             the fitting region, that part will be vignetted.  Setting      *
*             convolve = 1 will create derivative images for convolution     *
*             purposes.                                                      *
\****************************************************************************/

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

void convregion (struct derivs *df, float a[], struct image *d, 
    struct convpars *cpar);

struct derivs *create_work_arrays (struct fitpars *fpar, struct image *model,
    int convolve, struct convpars *cpar)
{

    struct derivs *new_obj_workarrays (struct fitpars *, struct image *model,
        int convolve, struct convpars *cpar);

    int i;
    struct derivs *df, *head;

    df = new_obj_workarrays (fpar, model, convolve, cpar);
    head = df;

    while (fpar->next != NULL) {
	fpar = fpar->next;
	df->next = new_obj_workarrays (fpar, model, convolve, cpar);
	df = df->next;
    };

    df->next = NULL;
    df = head;
    return (df);
}

/*****************************************************************************\
*  Create an entire work array for a single object, consisting of derivative  *
*  image space for both the classical parameters and high order parameters    *
\*****************************************************************************/

struct derivs *new_obj_workarrays (struct fitpars *fpar, struct image *model,
    int convolution, struct convpars *cpar)
{
    struct derivs *df;
    static int high = 0;
    static float *a;   /*  Need to keep track of centroid of a component    *\
                       \*  even when advancing through the high order terms */
    int i;

    if (high == 0) 
	a = fpar->a;

    /*----------------------------------------\
    |  Create arrays for this (sub)structure  |
    \----------------------------------------*/

    df = (struct derivs *) malloc ((size_t) (sizeof (struct derivs)));
    df->high = NULL;

    df->dpm = (float ***) malloc ((size_t)
	      (sizeof (float **) * (fpar->npars+1)));

    for (i = 0; i <= fpar->npars; i++) {
	if (!convolution) {

            if (i==0 || fpar->ia[i] > 0) {   
                /*  Allocate space for the entire fitting region  */

	        df->naxes[1] = model->naxes[1];
	        df->naxes[2] = model->naxes[2];
	        df->dpm[i] = matrix (1, df->naxes[2], 1, df->naxes[1]);
	    };
	} else if (convolution) {
            if (strncmp(fpar->objtype, "sky", 3)!=0) {
                if (i==0 || (i>0 && fpar->ia[i] > 0)) {

                    /*  Create space for convolution region only  */

                    convregion (df, a, model, cpar);
                    if (df->naxes[2] > 0 && df->naxes[1] > 0)
                        df->dpm[i] = matrix(1, df->naxes[2], 1, df->naxes[1]);
                };
            } else {
                df->naxes[1] = 0;
                df->naxes[2] = 0;
            };
	};
    };

    /*-------------------------------------------------------------------\
    |  If we're in the high order structured list, go down that list by  |
    |  recursively calling this subroutine.                              |
    \-------------------------------------------------------------------*/

    if (high >= 1 && fpar->next != NULL) 
        df->next = new_obj_workarrays (fpar->next, model, convolution, cpar);
    else if (high >= 1 && fpar->next == NULL)
	df->next = NULL;                 /*  End of high order list  */

    /*-------------------------------------------------------------\
    |  Now check to see if need to go into the higher order array  |
    \-------------------------------------------------------------*/

    if (fpar->high != NULL) {
	high += 1;
	df->high = new_obj_workarrays (fpar->high, model, convolution, cpar);
        high -= 1;
    };

    return (df);

}