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

#define OVERSAMPR 5.
#define peaksamp 50
#define dr 0.05 / peaksamp
#define EPS 1e-5

void trunc_pars (char objtype[], int ia[], struct fitpars *fpar,
    struct derivs *df, struct trunc_links *trlinks, struct trunc_val *trval);
struct fitpars *trunc_spiral (struct fitpars *r);
void highpars (struct fitpars *high, struct derivs *df,
    struct fitpars **bending, struct derivs **bdf, struct fitpars **four,
    struct derivs **fourdf, struct fitpars **c0, struct derivs **c0df,
    struct fitpars **r, struct derivs **rdf);
void drdblah (float x, float y, float *a, int *ia, float OcosPA,
    float OsinPA, struct fitpars *b, struct fitpars *c0, struct fitpars *f,
    struct fitpars *r, float *am, float *drd, float *bdrd, float *drdc,
    float *fdrd, float *rdrd);
void high_derivs (int ix, int iy, struct fitpars *h, struct derivs *hdf,
    float dfdr, float *drd, float dfdI, float *dIdfp);
void get_trunc_cons (char *objtype, int normtype, struct trunc_links *trlink,
    struct trunc_val *trval, float a[], float kap, float OcosPA,
    float OsinPA, float *trunc_norm, float *dfnorm);
void trunc_dfdr_calc (struct trunc_val *trval, struct trunc_links *trlink,
    float flux0, float *PP);
void trunc_derivs (int ix, int iy, float am, float flux,
    float dfdr, struct trunc_val *trval, struct trunc_links *tr);
void free_fitvecs (float **drd, struct fitpars **tilt, struct fitpars *b,
        float **bdrd, struct fitpars *f, float **fdrd, float **dIdfp,
        struct fitpars *r, float **rdrd, float **dfdrp, int normtype,
        struct trunc_links *trlinks, struct trunc_val *trval);
void writefits(char *phead, struct image *img, char *object, int add);

void nuker (char objtype[], float a[], int ia[], int normtype,
    struct fitpars *high, struct image *model, struct derivs *df)
{
    extern struct inpars input;

    unsigned long i, ix, iy, totsamp;
    float nukefunc (float ta[], float x, float y);

    float am, af, ta[NPARS+1], OcosPA, OsinPA, magfac, flux, B, D, GBA,
	a3zpt, rta4, rta4ta5, drdc, dfdr, c0v=2., flux0;
    float trunc_norm=1., dfnorm=0., normrad, PP;
    int iatmp[10], *iaptr;
    float atmp[10], *aptr;
    int n, j, os, mode, nsubsamp;
    float *drd, *bdrd, *fdrd, *rdrd, xf, yf, xc, yc, xp, yp, dfdI0, 
	  dfdIP, sumdfdI, step, r0a4, r0a4a5, D0;
    struct fitpars *b=NULL, *f=NULL, *c0=NULL, *r=NULL, *tilt=NULL,
	   *tiltptr=NULL;
    struct derivs *dptr, *bdf=NULL, *fdf=NULL, *c0df=NULL, *rdf=NULL;
    struct trunc_links trlink = {NULL, 1, 1, 1, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, 0., NULL, NULL, NULL, NULL, NULL, NULL},
        *trlptr;
    struct trunc_val trval = {0, 0., 0., 0., 0., 0., 0., 0., 0., 0., NULL},
        *trvptr;

/*  Ib=ta[3], rb=ta[4], alpha=ta[5], beta=ta[6], gamma=ta[7]  */

    os = input.sampfac;
    magfac = log(10.) / -2.5;

    a[5] = FMAX (a[5], 0.01);
    for (i=1; i <= NPARS; i++)
        ta[i] = a[i];

    ta[10] = -(ta[10])/ 180. * PI + PI/2.;  /* "Intuitive" orientation */

    /*---------------------------------------------------------\
    |  Set pointers to high order correction structures, then  |
    |  create arrays to hold dr/dpar values.                   |
    \---------------------------------------------------------*/

    highpars (high, df, &b, &bdf, &f, &fdf, &c0, &c0df, &r, &rdf);
    drd = vector (1, NPARS);

    if (b != NULL)
        bdrd = vector (1, b->npars);

    if (f != NULL)
        fdrd = vector (1, f->npars);

    if (r != NULL)
        rdrd = vector (1, r->npars);

    if (c0 != NULL)
        c0v += c0->a[1];

