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

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

/*******************************************************/
/* Meddis_network91                                    */
/* A multi-channel implementation of                   */
/* an inner hair cell model from (Meddis, 1991)        */
/*       By Aki Härmä 16. 12 1998  		       */
/*******************************************************/
void trans(double *sig, double fs, long int len, long int dim, 
	   double *OUT)
{
 long int t, j, ind;
 /* Constants */
 double A,B,g,y,l,x,r,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=5; B=300; g=2000; y=5.05; l=2500; 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;
  double yy, fs;
  
        sig = mxGetPr(prhs[0]);
      	len = (long int)mxGetM(prhs[0]);
	dim = (long int)mxGetN(prhs[0]);

        fs = mxGetScalar(prhs[1]); /* Sampling rate */
	
/* 	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);	
	return;
}