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

#define MINP 1e-10

/***************************************************************************\
*  Calculate truncation derivative parameters, i.e. the derivative of the   *
*  hyperbolic tangent transition function with respect to the inner         *
*  transition radius, ri, and the outer transition radius, ro.              *
\***************************************************************************/

int istrunc (char *objtype);

void high_derivs (int ix, int iy, struct fitpars *h, struct derivs *hdf,
    float dfdr, float *drd, float dfdI, float *dIdfp);

void trunc_derivs (int ix, int iy, float am, float flux, 
    float dfdr, struct trunc_val *trval, struct trunc_links *tr)
{
    void dP_drpars1 (float r, float ro, float del, int *ia, 
	struct trunc_val *trval, float *dP_drt, float *dP_ddel, int innout);

    void dP_drpars2 (float r, float ri, float ro, int *ia, struct trunc_val 
	*trval, float *dP_dri, float *dP_dro);

    float dP_drt, dP_dri, dP_dro, dP_ddel, dfdP, dPdam;

    extern struct inpars input;

    dfdP = flux / trval->P;

    if (tr->innout == 1)
	dfdP = -dfdP; 

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

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

    if (tr->lrtype == 1) {     /*  Softening length as free parameter  */

	if (tr->innout == 1)
            dP_drpars1 (am, tr->a[4] + tr->a[5], tr->a[5], tr->ia, trval, &dP_drt, 
	        &dP_ddel, tr->innout);
	else
            dP_drpars1 (am, tr->a[4], tr->a[5], tr->ia, trval, &dP_drt, 
	        &dP_ddel, tr->innout);

        if (tr->ia[4] == 1)
            tr->df->dpm[4][iy][ix] += (dfdP * dP_drt + -flux * trval->dfnorm) *
				      input.sampfac;
                                             
        if (tr->ia[5] == 1)
            tr->df->dpm[5][iy][ix] += (dfdP * dP_ddel) * input.sampfac; 

    } else {                /*  Softening radius as free parameter  */

        dP_drpars2 (am, tr->a[4], tr->a[5], tr->ia, trval, &dP_dri, &dP_dro);

        if (tr->ia[4] == 1)
            tr->df->dpm[4][iy][ix] += (dfdP * dP_dri) * input.sampfac;
                                             
        if (tr->ia[5] == 1)
            tr->df->dpm[5][iy][ix] += (dfdP * dP_dro + -flux * trval->dfnorm) *
				      input.sampfac;
    };

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

    if (tr->ia[9] != -1 && tr->ia[10] == 1)
        tr->df->dpm[10][iy][ix] += dfdr * tr->drd[10]; 

    if (tr->f != NULL)   /*  The last term is a dummy  */
        high_derivs (ix, iy, tr->f, tr->fdf, dfdr, tr->fdrd, 0., tr->fdrd);

    if (tr->b != NULL)   /*  The last term is a dummy  */
        high_derivs (ix, iy, tr->b, tr->bdf, dfdr, tr->bdrd, 0., tr->bdrd);

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

}


/*****************************************************************************\
*  Type 1 truncation:  Calculates dP_dri or dP_dro and dP_ddel.               *
*  Type 1 truncs have break radius and softening length as free parameters.   *
\*****************************************************************************/


void dP_drpars1 (float r, float ro, float del, int *ia, 
    struct trunc_val *trval, float *dP_drt, float *dP_ddel, int innout)
{
    float dSdro, dSddel;

    if (ia[4] == 1) {
        dSdro = -4.95 / del;
        *dP_drt = 0.5 * (1- trval->tanh * trval->tanh) * (dSdro * (1-r/ro) -
				r/ro/ro * trval->gtr);
    } else
	*dP_drt = 0.;

    if (ia[5] == 1) {
	if (innout == 0) {
            dSddel = 4.95 * ro / del / del;
            *dP_ddel = 0.5 * (1- trval->tanh * trval->tanh)* (dSddel * (1-r/ro));
	} else {
            dSddel = 4.95 * (ro - del) / del / del;
            *dP_ddel = 0.5 * (1- trval->tanh * trval->tanh)* (dSddel * (1-r/ro) -
				r/ro/ro * trval->gtr);
	};

    } else
	*dP_ddel = 0.;

}


/*****************************************************************************\
*  Type 2 truncation:  Calculates dP_dri and dP_dro.                          *
*  Type 2 trunc has softening radius (instead of length) as free parameter    *
\*****************************************************************************/


