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

void constraints (int mfit, float atry[], float a[], int ia[], 
    struct cons *constr, float *aorig)

{
    float hibound, lobound;
    int i, j;

    while (constr != NULL && constr->par != NULL) {

	i = constr->par[1];

	/****************************************\
	*  Strong coupling....                   *
	\****************************************/

	if (constr->op == 'o') {
	    for (j=2; j <= constr->ncomp; j++) {
		if (ia[constr->par[j]] > 0 &&
		    ia[constr->par[1]] == 1)

		    atry[constr->par[j]] = atry[constr->par[1]] + 
				       constr->cval[j];
	    };
	};

	if (constr->op == 'r') {
	    for (j=2; j <= constr->ncomp; j++) {
		if (ia[constr->par[j]] > 0 &&
		    ia[constr->par[1]] == 1)

		    atry[constr->par[j]] = atry[constr->par[1]] * 
				       constr->cval[j];
	    };
	};

	/****************************************\
	*  Constrain deviation from input value  *
	\****************************************/

	if (constr->op == 0) {
	    lobound = constr->cval[1];
	    hibound = constr->cval[2];

	    if (ia[i] == 1) {
	        if (atry[i] > aorig[i] + hibound)
	            atry[i] = aorig[i] + hibound;

	        if (atry[i] < aorig[i] + lobound)
		    atry[i] = aorig[i] + lobound;
	    };
	} else if (constr->op == 't') {

	/*****************************************\
	*  Constrain parameter by absolute range  *
	\*****************************************/

	    lobound = constr->cval[1];
	    hibound = constr->cval[2];

	    if (ia[i] == 1) {
	        if (atry[i] > hibound)
	            atry[i] = hibound;

	        if (atry[i] < lobound)
		    atry[i] = lobound;
	    };

	} else if (constr->op == '-') {

	    /************************************************************\
 	    *  Constrain difference between parameters of two different  *
	    *  components to be within a range hibound and lobound.      *
	    \************************************************************/

	    j = constr->par[2];

	    hibound = constr->cval[2];
	    lobound = constr->cval[1];

	    if (ia[i] > 0 && ia[j] > 0) {
	        if (atry[i] - atry[j] >= hibound) {

	            if (fabs(atry[i] - a[i]) <= fabs(atry[j] - a[j]))
		        atry[i] = atry[j] + hibound;
		    else
		        atry[j] = atry[i] - hibound; 
	        } else if (atry[i] - atry[j] <= lobound) {

		    if (fabs(atry[i] - a[i]) <= fabs(atry[j] - a[j]))
		        atry[i] = atry[j] + lobound;
		    else
		        atry[j] = atry[i] - lobound; 
	        };
	    };

	} else if (constr->op == '/') {

	    /************************************************\
	    *  Constrain ratio between 2 components to       *
	    *  be within a range hibound and lobound         *
	    \************************************************/

	    j = constr->par[2];

	    hibound = constr->cval[2];
	    lobound = constr->cval[1];

	    if (ia[i] > 0 && ia[j] > 0) {

	        if (atry[i] / atry[j] >= hibound) {
		    if (fabs(atry[i] - a[i]) <= fabs(atry[j] - a[j]))
		        atry[i] = fabs(atry[j]) * hibound;
		    else
		        atry[j] = fabs(atry[i]) / hibound;
	        } else  if (atry[i] / atry[j] <= lobound) {
		    if (fabs(atry[i] - a[i]) <= fabs(atry[j] - a[j]))
		        atry[i] = fabs(atry[j]) * lobound;
		    else
		        atry[j] = fabs(atry[i]) / lobound;
	        };
	    };
	};

	constr = constr->next;

    };
}