#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>

typedef struct sparse_mat_t
{
        double *a;
        int *col;
        int *start_row;
        int n;
        struct sparse_mat_t *l;
        struct sparse_mat_t *u;

} sparse_mat;


typedef struct contact_t
{
        int id;                     /* Identifies contact type */

        double xc;                  /* X coordinate of contact */
        double yc;                  /* Y coordinate of contact */
        double zc;                  /* Z coordinate of contact */
        double dc;                  /* Dendritic diameter at contact */

        double xp;                  /* Projected X coordinate */
        double yp;                  /* Projected Y coordinate */
        double zp;                  /* Projected Z coordinate */

        double sd;                  /* Shortest distance to dendrite (micron) */
        double pl;                  /* Measurement of physical length (micron) */
        double amp;                 /* Strength of contact */

        int xl;                     /* Left hand node */
        int xr;                     /* Right hand node */
        double frac;                /* Fraction of input to left hand node */

        struct contact_t *next;     /* Address of next contact */
} contact;


typedef struct soma_t
{
/*  Physical properties of soma */
        int nobs;                   /* No. of observations in somal specification */

        double *x;                  /* X-coords of defining point */
        double *y;                  /* Y-coords of defining point */
        double *z;                  /* Z-coords of defining point */
        double *d;                  /* Diameter of soma at point (x,y,z) */
        double p_len;               /* Length of soma */

/*  Static biophysical properties of soma */
        double cs;                  /* Somal membrane capacitance (mu F/cm^2) */
        double ga;                  /* Intracellular conductance of soma (mS/cm) */
        double gs;                  /* Membrane conductance of soma (mS/cm^2) */

/*  Contact information */
        contact *conlist;           /* List of contacts */
        int ncon;
} soma;


typedef struct branch_t
{
/*  Connectivity of branch */
        struct branch_t *parent;    /* Pointer to parent branch */
        struct branch_t *child;     /* Pointer to child branch */
        struct branch_t *peer;      /* Pointer to a peer branch */

/*  Physical properties of branch */
        int nobs;                   /* No. of observations in branch specification */
        double *x;                  /* X-coordinate of defining point */
        double *y;                  /* Y-coordinate of defining point */
        double *z;                  /* Z-coordinate of defining point */
        double *d;                  /* Diameter of dendrite */
        double p_len;               /* Length of branch */
        double *pl;                 /* Measurement of physical length (micron) */

/*  Biophysical properties of branch */
        double cm;                  /* Dendritic membrane capacitance (mu F/cm^2) */
        double ga;                  /* Intracellular conductance of dendrite (mS/cm) */
        double gm;                  /* Membrane conductance of dendrite (mS/cm^2) */

/*  Node information for spatial representation */
        int nodes;                  /* Total number nodes spanning branch */
        int junct;                  /* Junction node of the branch */
        int first;                  /* Internal node connected to junction */

/*  Contact information */
        contact *conlist;           /* List of contacts */
        int ncon;
} branch;



typedef struct dendrite_t
{
        branch *root;               /* Pointer to root branch of dendrite */
        double p_len;               /* Total length of dendrite */
} dendrite;


typedef struct neuron_t
{
        int ndend;                  /* Number of dendrites */
        dendrite *dendlist;         /* Pointer to an array of dendrites */
        soma *s;                    /* Soma structure */
} neuron;


/* Function type declarations */
neuron *Load_Sampled_Neuron(char *);
neuron *Load_Test_Neuron(char *);
void Destroy_Sampled_Neuron(neuron *);
void init_branch( branch *, int, double, double, double);
void BuildContactInfo(contact *, branch *, branch **);
void remove_branch( branch **, branch *);
void build_dendrite( branch **, branch *);
void clean_dendrite( branch *);
void destroy_dendrite( branch *);
int count_branches( branch *, branch *);
double branch_length( branch *, branch *);
int count_terminal_branches( branch *, branch *);
void enumerate_nodes( branch *, int *);
void cgs(int, sparse_mat *, double *, double *, double *, double);
void mat_vec_mult( sparse_mat *, double *, double *);
void mat_malloc( sparse_mat *, int, int);
void mat_free( sparse_mat *);
void Generate_Dendrite(branch *, int *);
void init_test_branch( branch *, int , double , double , double );
void build_test_dendrite( branch **, branch *);
void LU_Factor(sparse_mat *, int *);
void LU_Solve( sparse_mat *, double *, double *);
void Input_Current( branch *);
void Assign_Current( branch *, double *, double );
void count_branch ( branch *, int * );
void count_dendritic_length ( branch *, double *);
void phys_lengths( branch *, double * );

/* Global definitions */
#define             CS    1.0
#define             GS    0.091
#define             GA    14.286
#define             CM    1.0
#define             GM    0.091
#define         OUTPUT    "contactinfo.dat"
#define           TEND    10.0
#define             NT    1000
#define             DT    1.0
#define         NNODES    100

/* Parameters for exact solution */
#define            AMP    1.0e-3
#define            SIN    0.0e-3
#define              T    10.0
#define              M    1000

/* Parameters for TestCell1.d3 (unbranched dendrite) */
// #define             RS    0.002
// #define             RD    0.001
// #define              L    0.100
// #define              X    0.055

/* Parameters for TestCell3.d3 (Y-junction) */
// #define             RS    0.002
// #define             RD    0.000354487542
// #define              L    0.0500427736
// #define              X    0.039163520

/* Parameters for TestCell2.d3 (Y-junction) */
#define             RS    0.002
#define             RD    0.000354487542
#define              L    0.0667236981
#define              X    0.029735419

/* Parameters for TestCell4.d3 (large branched dendrite) */
// #define             RS    0.002
// #define             RD    0.000648741708
// #define              L    0.2256605
// #define              X    0.0135280571

/* Global Variables */
sparse_mat lhs, rhs;