void dP_drpars2 (float r, float ri, float ro, int *ia, 
    struct trunc_val *trval, float *dP_dri, float *dP_dro)
{
    float dSdri, dSdro, dGTRdri, dGTRdro;

    if (ia[4] == 1) {
        dSdri = -4.95 * ro / (ro - ri) / (ro - ri);
        dGTRdri = -dSdri;
        *dP_dri = 0.5 * (1- trval->tanh * trval->tanh) * (dGTRdri* r/ ro + 
		  dSdri);
    } else
        *dP_dri = 0.;

    if (ia[5] == 1) {
        dSdro = 4.95 * ri / (ro - ri) / (ro - ri);
        dGTRdro = -dSdro;
        *dP_dro = 0.5 * (1- trval->tanh * trval->tanh) * (dGTRdro* r/ ro + 
             dSdro + trval->gtr * -r / ro / ro);
    } else
        *dP_dro = 0.;

}


/*****************************************************************************\
*  Add the truncation derivative images from individual subcomponents,        *
*  in "fpar->high" into the main structure for the transition function.       *
\*****************************************************************************/

void assemble_trunc_df (struct fitpars *fpar, struct derivs *df)
{

    void add_trunc_derivs (struct fitpars *fpar, struct derivs *df, 
        struct fitpars *trfpar, struct derivs *trdf, int *comp);

    struct fitpars *fptr, *fptr2;
    struct derivs  *dptr, *dptr2;

    fptr = fpar;
    dptr = df;

    while (fptr != NULL) {
	if (istrunc (fptr->objtype)) {
	    add_trunc_derivs (fpar, df, fptr, dptr, fptr->outer);
	    add_trunc_derivs (fpar, df, fptr, dptr, fptr->inner);
	};

	fptr = fptr->next;
	dptr = dptr->next;
    };
}

/*****************************************************************************\
*  Add the truncation derivative images from individual subcomponents,        *
*  in "fpar->high" into the main structure for the transition function.       *
*									      *
*  fpar and df are the master structures, whereas trfpar and trdf are	      *
*  pointed at the truncation structure that we want to copy information       *
*  into from "df->dpm".                                                       *
\*****************************************************************************/

void add_trunc_derivs (struct fitpars *fpar, struct derivs *df, 
    struct fitpars *trfpar, struct derivs *trunc_df, int *comp)
{
    int i, j, k, ix, iy, done;
    struct fitpars *fptr, *fptr2;
    struct derivs *dptr, *dptr2, *trunc_df2;

    /*---------------------------------------------------\
    |  Look for components involved in the linking.      |
    \---------------------------------------------------*/

    for (i=1; i <= comp[0]; i++) {
	done = 0;

	fptr = fpar;
	dptr = df;

	/*  Step into the struct list to the right component  */

	for (j=1; j < comp[i]; j++) {
	    fptr = fptr->next;
	    dptr = dptr->next;
	};

	/*---------------------------------------------------------------\
        |  Once there, step into the fpar->high list.                    |
	|  Truncation parameters are stored in the high order struct     |
        |  of individual sub-components involved in linking.  To see     |
        |  why it's done this way, consider, for instance, when you      |
        |  link 2 sub-components f1(r) and f2(r).  In the product of     |
        |    F1(r) = f1(r) * P1(r) and  F2(r) = f2(r) * (1-P1(r))        |
        |  P1(r) shares the centroid of f1(r) because it is tapering     |
        |  the f1(r).  On the other hand, (1-P1(r)) shares the same      |
        |  centroid as f2(r).  Thus the derivative images have to be     |
        |  produced for f1*P1 and f2*(1-P1), individually, then          |
	|  convolved even though P1 and (1-P1) share the same set of     |
        |  trunc parameters.  The net derivative image for the trunc     |
        |  parameter is then assembled together below in trunc_df->dpm.  |
	\---------------------------------------------------------------*/
	       
	fptr = fptr->high;
	dptr = dptr->high;

	while (fptr != NULL && !done) {

	    /*---------------------------------------------------------------\
	    |  Search for the trunc component involved, then figure out if   |
	    |  the component is tapered on the outside or the inside         |
            |  because each subcomp in fptr->high can be one or both.  The   |
            |  way this is done is to see if the main trunc structure        |
	    |  component is the same as the number stored in fptr->tr_num,   |
	    |  which was set in "link_fpar_high.                             |
	    \---------------------------------------------------------------*/

	    if (istrunc (fptr->objtype) && fptr->tr_num == trfpar->compnum ) {

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

			for (k = 1; k <= fptr->npars; k++) {
		            if (fptr->ia[k] == 1)
		                trunc_df->dpm[k][iy][ix] += 
							dptr->dpm[k][iy][ix];
				    
			};

			/*-------------------------------------------\
			|  Now the fourier, diskiness/boxiness, or   |
			|  bending params.                           |
			\-------------------------------------------*/

			if (fptr->high != NULL) {
			    fptr2 = fptr->high;
			    dptr2 = dptr->high;
			    trunc_df2 = trunc_df->high;

			    while (fptr2 != NULL) {
				for (k = 1; k <= fptr2->npars; k++) {
		            	    if (fptr2->ia[k] == 1)  {
		                        trunc_df2->dpm[k][iy][ix] += 
						dptr2->dpm[k][iy][ix];
				    };
				};

				fptr2 = fptr2->next;
				dptr2 = dptr2->next;
				trunc_df2 = trunc_df2->next;
			    };
			};
		    };
		};
		done = 1;

	    };