    /*----------------------------------------------------------------\
    |  Figure out whether we're trying to fit a trunc.  If so, an     |
    |  outer truncation is innout = 1, whereas if there's an inner    |
    |  truncation, it's a innout = 0.  If there's both, then it's     |
    |  Type 3 ttype.  Also set pointers to trunc structures.          |
    |                                                                 |
    |  If the trunc parameters are modifying an inclined spiral       |
    |  light component, we have to remember what the inclination and  |
    |  sky PA are in "tilt", so we can tilt the truncation in         |
    |  the same way.  Otherwise, decoupling the two (truncation and   |
    |  light profile) makes things rather unintuitive for creating    |
    |  rings and such.                                                |
    \----------------------------------------------------------------*/

    trunc_pars (objtype, ia, high, df, &trlink, &trval);

    if (r != NULL)       /*  A spiral component.  */
            tilt = trunc_spiral (r);
    
    /*----------------------------------------------------------------\
    |  Pre-calculate some things we will frequently use later on but  |
    |  which needs to be calculated one time only.                    |
    |								      |
    |  a[3] is the galaxy surface brightness at break radius   	      |
    |  while ta[3] is the number count at that radius.                |
    \----------------------------------------------------------------*/

    GBA = ( ta[7] - ta[6] ) / ta[5];
    B = pow(2, -GBA);
    a3zpt = (ta[3] - model->muzpt)/-2.5;
    
    ta[3] = pow(10., a3zpt);

    /*-----------------------------\
    |  For the azimuthal profiles  |
    \-----------------------------*/

    OcosPA = cos (ta[10]);
    OsinPA = sin (ta[10]);

    /*----------------------------------------------------------------\
    |  Calculate truncation transition parameters and constants for   |
    |  derivative calculations.                                       |
    \----------------------------------------------------------------*/

    get_trunc_cons ("nuker", normtype, &trlink, &trval, ta, 0., OcosPA,
        OsinPA, &trunc_norm, &dfnorm);

    if (normtype == 3) {
        if (trlink.lrtype == 2 && trlink.innout == 0)
            normrad = trlink.a[5];       /*  Outer break radius  */
        else
            normrad = trlink.a[4];       /*  Break radius        */

        r0a4 = normrad/ta[4];
        r0a4a5 = pow(r0a4, ta[5]);
        D0 = 1. + r0a4a5;
    } else
	normrad = ta[4];

    /*------------------------------------------\
    |  Now fill the image and derivative grids  |
    \------------------------------------------*/

