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


/*****************************************************************************\
*  Extract the Fourier parameter values that the user wants to fit.           *
\*****************************************************************************/

void fourier_pars (char *string, struct fitpars *fpar)
{
    int npars, i, mode, pa, oldnpars;
    char dum[10];
    struct fitpars *f, *ft, *begin;

    sscanf (string, " %c%d)", &dum, &mode);

    /********************************************************\
    *  Check to see if the Fourier array is already active.  *
    *  First, check the fitpar tree for higher order terms   *
    \********************************************************/

    begin = fpar;

    f = fpar->high;
    while (f!= NULL && strncmp (f->objtype, "four", 4) != 0)
	f = f->next;

    if (f == NULL) {         /*  No Fourier terms, create new  */
	ft = (struct fitpars *) malloc ((size_t)(sizeof (struct fitpars)));
	ft->a = (float *) NULL;
	ft->ia = (int *) NULL;
	ft->sig = (float *) NULL;
	ft->inner = (int *) NULL;
	ft->outer = (int *) NULL;
	ft->high = (struct fitpars *) NULL;
	ft->next = (struct fitpars *) NULL;
	ft->compnum = 0;
	ft->npars = 0;
	ft->tr_num = 0;
	ft->normtype = 0;
	ft->outtype = 0;
	oldnpars = 0;
	sprintf (ft->objtype, "four");
        if (fpar->high == NULL)   /*  First in the list  */
	    fpar->high = ft;
        else{     /*  Not the first in the list; must advance downward  */
	    fpar = fpar->high;
	    while (fpar->next != NULL)
		fpar = fpar->next;
	    fpar->next = ft;
	};
	f = ft; 
    } else
	oldnpars = f->npars;

    /********************************************************************\
    *  Note that the number of parameters is always twice the largest    *
    *  mode being fitted, even if some of the modes are not being used.  *
    *  And since mode and PA come in pairs, fpar->npars is always twice  *
    *  the number of modes.                                              *
    \********************************************************************/

    npars = mode * 2;

    if (npars > f->npars) {
        f->a = (float *) realloc (f->a, (size_t)(sizeof (float) * (npars + OFFSET)));
        f->ia = (int *) realloc (f->ia, (size_t)(sizeof (int) * (npars + OFFSET)));
        f->sig = (float *) realloc (f->sig, (size_t)(sizeof (float) * (npars + OFFSET)));
        f->npars = npars;

        /*********************************************\
        *  Set unused mode a's to 0 and ia's to -1.   *
        \*********************************************/

        for (i=oldnpars+1; i <= npars; i++) {
	    f->a[i] = 0.;
	    f->ia[i] = -1;
	    f->sig[i] = 0.;
	};
    }

    f->a[0] = 0.;
    f->ia[0] = -1;
    f->sig[0] = 0.;

    mode = 2 * mode - 1;        /*  Mode index         */
    pa = mode + 1;		/*  Phase angle index  */

    /*****************************************************************\
    *  Read input.  Sometimes the user may forget to specify toggle   *
    *  mode.  Set them to 1 by default.                               *
    \*****************************************************************/
    
    f->ia[mode] = 1;
    f->ia[pa] = 1;
    sscanf (string, " %s %f %f %d %d", dum, &f->a[mode], &f->a[pa], 
	&f->ia[mode], &f->ia[pa]);

    if (f->a[mode] == 0. && f->ia[mode] == 1)
	f->a[mode] = 1.e-2;    /*  Initial guess of 1e-2 is for stability. *\
			       *   Also need to do this so galfit doesn't   *
			       \*  crash when the amplitude = 0.           */

    /*****************************************************************\
    *  Lastly, check to see if both a[i] and ia[i] are 0.  If so,     *
    *  that effectively means this parameter is disabled.  We should  *
    *  set ia[i] = -1 so this item is not listed as being fitted.     *
    *  On the other hand, we should reduce f->npars to make sure it   *
    *  accurately reflects the number of highest parameter.           *
    \*****************************************************************/

    oldnpars = 0;
    for (i=1; i <= f->npars; i+=2)  {
        if (f->a[i] == 0. && f->ia[i] == 0) {
            f->ia[i] = -1;
	    f->ia[i+1] = -1;
	}
	if (f->a[i] != 0. || f->ia[i] == 1)
	    oldnpars = i+1;
    };
    if (oldnpars < f->npars) {
        f->npars = oldnpars;
        f->a = (float *) realloc (f->a, (size_t)(sizeof (float) * (oldnpars + OFFSET)));
        f->ia = (int *) realloc (f->ia, (size_t)(sizeof (int) * (oldnpars + OFFSET)));
        f->sig = (float *) realloc (f->sig, (size_t)(sizeof (float) * (oldnpars + OFFSET)));
    };

    /*****************************************************\
    *  If everything's deleted, delete structure as well  *
    \*****************************************************/

    if (f->npars == 0) {
	fpar = begin;
	if (strncmp (fpar->high->objtype, "four", 4) == 0) {
	    free (fpar->high);
	    fpar->high = NULL;
	} else {
	    fpar = fpar->high;
	    while (strncmp (fpar->next->objtype, "four", 4) != 0) 
		fpar = fpar->next;

	    free (fpar->next);
	    fpar->next = NULL;
	};
    };
}