	    fptr = fptr->next;
	    dptr = dptr->next;
	};
    };
}


/*****************************************************************************\
*  Calculates df/dr, i.e. dfdam, and df/dI * P, where P is the hyperbolic     *
*  tangent function.  Truncation profiles are multiplied by a hyperbolic      *
*  tangent.                                                                   *
*                                                                             *
*  Note that trunc_out (outer taper) and trunc_in (inner taper) are actually  *
*  two independent "components", and each one has its own inner (a[3])        *
*  transition radius and softening length (a[4]).  Don't get trunc_out and    *
*  trunc_in confused with a[3] or a[4].                                       *
\*****************************************************************************/

void trunc_dfdr_calc (struct trunc_val *trval, struct trunc_links *trlink, 
    float flux0, float *PP)
{

    struct trunc_links *trlptr;
    struct trunc_val *trvptr;

    trlptr = trlink;
    trvptr = trval;

    while (trlptr != NULL && trlptr->a != NULL) {
        if (trlptr->innout == 1) {             		/*  Outer truncation  */

	    if (trlptr->lrtype == 1) {   /*  Softening *length* as free param.  */

                trvptr->tanh = tanh(trvptr->gtr * trvptr->am/(trlptr->a[4] +
			     trlptr->a[5]) + trvptr->s);
                trvptr->P = FMAX (1 - 0.5*(trvptr->tanh + 1), MINP); 

		/*************************************************************\
		*  dfdam is only completely partially here, and fully below.  *
		\*************************************************************/

                trvptr->dfdam = 0.5 * (1-trvptr->tanh * trvptr->tanh)*trvptr->gtr/
		               (trlptr->a[4] + trlptr->a[5]);

	    } else {                    /*  Softening *radius* as free param.  */

                trvptr->tanh = tanh(trvptr->gtr * trvptr->am/
			      trlptr->a[5] + trvptr->s);
                trvptr->P = FMAX (1 - 0.5*(trvptr->tanh + 1), MINP);
                trvptr->dfdam = 0.5 * (1-trvptr->tanh * trvptr->tanh)*trvptr->gtr / 
		               trlptr->a[5];
	    };
	};

        if (trlptr->innout == 0) {  	 	      /*  Inner truncation  */

	    if (trlptr->lrtype == 1) {   /*  Softening *length* as free param.  */

                trvptr->tanh = tanh(trvptr->gtr * trvptr->am/trlptr->a[4] + 
			      trvptr->s);
                trvptr->P = FMAX (0.5*(trvptr->tanh + 1), MINP);

                trvptr->dfdam = 0.5* (1-trvptr->tanh * trvptr->tanh) * trvptr->gtr/
		               trlptr->a[4];

	    } else {                    /*  Softening *radius* as free param.  */

                trvptr->tanh = tanh(trvptr->gtr * trvptr->am/trlptr->a[5] + 
			      trvptr->s);
                trvptr->P = FMAX (0.5*(trvptr->tanh + 1), MINP);
                trvptr->dfdam = 0.5 * (1-trvptr->tanh * trvptr->tanh) * 
		               trvptr->gtr / trlptr->a[5];
	    };
	};

	*PP = *PP * trvptr->P;

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

    trlptr = trlink;
    trvptr = trval;

    while (trlptr != NULL && trlptr->a != NULL) {

	trvptr->dfdam = trvptr->dfdam * flux0 * *PP / trvptr->P;

	if (trlptr->innout == 1) 
	    trvptr->dfdam = -trvptr->dfdam; 

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

}