int main( int argc, char **argv )
{
    int k, id, start, nodes, nc, i, in, FirstNode;
    int counter, nb;
    double *v, *x, max, *eta, as, gama, *chi, xold, xnew,
           frac, arg, sum, tmp, vs, pi, dt, tnow, tout,
           length, h;
    neuron *n;
    extern sparse_mat lhs, rhs;
    branch *bnow;
    FILE *fp;

/*  Compute ancillary variables */
    pi = 4.0*atan(1.0);
    frac = X/L;

/*  Load sampled neuron */
    if ( argc != 2 && argc != 3 ) {
        printf("\n Invoke program with load <input>\n");
        return(1);
    } else if ( argc == 2 ) {
        n = Load_Sampled_Neuron( argv[1] );
        if ( !n ) {
            printf("\n Failed to find sampled neuron\n");
            return(1);
        }
    } else if ( strcmp(argv[1], "test") == 0  || strcmp(argv[1], "Test") == 0 ) {
        n = Load_Test_Neuron( argv[2] );
        nb = 0;
        length = 0.0;

        for(i=0; i<n->ndend; i++) count_branch( n->dendlist[i].root, &nb );
        for(i=0; i<n->ndend; i++) count_dendritic_length( n->dendlist[i].root, &length );

        h = length/((double) NNODES - nb);
        for(i=0; i<n->ndend; i++) phys_lengths( n->dendlist[i].root, &h );

        if ( !n ) {
            printf("\n Failed to find sampled neuron\n");
            return(1);
        }
        as = 4.0*pi*RS*RS;
        gama = 0.5*as/(pi*RD*L);
        eta = (double *) malloc( (M+1)*sizeof(double) );
        chi = (double *) malloc( (M+1)*sizeof(double) );
        eta[0] = GM/CM;
        chi[0] = 0.5*(AMP+SIN)/(pi*CM*RD*L*(1.0+gama));
        chi[0] /= eta[0];
        for ( k=1 ; k<=M ; k++ ) {
            xnew = arg = pi*((double) k );
            do {
                xold = xnew;
                xnew = arg-atan(gama*xold);
            } while ( fabs(xold-xnew) > 5.e-7 );
            tmp = cos(xnew);
            eta[k] = (GM+0.5*RD*GA*xnew*xnew/(L*L))/CM;
            chi[k] = tmp*(AMP*cos(xnew*(1.0-frac))+SIN*tmp)/
                         (pi*CM*RD*L*(1.0+gama*tmp*tmp));
            chi[k] /= eta[k];
        }
    } else {
        printf("\n Failed to find sampled neuron\n");
        return(1);
    }

/*  Enumerate Nodes */

    FirstNode = 0;
    for ( k=0 ; k<n->ndend ; k++ ) enumerate_nodes( n->dendlist[k].root, &FirstNode );
    for ( k=0 ; k<n->ndend ; k++ ) n->dendlist[k].root->junct = FirstNode;
    printf("Number of nodes is %d\n", FirstNode+1);
    getchar( );

/* Construct Sparse Matrices */

    counter = 0;
    nodes = FirstNode+1;
    mat_malloc( &lhs, nodes, 3*nodes-2 );
    mat_malloc( &rhs, nodes, 3*nodes-2 );
    lhs.n = rhs.n = nodes;
    lhs.a[3*nodes-3] = rhs.a[3*nodes-3] = 0.0;
    lhs.start_row[0] = rhs.start_row[0] = 0;

    for ( k=0 ; k<n->ndend ; k++ ) {
        bnow = n->dendlist[k].root;
        Generate_Dendrite( bnow, &counter);
    }

    for ( k=0 ; k<n->ndend ; k++ ) {
        bnow = n->dendlist[k].root;
        lhs.a[counter] = 0.5*bnow->d[1]*bnow->pl[1];
        rhs.a[counter] = -(bnow->d[0])*(bnow->d[1])/bnow->pl[1];
        lhs.col[counter] = rhs.col[counter] = bnow->first;
        lhs.a[3*nodes-3] += 1.5*(bnow->d[0])*bnow->pl[1];
        rhs.a[3*nodes-3] += (bnow->d[0])*(bnow->d[1])/bnow->pl[1];
        counter++;
    }
    lhs.col[3*nodes-3] = rhs.col[3*nodes-3] = nodes-1;
    lhs.start_row[nodes] = rhs.start_row[nodes] = 3*nodes-2;

    for( i=0; i<n->ndend ; i++ ) Input_Current(n->dendlist[i].root);

    fp = fopen("out.res", "w");
    fclose(fp);

    dt = 1.0/((double) NT);
    for ( i=0 ; i<3*nodes-2 ; i++ ) {
        rhs.a[i] *= GA;
        rhs.a[i] += GM*lhs.a[i];
        lhs.a[i] *= CM;
        rhs.a[i] *= 0.5*dt;
        lhs.a[i] += rhs.a[i];
        rhs.a[i] = lhs.a[i] - 2.0*rhs.a[i];
    }

/*  Add capacitive term of soma */

    lhs.a[3*nodes-3] += (4.0*as/pi)*(CS+0.5*GS*dt);
    rhs.a[3*nodes-3] += (4.0*as/pi)*(CS-0.5*GS*dt);

    v = (double *) malloc( (nodes)*sizeof(double) );
    x = (double *) malloc( (nodes)*sizeof(double) );
    for ( i=0 ; i<nodes ; i++ ) v[i] = x[i] = 0.0;
    start = 1;
    LU_Factor(&lhs, &start);

/*  Initialise temporal integration */

    tnow = 0.0;
    tout = DT;
    while ( tnow < TEND ) {
        tnow += dt;
        mat_vec_mult(&rhs,v,x);
        x[nodes-1] -= 4.0*dt*SIN/pi;
        tmp = 4.0*dt/pi;
        for ( i=0 ; i<nodes ; i++ ) v[i] = x[i];
        if ( tnow < T ) {
            for ( i=0 ; i<n->ndend ; i++ ) Assign_Current(n->dendlist[i].root, x, tmp);
        }
        LU_Solve( &lhs, v, x );

        if ( tnow > tout && tnow <= T ) {
            printf("\rReached time %5.1lf ms\t", tout);
            for ( vs=0.0,i=M ; i>=0 ; i-- ) {
                arg = tnow*eta[i];
                if ( arg > 20.0 ) {
                    tmp = 1.0;
                } else {
                    tmp = (1.0-exp(-arg));
                }
                vs -= tmp*chi[i];
            }
            printf("\nNumerical Voltage %12.6lf mV\n",v[nodes-1]);
            printf("Exact Voltage %12.6lf mV\n", vs);
            tout += DT;
            fp = fopen("out.res", "a");
            fprintf(fp,"Error %12.6lf% \n",fabs((v[nodes-1]-vs)/(0.01*vs)));
            fclose(fp);
        }
        if ( tnow > tout && tnow > T ) {
            printf("\rReached time %5.1lf ms\t", tout);
            for ( vs=0.0,i=M ; i>=0 ; i-- ) {
                arg = (tnow-T)*eta[i];
                if ( arg < 20.0 ) {
                    tmp = exp(-arg);
                    arg = tnow*eta[i];
                    if ( arg < 20.0 ) tmp -= exp(-arg);
                    vs -= tmp*chi[i];
                }
            }
            printf("\nNumerical Voltage %12.6lf mV\n",v[nodes-1]);
            printf("Exact Voltage %12.6lf mV\n",vs);
            tout += DT;
            fp = fopen("out.res", "a");
            fprintf(fp,"Error %12.6lf% \n",fabs((v[nodes-1]-vs)/(0.01*vs)));
            fclose(fp);
        }
    }

    Destroy_Sampled_Neuron( n );
    return(0);
}

 /******************************************************
        Function to constuct sparse matrices
  ******************************************************/

