/*=================================================================
 * Modele BaranyiRoberts en C et double precision (VERSION log neperien)
 * Pour la
 * Simulation de strategie optimale, version C
 * C'est la version log du modèle car elle donne de meilleurs résultats que la version exp utilisée dans les autres tâches
 * Voir BaranyiRoberts pour les E/S
 * COMPILATION sous matlab par: 
 * mex  /home/abouvier/MONFILTREX/v2013-1.0/SRC/SIMULATION_STRATEGIE/C/CLnBaranyiRoberts.c

 *============================================================*/
#include <math.h>
#include "mex.h"
double exp( double arg );
double log( double arg );

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

void  CLnBaranyiRobertsDP(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 lnN0= (double) param[2*N+Ncur];
    double lnNmax=(double) param[3*N+Ncur];
    double At, A0t, Bt, t, rep0, kk, e4;
    int ind, i;
    
    for ( i=0; i< ntemps; i++) {
      t=temps[i];
      A0t = (exp(-mumax*t))+ (exp(-mumax * lambda)) -
	           exp(-mumax*t - mumax * lambda);
    if (A0t <= 0) {
      A0t=10.e-40;
    }
    At= t+log(A0t)/mumax;
    e4= exp( mumax * At);
    Bt= 1 + (e4 -1)/exp(lnNmax-lnN0);
    rep0=lnN0+  (mumax*At) - log(Bt);
    ind=N*i+Ncur;
     rep[ind]=rep0;
 /* Ramener dans les bornes les réponses invalides */
/*     mexPrintf("valmin %g valmax %g\n", *valmin, *valmax); */
/* VOIR
 * if (rep[ind] < (*valmin)) {
      rep[ind]=(double)(*valmin);
    }
     
    if (rep[ind] > (*valmax)) {
      rep[ind]=(double)(*valmax);
    }
 
*FIN VOIR */     

    } /* fin i */
    
    
}
/* ============================================================== */
void CLnBaranyiRoberts(double *param, double *temps,
 double *valmin, double *valmax,
 double *rep, int ntemps, int N)
{

    /* boucle sur l'appel à LnBaranyiRobertsDP */
    int u;
    for (u=0; u<N; u++) {
      CLnBaranyiRobertsDP( param,  temps,valmin, valmax, rep, ntemps, N, u);
    }

}

/* ================================================== */
/* The gateway function:
INPUT
nlhs: le nombre de sorties (la fonction a pls sorties quand
 l'utilisateur affecte une structure à pls composants
 à son appel a CLnBaranyiRoberts)
plhs: les sorties
nrhs: le nombre d'arguments d'entrée
prhs: les arguments d'entrée
 */
void mexFunction( int nlhs, mxArray *plhs[],
                  int nrhs, const mxArray *prhs[] )
{
  double *rep, *param, *temps, *valmin, *valmax;
  int N, ntemps;

  
  /* VERIF AU CAS OU, a oter  pour gagner du temps */
  /* Check for proper number of arguments. */
  /* 5 input arguments: les parametres (N,4), et les temps (1,ntemps):

  if(nrhs!=5) {
    mexErrMsgIdAndTxt( "LnBaranyiRoberts:invalidNumInputs",
            "5 input required.");
  } 
 */  

/* The first input argument param is Nx4 */
  N = (int) mxGetM(prhs[0]);
 
  /* VERIF AU CAS OU, a oter  pour gagner du temps */
 /*  param must have 4 columns */
  /*
if(mxGetN(prhs[0])!=4) {
    mexErrMsgIdAndTxt("LnBaranyiRoberts:notRowVector",
                      "Input param must have 4 columns.");
}
  */

/* create a pointer to the  input   */
param = mxGetPr(prhs[0]);

/* Second input argument: times is 1 x ntemps */
ntemps = (int) mxGetN(prhs[1]);

temps=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]);

/* create the output matrix
 plhs[0] = mxCreateDoubleMatrix((mwSize)N,(mwSize)ntemps,mxREAL);
rep = mxGetPr(plhs[0]);
 */
/* call the computational routine */
 CLnBaranyiRoberts( param, temps,valmin, valmax,  rep,(int)ntemps, (int)N);

} 
