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

/***********************************************************************\
*  Calculate hybrid 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.          *
\***********************************************************************/

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

void hybderivs (int ix, int iy, float am, float flux, float flux0,
    struct hybval hybc, float dfdr, struct hyblinks *h)
{
    void dP_drpars1 (float r, float ri, float del, int *ia, struct hybval hybc, 
        float *dP_dro, float *dP_ddel);

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

    float dP_dri, dP_dro, dP_ddel, dfdP, dPdam;

    /*--------------------------------------------------------------------\
    |  Note that  hout->df has opposite sign as hin->df.  That is taken   |
    |  care of in the passing of parameter, which changes sign depending  |
    |  on which derivative (hout or hin) are being calculated.            |
    \--------------------------------------------------------------------*/

    if (h->type == 1) {

        dP_drpars1 (am, h->a[3], h->a[4], h->ia, hybc, &dP_dro, &dP_ddel);
    
        if (h->ia[3] == 1)
            h->df->dpm[3][iy][ix] += (flux0 * dP_dro + flux * hybc.dfnorm);
                                             
        if (h->ia[4] == 1)
            h->df->dpm[4][iy][ix] += (flux0 * dP_ddel + flux * hybc.dfnorm); 

    } else {

        dP_drpars2 (am, h->a[3], h->a[4], h->ia, hybc, &dP_dri, &dP_dro);

        if (h->ia[3] == 1)
            h->df->dpm[3][iy][ix] += (flux0 * dP_dri);
                                             
        if (h->ia[4] == 1)
            h->df->dpm[4][iy][ix] += (flux0 * dP_dro + flux * hybc.dfnorm);

    };

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

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

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

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

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

}


/*****************************************************************************\
*  Type 1 hybrid:  Calculates dP_dri and dP_ddel.                             *
*  Type 1 hybrids have an outer truncation and no inner truncation.           *
\*****************************************************************************/


void dP_drpars1 (float r, float ro, float del, int *ia, struct hybval hybc, 
    float *dP_dro, float *dP_ddel)
{
    float dSdro, dSddel;

    if (ia[3] == 1) {
        dSdro = -4.95 / del;
        *dP_dro = 0.5 * (1- hybc.tanh * hybc.tanh) * (dSdro -
				r/ro/ro * hybc.gtr - dSdro * r/ro);
    } else
	*dP_dro = 0.;

    if (ia[4] == 1) {
        dSddel = 4.95 * ro / del / del;
        *dP_ddel = 0.5 * (1- hybc.tanh * hybc.tanh) * (dSddel * (1-r/ro));

    } else
	*dP_ddel = 0.;

}