void Generate_Dendrite( branch *bnow, int *counter)
{
    int i, k;
    extern sparse_mat lhs, rhs;
    branch *btmp;
    double SumL, SumR;

/* Step 1 - Recurse to the end of the dendrite */

    if ( bnow->child != NULL ) Generate_Dendrite( bnow->child, counter);
    if ( bnow->peer != NULL ) Generate_Dendrite( bnow->peer, counter);

    for ( k=bnow->first-bnow->nobs+2,i=bnow->nobs-1 ; i>0 ; i--,k++ ) {

/* Step 2 - Fill in matrix entries for terminal points */

        if ( bnow->child == NULL && i == bnow->nobs - 1 ) {

            lhs.a[*counter] = 1.5*(bnow->d[i])*(bnow->pl[i]-bnow->pl[i-1]);
            rhs.a[*counter] = (bnow->d[i-1]*bnow->d[i])/(bnow->pl[i]-bnow->pl[i-1]);
            lhs.col[*counter] = rhs.col[*counter] = k;
            (*counter)++;

            lhs.a[*counter] = 0.5*(bnow->d[i-1])*(bnow->pl[i] - bnow->pl[i-1]);
            rhs.a[*counter] = -(bnow->d[i-1]*bnow->d[i])/(bnow->pl[i] - bnow->pl[i-1]);

            if ( k == bnow->first ) {
                lhs.col[*counter] = rhs.col[*counter] = bnow->junct;
            } else {
                lhs.col[*counter] = rhs.col[*counter] = k + 1;
            }

            (*counter)++;
            lhs.start_row[k+1] = rhs.start_row[k+1] = *counter;

/* Step 3 - Fill in matrix entries for branch points */

        } else if ( bnow->child != NULL && i == bnow->nobs - 1 ) {

            btmp = bnow->child;
            SumR = SumL = 0.0;

            while ( btmp != NULL ) {

                lhs.a[*counter] = 0.5*btmp->d[1]*btmp->pl[1];
                rhs.a[*counter] = -(btmp->d[0]*btmp->d[1])/btmp->pl[1];
                lhs.col[*counter] = rhs.col[*counter] = btmp->first;
                (*counter)++;
                SumL += 1.5*btmp->d[0]*btmp->pl[1];
                SumR += (btmp->d[0]*btmp->d[1])/btmp->pl[1];
                btmp = btmp->peer;
            }

            lhs.a[*counter] = SumL+1.5*bnow->d[i]*(bnow->pl[i]-bnow->pl[i-1]);
            rhs.a[*counter] = SumR+(bnow->d[i-1]*bnow->d[i])/(bnow->pl[i]-bnow->pl[i-1]);
            lhs.col[*counter] = rhs.col[*counter] = k;
            (*counter)++;

            lhs.a[*counter] = 0.5*(bnow->d[i-1])*(bnow->pl[i]-bnow->pl[i-1]);
            rhs.a[*counter] = -(bnow->d[i-1]*bnow->d[i])/(bnow->pl[i] - bnow->pl[i-1]);

            if ( k == bnow->first ) {
                lhs.col[*counter] = rhs.col[*counter] = bnow->junct;
            } else {
                lhs.col[*counter] = rhs.col[*counter] = k + 1;
            }

            (*counter)++;
            lhs.start_row[k+1] = rhs.start_row[k+1] = *counter;

        } else {

    /* Step 4 - Fill in matrix entries for internal point */

            lhs.a[*counter] = 0.5*(bnow->d[i+1])*(bnow->pl[i+1] - bnow->pl[i]);
            rhs.a[*counter] = -(bnow->d[i]*bnow->d[i+1])/(bnow->pl[i+1] - bnow->pl[i]);
            lhs.col[*counter] = rhs.col[*counter] = k - 1;
            (*counter)++;

            lhs.a[*counter] = 1.5*(bnow->d[i])*(bnow->pl[i+1] - bnow->pl[i-1]);
            rhs.a[*counter] = (bnow->d[i-1]*bnow->d[i])/(bnow->pl[i] - bnow->pl[i-1])
                            + (bnow->d[i]*bnow->d[i+1])/(bnow->pl[i+1] - bnow->pl[i]);
            lhs.col[*counter] = rhs.col[*counter] = k;
            (*counter)++;

            lhs.a[*counter] = 0.5*(bnow->d[i-1])*(bnow->pl[i] - bnow->pl[i-1]);
            rhs.a[*counter] = -(bnow->d[i-1]*bnow->d[i])/(bnow->pl[i] - bnow->pl[i-1]);
            lhs.col[*counter] = rhs.col[*counter] = k + 1;

            if ( k == bnow->first ) {
                lhs.col[*counter] = rhs.col[*counter] = bnow->junct;
            } else {
                lhs.col[*counter] = rhs.col[*counter] = k + 1;
            }

            (*counter)++;

            lhs.start_row[k+1] = rhs.start_row[k+1] = *counter;
        }
    }
    return;
}


