/*=================================================================
 * FUNCTION
 * Compute Baranyi- Roberts model,
 * INPUT
 * param: matrix (N x nf). Values of the 'nf' model parameters.
 *   (N=LHS, with the Sobol method N=1, with the Tornsey method)
 * temps : matrix (1, ntimes). Times
 * valmin, valmax: scalars. Lower and upper bounds of admissible model values.
 * INPUT/OUTPUT
 * rep: matrix (N x ntimes) allocated before the call. 
 * Model values. 
 * COMPILATION in Matlab session by:
 * mex  <FILTREX pathname>/SRC/SIMULATION_STRATEGIE/C/CBaranyiRoberts.c
 * (or by file 'compil.m')
 *  CALLED BY: calcStrategieOptSobolDP (Sobol method), deriv1p (Tornsey method)
 *  in the task: "simulation of an optimal sampling". *
 *============================================================*/
#include <math.h>
#include "mex.h"
double exp( double arg );


/* ============================================================== */

void  CBaranyiRobertsDP(double *param, double * temps,
        double *valmin, double *valmax, double *rep,
        int ntemps, int N, int Ncur)
{
    
    double mumax=(double) param[Ncur];
    double lambda=(double) param[N+Ncur];
    double N0= (double) param[2*N+Ncur];
    double Nmax=(double) param[3*N+Ncur];
    double e1, B, expmumaxA, ti;
    int i, ind;
    
     
    for ( i=0; i< ntemps; i++) {
        ti=(double)temps[i];
        e1= mumax *  (ti-lambda);
        if (e1>100) {
            e1=100;
        }
        expmumaxA=1.0 + exp(e1)-exp(-mumax * lambda);
        B=1+(expmumaxA-1) * (N0/Nmax);
        ind=N*i+Ncur;
        rep[ind]=(double)(N0 * (expmumaxA/B));

        /* Put inside the bounds the invalid model output */
    if (rep[ind] < (*valmin)) {
      rep[ind]=(double)(*valmin);
    }
    if (rep[ind] > (*valmax)) {
      rep[ind]=(double)(*valmax);
    }
   
        
    } /* fin i */
    
    
}
/* ============================================================== */
void CBaranyiRoberts(double *param, double *temps,
        double *valmin, double *valmax,
        double *rep, int ntemps, int N)
{
    
    
    /* Loop on the call to  BaranyiRobertsDP */
    int u;
    for (u=0; u<N; u++) {
                
        CBaranyiRobertsDP( param,  temps, valmin, valmax, rep, ntemps, N, u);
    }
        
}

/* ================================================== */
/* The gateway function:
 * INPUT
 * nlhs: number of components (la fonction a pls sorties quand
 * l'utilisateur affecte une structure à pls composants
 * à son appel a CBaranyiRoberts)
 * plhs: output
 * nrhs: number of input parameters
 * prhs: nput parameters
 */
void mexFunction( int nlhs, mxArray *plhs[],
        int nrhs, const mxArray *prhs[] )
{
    double  *param, *temps, *rep, *valmin, *valmax;
    int N, ntemps;
    
     
    /* The first input argument param is Nxnf */
    N = (int) mxGetM(prhs[0]);
    
    /* Check : param must have 4 columns */
    /* Check removed in the distribution release
if(mxGetN(prhs[0])!=4) {
    mexErrMsgIdAndTxt("BaranyiRobertsDP:notRowVector",
                      "Input param must have 4 columns.");
}
     */
    
    /* create a pointer to the  input   */   
    param = (double *)mxGetPr(prhs[0]);
   
    /* Second input argument: times is 1 x ntimes (or network npoints when  D-optima method)*/
    ntemps = (int) mxGetN(prhs[1]);
    temps=(double *)mxGetPr(prhs[1]);
     
    /* Third input argument:valmin */
    valmin=(double *)mxGetPr(prhs[2]);
    
    /* Fourth input argument:valmax */
    valmax=(double *)mxGetPr(prhs[3]);
        
    /* 5th input argument: result matrix */
    rep = (double *)mxGetPr(prhs[4]);
    /* Check removed in the distribution release
  
 if(mxGetM(prhs[4]) != N) {
  
   mexErrMsgIdAndTxt("BaranyiRobertsDP: Y should have N rows", "Input error");
 }
 if(mxGetN(prhs[4]) != ntemps) {
   mexErrMsgIdAndTxt("BaranyiRobertsDP: Y should have ntimes columns", "Input error");
 }
     */
 
    /* call the computational routine */
    CBaranyiRoberts( param, temps, valmin, valmax, rep,(int)ntemps, (int)N);
        
}