/*****************************************************************************\
*  Type 2 hybrid:  Calculates dP_dri and dP_dro.                              *
*  Type 2 hybrids have an inner truncation (e.g. a ring) and no outer         *
*  truncation.                                                                *
\*****************************************************************************/


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

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

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

}


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

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

    void add_hybderivs (struct fitpars *fpar, struct derivs *df, 
        struct fitpars *hfpar, struct derivs *hdf, int *comp);

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

    fptr = fpar;
    dptr = df;

    while (fptr != NULL) {
	if (strncmp (fptr->objtype, "hybrid", 6) == 0) {
	    add_hybderivs (fpar, df, fptr, dptr, fptr->outer);
	    add_hybderivs (fpar, df, fptr, dptr, fptr->inner);
	};

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

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

void add_hybderivs (struct fitpars *fpar, struct derivs *df, 
    struct fitpars *hfpar, struct derivs *hybdf, int *comp)
{
    int i, j, k, ix, iy, done;
    struct fitpars *fptr, *fptr2;
    struct derivs *dptr, *dptr2, *hybdf2;

    /*---------------------------------------------------\
    |  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.                   |
	|  Hybrid 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    |
        |  hybrid parameters.  The net derivative image for the hybrid  |
        |  parameter is then assembled together below in hybdf->dpm.    |
	\--------------------------------------------------------------*/
	       
	fptr = fptr->high;
	dptr = dptr->high;

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

	    /*---------------------------------------------------------------\
	    |  Search for the hybrid 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 hybrid structure       |
	    |  component is the same as the number stored in fptr->hybnum,   |
	    |  which was set in "link_fpar_high.                             |
	    \---------------------------------------------------------------*/

	    if (strncmp (fptr->objtype, "hybrid", 6) == 0 && 
		fptr->hybnum == hfpar->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)
		                hybdf->dpm[k][iy][ix] += dptr->dpm[k][iy][ix];
			};

			/*-------------------------------------------\
			|  Now the fourier, diskyness/boxyness, or   |
			|  bending params.                           |
			\-------------------------------------------*/

			if (fptr->high != NULL) {
			    fptr2 = fptr->high;
			    dptr2 = dptr->high;
			    hybdf2 = hybdf->high;

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

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

	    };

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


/*****************************************************************************\
*  Calculates df/dr and df/dI * P, where P is the hyperbolic tangent          *
*  function.  Hybrid profile are multiplied by a hyperbolic tangent.          *
*                                                                             *
*  Note that hout (outer taper) and hin (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 hout and hin     *
*  confused with a[3] or a[4].                                                *
\*****************************************************************************/

void hyb_dfdx_calc (float *dfdIP, float *dfdr_P1, float *dfdr_P2, int hybrid, 
    struct hybval *hybc1, struct hybval *hybc2, float hout_am, float hin_am,
    struct hyblinks *hout, struct hyblinks *hin, float dfdI, float flux0, 
    float hybnorm)
{

    float P1, P2, dP1dam, dP2dam;


    if (hybrid == 1) {             		/*  Outer truncation  */

	if (hout->type == 1) {      /*  Break radius and softening *length*  */

            hybc1->tanh = tanh(hybc1->gtr * hout_am/hout->a[3] + hybc1->s);
            P1 = 0.5*(hybc1->tanh + 1); 
            *dfdIP = dfdI * (1-P1);

            dP1dam = 0.5 * (1-hybc1->tanh * hybc1->tanh)*hybc1->gtr/hout->a[3];
            *dfdr_P1 = -flux0 * dP1dam;

	} else {                    /*  Break radius and softening *radius*  */

            hybc1->tanh = tanh(hybc1->gtr * hout_am/hout->a[4] + hybc1->s);
            P1 = 0.5*(hybc1->tanh + 1);
            *dfdIP = dfdI * (1-P1);

            dP1dam = 0.5 * (1-hybc1->tanh * hybc1->tanh)*hybc1->gtr / hout->a[4];
            *dfdr_P1 = -flux0 * dP1dam;
	};

    } else if (hybrid == 2) {      		/*  Inner truncation  */

	if (hin->type == 1) {	    /*  Break radius and softening *length*  */

            hybc2->tanh = tanh(hybc2->gtr * hin_am/hin->a[3] + hybc2->s);
            P2 = 0.5*(hybc2->tanh + 1);
            *dfdIP = dfdI * P2 / hybnorm;

            dP2dam = 0.5* (1-hybc2->tanh * hybc2->tanh) * hybc2->gtr/hin->a[3];
            *dfdr_P2 = flux0 * dP2dam;

	} else {                    /*  Break radius and softening *radius*  */

            hybc2->tanh = tanh(hybc2->gtr * hin_am/hin->a[4] + hybc2->s);
            P2 = 0.5*(hybc2->tanh + 1);
            *dfdIP = dfdI * P2 / hybnorm;

            dP2dam = 0.5 * (1-hybc2->tanh * hybc2->tanh) * hybc2->gtr / hin->a[4];
            *dfdr_P2 = flux0 * dP2dam;

	};

    } else if (hybrid == 3) { 		     /* Both inner & outer trans.  */

	if (hout->type == 1) {

            hybc1->tanh = tanh(hybc1->gtr * hout_am/hout->a[3] + hybc1->s);
            P1 = 0.5*(hybc1->tanh+1);
            dP1dam = 0.5 * (1-hybc1->tanh * hybc1->tanh)*hybc1->gtr/hout->a[3];
            *dfdr_P1 = -flux0 * dP1dam;

	} else {

            hybc1->tanh = tanh(hybc1->gtr * hout_am/hout->a[4] + hybc1->s);
            P1 = 0.5*(hybc1->tanh+1);
            dP1dam = 0.5 * (1-hybc1->tanh * hybc1->tanh)*hybc1->gtr / hout->a[4];
            *dfdr_P1 = -flux0 * dP1dam;

	}

	if (hin->type == 1) {

            hybc2->tanh = tanh(hybc2->gtr * hin_am/hin->a[3] + hybc2->s);
            P2 = 0.5*(hybc2->tanh+1);
            dP2dam = 0.5 * (1-hybc2->tanh * hybc2->tanh)*hybc2->gtr /hin->a[3];
            *dfdr_P2 = flux0 * dP2dam;

	} else {

            hybc2->tanh = tanh(hybc2->gtr * hin_am/hin->a[4] + hybc2->s);
            P2 = 0.5*(hybc2->tanh+1);
            dP2dam = 0.5 * (1-hybc2->tanh * hybc2->tanh)*hybc2->gtr / hin->a[4];
            *dfdr_P2 = flux0 * dP2dam;

	};

        *dfdIP = dfdI * (P2 - P1) / hybnorm;
    };

}