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

/***************************************************************************\
*  Calculate truncation amplitude and gradient parameters, i.e. the         *
*  hyperbolic tangent transition function and the its derivative with       *
*  respect to the radial coordinate.                                        *
\***************************************************************************/

void trunc_tanh_calc (float *P1, float *dP1dam, float *P2, float *dP2dam, 
    int ttype, struct trunc_val *trunc_c1, struct trunc_val *trunc_c2, 
    float trunc_out_am, float trunc_in_am, struct trunc_links *trunc_out, 
    struct trunc_links *trunc_in)
{

    if (ttype == 1 || ttype == 3) {           		/*  Outer truncation  */

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

            trunc_c1->tanh = tanh(trunc_c1->gtr * trunc_out_am/(trunc_out->a[3]+
		trunc_out->a[4]) + trunc_c1->s);
            *P1 = 0.5*(trunc_c1->tanh + 1); 
	    *P1 = (1 - *P1);

            *dP1dam = 0.5 * (1-trunc_c1->tanh * trunc_c1->tanh)*trunc_c1->gtr/
		     (trunc_out->a[3] + trunc_out->a[4]);
	    *dP1dam = -*dP1dam;
	} else {                    /*  Softening *radius* as free param.  */

            trunc_c1->tanh = tanh(trunc_c1->gtr * trunc_out_am/
			     trunc_out->a[4] + trunc_c1->s);
            *P1 = 0.5*(trunc_c1->tanh + 1);
	    *P1 = (1 - *P1);

            *dP1dam = 0.5 * (1-trunc_c1->tanh * trunc_c1->tanh)*trunc_c1->gtr / 
		     trunc_out->a[4];
	    *dP1dam = -*dP1dam;
	};
    };

    if (ttype == 2 || ttype == 3) {     		/*  Inner truncation  */

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

            trunc_c2->tanh = tanh(trunc_c2->gtr * trunc_in_am/trunc_in->a[3] + 
			     trunc_c2->s);
            *P2 = 0.5*(trunc_c2->tanh + 1);

            *dP2dam = 0.5* (1-trunc_c2->tanh * trunc_c2->tanh) * trunc_c2->gtr/
		     trunc_in->a[3];
	} else {                    /*  Softening *radius* as free param.  */

            trunc_c2->tanh = tanh(trunc_c2->gtr * trunc_in_am/trunc_in->a[4] + 
			     trunc_c2->s);
            *P2 = 0.5*(trunc_c2->tanh + 1);

            *dP2dam = 0.5 * (1-trunc_c2->tanh * trunc_c2->tanh) * 
		     trunc_c2->gtr / trunc_in->a[4];
	};
    };


/* 

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

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

            trunc_c1->tanh = tanh(trunc_c1->gtr * trunc_out_am/
			     (trunc_out->a[3] + trunc_out->a[4]) + trunc_c1->s);
            *P1 = 0.5*(trunc_c1->tanh+1);
            *dP1dam = 0.5 * (1-trunc_c1->tanh * trunc_c1->tanh)*trunc_c1->gtr/
		     trunc_out->a[3];
	    *dP1dam = -*dP1dam;
	} else {                      /*  Softening *radius* as free param.  

            trunc_c1->tanh = tanh(trunc_c1->gtr * trunc_out_am/
			     trunc_out->a[4] + trunc_c1->s);
            *P1 = 0.5*(trunc_c1->tanh+1);
            *dP1dam = 0.5 * (1-trunc_c1->tanh * trunc_c1->tanh)*trunc_c1->gtr / 
		     trunc_out->a[4];
	    *dP1dam = -*dP1dam;
	}

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

            trunc_c2->tanh = tanh(trunc_c2->gtr * trunc_in_am/trunc_in->a[3] + 
			     trunc_c2->s);
            *P2 = 0.5*(trunc_c2->tanh+1);
            *dP2dam = 0.5 * (1-trunc_c2->tanh * trunc_c2->tanh)*trunc_c2->gtr /
		     trunc_in->a[3];
	} else {                      /*  Softening *radius* as free param. 

            trunc_c2->tanh = tanh(trunc_c2->gtr * trunc_in_am/trunc_in->a[4] + 
			     trunc_c2->s);
            *P2 = 0.5*(trunc_c2->tanh+1);
            *dP2dam = 0.5 * (1-trunc_c2->tanh * trunc_c2->tanh)*trunc_c2->gtr / 
		     trunc_in->a[4];
	};
    };

*/

}