typedef struct cond_t {
        int n;               /* No. timesteps in conductance profile */
        double dt;           /* Integration time step (msecs) */
        double *g;           /* Conductance profile for prescribed dt (mu a) */
        double tol;          /* Determines off-threshold for synapse */
        double tau;          /* Rise time for synapse (msecs) */
        double gmax;         /* Peak conductance per synapse (mu) */
} cond;

typedef struct synapse_t
{
        cond *cp;            /* Pointer to profile of synaptic conductance */
        int ind;             /* Indicator of status of synapse */
        int max;             /* Maximum value of status counter */
        double gold;         /* Previous synaptic conductance */
        double vsyn;         /* Reversal potential of synapse (mV) */
} synapse;

void cgs(int, sparse_mat *, double *, double *, double *, double),
     FreeConductanceProfile( cond *),
     UpdateSynapticConductance( synapse *syn);
cond *ConductanceProfile( double, double, double, double);
synapse *SynapticProfile( int, int, double, double, double, cond * );
double c_in(synapse *, long int);

/* Declaration of HH coefficient functions */
double alfa_h( double );
double alfa_m( double );
double alfa_n( double );
double beta_h( double );
double beta_m( double );
double beta_n( double );

double temperature_correction;


 /********************************************
        Computes a conductance profile
  ********************************************/
cond *ConductanceProfile( double dt,   /* Integration time step */
                          double tau,  /* Synaptic time constant */
                          double tol,  /* Determines off-threshold for synapse */
                          double gmax  /* Maximum conductance (mS) */ )
{
    int k;
    cond *con;
    double tmp, told, tnew;

    con = (cond *) malloc( sizeof(cond) );
    con->dt = dt;
    con->tau = tau;
    con->tol = tol;
    con->gmax = gmax;

/* Iterate to find duration of pulse */
    tmp = 1.0-log(tol);
    tnew = tmp;
    do {
        told = tnew;
        tnew = tmp+log(told);
    } while ( fabs(tnew-told)>5.e-11 );
    con->n = ((int) ceil(tau*tnew/dt));
    con->g = (double *) malloc( (con->n)*sizeof(double) );
    con->g[0] = 0.0;
    for ( k=1 ; k<(con->n) ; k++ ) {
        tmp = dt*((double) k)/tau;
        con->g[k] = gmax*tmp*exp(1.0-tmp);
    }
    return con;
}


 /********************************************
     Function to free a conductance profile
  ********************************************/
void FreeConductanceProfile( cond *con)
{
    free(con->g);
    free(con);
    return;
}


 /*******************************************************
    Function to initialise and update synaptic activity
  *******************************************************/
void UpdateSynapticConductance( synapse *syn)
{
    extern unsigned long int ix, iy, iz;
    static double RealRate;
    static int start=1;
    int j, k;
    double deadtime, dt, interval;

/*  Step 1. - Initialise synapse */
    if ( start ) {
        if ( syn->cp ) {
            dt = (syn->cp->dt);
            syn->max = syn->cp->n;
            syn->gold = 0.0;
            deadtime = ((double) syn->cp->n)*dt;
            RealRate = (1000.0/RATE)-deadtime;
            if ( RealRate <= 0.0 ) {
                printf("\nImpossible firing rate requested\n");
                fflush(stdout);
                return NULL;
            }
            interval = -RATE*log(ran( &ix, &iy, &iz))/dt;
            if ( fmod(interval,1.0) <= 0.5 ) {
               syn->ind = -((int) floor(interval));
            } else {
               syn->ind = -((int) ceil(interval));
            }
        } else {
            printf("\nSynapse has no conductance profile!\n");
            return NULL;
        }
        start =0;
    } else {
        (syn->ind)++;
        if ( syn->ind == syn->max ) {
            dt = (syn->cp->dt);
            syn->gold = 0.0;
            interval = -RealRate*log(ran( &ix, &iy, &iz))/dt;
            if ( fmod(interval,1.0) <= 0.5 ) {
               syn->ind = -((int) floor(interval));
            } else {
               syn->ind = -((int) ceil(interval));
            }
        }
    }
    return;
}



 /********************************************************************
                    ALPHA for ACTIVATION OF SODIUM
  *******************************************************************/
double alfa_m( double volt )
{
    extern double temperature_correction;
    double tmp;

    tmp = 2.5-0.1*volt;
    if ( fabs(tmp)<0.001 ) {
        tmp = 1.0/(((tmp/24.0+1.0/6.0)*tmp+0.5)*tmp+1.0);
    } else {
        tmp = tmp/(exp(tmp)-1.0);
    }
    return tmp*temperature_correction;
}


 /********************************************************************
                    BETA for ACTIVATION OF SODIUM
  ********************************************************************/
double beta_m( double volt )
{
    extern double temperature_correction;
    double tmp;

    tmp = volt/18.0;
    return 4.0*temperature_correction*exp(-tmp);
}


 /********************************************************************
                    ALPHA for INACTIVATION OF SODIUM
  ********************************************************************/
double alfa_h( double volt )
{
    extern double temperature_correction;
    double tmp;

    tmp = 0.05*volt;
    return 0.07*temperature_correction*exp(-tmp);
}


 /********************************************************************
                    BETA for INACTIVATION OF SODIUM
  ********************************************************************/
