#include <stdio.h>
#include <math.h>
#include "mex.h"

double  *NFArray(int size)
{
	double	*p;
	p = (double *)mxCalloc(sizeof(*p),size);
	return p;
}

/*******************************************************/
/* Meddis_network91M                                   */
/* A multi-channel implementation of                   */
/* an inner hair cell model from (Meddis, 1991)        */
/* Free parametrization                                */
/*                                                     */
/*       By Aki Härmä 6. 5 1999  		       */
/*******************************************************/
void trans(double *sig, double fs, long int len, long int dim, 
	   double *OUT, double A, double B, double g, double y,
	   double l, double r, double x)
{
 long int t, j, ind;
 /* Constants */
 double m,h,dt,gdt,ydt,ldt,rdt,xdt,hdt; 
 /* States */
 double *c, *q, *w;
 /* Variables */
 double k, stA, replenish, eject, loss, reuptake, reprocess;
 /* From (Meddis, 1990) */
 /* A=2; B=5000; g=200; y=1.05; l=3500; x=66.31; r=6580; */
 m=1.0; h=50000;
 /* Allocate and initialize */
 c=NFArray(dim);  q=NFArray(dim);  w=NFArray(dim); 
 dt=1/fs; gdt = g*dt; ydt = y*dt; ldt = l*dt; rdt = r*dt;
 xdt = x*dt; hdt = h*dt;
 k=g*A/(A+B); c[0]=k*y*m/(y*(l+r)+k*l);q[0]=c[0]*(l+r)/k;w[0]=c[0]*r/x;
 for (j=1;j<dim;j++){c[j]=c[0];q[j]=q[0];w[j]=w[0];}

 for (t=0;t<len;t++)
   {
   for (j=0;j<dim;j++)
     {
     ind=t+j*len;
     stA=(sig[ind]+A>0)?sig[ind]+A:0.0;
     k=gdt*stA/(stA+B);
     replenish=(q[j]<m)?ydt*(m-q[j]):0.0;
     eject=k*q[j];
     loss=ldt*c[j];
     reuptake=rdt*c[j];
     reprocess=xdt*w[j];
     q[j]+=replenish-eject+reprocess;
     c[j]+=eject-loss-reuptake;
     w[j]+=reuptake-reprocess;
     OUT[ind]=c[j];
     }
   }
return;
}

/************************  MEXFUNCTION *****************************/
void mexFunction(
	int		nlhs,
	mxArray	*plhs[],
	int		nrhs,
	const mxArray    *prhs[]
	)
{
  long int len, dim, m, n;
  double *sig, *OUT, *params;
  double yy, fs;
  double A,B,g,y,l,r,x;
  
        sig = mxGetPr(prhs[0]);
      	len = (long int)mxGetM(prhs[0]);
	dim = (long int)mxGetN(prhs[0]);

        fs = mxGetScalar(prhs[1]); /* Sampling rate */
   if (nrhs==3) /* Custom parameters */ 
     {
	params=mxGetPr(prhs[2]); /* Parameters */	
	A=params[0];B=params[1];g=params[2];y=params[3];
	l=params[4];r=params[5];x=params[6];
     }
   else /* Use medium spont. rate model from (Meddis&Hewitt, 1991) */
     {
        A=5; B=300; g=2000; y=5.05; l=2500; x=66.31; r=6580; 
     }
/* 	mexPrintf("Dim: %d   Len: %d",dim,len); */

       plhs[0] = mxCreateDoubleMatrix(len, dim, mxREAL);
	OUT = mxGetPr(plhs[0]);
	if(len<dim) mexPrintf("More channels than samples?");
	else trans(sig,fs,len,dim,OUT,A,B,g,y,l,r,x);	
	return;
}