/*****************************************************************************\
*  Calculate truncation transition parameters and constants for derivative    *
*  calculations.  The parameters calculated here don't change with radial     *
*  position of the pixels, thus only need to be performed once per model.     *
*									      *
*  Note that:								      *
*    trlink->lrtype = 1  has break radius and softening *length* as params.   *
*    trlink->lrtype = 2  has break radius and softening *radius* as params.   *
*									      *
\*****************************************************************************/

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

#define EPS 1e-5
float bessk1 (float x);
float bessk0 (float x);

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 normfacs (char *objtype, int normtype, float *trunc_norm, 
	float *dfnorm, struct trunc_links *trlink, struct trunc_val *trval, 
	float a[], float kap);

    float pa;

    /*--------------------------------------------------------------\
    |  Calculate the parameters needed for hyperbolic tangent:      |
    \--------------------------------------------------------------*/

    normfacs (objtype, normtype, trunc_norm, dfnorm, trlink, trval, a, kap);

    while (trlink != NULL && trlink->objtype != NULL) {                 

	if (trlink->lrtype == 1) {

	    if (trlink->innout == 1)             /*  Outer truncation  */
                trval->s = 2.65 - 4.95 * (trlink->a[4] + trlink->a[5]) /
		              trlink->a[5];

	    else if (trlink->innout == 0)        /*  Inner truncation  */
                trval->s = 2.65 - 4.95 * trlink->a[4] / trlink->a[5];

	} else if (trlink->lrtype == 2) {    

            trval->s = 2.65 - 4.95 * trlink->a[5] /
		          (trlink->a[5] - trlink->a[4]);
	};

        trval->gtr = 2 - trval->s;

        /*--------------------------------------------------------------------\
        |  If the truncation doesn't have independent q and PA, then the      |
        |  default is to assume that they are the same as the Sersic comp..   |
        \--------------------------------------------------------------------*/

	if (trlink->ia[1] == -1) {
            trlink->a[1] = a[1];  /* This is still needed for convolution box */
            trlink->a[2] = a[2];
        };

	/*------------------------------------------\
        |  PA and axis ratio are free parameters    |
	\*------------------------------------------/

	if (trlink->ia[9] != -1) {      
            pa = -trlink->a[10] / 180. * PI + PI/2.;
            trunc_c1->sinPA = sin (pa);
            trunc_c1->cosPA = cos (pa);
	} else {
	    trlink->a[9] = a[9];
	    trlink->a[10] = a[10];
            trunc_c1->sinPA = OsinPA;
            trunc_c1->cosPA = OcosPA;
	};

        /*--------------------------------------------------------------------\
        |  If the truncation doesn't have independent q and PA, then the      |
        |  default is to assume that they are the same as the Sersic comp..   |
        \--------------------------------------------------------------------*/

	if (trlink->ia[1] == -1) {
            trlink->a[1] = a[1];  /* This is still needed for convolution box  */
            trlink->a[2] = a[2];
	};

        /* PA and axis ratio are free parameters */

	if (trlink->ia[9] != -1) {
            pa = -trlink->a[10] / 180. * PI + PI/2.;
            trval->sinPA = sin (pa);
            trval->cosPA = cos (pa);
	} else {
	    trlink->a[9] = a[9];
	    trlink->a[10] = a[10];
            trval->sinPA = OsinPA;
            trval->cosPA = OcosPA;
	};

	trlink = trlink->next;
	trval = trval->next;
    };
}


/****************************************************************************\
*  Calculate the normalization factor for truncation, depending on whether   *
*  to normalize at the galaxy center, effective radius, outer radius, mean   *
*  radius, or what not.  This is specified in 'normtype:'                    *
*									     *
*    0 = total mag,  1 = central SB,  2 = effective radius                   *
*    3 = break radius,  4 = mean of inner and outer break radius             *
*								             *
*  dfnorm is the derivative of the trunc normalization with respect to       *
*  the radius parameter that's used as a free parameter (e.g. df/dro,        *
*  df/dre, etc.), divided by fnorm.  For instance, given the definition      *
*  that:               							     *
*									     *
*       flux (rp) = flux0 (r=rp) * P / fnorm                                 *
*                 = a[3] * P,              			             *
*								             *
*	"dfnorm" =  dfnorm/drp / fnorm					     *
*   									     *
*  Note:  be careful about the sign in dfnorm!		                     *
*   									     *
\****************************************************************************/

void normfacs (char *objtype, int normtype, float *trunc_norm, float *dfnorm, 
    struct trunc_links *trlink, struct trunc_val *trval, float a[], float kap)
{
    float invp5, hh, rr, K1, K0, rr2, rdiff, GBA, r0a4, D0;

