/*=================================================================
 * Modele Leguerinel en C et double precision
 * Pour la
 * Simulation de strategie optimale, version C
 * Voir BaranyiRoberts pour les E/S
 * COMPILATION sous matlab par: 
 * mex  /home/abouvier/MONFILTREX/v2013-1.0/SRC/SIMULATION_STRATEGIE/C/CLeguerinel.c

 *============================================================*/
#include <math.h>
#include "mex.h"
double pow( double base, double exponent );
double log10( double arg );

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

void  CLeguerinelDP(double *param, double * temps, 
        double *valmin, double *valmax, double *rep,
        int ntemps, int N, int Ncur)
{

double delta=(double) param[Ncur];
double p=(double) param[N+Ncur];
double K= (double) param[2*N+Ncur];
double pHopt= (double) param[3*N+Ncur];
double Awopt= (double) param[4*N+Ncur];
double ZpH = (double) param[5*N+Ncur];
double ZAw = (double) param[6*N+Ncur];
double pH  = (double) param[7*N+Ncur];
double Aw  = (double) param[8*N+Ncur];
 
 int i, ind, u;
 double aux1, aux2, aux3,H, ti;

 for ( i=0; i< ntemps; i++) {
   ti=temps[i];
   aux1 = (pH -pHopt) / ZpH;
   aux1 = aux1 * aux1;
   aux2 = (Aw - Awopt) / ZAw;
   aux2 = aux2 * aux2;
   H = pow(10.0, (aux1+aux2));
   aux3 = ti / (delta * H);

   ind=N*i+Ncur;
   rep[ind]= K -  pow(aux3, p);

 /* Ramener dans les bornes les réponses invalides */
   /*ATT: ce qu'on a calculé, c'est log10(Nt) */
   /*  Nt doit etre >0 => log10(Nt) >1 */

    if (rep[ind] < (*valmin)) {
      rep[ind]=(double)(*valmin);
    }
    if (rep[ind] > (*valmax)) {
      rep[ind]=(double)(*valmax);
    }
     
 } /* fin i */
}
/* ============================================================== */
void CLeguerinel(double *param,  double *temps, 
 double *valmin, double *valmax,
  double *rep,int ntemps,
  int N)
{

    /* boucle sur l'appel à LeguerinelDP */
    int u;
    for (u=0; u<N; u++) {
        CLeguerinelDP( 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 CLeguerinel)
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;

  

/* The first input argument param is Nxnf */
  N = (int) mxGetM(prhs[0]);
   
/* 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 */
 CLeguerinel( param,  temps, valmin, valmax, rep,(int)ntemps, (int)N);

} 
