#include <stdio.h>
#include <string.h>
#include "structs.h"

int par2num (char *par);

void strong_links (struct fitpars *fpar, struct derivs *df, 
    struct cons *constr)
{
    struct fitpars *fptr1, *fptr2;
    struct derivs *dptr1, *dptr2;
    int i, j, m, n, parnum;

    while (constr != NULL) {

	if (constr->op == 'o' || constr->op == 'r') {
            parnum = par2num (constr->parname);

	    fptr1 = fpar;
	    dptr1 = df;

	    /******************************************************\
	    *  fptr1 and dptr1 point at the group lead component.  *
	    \******************************************************/

	    for (i=1; i < constr->comp[1]; i++) {
		if (fptr1->next != NULL && dptr1->next != NULL) {
		    fptr1 = fptr1->next;
		    dptr1 = dptr1->next;
		} else
		    constr->comp[1] = -1;
	    };

	    if (strncmp(constr->parname, "f", 1) == 0) {

		fptr1 = fptr1->high;
		dptr1 = dptr1->high;
		while (fptr1 != NULL && dptr1 != NULL &&
		       strncmp (fptr1->objtype, constr->parname, 1) != 0) {

		    fptr1 = fptr1->next;
		    dptr1 = dptr1->next;
		};
	    };

	    /****************************************************************\
	    *  fptr2 and dptr2 are trailing components.  Cycle them through  *
	    *  the constraint component list constr->comp[2]...[n]           *
	    \****************************************************************/

	    for (i=2; i <= constr->ncomp; i++) {
	        fptr2 = fpar;
	        dptr2 = df;

	        for (j=1; j < constr->comp[i]; j++) {
		    if (fptr1->next != NULL && dptr2->next != NULL) {
		        fptr2 = fptr2->next;
		        dptr2 = dptr2->next;
		    } else
			constr->comp[i] = -1;
		};
		
	        if (strncmp(constr->parname, "f", 1) == 0) {

		    fptr2 = fptr2->high;
		    dptr2 = dptr2->high;
		    while (fptr2 != NULL && 
		           strncmp (fptr2->objtype, constr->parname, 1) != 0) {

		        fptr2 = fptr2->next;
		        dptr2 = dptr2->next;
		    };
	        };

		for ( n=1; n <= dptr1->naxes[2]; n++) {
		    for ( m=1; m <= dptr1->naxes[1]; m++) 
		        if (fptr1->ia[parnum] > 0 && fptr2->ia[parnum] > 0)
		            dptr1->dpm[parnum][n][m] += dptr2->dpm[parnum][n][m];
		};

	    };

	};

        constr = constr->next;
    };

}