    /*-----------------\
    |  Sersic profile  |
    \-----------------*/

    if (strncmp (objtype, "sersic", 6) == 0) {
        invp5 = 1./a[5];

	if (normtype == 1) {
            *trunc_norm = exp(kap);
	    *dfnorm = 0.;
	};

	if (trlink->lrtype == 1) {  /* Truncation length being free par  */
	    if (normtype == 3)  {   /* Norm at break radius  */

                rr = pow(trlink->a[4]/a[4], invp5);

                *trunc_norm = exp(-kap*(rr-1.));
                *dfnorm = -kap*invp5 / a[4] * pow(trlink->a[4]/a[4],
                              invp5-1);
	    };
	} else {                    /* Truncation radius being free par  */
            if (normtype == 3) {

                rr = pow(trlink->a[5]/a[4], invp5);

                *trunc_norm = exp(-kap*(rr-1.));
                *dfnorm = -kap*invp5 / a[4] * pow(trlink->a[5]/a[4],
                          invp5-1);
	    };
	};
    };

    /*----------------------\
    |  Exponential profile  |
    \----------------------*/

    if (strncmp (objtype, "expdisk", 7) == 0) {
        if (normtype == 1) {
            *trunc_norm = 1.;
            *dfnorm = 0.;
        };

        if (trlink->lrtype == 1) {  /* Truncation length being free par  */
            if (normtype == 3)  {   /* Norm at break radius  */

                *trunc_norm = exp(- trlink->a[4]/a[4]);
                *dfnorm = -1. / a[4];
            };
        } else {                    /* Truncation radius being free par  */
            if (normtype == 3) {

                *trunc_norm = exp(- trlink->a[5]/a[4]);
                *dfnorm = -1. / a[4];
            };
        };
    };

    /*-------------------\
    |  Gaussian profile  |
    \-------------------*/

    if (strncmp (objtype, "gaussian", 8) == 0) {
        if (normtype == 1) {
            *trunc_norm = 1.;
            *dfnorm = 0.;
        };

        if (trlink->lrtype == 1) {  /* Truncation length being free par  */
            if (normtype == 3)  {   /* Norm at break radius  */

                rr = pow(trlink->a[4]/a[4], 2);
                *trunc_norm = exp(-rr / 2.);
                *dfnorm = -trlink->a[4] / a[4] / a[4];

            };
        } else {                    /* Truncation radius being free par  */
            if (normtype == 3) {

                rr = pow(trlink->a[5]/a[4], 2);

                *trunc_norm = exp(-rr / 2.);
                *dfnorm = -trlink->a[5] / a[4] / a[4];
            };
        };
    };

    /*-----------------\
    |  Moffat profile  |
    \-----------------*/

    if (strncmp (objtype, "moffat", 6) == 0) {
        if (normtype == 1) {
            *trunc_norm = 1.;
            *dfnorm = 0.;
        };

        if (trlink->lrtype == 1) {  /* Truncation length being free par  */
            if (normtype == 3)  {   /* Norm at break radius  */

                rr = trlink->a[4]/a[4];
                *trunc_norm = 1./ pow(1 + rr*rr, a[5]);
                *dfnorm = -2 * a[5] * trlink->a[4] / a[4] / a[4] / (1 + rr*rr);
            };
        } else {                    /* Truncation radius being free par  */
            if (normtype == 3) {

                rr = trlink->a[5]/a[4];
                *trunc_norm = 1./ pow(1 + rr*rr, a[5]);
                *dfnorm = -2 * a[5] * trlink->a[5] / a[4] / a[4] / (1 + rr*rr);
            };
        };
    };

    /*-----------------\
    |  Ferrer profile  |
    \-----------------*/

    if (strncmp (objtype, "ferrer", 6) == 0) {
        if (normtype == 1) {
            *trunc_norm = 1.;
            *dfnorm = 0.;
        };

        if (trlink->lrtype == 1) {  /* Truncation length being free par  */
            if (normtype == 3)  {   /* Norm at break radius  */

                rr = trlink->a[4]/a[4];
                rr2 = 1 - pow(rr, 2 - a[6]);

                *trunc_norm = pow(rr2, a[5]);
                *dfnorm = -a[5]*(2-a[6])/rr2* pow(rr, 2 - a[6])/trlink->a[4];
            };
        } else {                    /* Truncation radius being free par  */
            if (normtype == 3) {

                rr = trlink->a[5]/a[4];
                rr2 = 1 - pow(rr, 2 - a[6]);

                *trunc_norm = pow(rr2, a[5]);
                *dfnorm = -a[5]*(2-a[6])/rr2* pow(rr, 2 - a[6])/trlink->a[5];
            };
        };
    };