neuron *Load_Sampled_Neuron(char *filename)
{
    int j, k, ncon, n, id, connected, ignored;
    double tmp, piby2, xold, xnew, yold, ynew, zold, znew, diam,
           xl, xr, dl, dr, px, py, pz, min;
    neuron *cell;
    soma *s;
    contact *oldcon, *newcon, *firstcon;
    branch *bold, *bnew, *first_branch, *bopt;
    char temp[100];
    FILE *input;

/*  STEP 1. - Open neuron data file */
    printf("\nOpening file %s\n",filename);
    if ( (input=fopen(filename,"r"))==NULL ) return NULL;

/*  STEP 2. - Get memory for neuron structure */
    cell = (neuron *) malloc( sizeof(neuron) );
    s = cell->s = (soma *) malloc( sizeof(soma) );

/*  STEP 3, - Get information about soma */
    fscanf(input,"%s",temp);
    fscanf(input,"%d",&n);
    printf("Identified a %s defined by %d data\n", temp, n);

/*  STEP 4. - Initialise soma structure */
    s->x = (double *) malloc( n*sizeof(double) );
    s->y = (double *) malloc( n*sizeof(double) );
    s->z = (double *) malloc( n*sizeof(double) );
    s->d = (double *) malloc( n*sizeof(double) );
    s->gs = GS;
    s->ga = GA;
    s->cs = CS;
    s->conlist = NULL;
    s->ncon = 0;
    printf("%s initialised\n", temp);

/*  STEP 5. - Get soma morphological data */
    fscanf(input,"%lf %lf %lf %lf",&xold, &yold, &zold, &diam);
    s->x[0] = xold;
    s->y[0] = yold;
    s->z[0] = zold;
    s->d[0] = diam;
    for ( k=j=1,s->p_len=0.0 ; k<n ; k++ ) {
        fscanf(input,"%lf %lf %lf %lf",&xnew, &ynew, &znew, &diam);
        tmp = pow(xnew-xold,2)+pow(ynew-yold,2)+pow(znew-zold,2);
        if ( tmp > 0.01 ) {
            s->x[j] = xold = xnew;
            s->y[j] = yold = ynew;
            s->z[j] = zold = znew;
            s->d[j] = diam;
            j++;
            s->p_len += sqrt(tmp);
        }
    }
    s->nobs = j;
    printf("Finished soma\n");

/*  STEP 6. - Get branch and contact data */
    firstcon = oldcon = NULL;
    bold = NULL;
    while  ( fscanf(input,"%s", temp)!=EOF ) {
        if ( strcmp(temp, "Branch") == 0 || strcmp(temp, "branch") == 0 ) {
            fscanf(input, "%d", &n);
            printf("Found a branch defined by %d nodes\n", n);
            bnew = (branch *) malloc( sizeof(branch) );
            if ( bold != NULL) {
                bold->child = bnew;
            } else {
                first_branch = bnew;
            }
            bnew->parent = bold;
            bnew->peer = NULL;
            bnew->child = NULL;

/*  STEP 7. - Initialise branch */
            init_branch( bnew, n, CM, GM, GA);
            fscanf(input,"%lf %lf %lf %lf", &xold, &yold, &zold, &diam);
            bnew->x[0] = xold;
            bnew->y[0] = yold;
            bnew->z[0] = zold;
            bnew->d[0] = diam;
            for ( bnew->p_len=bnew->pl[0]=0.0,k=j=1 ; k<n ; k++ ) {
                fscanf(input,"%lf %lf %lf %lf", &xnew, &ynew, &znew, &diam);
                tmp = pow(xnew-xold,2)+pow(ynew-yold,2)+pow(znew-zold,2);
                if ( tmp > 0.01 ) {
                    bnew->p_len += sqrt(tmp);
                    bnew->pl[j] = bnew->p_len;
                    bnew->x[j] = xold = xnew;
                    bnew->y[j] = yold = ynew;
                    bnew->z[j] = zold = znew;
                    bnew->d[j] = diam;
                    j++;
                }
            }
            bnew->nobs = j;
            bold = bnew;
        } else if ( strcmp(temp, "Marker") == 0 || strcmp(temp, "marker") == 0 ) {

/*  STEP 8. - Initialise marker */
            fscanf(input, "%d %d", &id, &n);
            printf("Found and initialised %d branch contacts of type %d\n", n, id);
            for ( k=0 ; k<n ; k++ ) {
                newcon = (contact *) malloc( sizeof(contact) );
                newcon->next = NULL;
                newcon->sd = NULL;
                newcon->id = id;
                if ( oldcon != NULL ) {
                    oldcon->next = newcon;
                } else {
                    firstcon = newcon;
                }
                fscanf(input,"%lf %lf %lf %lf", &newcon->xc, &newcon->yc, &newcon->zc, &newcon->dc );
                oldcon = newcon;
            }
        } else {
            printf("Unknown block type %s!\n", temp);
            return NULL;
        }
    }
    fclose(input);

/*  STEP 9. - Associate contacts with branches and soma */
    ignored = 0;
    while ( firstcon ) {
        bold = first_branch;
        bopt = NULL;
        while ( bold ) {
            BuildContactInfo( firstcon, bold, &bopt);
            bold = bold->child;
        }
        newcon = firstcon->next;
        if ( firstcon->sd > 4.0 ) {

/*  STEP 10a. - Check for proximity to soma */
            px = firstcon->xc;
            py = firstcon->yc;
            pz = firstcon->zc;

/*  STEP 1. - First stage is different from others */
            xnew = s->x[0]; ynew = s->y[0]; znew = s->z[0];
            firstcon->sd = min = sqrt(pow(xnew-px,2)+pow(ynew-py,2)+pow(znew-pz,2))-(s->d[0]);

/*  STEP 2. - Second stage compares points and projected points */
            for ( k=1 ; k<s->nobs ; k++ ) {
                xnew = s->x[k]; ynew = s->y[k]; znew = s->z[k];
                min = sqrt(pow(xnew-px,2)+pow(ynew-py,2)+pow(znew-pz,2))-(s->d[k]);
                if ( min < firstcon->sd ) firstcon->sd = min;
            }
            if ( firstcon->sd < 4.0 ) {
                oldcon = s->conlist;
                if ( oldcon ) {
                    while ( oldcon->next ) oldcon = oldcon->next;
                    oldcon->next = firstcon;
                } else {
                    s->conlist = firstcon;
                }
                firstcon->next = NULL;
                s->ncon++;
            } else {
                free(firstcon);
                ignored++;
            }
        } else {

/*  STEP 3. - Add contact to conlist of optimal branch */
            oldcon = bopt->conlist;
            if ( oldcon ) {
                while ( oldcon->next ) oldcon = oldcon->next;
                oldcon->next = firstcon;
            } else {
                bopt->conlist = firstcon;
            }
            firstcon->next = NULL;
            bopt->ncon++;
        }
        firstcon = newcon;
    }
    printf("\nIgnored %d contacts", ignored);

/*  STEP 11. - Count dendritic branches at soma */
    bold = first_branch;
    n = 0;
    while ( bold ) {
        bnew = first_branch;
        do {
            k = bnew->nobs-1;
            tmp = pow(bold->x[0]-bnew->x[k],2)+
                  pow(bold->y[0]-bnew->y[k],2)+
                  pow(bold->z[0]-bnew->z[k],2);
            connected = ( tmp < 0.01 );
            bnew = bnew->child;
        } while ( bnew && !connected );
        if ( !connected ) n++;
        bold = bold->child;
    }
    cell->ndend = n;
    printf("\n\nTree contains %d individual dendrite(s) ...\n", n);

/* STEP 12. - Identify somal dendrites but extract nothing */
    cell->dendlist = (dendrite *) malloc( (cell->ndend)*sizeof(dendrite) );
    bold = first_branch;
    n = 0;
    while ( n < cell->ndend ) {
        bnew = first_branch;
        do {
            k = bnew->nobs-1;
            tmp = pow(bold->x[0]-bnew->x[k],2)+
                  pow(bold->y[0]-bnew->y[k],2)+
                  pow(bold->z[0]-bnew->z[k],2);
            connected = ( tmp < 0.01 );
            bnew = bnew->child;
        } while ( bnew && !connected );
        if ( !connected ) {
            cell->dendlist[n].root = bold;
            n++;
        }
        bold = bold->child;
    }

/*  STEP 13. - Extract root of each dendrite from dendrite list */
    for ( k=0 ; k<cell->ndend ; k++ ) {
        bold = cell->dendlist[k].root;
        remove_branch( &first_branch, bold);
    }

/*  STEP 14. - Build each dendrite from its root branch */
    for ( k=0 ; k<cell->ndend ; k++ ) {
        build_dendrite( &first_branch, cell->dendlist[k].root);
        clean_dendrite( cell->dendlist[k].root);
    }
    if ( first_branch != NULL ) printf("\nWarning: Unconnected branch segments still exist\n");
    return cell;
}

