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

#define OVERSAMPR 10.
#define MINR 1.e-15

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 **dI0dfp,
        struct fitpars *r, float **rdrd, float **dI0drp, int normtype,
        struct trunc_links *trlinks, struct trunc_val *trval);

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

    unsigned long i, totsamp;
    int j, ix, iy, nsubsamp, os, mode, n, trunc=0;
    float am, ta[NPARS+1], step, xp, yp, xf, yf, xc, yc, rrad, dfdA, flux,
	  flux0, dfdr, OcosPA, OsinPA, c=2., err, rro, r0ro, rro2mb, 
	  r0ro2mb, r0r2, pre_dfdal, sumdfdA, zpt, magfac, 
	  drdc, twomb;
    float *drd, *bdrd, *fdrd, *rdrd;
    float hdrdc, trunc_norm=1., dfnorm=0., normrad, PP;
    int iatmp[10], *iaptr;
    float atmp[10], *aptr;
    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;

    os = input.sampfac;
    magfac = log(10.) / -2.5;
    
    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)
        c += 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);

    twomb = 2-ta[6];

    zpt = (ta[3] - model->muzpt)/-2.5;
    ta[3] = pow (10., zpt);

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

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

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

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

    if (normtype == 1)
        normrad = 0.;                   /*  Central surface brightness  */

    if (normtype == 2 || normtype == 3) {

        if (normtype == 2)
            normrad = ta[4];                 /*  Scale radius        */

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

	r0ro = normrad / ta[4];
	r0ro2mb = pow(r0ro, twomb);
	r0r2 = 1-r0ro2mb;
    };

    /*------------------------------------------\
    |  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 (am < OVERSAMPR) {            
                if (ta[4] * ta[9] <= 1. && am <= os) {
	    
                    /* If scalelength is small, oversample by a lot more.  */
	    
                    nsubsamp = IMIN(100, (int)(1./fabs(ta[4] * ta[9]) * 
			       40 / 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;
	    
            step = 1./nsubsamp;
	    xf = ix - 0.5 + 1./(2 * nsubsamp);
	    yf = iy - 0.5 + 1./(2 * nsubsamp);
	    
	    sumdfdA = 0.;

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

	    /*  Truncate the fit at the tidal radius  */	    

	    xc = 0.;
	    while (xc <= nsubsamp - 1 && am < ta[4]) {
	        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  |
                    \-----------------------*/

	            rro = am / ta[4];
		    rro2mb = pow(rro, twomb);
		    rrad =  1 - rro2mb;

		    if (rrad < 1e-5) {
			if (rrad < 0.)
			    rrad = 0.;
		   	else
			    rrad = 1e-5;
		    };

		    dfdA = pow (rrad, ta[5]) / totsamp / trunc_norm;
	            flux0 = ta[3] * dfdA;

		    pre_dfdal = flux0 * ta[5] / rrad; 
	            dfdr = -pre_dfdal * twomb * rro2mb / am;

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

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

                    flux = ta[3] * dfdA;

                    /*--------------------------------------------------------\
                    |  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] += pre_dfdal* twomb * rro2mb / 
					      ta[4] * os;

                        if (normtype == 3)
                            df->dpm[4][iy][ix] += flux * dfnorm * normrad /
				                    ta[4] * os; 
		    };

	            if (ia[5] > 0) {
	                df->dpm[5][iy][ix] += flux * log(rrad);

                        if (normtype == 3)
                            df->dpm[5][iy][ix] += -flux * log(r0r2);
		    };
		    
	            if (ia[6] > 0) {
	                df->dpm[6][iy][ix] += pre_dfdal * rro2mb * 
					      log(rro);

                        if (normtype == 3)
	                    df->dpm[6][iy][ix] += -flux * ta[5] / r0r2 * 
						  r0ro2mb * log(r0ro);
		    };
		    
                    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;
                    };

	            sumdfdA += dfdA;
		    yc = yc + 1.;
	        };	    
	        xc = xc + 1.;
	    };	    
		    
	    flux = sumdfdA * ta[3];
	    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);

}