#include "mex.h"
#include "of_time_interp.c"
/*--------------------------------------------------------------------*/
/* LININTERP   */
/*--------------------------------------------------------------------*/



/* the mex function -----------------------------------------------*/
void mexFunction (int nlhs, mxArray *plhs[],
                  int nrhs, const mxArray *prhs[]) {
char           *buf ;


  if (nrhs==0)  {
  double *a;
  plhs[0] = mxCreateDoubleMatrix (1, 1, mxREAL);
  a = mxGetPr (plhs[0]);   *a=5.0001;
  return; 
  }

buf = mxArrayToString( prhs[0] );

/*----------------------------------------------------------- cvs */
if (!strcmp("cvs",buf))  {

  plhs[0]=mxCreateString("$Revision: 1.18 $  $Date: 2006/01/16 17:31:29 $");

}
/*----------------------------------------------------------- LinInterp */
else if (!strcmp("lininterp",buf))  {

  /* out = of_time ('lininterp',table,val,last) */

  double    *out, *table, *val, *last;
  int       Nc, M, Npoints;


/*
table=[   1     1     2
          2     2     4
          3     3     6];
val=2.5;
last=[2 1 .5]; % index starts at 0
out = of_time ('lininterp',table,val,last)
last

of_time('lininterp',[0 0;1 1],linspace(0,1,10)',[1 0 0])
out=of_time('lininterp',[0;1;2],linspace(0,2,10)',[0 0 0])


*/

  if (nrhs!=4)  mexErrMsgTxt ("4 inputs required.");

  /* inputs */
  M     = mxGetM(prhs[1]);   /* number of points */
  Nc    = mxGetN(prhs[1])-1; /* number of curves */
  table = mxGetPr(prhs[1]);  /* [xi yi ...] */

  if ((mxGetM(prhs[2])==0) || (mxGetN(prhs[2])==0)) {
   plhs[0] = mxCreateDoubleMatrix (0, 0, mxREAL); return;
  }

  val   = mxGetPr(prhs[2]);  /*  xs         */
  last  = mxGetPr(prhs[3]);  /* [i1 xi si  ] */
  Npoints= mxGetM(prhs[2]);
  if (mxGetN(prhs[2])>Npoints) Npoints=mxGetN(prhs[2]);

  if (nlhs>0) {
   if (Nc>0) { plhs[0] = mxCreateDoubleMatrix (Npoints, Nc, mxREAL);
   } else {    plhs[0] = mxCreateDoubleMatrix (Npoints, 1, mxREAL);  }
   out= mxGetPr(plhs[0]);
  } else {out=NULL;}

  of_time_LinInterp(table,val,last,out,M,Nc,Npoints,nlhs);


}
/*----------------------------------------------------------- storelaststep */
else if (!strcmp("storelaststep",buf))  {
  /*
   of_time('storelaststep',uva,u,v,a);
  */
double     *uva, *u, *v, *a;
int        N;

N=mxGetM(prhs[1]); uva=mxGetPr(prhs[1]);
u=mxGetPr(prhs[2]); memcpy(uva,u,N*sizeof(double)); 
v=mxGetPr(prhs[3]); memcpy(uva+N,v,N*sizeof(double)); 
a=mxGetPr(prhs[4]); memcpy(uva+2*N,a,N*sizeof(double)); 

}
/*----------------------------------------------------------- newmarkinterp */
else if (!strcmp("newmarkinterp",buf))  {
  /* of_time ('newmarkinterp', out, beta,gamma,uva,a1, t0,t1 )
     %Geradin p.371 equation 7.3.9   */
double     *beta, *gamma, *uva, *a1,*t0, *t1, *def, *v, *a, *cur, *tout, dt;
int         j1, j2, M, *ind, Nind, Ntout;

/* inputs */
if (nrhs!=8)  mexErrMsgTxt ("8 inputs required.");
beta = mxGetPr(prhs[2]); gamma = mxGetPr(prhs[3]);
uva = mxGetPr(prhs[4]); a1 = mxGetPr(prhs[5]);
t0 = mxGetPr(prhs[6]);  t1 = mxGetPr(prhs[7]); 

/* out structure (input and output) */
cur = mxGetPr (mxGetField(prhs[1],0,"cur"));    /* [last_ind  last_t] */
tout = mxGetPr (mxGetField(prhs[1],0,"data"));   
 Ntout=mxGetM(mxGetField(prhs[1],0,"data"));
if (!mxIsInt32(mxGetField(prhs[1],0,"OutInd"))) mexErrMsgTxt("ind must be int32");
ind = mxGetData(mxGetField(prhs[1],0,"OutInd"));

def = mxGetPr (mxGetField(prhs[1],0,"def"));  
/*
v   = mxGetPr (mxGetField(prhs[1],0,"v"));  
a   = mxGetPr (mxGetField(prhs[1],0,"a"));  
*/
Nind=mxGetN(mxGetField(prhs[1],0,"OutInd"));   /* length of ind        */
if (mxGetM(mxGetField(prhs[1],0,"OutInd"))>Nind)
  Nind=mxGetM(mxGetField(prhs[1],0,"OutInd"));
M=mxGetM(prhs[4]);   /* Nddl        */

 j1=(int)cur[0]; 
 if (j1>Ntout || tout[j1]>t1[0]) j1=0;

while(j1<=Ntout) {
 if (tout[j1]>=t0[0] && tout[j1]<=t1[0]) {
  cur[0]=(double)(j1+1); cur[1]=tout[j1];   dt=tout[j1]-t0[0];
  /* mexPrintf(" j1=%i [%g : %g : %g ]\n",j1,t0[0],tout[j1],t1[0]); */

  for (j2=0; j2<Nind; j2++) { 
   def[j2+Nind*j1]= uva[ind[j2]-1] + dt*uva[ind[j2]-1+M] 
     + dt*dt*(.5-beta[0])*uva[ind[j2]-1+2*M] + dt*dt*beta[0]*a1[ind[j2]-1];
   /*
   v[j2+Nind*j1]=uva[ind[j2]-1+M]+(1-gamma[0])*dt*uva[ind[j2]-1+2*M]
     +gamma[0]*dt*a1[ind[j2]-1];
   a[j2+Nind*j1]= a1[ind[j2]-1];
   */
  }
 }
 if (tout[j1]>t1[0])  {  break; }
 j1++;
}


}
/*----------------------------------------------------------- end command */
} /* end mexFunction */    