/****************************************************************

                Function to initialise a BRANCH

****************************************************************/

void init_branch( branch *b, int n, double cm, double gm, double ga)
{

/*  Allocate memory for spatial orientation of dendrite */
    b->x = (double *) malloc( n*sizeof(double) );
    b->y = (double *) malloc( n*sizeof(double) );
    b->z = (double *) malloc( n*sizeof(double) );

/*  Allocate memory for branch geometry */
    b->d = (double *) malloc( n*sizeof(double) );
    b->pl = (double *) malloc( n*sizeof(double) );

/*  Set parameter values */
    b->cm = cm;
    b->gm = gm;
    b->ga = ga;

/*  Initialise node information */
    b->nodes = 0;
    b->first = 0;
    b->junct = 0;

/*  Initialise contact information */
    b->conlist = NULL;
    b->ncon = 0;
    return;
}

/**************************************************************

                Function to initialise a TEST BRANCH

***************************************************************/

void init_test_branch( branch *b, int n, double cm, double gm,double ga)
{

/*  Allocate memory for spatial orientation of dendrite */
    b->x = (double *) malloc( 2*sizeof(double) );
    b->y = (double *) malloc( 2*sizeof(double) );
    b->z = (double *) malloc( 2*sizeof(double) );

/*  Allocate memory for branch geometry */
    b->d = (double *) malloc( n*sizeof(double) );
    b->pl = (double *) malloc( n*sizeof(double) );

/*  Set parameter values */
    b->cm = cm;
    b->gm = gm;
    b->ga = ga;

/*  Initialise node information */
    b->nodes = 0;
    b->first = 0;
    b->junct = 0;

/*  Initialise contact information */
    b->conlist = NULL;
    b->ncon = 0;
    return;
}

/***************************************************************

             Function to build CONTACT information

***************************************************************/

void BuildContactInfo(contact *con, branch *b, branch **bopt)
{
    int k;
    double px, py, pz, tmp, xold, xnew, yold, ynew, zold, znew,
           numer, denom, xmin, ymin, zmin, min;

    px = con->xc;
    py = con->yc;
    pz = con->zc;

/*  STEP 1. - First stage is different from others */
    xnew = b->x[0]; ynew = b->y[0]; znew = b->z[0];
    min = sqrt(pow(xnew-px,2)+pow(ynew-py,2)+pow(znew-pz,2));
    if ( con->sd == NULL || ( con->sd != NULL && min < con->sd ) ) {
        con->sd = min;
        con->xp = xnew; con->yp = ynew; con->zp = znew;
        con->pl = 0.0;
        *bopt = b;
    }

/*  STEP 2. - Second stage compares points and projected points */
    for ( k=1 ; k<b->nobs ; k++ ) {

        xold = xnew; yold = ynew; zold = znew;
        xnew = b->x[k]; ynew = b->y[k]; znew = b->z[k];
        numer = (xnew-xold)*(px-xold)+(ynew-yold)*(py-yold)+(znew-zold)*(pz-zold);
        denom = pow(xnew-xold,2)+pow(ynew-yold,2)+pow(znew-zold,2);

/*  STEP 2a. - Project onto branch */
        if ( 0.0 <= numer && numer <= denom ) {
            tmp = numer/denom;
            xmin = (1.0-tmp)*xold+tmp*xnew;
            ymin = (1.0-tmp)*yold+tmp*ynew;
            zmin = (1.0-tmp)*zold+tmp*znew;
            min = sqrt(pow(xmin-px,2)+pow(ymin-py,2)+pow(zmin-pz,2));
            if ( !(con->sd) || ( con->sd && min < con->sd ) ) {
                con->sd = min;
                con->xp = xmin; con->yp = ymin; con->zp = zmin;
                con->pl = (1.0-tmp)*b->pl[k-1]+tmp*b->pl[k];
                *bopt = b;
            }
        }

/*  STEP 2b. - Check proximity to points of branch */
        min = sqrt(pow(xnew-px,2)+pow(ynew-py,2)+pow(znew-pz,2));
        if ( !(con->sd) || ( con->sd && min < con->sd ) ) {
            con->sd = min;
            con->xp = xnew; con->yp = ynew; con->zp = znew;
            con->pl = b->pl[k];
            *bopt = b;
        }
    }
    return;
}

/***************************************************************

          Function to remove a branch from a branch list

****************************************************************/

void remove_branch(branch **head, branch *b)
{
    if ( *head == NULL || b == NULL ) return;
    if ( *head == b ) {
        *head = b->child;
        if ( *head != NULL )  (*head)->parent = NULL;
    } else {
        b->parent->child = b->child;
        if ( b->child != NULL ) b->child->parent = b->parent;
    }
    b->parent = NULL;
    b->child = NULL;
    return;
}

/***************************************************************

          Function to build a dendrite from its root

***************************************************************/

void build_dendrite( branch **head, branch *root)
{
    int k;
    double tmp;
    branch *bnow, *bnext, *btmp;

    bnow = *head;
    while ( bnow != NULL ) {

/* Store bnow's child in case it's corrupted */
        bnext = bnow->child;

/* Search if proximal end of bnow is connected to distal end of root */
        k = (root->nobs)-1;
        tmp = pow(bnow->x[0]-root->x[k],2)+pow(bnow->y[0]-root->y[k],2)+
              pow(bnow->z[0]-root->z[k],2);
        if ( tmp <= 0.01 ) {

/* Take bnow out of the branch list */
            remove_branch( head, bnow);

/* Connect bnow to the root as the child or a peer of the child.
   Initialise childs' children and peers to NULL as default */
            bnow->child = NULL;
            bnow->peer = NULL;
            bnow->parent = root;

/* Inform root about its child if it's the first child, or add
   new child to first child's peer list */
            if ( root->child != NULL ) {
                btmp = root->child;
                while ( btmp->peer != NULL ) btmp = btmp->peer;
                btmp->peer = bnow;
            } else {
                root->child = bnow;
            }
        }

/* Initialise bnow to next branch in list */
        bnow = bnext;
    }

/* Iterate through remaining tree */
    if ( root->child ) build_dendrite( head, root->child);
    if ( root->peer ) build_dendrite( head, root->peer);
    return;
}

