#include <stdio.h>
#include <math.h>
#include <string.h>
#include <strings.h>
#include "structs.h"

/* #define CSQ(i) (sqrt(chi2nu * covar[(i)][(i)]))  */


/**********************************************************************\
*  The parameter uncertainty SIG(i) is scaled by chi2/nu.              *
\**********************************************************************/

/* #define SIG(i) (sqrt(chi2nu) * sig[i]) */

#define SIG(i) (sig[i])      /*  sqrt(chi2nu) is modified in 'assign_err.c'  */

void get_fparc (char *par, float a, int ia, int natspc);
void get_eparc (char *par, float a, int ia, int natspc);
int istrunc (char *objtype);

void fprintpar (int lowx, int lowy, struct fitpars *fpar, FILE *ftype, 
	char newout[], double *sig, float *coerr, double chisq, int ndof)
{
    extern struct inpars input;
    extern float xskycent, yskycent;
    struct fitpars *h;
    float x, y, chi2nu;
    char *cptr, xc[80], yc[80], parc[80], par2c[80], sigx[80], sigy[80];
    int noff = 0, noff2, c, i, j, ibreak, natspc;

    chi2nu = chisq / ndof;

    fprintf (ftype, "-----------------------------------------------------------------------------\n\n");
    fprintf (ftype, "Input image     : %s%s \n", input.data, input.imgsect);
    fprintf (ftype, "Init. par. file : %s \n", input.initparfile);
    fprintf (ftype, "Restart file    : %s \n", newout);
    fprintf (ftype, "Output image    : %s \n\n", input.output);

    while (fpar != NULL) {
        x = fpar->a[1] + lowx - 1.;
	y = fpar->a[2] + lowy - 1.;
        
        get_fparc (xc, x, fpar->ia[1], 8);
        get_fparc (yc, y, fpar->ia[2], 8);

        /*------------------------------------\
        |  Most radial profile functions....  |
        \------------------------------------*/

        if (strncmp(fpar->objtype, "sersic", 6)==0 || 
            strncmp(fpar->objtype, "moffat", 6)==0 ||
	    strncmp(fpar->objtype, "devauc", 6)==0 ||
            strncmp(fpar->objtype, "expdisk", 7)==0 ||
	    strncmp(fpar->objtype, "edgedisk", 8)==0 ||
	    strncmp(fpar->objtype, "gaussian", 8)==0 ||
	    strncmp(fpar->objtype, "psf", 3)==0 ||
	    strncmp(fpar->objtype, "nuker", 5)==0 ||
	    strncmp(fpar->objtype, "king", 4)==0 || 
	    strncmp(fpar->objtype, "ferrer", 6)==0) {

	    if (strncmp(fpar->objtype, "nuker", 5)==0)
		natspc = 5;
	    else if (strncmp(fpar->objtype, "king", 4)==0 || 
		strncmp(fpar->objtype, "ferrer", 6)==0)
		natspc = 6;
	    else
		natspc = 7;

            fprintf (ftype, " %-10s: (%s, %s)", fpar->objtype, xc, yc);
	    for (i=3; i <= 10; i++) {
		if (fpar->ia[i] > -1) {
                    get_fparc (parc, fpar->a[i], fpar->ia[i], natspc);
                    fprintf (ftype, " %s", parc);
		};
	    };
	    fprintf (ftype, "\n");

	    natspc = 7;
	    get_fparc (xc, SIG(1+noff), fpar->ia[1], natspc-1);
	    get_fparc (yc, SIG(2+noff), fpar->ia[2], natspc-1);

            fprintf (ftype, " %-13s (%s,   %s)", "", xc, yc);

	    if (strncmp(fpar->objtype, "nuker", 5)==0)
		natspc = 5;
	    else if (strncmp(fpar->objtype, "king", 4)==0 ||
		strncmp(fpar->objtype, "ferrer", 6)==0)
		natspc = 6;

	    for (i=3; i <= 10; i++) {
		if (fpar->ia[i] > -1) {
                    get_fparc (parc, SIG(i+noff), fpar->ia[i], natspc);
                    fprintf (ftype, " %s", parc);
		};
	    };
	    fprintf (ftype, "\n");
	};

        /*------\
        |  Sky  |
        \------*/

	if (strncmp(fpar->objtype, "sky", 3) ==0) {

	    x = xskycent + lowx - 1;
	    y = yskycent + lowy - 1;

	    fprintf (ftype, " sky       : [%6.2f, %6.2f]", x, y);

            for (i=1; i <= 3; i++) {
                if (fabs(fpar->a[i]) > 1.)
                    get_fparc (parc, fpar->a[i], fpar->ia[i], 7);

	    	if (fabs(fpar->a[i]) < 1.)
                    get_eparc (parc, fpar->a[i], fpar->ia[i], 7);

		fprintf (ftype, "  %s", parc);
            };
            fprintf (ftype, "\n");

	    fprintf (ftype, "                              ");
            for (i=1; i <= 3; i++) {
                if (fabs(fpar->a[i]) > 1.)
                    get_fparc (parc, SIG(i+noff), fpar->ia[i], 6);

	    	if (fabs(fpar->a[i]) < 1.)
                    get_eparc (parc, SIG(i+noff), fpar->ia[i], 6);

		fprintf (ftype, "  %s", parc);
            };
            fprintf (ftype, "\n");
	};

        /*-----------------------\
        |  Truncation functions  |
        \-----------------------*/

        if (istrunc(fpar->objtype)) {

            if ((cptr = index (fpar->objtype, '2')) != NULL) {
                if (rindex (cptr, 'b') != NULL)
                    fprintf (ftype, " %s : ", fpar->objtype);
                else
                    fprintf (ftype, " %s   : ", fpar->objtype);
            } else if ((cptr = index (fpar->objtype, '-')) != NULL)
                fprintf (ftype, " %s  : ", fpar->objtype);
            else
                fprintf (ftype, " %s   : ", fpar->objtype);

            if (fpar->inner[0] != 0) {
                fprintf (ftype, "[%d", fpar->inner[1]);

                for (i=2; i <= fpar->inner[0]; i++)
                    fprintf (ftype, ", %d", fpar->inner[i]);

                fprintf (ftype, "] ");
            } else
                fprintf (ftype, "[0] ");

            if (fpar->outer[0] != 0) {
                fprintf (ftype, "[%d", fpar->outer[1]);

                for (i=2; i <= fpar->outer[0]; i++)
                    fprintf (ftype, ", %d", fpar->outer[i]);

                fprintf (ftype, "] ");
            } else
                fprintf (ftype, "[0] ");

	    /**********************************\
	    *  (x, y), q, PA free parameters?  *
	    \**********************************/

            if (fpar->ia[1] != -1 || fpar->ia[2] != -1)
                fprintf (ftype, " (%s, %s)", xc, yc);

            for (i=3; i <= 10; i++) {
                if (fpar->ia[i] > -1) {
                    get_fparc (parc, fpar->a[i], fpar->ia[i], 7);
                    fprintf (ftype, " %s", parc);
                };
            };
            fprintf (ftype, "\n");

            if (fpar->ia[1] != -1 || fpar->ia[2] != -1) {
                get_fparc (xc, SIG(1+noff), fpar->ia[1], 6);
                get_fparc (yc, SIG(2+noff), fpar->ia[2], 6);
                fprintf (ftype, " %-22s (%s,   %s)", "", xc, yc);
	    };

	    if (fpar->ia[1] == -1 || fpar->ia[2] == -1)
                fprintf (ftype, "                     ");

            for (i=3; i <= 10; i++) {
                if (fpar->ia[i] > -1) {
                    get_fparc (parc, SIG(i+noff), fpar->ia[i], 7);
                    fprintf (ftype, " %s", parc);
                };
            };
            fprintf (ftype, "\n");
        };

	noff += fpar->npars;   /*  Advance counter  */

        /*------------------------\
        |  High order parameters  |
        \------------------------*/

        if (fpar->high != NULL) {
            h = fpar->high;

	    noff2 = noff;

            while (h != NULL) {
                c = 0;

                /*---------------------\
                |  Rotation functions  |
                \---------------------*/

                if (strncmp (h->objtype, "power", 5) == 0 ||
		    strncmp (h->objtype, "log", 3) == 0 ) {
                    fprintf (ftype, "   %6s  :", h->objtype);

                    for (i=1; i <= 10; i++) {
                        if (h->ia[i] > -1) {
                            get_fparc (parc, h->a[i], h->ia[i], 7);
                            fprintf (ftype, " %s", parc);
			};
		    };
		    fprintf (ftype, "\n");

                    fprintf (ftype, "           :");
                    for (i=1; i <= 10; i++) {
                        if (h->ia[i] > -1) {
                            get_fparc (parc, SIG(i+noff2), h->ia[i], 7);
                            fprintf (ftype, " %s", parc);
			};
		    };
		    fprintf (ftype, "\n");
		};

                /*----------------\
                |  Bending modes  |
                \----------------*/

                if (strncmp (h->objtype, "bend", 4) == 0) {
		    ibreak = 1;
                    for (i=1; i <= h->npars; i++) {
                        if (!(h->a[i] == 0 && h->ia[i] == 0) &&
                            h->ia[i] !=-1) {

                            if (c == 0) {
                                fprintf (ftype, "   bend    : ");
			    };

                            get_eparc (parc, h->a[i], h->ia[i], 8);
                            fprintf (ftype, "(%d: %s)  ", i, parc);
                                    
                            c++;
                        };

                        if (c == 4) {
                            c = 0;

                            fprintf (ftype, "\n");
                            fprintf (ftype, "             ");
			    for (j=ibreak; j <= i; j++) {
                        	if (!(h->a[j] == 0 && h->ia[j] == 0) &&
                            	    h->ia[j] !=-1) { 

                            	    get_fparc (parc, SIG(j+noff2), h->ia[j], 8);
                                    fprintf (ftype, "(%d: %s)  ", j, parc);
				};
			    };

                            fprintf (ftype, "\n");
			    ibreak = i+1; 
                        };

                    };
                    if (c%4 != 0) {
                        fprintf (ftype, "\n");

                        fprintf (ftype, "             ");
			for (j=ibreak; j < i; j++) {
                            if (!(h->a[j] == 0 && h->ia[j] == 0) &&
                                h->ia[j] !=-1) {

                            	get_fparc (parc, SIG(j+noff2), h->ia[j], 8);
                                fprintf (ftype, "(%d:%s)  ", j, parc);
			    };
			};

                        fprintf (ftype, "\n");

		    }
                };

                /*----------------\
                |  Fourier modes  |
                \----------------*/

                if (strncmp (h->objtype, "four", 4) == 0) {
		    ibreak = 1;
                    for (i=1; i <= h->npars; i=i+2) {
                        if (!(h->a[i] == 0 && h->ia[i] == 0) &&
                            h->ia[i] !=-1) {

                            if (c == 0) {
                                fprintf (ftype, "   fourier : ");
			    };

                            get_fparc (parc, h->a[i], h->ia[i], 6);
                            get_fparc (par2c, h->a[i+1], h->ia[i+1], 7);

                            fprintf (ftype, "(%d:%s, %s)  ", i/2+1, 
				parc, par2c);
                                    
                            c++;
                        };

                        if (c == 3) {
                            c = 0;

                            fprintf (ftype, "\n");
                            fprintf (ftype, "             ");
			    for (j=ibreak; j <= i; j+=2) {
                        	if (!(h->a[j] == 0 && h->ia[j] == 0) &&
                            	    h->ia[j] !=-1) {

                            	    get_fparc (parc, SIG(j+noff2), h->ia[j], 6);
                            	    get_fparc (par2c, SIG(j+1+noff2), 
					h->ia[j+1], 7);

                                    fprintf (ftype, "(%d:%s, %s)  ", 
                                        j/2+1, parc, par2c);
				};
			    };

                            fprintf (ftype, "\n");
			    ibreak = i+2;
                        };


                    };
                    if (c%3 != 0) {
                        fprintf (ftype, "\n");

                        fprintf (ftype, "             ");
			for (j=ibreak; j < i; j+=2) {
                            if (!(h->a[j] == 0 && h->ia[j] == 0) &&
                                h->ia[j] !=-1) {

                            	get_fparc (parc, SIG(j+noff2), h->ia[j], 6);
                            	get_fparc (par2c, SIG(j+1+noff2), 
				    h->ia[j+1], 7);

                                fprintf (ftype, "(%d:%s, %s)  ", j/2+1,
                                    parc, par2c);
			    };
			};

                        fprintf (ftype, "\n");

		    }
                };

                /*---------------------\
                |  Diskyness/boxyness  |
                \---------------------*/

                if (strncmp (h->objtype, "C0", 2) == 0 ) {
                    get_fparc (parc, h->a[1], h->ia[1], 6);
                    fprintf (ftype, "     c0    :  %s \n", parc);

                    get_fparc (parc, SIG(1+noff2), h->ia[1], 6);
		    fprintf (ftype, "           :  %s \n", parc);
		};

		if (!istrunc(h->objtype))
		    noff += h->npars;
                h = h->next;
            };
        };

	fpar = fpar->next;
    };

    fprintf (ftype, " Chi^2 = %.5f,  ndof = %d\n", chisq, ndof);
    fprintf (ftype, " Chi^2/nu = %.3f \n\n", chisq / ndof);
    fprintf (ftype, "-----------------------------------------------------------------------------\n");

    return;
}