    /*---------------\
    |  King profile  |
    \---------------*/

    if (strncmp (objtype, "king", 4) == 0) {
        if (normtype == 1) {
            *trunc_norm = 1.;
            *dfnorm = 0.;
        };

        if (normtype == 2) {
            rr2 = pow (a[5]/a[4], 2.);
            rdiff = 1./pow(2, 1./a[6]) - 1./pow(1+rr2, 1./a[6]);

            *trunc_norm = pow (rdiff, a[6]);
            *dfnorm = 0.;
        };

        if (trlink->lrtype == 1) {  /* Truncation length being free par  */
            if (normtype == 3)  {   /* Norm at break radius  */

                rr = pow (trlink->a[4] / a[4], 2.);
                rr2 = pow (a[5]/a[4], 2.);
                rdiff = 1./pow(1+rr, 1./a[6]) - 1./pow(1+rr2, 1./a[6]);

                *trunc_norm = pow (rdiff, a[6]);
                *dfnorm = -2 / rdiff * pow (1+rr, -1./a[6] - 1.) * rr /
                                trlink->a[4];
            };
        } else {                    /* Truncation radius being free par  */
            if (normtype == 3) {

                rr = pow (trlink->a[5] / a[4], 2.);
                rr2 = pow (a[5]/a[4], 2.);
                rdiff = 1./pow(1+rr, 1./a[6]) - 1./pow(1+rr2, 1./a[6]);

                *trunc_norm = pow (rdiff, a[6]);
                *dfnorm = -2 / rdiff * pow (1+rr, -1./a[6] - 1.) * rr /
                                trlink->a[5];
            };
        };
    };

    /*----------------\
    |  Nuker profile  |
    \----------------*/

    if (strncmp (objtype, "nuker", 5) == 0) {
        if (normtype >= 0 && normtype <= 2) {
            *trunc_norm = 1.;
            *dfnorm = 0.;
        };

        if (trlink->lrtype == 1) {  /* Truncation length being free par  */
            if (normtype == 3)  {   /* Norm at break radius  */

                GBA = (a[7]-a[6])/a[5];
                r0a4 = trlink->a[4] / a[4];
                D0 = 1 + pow(r0a4, a[5]);

                *trunc_norm = pow (2, -GBA) * pow (r0a4, -a[7]) *
                           pow (D0, GBA);
                *dfnorm = 1. / trlink->a[4] * (-a[7] + GBA / D0 *
                          pow(r0a4, a[5]) * a[5]);
            };
        } else {                    /* Truncation radius being free par  */
            if (normtype == 3) {

                GBA = (a[7]-a[6])/a[5];
                r0a4 = trlink->a[5] / a[4];
                D0 = 1 + pow(r0a4, a[5]);

                *trunc_norm = pow (2, -GBA) * pow (r0a4, -a[7]) *
                           pow (D0, GBA);
                *dfnorm = 1. / trlink->a[5] * (-a[7] + GBA / D0 *
                          pow(r0a4, a[5]) * a[5]);
            };
        };
    };

    /*-----------------------\
    |  Edge-on disk profile  |
    \-----------------------*/

    if (strncmp (objtype, "edgedisk", 8) == 0) {
        if (normtype == 0 || normtype == 1) {
            *trunc_norm = 1.;
            *dfnorm = 0.;
        };

        if (normtype == 2) {

	    if (strncmp (trlink->objtype, "height", 6) == 0) {
                *trunc_norm = 1./pow(cosh(1), 2.);
                *dfnorm = 0.;
            };

	    if (strncmp (trlink->objtype, "length", 6) == 0) {
                *trunc_norm = bessk1 (1.);
                *dfnorm =  0.;
            };
        };

	if (normtype == 3) {

	    if (strncmp (trlink->objtype, "height", 6) == 0) {

                if (trlink->lrtype == 1)        /* Truncation length being free par  */
                    hh = trlink->a[4] / a[4];
                else if (trlink->lrtype == 2)   /* Truncation radius being free par  */
                    hh = trlink->a[5] / a[4];

                *trunc_norm = 1./pow(cosh(hh), 2.);
                *dfnorm = -sinh(hh)/cosh(hh) * 2 / a[4];
            };


	    if (strncmp (trlink->objtype, "length", 6) == 0) {

                if (trlink->lrtype == 1)        /* Truncation length being free par  */
                    rr = trlink->a[4] / a[5];
                else if (trlink->lrtype == 2)   /* Truncation radius being free par  */
                    rr = trlink->a[5] / a[5];

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

                *trunc_norm = rr * K1;
                *dfnorm =  (K1 + rr * (-K0 - K1/rr))/ K1/rr/a[5];
            };
        };
    };

    trval->dfnorm = *dfnorm;
}