/***************************************************************

          Function to remove peerless children

***************************************************************/
void clean_dendrite( branch *root)
{
    int k, np, nc, mem, n;
    double tmp, sarea;
    contact *con;
    branch *btmp, *brem;

/*  Iterate through remaining tree */
    if ( root->child != NULL ) clean_dendrite( root->child );
    if ( root->peer != NULL ) clean_dendrite( root->peer );

/*  Extend original parent limb */
    brem = root->child;
    if ( brem != NULL && brem->peer == NULL ) {
        root->child = brem->child;
        if ( brem->child ) {
            brem->child->parent = root;
            btmp = brem->child->peer;
            while ( btmp ) {
                btmp->parent = root;
                btmp = btmp->peer;
            }
        }
//        root->nodes += (brem->nodes)-1;
        np = root->nobs;
        nc = brem->nobs;
        mem = np+nc-1;
        root->nobs = mem;
        root->x = (double *) realloc( (void *)root->x, mem*sizeof(double) );
        for ( k=np ; k<mem ; k++ ) root->x[k] = brem->x[k-np+1];
        root->y = (double *) realloc( (void *)root->y, mem*sizeof(double) );
        for ( k=np ; k<mem ; k++ ) root->y[k] = brem->y[k-np+1];
        root->z = (double *) realloc( (void *)root->z, mem*sizeof(double) );
        for ( k=np ; k<mem ; k++ ) root->z[k] = brem->z[k-np+1];
        root->d = (double *) realloc( (void *)root->d, mem*sizeof(double) );
        for ( k=np ; k<mem ; k++ ) root->d[k] = brem->d[k-np+1];
        root->pl = (double *) realloc( (void *)root->pl, mem*sizeof(double) );
        for ( k=np ; k<mem ; k++ ) root->pl[k] = root->p_len+brem->pl[k-np+1];
        root->ncon += brem->ncon;
        con = root->conlist;
        if ( con ) {
            while ( con->next ) con = con->next;
            con->next = brem->conlist;
        } else {
            root->conlist = brem->conlist;
        }
        brem->conlist = NULL;
        free(brem->x);
        free(brem->y);
        free(brem->z);
        free(brem->d);
        free(brem->pl);
        free ( brem );
    }
    return;
}

/***************************************************************

                Function to destroy a NEURON

****************************************************************/

void Destroy_Sampled_Neuron(neuron *cell)
{
    int i;
    contact *prevcon, *nextcon;

/*  Free Soma */

    if ( cell->s != NULL) {
        if ( cell->s->x ) free ( cell->s->x );
        if ( cell->s->y ) free ( cell->s->y );
        if ( cell->s->z ) free ( cell->s->z );
        if ( cell->s->d ) free ( cell->s->d );
        prevcon = cell->s->conlist;
        while ( prevcon ) {
            nextcon = prevcon->next;
            free ( prevcon );
            prevcon = nextcon;
        }
        free ( cell->s );
    }
    for ( i=0 ; i<cell->ndend ; i++ ) destroy_dendrite( cell->dendlist[i].root );
    free(cell);
    return;
}

/***************************************************************

               Function to destroy DENDRITE

****************************************************************/

void destroy_dendrite( branch *b )
{
    int i;
    contact *prevcon, *nextcon;

    if ( b->child ) destroy_dendrite(b->child);
    if ( b->peer ) destroy_dendrite(b->peer);
    free(b->x);
    free(b->y);
    free(b->z);
    free(b->d);
    free(b->pl);
    prevcon = b->conlist;
    while ( prevcon ) {
        nextcon = prevcon->next;
        free ( prevcon );
        prevcon = nextcon;
    }
    free ( b );
    return;
}

/***************************************************************

             Function to count number of branches

***************************************************************/

int count_branches( branch *bstart, branch *bnow)
{
    static int n;

    if ( bstart == bnow ) n = 0;
    if ( bnow ) {
        if ( bnow->child ) count_branches(bstart, bnow->child);
        if ( bnow->peer ) count_branches(bstart, bnow->peer);
        n++;
    }
    return n;
}

/***************************************************************

             Function to find length of dendrite from
                    current branch to tips.

****************************************************************/

double branch_length( branch *bstart, branch *bnow)
{
    static double length;

    if ( bstart == bnow ) length = 0.0;
    if ( bnow ) {
        if ( bnow->child ) branch_length(bstart, bnow->child);
        if ( bnow->peer ) branch_length(bstart, bnow->peer);
        length += bnow->p_len;
    }
    return length;
}

/****************************************************************

        Function to count number of terminal branches

****************************************************************/

int count_terminal_branches( branch *bstart, branch *bnow)
{
    static int n;

    if ( bstart == bnow ) n = 0;
    if ( bnow ) {
        if ( bnow->child ) count_terminal_branches(bstart, bnow->child);
        if ( bnow->peer ) count_terminal_branches(bstart, bnow->peer);
        if ( !bnow->child ) n++;
    }
    return n;
}

/**********************************************************

        Function to enumerate the nodes on a dendrite

**********************************************************/

void enumerate_nodes(branch *bnow, int *FirstNode )
{
    branch *btmp;

    if ( (bnow->child) != NULL ) enumerate_nodes( bnow->child, FirstNode );
    if ( (bnow->peer) != NULL ) enumerate_nodes( bnow->peer, FirstNode );

    if ( bnow->child != NULL ) {
        btmp = bnow->child;
        while( btmp != NULL ){
            btmp->junct = *FirstNode;
            btmp = btmp->peer;
        }
    }

    bnow->first = *FirstNode + bnow->nobs - 2;
    *FirstNode += (bnow->nobs)-1;
    return;
}


/*************************************************************

    Multiplies sparse matrix a[ ][ ] with vector v[ ]

*************************************************************/

void mat_vec_mult( sparse_mat *a, double *v , double *b)
{
    int i, j, ntop, n;

    n = a->n;
    for ( j=0 ; j<n ; j++) {
        ntop = a->start_row[j+1];
        for( b[j]=0.0,i=(a->start_row[j]) ; i<ntop ; i++ ) {
            b[j] += (a->a[i])*v[a->col[i]];
        }
    }
    return;
}


/*******************************************************************

    Allocates memory to a sparse matrix - returns 0 if successful

*******************************************************************/

void mat_malloc( sparse_mat *a, int n, int w)
{
    a->a = (double *) malloc( w*sizeof(double) );
    a->col = (int *) malloc( w*sizeof(int) );
    a->start_row = (int *) malloc( (n+1)*sizeof(int) );
    a->n = n;
    a->l = malloc(sizeof(sparse_mat));
    a->u = malloc(sizeof(sparse_mat));
    a->l->a = (double *) malloc( (2*n-1)*sizeof(double) );
    a->l->col = (int *) malloc( (2*n-1)*sizeof(int) );
    a->l->start_row = (int *) malloc( (n+1)*sizeof(int) );
    a->l->n = n;
    a->u->a = (double *) malloc( (2*n-1)*sizeof(double) );
    a->u->col = (int *) malloc( (2*n-1)*sizeof(int) );
    a->u->start_row = (int *) malloc( (n+1)*sizeof(int) );
    a->u->n = n;

    if ( !a->start_row ) exit(1);

}