double beta_h( double volt )
{
    extern double temperature_correction;
    double tmp;

    tmp = 3.0-0.1*volt;
    return temperature_correction/(exp(tmp)+1.0);
}


 /********************************************************************
                    ALPHA for ACTIVATION OF POTASSIUM
  ********************************************************************/
double alfa_n( double volt )
{
    extern double temperature_correction;
    double tmp;

    tmp = 1.0-0.1*volt;
    if ( fabs(tmp)<0.001 ) {
        tmp = 0.1/(((tmp/24.0+1.0/6.0)*tmp+0.5)*tmp+1.0);
    } else {
        tmp = 0.1*tmp/(exp(tmp)-1.0);
    }
    return tmp*temperature_correction;
}


 /********************************************************************
                    BETA for ACTIVATION OF POTASSIUM
  ********************************************************************/
double beta_n( double volt )
{
    extern double temperature_correction;
    double tmp;

    tmp = 0.0125*volt;
    return 0.125*temperature_correction*exp(-tmp);
}



void cgs(int getmem, sparse_mat *a, double *b, double *x0, double *x, double tol)
{
    long int i, k, n, finish;
    static int start=1;
    double rho_old, rho_new, alpha, beta, sigma, err;
    static double *p, *q, *u, *v, *r, *rb, *y;

/* Step 1 - Check memory status */
    n = a->n;
    if ( start ) {
        r = (double *) malloc( n*sizeof(double) );
        rb = (double *) malloc( n*sizeof(double) );
        p = (double *) malloc( n*sizeof(double) );
        q = (double *) malloc( n*sizeof(double) );
        u = (double *) malloc( n*sizeof(double) );
        v = (double *) malloc( n*sizeof(double) );
        y = (double *) malloc( n*sizeof(double) );
        start = 0;
    }

/* Step 2 - Initialise residual, p[ ] and q[ ] */
    mat_vec_mult( a, x0, r);
    for ( rho_old=0.0,i=0 ; i<n ; i++ ) {
        r[i] = b[i]-r[i];
        rho_old += r[i]*r[i];
        rb[i] = r[i];
        p[i] = 0.0;
        q[i] = 0.0;
    }
    if ( rho_old<tol*((double) n) ) {
        for ( i=0 ; i<n ; i++ ) x[i] = x0[i];
        return;
    }
    rho_old = 1.0;
    finish = 0;

/* The main loop */
    while ( !finish ) {

/* Compute scale parameter for solution update */
        for ( rho_new=0.0,i=0 ; i<n ; i++ ) rho_new += r[i]*rb[i];
        beta = rho_new/rho_old;

/* Update u[ ] and p[ ] */
        for ( i=0 ; i<n ; i++ ) {
            u[i] = r[i]+beta*q[i];
            p[i] = u[i]+beta*(q[i]+beta*p[i]);
        }

/* Update v[ ] and compute sigma */
        mat_vec_mult( a, p, v);
        for ( sigma=0.0,i=0 ; i<n ; i++ ) sigma += rb[i]*v[i];

/* Compute alpha and update q[ ], v[ ] and x[ ] */
        alpha = rho_new/sigma;
        for ( i=0 ; i<n ; i++ ) {
            q[i] = u[i]-alpha*v[i];
            v[i] = alpha*(u[i]+q[i]);
            x[i] += v[i];
        }

 /* Update r[ ] and estimate error */
        mat_vec_mult( a, v, y);
        for ( err=0.0,i=0 ; i<n ; i++ ) {
            r[i] -= y[i];
            err += r[i]*r[i];
        }
        rho_old = rho_new;
        if ( err<tol*((double) n) ) finish = 1;
    }

/* Check memory status */
    if ( getmem<=0 ) start = 1;
    if ( start ) {
        free(r);
        free(rb);
        free(p);
        free(q);
        free(u);
        free(v);
        free(y);
    }
    return;
}

/**********************************************************
                Returns Gaussian deviate
  **********************************************************/
double normal( double mean, double sigma)
{
    static int start=1;
    long int n;
    static double g1, g2;
    double v1, v2, w, ran(long int *);

    if ( start ) {
        n = 1;
        do {
            v1 = 2.0*ran(&n)-1.0;
            v2 = 2.0*ran(&n)-1.0;
            w = v1*v1+v2*v2;
        } while ( w==0.0 || w>=1.0 );
        w = log(w)/w;
        w = sqrt(-w-w);
        g1 = v1*w;
        g2 = v2*w;
        start = !start;
        return (mean+sigma*g1);
    } else {
        start = !start;
        return (mean+sigma*g2);
    }
}

 /**********************************************************
         Order integers in ascending order by heapsort
  **********************************************************/
void heapsort( long int n, long int *x)
{
    int finish;
    long int i, ir, j, k, tmp;

    if ( n<2 ) return;
    k = n/2;
    ir = n-1;
    finish = 0;
    while ( !finish ) {
        if ( k>0 ) {
            tmp = x[--k];
        } else {
            tmp = x[ir];
            x[ir] = x[0];
            if ( --ir==0 ) {
                x[0] = tmp;
                finish = 1;
            }
        }
        i = k;
        j = 2*k+1;
        while ( j<=ir ) {
            if ( j<ir && x[j]<x[j+1] ) j++;
            if ( tmp<x[j] ) {
                x[i] = x[j];
                i = j;
                j = 2*j+1;
            } else {
                j = ir+1;
            }
        }
        x[i] = tmp;
    }
    return;
}