    for (iy=1; iy <= df->naxes[2]; iy++) {
        for (ix=1; ix <= df->naxes[1]; ix++) {

            drdblah ((float)ix, (float)iy, ta, ia, OcosPA,
                OsinPA, b, c0, f, r, &am, drd, bdrd, &drdc, fdrd, rdrd);

            /*-----------------------------------------------------\
            |  If we're near the center of the galaxy we'll        |
            |  subdivide the sampling.                             |
            \-----------------------------------------------------*/

            if (fabs(am) < OVERSAMPR) {
                if (ta[4] * ta[9] <= 1. && fabs(am) <= os) {

                    /* If scalelength is small, oversample by a lot more.  */

                    nsubsamp = IMIN(100, (int)(1./fabs(ta[4]*ta[9])*2*
                               OVERSAMPR/os));
                } else  if (fabs(am) <= 4. * os)
                    nsubsamp = 2 * OVERSAMPR / os;
                else {
                    /* If scalelength is large, bigger step */

                    nsubsamp = IMIN(100, NINT(2*OVERSAMPR/fabs(am)/os));
                };
            } else
                nsubsamp = 1;

            nsubsamp = IMAX (1, nsubsamp);
            totsamp = nsubsamp * nsubsamp;

            xf = ix - 0.5 + 1./(2 * nsubsamp);
            yf = iy - 0.5 + 1./(2 * nsubsamp);
            step = 1./nsubsamp;

	    sumdfdI = 0.;

            /*---------------------------\
            |  Now loop for subsampling  |
            \---------------------------*/

            xc = 0.;
            while (xc <= nsubsamp - 1) {
                yc = 0.;
                xp = xf + xc * step;
                while (yc <= nsubsamp - 1) {
                    yp = yf + yc * step;

                    if (nsubsamp > 1)
                        drdblah (xp, yp, ta, ia, OcosPA, OsinPA, b, c0,
                            f, r, &am, drd, bdrd, &drdc, fdrd, rdrd);

                    /*-------------------------------------------------------\
                    |  If the truncation has its own PA and axis ratio,      |
                    |  calculate the truncation shape separately from        |
                    |  the light profile model.  Note that if the light      |
                    |  profile is a tilted spiral component, the truncation  |
                    |  is also tilted in the same way, given by parameters   |
                    |  stored in "tilt."                                     |
                    \-------------------------------------------------------*/

                    trlptr = &trlink;
                    trvptr = &trval;
                    while (trlptr != NULL && trlptr->objtype != NULL) {
                        if (trlptr->tilt == 1)
                            tiltptr = tilt;
                        else
                            tiltptr = NULL;

                        /*--------------------------------------------------\
                        |  If the truncation type is "radial2", i.e. inner  |
                        |  and outer break radius, and if it's an inner     |
                        |  truncation, then the radius scale paramter is    |
                        |  the outer break radius.                          |
                        \--------------------------------------------------*/

                        if (trlptr->innout == 0 && trlptr->lrtype == 2) {

                            for (i=1; i<=10; i++) {
                                atmp[i] = trlptr->a[i];
                                iatmp[i] = trlptr->ia[i];
                            };
                            atmp[4] = trlptr->a[5];
                            atmp[5] = trlptr->a[4];
                            iatmp[4] = trlptr->ia[5];
                            iatmp[5] = trlptr->ia[4];
                            aptr = atmp;
                            iaptr = iatmp;

                        } else {
                            aptr = trlptr->a;
                            iaptr = trlptr->ia;
                        };

                        drdblah (xp, yp, aptr, iaptr,
                            trvptr->cosPA, trvptr->sinPA,
                            trlptr->b, trlptr->c0, trlptr->f, tiltptr,
                            &trvptr->am, trlptr->drd, trlptr->bdrd,
                            &trlptr->drdc, trlptr->fdrd, NULL);

                        trlptr = trlptr->next;
                        trvptr = trvptr->next;
                    };

                    /*-----------------------\
                    |  Classical parameters  |
                    \-----------------------*/

            	    rta4 = FMAX(am/ta[4], EPS);
	    	    rta4ta5 = pow(rta4, ta[5]);
	    	    D = 1. + rta4ta5;

		    dfdI0 = B * pow(rta4, -ta[7]) * pow(D, GBA) / totsamp / 
			    trunc_norm;
	            flux0= ta[3] * dfdI0;
                    dfdr = flux0 / (am+EPS) * (-ta[7] + GBA / D * rta4ta5 * 
			   ta[5]);

                    if (trlink.objtype != NULL && normtype > 0) {
                        PP = 1;
                        trunc_dfdr_calc (&trval, &trlink, flux0, &PP);
                        dfdI0 = dfdI0 * PP;
                        dfdr = dfdr * PP;
                    };

                    /*--------------------------------------------------------\
                    |  Now put it all together with the normal parameters.    |
                    |  Truncation P included here                             |
                    \--------------------------------------------------------*/

                    flux = ta[3] * dfdI0;

                    /*--------------------------------------------------------\
                    |  What does dfdr operate on when there's a truncation?   |
                    |  See "sersic.c" for more details.                       |
                    \--------------------------------------------------------*/

            	    if (ia[1] > 0) {
			df->dpm[1][iy][ix] += dfdr * drd[1];

                        trlptr = &trlink;
                        trvptr = &trval;
                        while (trlptr != NULL && trlptr->drd != NULL) {
                            if (trlptr->ia[1] == 3)
                                df->dpm[1][iy][ix] += trvptr->dfdam * trlptr->drd[1];
                            trlptr = trlptr->next;
                            trvptr = trvptr->next;
                        };
		    };

            	    if (ia[2] > 0) {
			df->dpm[2][iy][ix] += dfdr * drd[2];

                        trlptr = &trlink;
                        trvptr = &trval;
                        while (trlptr != NULL && trlptr->drd != NULL) {
                            if (trlptr->ia[2] == 3)
                                df->dpm[2][iy][ix] += trvptr->dfdam * trlptr->drd[2];
                            trlptr = trlptr->next;
                            trvptr = trvptr->next;
                        };
		    }

		    if (ia[3] > 0) df->dpm[3][iy][ix] += flux * magfac;

	    	    if (ia[4] > 0) {
			df->dpm[4][iy][ix] += 
	       		    flux / ta[4] * (ta[7] - GBA/D * ta[5] * rta4ta5) * os;

                        if (normtype == 3)
			    df->dpm[4][iy][ix] -= flux / ta[4] * (ta[7] - 
				GBA/D0 * ta[5] * r0a4a5) * os;  
		    };

		    if (ia[5] > 0) {
			df->dpm[5][iy][ix] += 
			    flux * GBA/ta[5] * ( log (2.) - log (D) +
			    ta[5]/D * rta4ta5 * log (rta4) );

                        if (normtype == 3)
			    df->dpm[5][iy][ix] -= flux * GBA/ta[5] * 
				( log (2.) - log (D0) + ta[5]/D0 * r0a4a5 * 
				log (r0a4) );
		    };

	    	    if (ia[6] > 0) {
			df->dpm[6][iy][ix] += 
			    flux/ta[5] * (log (2.) - log(D));

                        if (normtype == 3)
			    df->dpm[6][iy][ix] -= 
			        flux/ta[5] * (log (2.) - log(D0));
		    };

		    if (ia[7] > 0) {
			df->dpm[7][iy][ix] += 
			    flux* (-log (2.)/ta[5] - log(rta4) + log(D)/ta[5]);

                        if (normtype == 3)
			    df->dpm[7][iy][ix] -= flux* (-log (2.)/ta[5] 
			        - log(r0a4) + log(D0)/ta[5]);
		    };

            	    if (ia[9] > 0) {
			df->dpm[9][iy][ix] += dfdr * drd[9];

                        trlptr = &trlink;
                        trvptr = &trval;
                        while (trlptr != NULL && trlptr->drd != NULL) {
                            if (trlptr->ia[9] == -1)
                                df->dpm[9][iy][ix] += trvptr->dfdam * drd[9];

                            trlptr = trlptr->next;
                            trvptr = trvptr->next;
                        };
		    };

		    if (ia[10] > 0) {
			df->dpm[10][iy][ix] += dfdr * drd[10];

                        trlptr = &trlink;
                        trvptr = &trval;
                        while (trlptr != NULL && trlptr->drd != NULL) {
                            if (trlptr->ia[10] == -1)
                                df->dpm[10][iy][ix] += trvptr->dfdam * drd[10];

                            trlptr = trlptr->next;
                            trvptr = trvptr->next;
                        };
		    };

                    /*-------------------------------------------\
                    |  Traditional diskyness/boxyness parameter  |
                    \-------------------------------------------*/

            	    if (c0 != NULL && c0->ia[1] > 0)
                	c0df->dpm[1][iy][ix] += (dfdr * drdc);

                    /*-------------------------------------------\
                    |  Bending, Fourier and PA rotational terms  |
                    \-------------------------------------------*/

                    if (b != NULL)
                        high_derivs (ix, iy, b, bdf, dfdr, bdrd, 0., NULL);

                    if (f != NULL)
                        high_derivs (ix, iy, f, fdf, dfdr, fdrd, 0., NULL);

                    if (r != NULL)
                        high_derivs (ix, iy, r, rdf, dfdr, rdrd, 0., NULL);

                    /*-------------------------------------------------------\
                    |  Hybrid or meta-profile terms.  Note that trlink       |
                    |  outer taper and inner taper correspond to two         |
                    |  different, and independent components so that         |
                    |  trlink->a for inner and trlink->a for outer are       |
                    |  mutually exclusive parameters.   Note also that a[4]  |
                    |  a[5] do not "mix" between the two.                    |
                    \-------------------------------------------------------*/

                    trlptr = &trlink;
                    trvptr = &trval;
                    while (trlptr != NULL && trlptr->a != NULL) {
                        trunc_derivs (ix,iy, trvptr->am, flux, trvptr->dfdam, trvptr,
                            trlptr);

                        trlptr = trlptr->next;
                        trvptr = trvptr->next;
                    };

 		    sumdfdI += dfdI0;
		    yc = yc + 1;
        	};
		xc = xc + 1;
    	    };

            flux = ta[3] * sumdfdI;
            df->dpm[0][iy][ix] = flux;               /* ... for convlution   */

            trlptr = &trlink;
            while (trlptr != NULL && trlptr->df != NULL) {
                if (normtype > 0)
                    trlptr->df->dpm[0][iy][ix] = 0.; /* Trunc. are not real components */

                trlptr = trlptr->next;
            };
	};
    };

   free_fitvecs (&drd, &tilt, b, &bdrd, f, &fdrd, NULL, r,
        &rdrd, NULL, normtype, &trlink, &trval);

}

#if DEBUG
    sprintf (model->name, "model-0.fits");  /* Output image */
    writefits ("test.fits", model, "", 0);
#endif

/*

float nukefunc (float ta[], float x, float y)
{
    float nuke, r, c, xd, yd, Ib;

    r = rdist (x, y, ta);
    nuke = pow (2, ( ta[6] - ta[7])/ta[5]) * ta[3] * pow((ta[4] / r), ta[7]) *
                pow(1 + pow((r/ta[4]), ta[5]), (ta[7]-ta[6])/ta[5]);
    return (nuke); 
}

*/