/* De-allocates memory of a sparse matrix */
void mat_free( sparse_mat *a)
{
    free(a->a);
    free(a->col);
    free(a->start_row);
    free(a);
}



 /***************************************************************
            Function To Load A Test Neuron
 ***************************************************************/
neuron *Load_Test_Neuron(char *filename)
{
    int j, k, ncon, n, id, connected, ignored;
    double tmp, piby2, xl, xr, dl, dr, px, py, pz, min, radius, dx;
    neuron *cell;
    contact *newcon;
    branch *bold, *bnew, *first_branch, *bopt;
    char temp[100];
    FILE *input;

/*  STEP 1. - Open neuron data file */
    printf("\nOpening file %s\n",filename);
    if ( (input=fopen(filename,"r"))==NULL ) return NULL;

/*  STEP 2. - Get memory for neuron structure */
    cell = (neuron *) malloc( sizeof(neuron) );

/*  STEP 3. - Get branch and contact data */
    bold = NULL;
    while  ( fscanf(input,"%s", temp)!=EOF ) {
        if ( strcmp(temp, "Branch") == 0 || strcmp(temp, "branch") == 0 ) {
            fscanf(input, "%d", &n);
            printf("Found a branch defined by %d nodes\n", n);
            bnew = (branch *) malloc( sizeof(branch) );
            if ( bold != NULL) {
                bold->child = bnew;
            } else {
                first_branch = bnew;
            }
            bnew->parent = bold;
            bnew->peer = NULL;
            bnew->child = NULL;

/*  STEP 4. - Initialise branch */
            init_test_branch( bnew, n, CM, GM, GA);
            fscanf(input,"%lf %lf %lf", &(bnew->x[0]), &(bnew->y[0]), &(bnew->z[0]));
            fscanf(input,"%lf %lf %lf", &(bnew->x[1]), &(bnew->y[1]), &(bnew->z[1]));
            fscanf(input,"%lf %lf", &(bnew->p_len), &radius);
            dx = bnew->p_len / ((double) n-1 );
            for ( j=0; j<n; j++) {
                bnew->pl[j] = dx*((double) j );
                bnew->d[j] = 2.0*radius;
            }
            bnew->nobs = n;
            bold = bnew;
        } else if ( strcmp(temp, "Marker") == 0 || strcmp(temp, "marker") == 0 ) {

/*  STEP 5. - Initialise marker */
            printf("Found and initialised a branch contact\n");
            newcon = (contact *) malloc( sizeof(contact) );
            newcon->next = NULL;
            newcon->sd = NULL;
            fscanf(input,"%lf %lf", &newcon->pl, &newcon->amp );
            bnew->conlist = newcon;
        }
    }
    fclose(input);

/*  STEP 6. - Count dendritic branches at soma */
    bold = first_branch;
    n = 0;
    while ( bold ) {
        bnew = first_branch;
        do {
            tmp = pow(bold->x[0]-bnew->x[1],2)+
                  pow(bold->y[0]-bnew->y[1],2)+
                  pow(bold->z[0]-bnew->z[1],2);
            connected = ( tmp < 0.01 );
            bnew = bnew->child;
        } while ( bnew && !connected );
        if ( !connected ) n++;
        bold = bold->child;
    }
    cell->ndend = n;
    printf("\n\nTree contains %d individual dendrite(s) ...\n", n);

/* STEP 7. - Identify somal dendrites but extract nothing */
    cell->dendlist = (dendrite *) malloc( (cell->ndend)*sizeof(dendrite) );
    bold = first_branch;
    n = 0;
    while ( n < cell->ndend ) {
        bnew = first_branch;
        do {
            tmp = pow(bold->x[0]-bnew->x[1],2)+
                  pow(bold->y[0]-bnew->y[1],2)+
                  pow(bold->z[0]-bnew->z[1],2);
            connected = ( tmp < 0.01 );
            bnew = bnew->child;
        } while ( bnew && !connected );
        if ( !connected ) {
            cell->dendlist[n].root = bold;
            n++;
        }
        bold = bold->child;
    }

/*  STEP 8. - Extract root of each dendrite from dendrite list */
    for ( k=0 ; k<cell->ndend ; k++ ) {
        bold = cell->dendlist[k].root;
        remove_branch( &first_branch, bold);
    }

/*  STEP 9. - Build each test dendrite from its root branch */
    for ( k=0 ; k<cell->ndend ; k++ ) {
        build_test_dendrite( &first_branch, cell->dendlist[k].root);
    }
    if ( first_branch != NULL ) printf("\nWarning: Unconnected branch segments still exist\n");
    return cell;
}



/**************************************************************

        Function to build a test dendrite from its root

***************************************************************/
void build_test_dendrite( branch **head, branch *root)
{
    double tmp;
    branch *bnow, *bnext, *btmp;

    bnow = *head;
    while ( bnow != NULL ) {

/* Store bnow's child in case it's corrupted */
        bnext = bnow->child;

/* Search if proximal end of bnow is connected to distal end of root */
        tmp = pow(bnow->x[0]-root->x[1],2)+pow(bnow->y[0]-root->y[1],2)+
              pow(bnow->z[0]-root->z[1],2);
       if ( tmp <= 0.01 ) {

/* Take bnow out of the branch list */
            remove_branch( head, bnow);

/* Connect bnow to the root as the child or a peer of the child.
   Initialise childs' children and peers to NULL as default */
            bnow->child = NULL;
            bnow->peer = NULL;
            bnow->parent = root;

/* Inform root about its child if it's the first child, or add
   new child to first child's peer list */
            if ( root->child != NULL ) {
                btmp = root->child;
                while ( btmp->peer != NULL ) btmp = btmp->peer;
                btmp->peer = bnow;
            } else {
                root->child = bnow;
            }
        }

/* Initialise bnow to next branch in list */
        bnow = bnext;
    }

/* Iterate through remaining tree */
    if ( root->child ) build_test_dendrite( head, root->child);
    if ( root->peer ) build_test_dendrite( head, root->peer);
    return;
}

 /***************************************************************
            Function To Factorise A Sparse Matrix
 ***************************************************************/
