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

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

float bessk1 (float x);
float bessk0 (float x);

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 edgedisk (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; 
    int j, nsubsamp, os, sgnx, sgnh, sgn; 
    int trunc=0, ttype=0, httype=0, lttype=0;
    float xp, yp, xd, yd, ta[NPARS+1], cosPA, sinPA, sumdfdI, xf, yf,
	  xc, yc, flux, dfdI, step, err, magfac, x, h, xrs, yhs, dxdx0, 
	  dxdy0, dhdx0, dhdy0, dfdx, dfdh, K0, K1, dK1dx, dK1drs,
	  dxdtheta, dhdtheta, hd, rd;
    float trunc_norm=1., flux0, normrad, dfnorm, PP;
    float *bdrd, *fdrd, *rdrd;
    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;
    float dum, *dumptr=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;

    /*-------------------------------------------------------------------\
    |  Edge on disk can have truncations either in the height or length  |
    |  direction, each of which can have an inner & outer truncation.    |
    |  Hence, need twice the storage space.                              |
    \-------------------------------------------------------------------*/

    os = input.sampfac;

    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);

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

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

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

    /*----------------------------------------------------------------\
    |  Pre-calculate some things we will frequently use later on but  |
    |  which needs to be calculated one time only.                    |
    |								      |
    |  a[3] is the galaxy total brightness while ta[3] is the count   |
    |  rate at the disk scalelength.      		      	      |
    \----------------------------------------------------------------*/


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

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

    cosPA = cos (ta[10]);
    sinPA = sin (ta[10]);

    /*------------------------------------------------------------\
    |  Take care of truncation parameters and pointers.           |
    \------------------------------------------------------------*/
   
    trunc_pars (objtype, ia, high, df, &trlink, &trval);

    if (r != NULL)       /*  A spiral component.  */
        tilt = trunc_spiral (r);

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

    get_trunc_cons ("edgedisk", normtype, &trlink, &trval, ta, 0., cosPA,
        sinPA, &trunc_norm, &dfnorm);

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

    if (normtype == 2) {
	if (strncmp (trlink.objtype, "height", 6) == 0)
            normrad = ta[4];            /*  Scale height        */
	else if (strncmp (trlink.objtype, "length", 6) == 0)
	    normrad = ta[5];            /*  Scale length        */
    };

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

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

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

	    xd = ix - ta[1];
	    yd = iy - ta[2];

	    x = xd * cosPA - yd * sinPA;
	    h = xd * sinPA + yd * cosPA;

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

	    nsubsamp = 1;
	    yhs = fabs(h) / ta[4];
            xrs = fabs(x) / ta[5];

            if (h < OVERSAMPR) {
                if (ta[4] <= 1 && FMIN(fabs(h), fabs(x)) <= os ) {

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

                    nsubsamp = IMIN(100, (int)(20./FMIN(fabs(h), fabs(x))/os)); 
                } else  if (FMIN(fabs(h), fabs(x)) <= 3. * os)
                    nsubsamp = 5 / os;
                else {
                    /* If scalelength is large, bigger step */

                    nsubsamp = IMIN(100, NINT(10./FMIN(fabs(h), fabs(x))/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 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  |
                    \-----------------------*/

		    xd = xp - ta[1];
		    yd = yp - ta[2];

		    x = xd * cosPA - yd * sinPA;
		    h = xd * sinPA + yd * cosPA;

		    hd = fabs(h);
		    rd = fabs(x);
		    yhs = hd / ta[4];
	            xrs = rd / ta[5];

		    sgnx = (x < 0. ? -1 : 1);
		    sgnh = (h < 0. ? -1 : 1);

		    K1 = bessk1 (xrs);
		    K0 = bessk0 (xrs);

		    if (xrs != 0.)
	                dfdI = xrs * K1 * 1./cosh (yhs) * 1./cosh (yhs) / 
			   totsamp / trunc_norm;
		    else
			dfdI = 1./cosh (yhs) * 1./cosh (yhs) / 
			   totsamp / trunc_norm;

		    if (xrs != 0.)
		        dK1dx = sgnx * (-K0 - K1 / xrs) / ta[5];
		    else
			dK1dx = 0.;

		    flux0 = ta[3] * dfdI;

		    if (x != 0.) 
			dfdx = flux0 * (sgnx/fabs(x) + dK1dx / K1);
		    else 
			dfdx = 0.;

		    if (h== 0.) 
			dfdh = 0.;
		    else
		        dfdh = -2 * flux0 * tanh(yhs) / ta[4] * sgnh;

                    if (trlink.objtype != NULL && normtype > 0) {

			/*  Assign radial coordinate for truncation  */

                        trlptr = &trlink;
                        trvptr = &trval;
                        while (trvptr != NULL && trlptr->drd != NULL) {
			    if (strncmp (trlptr->objtype, "height", 6) == 0)
				trvptr->am = hd;
			    else if (strncmp (trlptr->objtype, "length", 6) == 0)
				trvptr->am = rd;

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

                        PP = 1;
                        trunc_dfdr_calc (&trval, &trlink, flux0, &PP);
                        dfdI = dfdI * PP;
                        dfdx = dfdx * PP;
                        dfdh = dfdh * PP;
                    };

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

		    flux = ta[3] * dfdI;

                    if (ia[1] > 0) {
			dxdx0 = -cosPA;
			dhdx0 = -sinPA;
			df->dpm[1][iy][ix] += (dfdx*dxdx0 + dfdh*dhdx0); 

                        trlptr = &trlink;
                        trvptr = &trval;
                        while (trlptr != NULL && trlptr->drd != NULL) {
                            if (trlptr->ia[1] == 3) {
				if (strncmp (trlptr->objtype, "height", 6) == 0) 
                                    df->dpm[1][iy][ix] += sgnh * trvptr->dfdam * dhdx0;
				else if (strncmp (trlptr->objtype, "length", 6) == 0) 
                                    df->dpm[1][iy][ix] += sgnx * trvptr->dfdam * dxdx0;
			    };
                            trlptr = trlptr->next;
                            trvptr = trvptr->next;
                        };
		    };

                    if (ia[2] > 0) {
			dxdy0 = sinPA;
			dhdy0 = -cosPA; 
			df->dpm[2][iy][ix] += (dfdx*dxdy0 + dfdh*dhdy0);

                        trlptr = &trlink;
                        trvptr = &trval;
                        while (trlptr != NULL && trlptr->drd != NULL) {
                            if (trlptr->ia[1] == 3) {
                                if (strncmp (trlptr->objtype, "height", 6) == 0) 
                                    df->dpm[2][iy][ix] += sgnh * trvptr->dfdam * dhdy0;
                                else if (strncmp (trlptr->objtype, "length", 6) == 0)
                                    df->dpm[2][iy][ix] += sgnx * trvptr->dfdam * dxdy0;
                            };
                            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] += (2 * flux*tanh (yhs) * 
							yhs/ta[4]) * os;
                        if (normtype == 3 && strncmp (trlink.objtype, "height", 6) == 0)
		 	    df->dpm[4][iy][ix] += flux * dfnorm * 
					          normrad / ta[4] * os;
		    };

		    if (ia[5] > 0) {
			if (xrs != 0.) {
			    dK1drs = - xrs * (-K0 - K1 / xrs) / ta[5];
		            df->dpm[5][iy][ix] += flux * (-1./ ta[5] + dK1drs/K1) * os;
			} else
		            df->dpm[5][iy][ix] += flux * (-1./ ta[5]) * os;

                        if (normtype == 3 && strncmp (trlink.objtype, "length", 6) == 0)
			        df->dpm[5][iy][ix] += flux * dfnorm * 
					              normrad / ta[5] * os;
		    };

		    if (ia[10] > 0) {
			dxdtheta = (-xd * sinPA - yd * cosPA) * -PI / 180.;
			dhdtheta = (xd * cosPA - yd * sinPA) * -PI / 180.;
                        df->dpm[10][iy][ix] += (dfdx*dxdtheta + dfdh*dhdtheta);

                        trlptr = &trlink;
                        trvptr = &trval;
                        while (trlptr != NULL && trlptr->drd != NULL) {

                            if (trlptr->ia[10] == -1) {
                                if (strncmp (trlptr->objtype, "height", 6) == 0)
                                    df->dpm[10][iy][ix] += sgnh * trvptr->dfdam * dhdtheta;
                                else if (strncmp (trlptr->objtype, "length", 6) == 0)
                                    df->dpm[10][iy][ix] += sgnx * trvptr->dfdam * dxdtheta;
			    };

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

                    /*-------------------------------------------------------\
                    |  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 += dfdI;
		    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 (&dumptr, &tilt, b, &bdrd, f, &fdrd, NULL, r,
        &rdrd, NULL, normtype, &trlink, &trval);
}


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