void LU_Factor(sparse_mat *m, int *start)
{
    double tmp, sum;
    int i, j, k, r, n, cl, cu, col, row;

/*  Step 1. - Identify matrix dimension */
    n = m->n;

/*  Step 2. - Fill column vectors for triangular matrices */
    if ( *start ) {
        cl = cu = 0;
        for ( i=k=0 ; i<n ; i++ ) {
            m->l->start_row[i] = cl;
            m->u->start_row[i] = cu;
            while ( m->col[k] < i ) m->l->col[cl++] = m->col[k++];
            m->l->col[cl++] = m->col[k];
            m->u->col[cu++] = m->col[k++];
            while ( k < m->start_row[i+1] ) m->u->col[cu++] = m->col[k++];
        }
        m->l->start_row[n] = cl;
        m->u->start_row[n] = cu;
        *start = 0;
    }

/*  Step 3. - Fill remaining entries of L and U row by row */
    cl = cu = 0;
    for ( i=0 ; i<n ; i++ ) {
        for ( k=m->start_row[i] ; k < m->start_row[i+1] ; k++ ) {
            if ( m->col[k] < i ) {
                sum = m->a[k];
                for ( j=m->l->start_row[i] ; j<cl ; j++ ) {
                    col = m->l->col[j];
                    row = m->u->start_row[col];
                    while ( m->u->col[row] < m->col[k] && row < m->u->start_row[col+1] ) row++;
                    if ( m->u->col[row] == m->col[k] ) sum -= (m->l->a[j])*(m->u->a[row]);
                }
                row = m->u->start_row[m->l->col[cl]];
                while ( m->u->col[row] < m->col[k] ) row++;
                m->l->a[cl++] = sum/(m->u->a[row]);
            } else if ( m->col[k] == i ) {
                m->l->a[cl++] = 1.0;
                sum = m->a[k];
                for ( j=m->l->start_row[i] ; j<m->l->start_row[i+1]-1; j++ ) {
                    col = m->l->col[j];
                    row = m->u->start_row[col];
                    while ( m->u->col[row] < i && row < m->u->start_row[col+1] ) row++;
                    if ( m->u->col[row] == i ) sum -= (m->l->a[j])*(m->u->a[row]);
                }
                m->u->a[cu++] = sum;
            } else {
                sum = m->a[k];
                for ( j=m->l->start_row[i] ; j<m->l->start_row[i+1]-1; j++ ) {
                    col = m->l->col[j];
                    row = m->u->start_row[col];
                    while ( m->u->col[row] < m->col[k] && row < m->u->start_row[col+1] ) row++;
                    if ( m->u->col[row] == m->col[k] ) sum -= (m->l->a[j])*(m->u->a[row]);
                }
                m->u->a[cu++] = sum;
            }
        }
    }
    return;
}

/*******************************************************************

          Function to Solve the matrix problem

*******************************************************************/

void LU_Solve(sparse_mat *m, double *x, double *b )
{
    int i,j;
    double *z;

    z = (double *) malloc( (m->n)*sizeof(double) );

    for ( i=0 ; i<m->n ; i++ ) {
        z[i] = b[i];
        for (j=m->l->start_row[i];j<m->l->start_row[i+1]-1;j++)
        { z[i] -= (m->l->a[j])*(z[(m->l->col[j])]); }
        z[i] /= m->l->a[m->l->start_row[i+1]-1];
    }

    for ( i = (m->n) - 1 ; i>=0 ; i-- ) {
        x[i] = z[i];
        for (j=m->u->start_row[i]+1;j<m->u->start_row[i+1];j++)
        { x[i] -= (m->u->a[j])*(x[m->u->col[j]]); }
        x[i] /= m->u->a[m->u->start_row[i]];
    }
    free(z);
    return;
}



 /*************************************************************
            Function to input current to dendrite
 **************************************************************/
void Input_Current( branch *bnow )
{
    int k;
    double tmp;

    if ( bnow->child != NULL ) Input_Current( bnow->child );
    if ( bnow->peer != NULL ) Input_Current( bnow->peer );

    if ( bnow->conlist != NULL ) {
        if ( bnow->conlist->pl <= bnow->pl[1] ) {
            bnow->conlist->xl = bnow->junct;
            bnow->conlist->xr = bnow->first;
            bnow->conlist->frac = (1.0-(bnow->conlist->pl/bnow->pl[1]));
        } else {
            k = 1;
            while ( bnow->conlist->pl > bnow->pl[k] ) k++;
            bnow->conlist->xl = bnow->first-k+2;
            bnow->conlist->xr = bnow->first-k+1;
            bnow->conlist->frac = (bnow->pl[k]-bnow->conlist->pl)/(bnow->pl[k]-bnow->pl[k-1]);
        }
    }
    return;
}

 /************************************************************
                 Function to assign current
 **************************************************************/
void Assign_Current(branch *bnow, double *x, double fac )
{
    if (bnow->child != NULL ) Assign_Current(bnow->child, x, fac );
    if (bnow->peer != NULL ) Assign_Current(bnow->peer, x, fac );

    if ( bnow->conlist != NULL ) {
        x[bnow->conlist->xl] -= fac*(bnow->conlist->frac)*(bnow->conlist->amp);
        x[bnow->conlist->xr] -= fac*(1.0-(bnow->conlist->frac))*(bnow->conlist->amp);
    }
    return;
}

 /*********************************************************************
                    Function to count number of branches
 *********************************************************************/

void count_branch ( branch *b, int *nb )
{
    (*nb)++;
    if ( b->child != NULL ) count_branch ( b->child, nb );
    if ( b->peer != NULL ) count_branch ( b->peer, nb );
    return;
}

 /*********************************************************************
            Function to count total dendritic length
 *********************************************************************/

void count_dendritic_length ( branch *b, double *length )
{
    (*length) += b->p_len;
    if ( b->child != NULL ) count_dendritic_length ( b->child, length );
    if ( b->peer != NULL ) count_dendritic_length ( b->peer, length );
    return;
}

 /********************************************************************
            Function to re-compute physical lengths
 ********************************************************************/

void phys_lengths( branch *b, double *h )
{
    int i,n,k;
    double tmp, hnow, *dtmp, *ptmp, *ptr;


    n = ((int) ceil((b->p_len)/(*h))) + 1;
    hnow = b->p_len/((double) n-1);
    dtmp = (double *) malloc( n*sizeof(double));
    ptmp = (double *) malloc( n*sizeof(double));
    for ( i=0; i<(n-1); i++ ) ptmp[i] = ((double) i)*hnow;
    ptmp[n-1] = b->p_len;

    dtmp[0] = b->d[0];
    for ( i=1 ; i<(n-1) ; i++ ) {
        k = 0;
        while ( ptmp[i] > b->pl[k] ) k++;
        dtmp[i] = b->d[k-1] + (b->d[k] - b->d[k-1])*(ptmp[i] - b->pl[k-1])/(b->pl[k]-b->pl[k-1]);
    }
    dtmp[n-1] = b->d[b->nobs-1];

    ptr = b->d;
    b->d = dtmp;
    free(ptr);
    ptr = b->pl;
    b->pl = ptmp;
    free(ptr);

    b->nobs = n;

    if (b->child != NULL) phys_lengths( b->child, h);
    if (b->peer != NULL) phys_lengths( b->peer, h);
}