#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>

#include "mex.h"
#define PI 3.14159265

typedef struct
{
	double xi;
   double yi;
   double zi;
   double ri;
}Boule;



void mexFunction(int nlhs,mxArray *plhs[],int nrhs,const mxArray *prhs[])
{
   if ((nrhs==7 || nrhs==8 || nrhs ==9) && (nlhs==4))
   {

	   double Xmax,Ymax,Zmax,res,deltat,N2,C,mu,sigma2,temp,K;
      int nbpasX,nbpasY,nbpasZ,N;
      double law;
      int i,j,k,l;
      double umax,H1;
      double *u,*W,*tmp,***Qr,*xt,*yt,*zt;
      Boule *b;
      time_t t;

      int **BX,**BY,**BZ;
	   double *qrbis;
      int law2;

      int Arrondi(double x);    /*Arrondi un reel à l'entier le plus proches*/
      int TirageAleatoire(int n); /*Tirage uniforme entre 0 et n-1*/
      double NormalBoxMuller(double mean,double s); /*Tirage loi normale*/



   for (i=0;i<nrhs;i++)
 	   if (mxGetN(prhs[0])*mxGetM(prhs[0])!=1)
		{    mexErrMsgTxt("All inputs must be scalar");
      }

/*Affectation des arguments de la ligne de commandes*/
      Xmax=mxGetScalar(prhs[0]);
      printf("Xmax=%f\n",Xmax);
      Ymax=mxGetScalar(prhs[1]);
      printf("Ymax=%f\n",Ymax);
      Zmax=mxGetScalar(prhs[2]);

      printf("Zmax=%f\n",Zmax);
      res=mxGetScalar(prhs[3]);
      printf("res=%f\n",res);
      deltat=mxGetScalar(prhs[4]);
      printf("deltat=%f\n",deltat);
      law=mxGetScalar(prhs[5]);
      printf("law=%f\n",law);
      printf ("argc=%d\n",nrhs);

      law2=Arrondi(law);
      printf("law2=%d\n",law2);

/* choix aleatoire d'un germe pour la fonction rand avec l'heure calendaire */
      srand((unsigned) time(&t));

      switch(nrhs)
      {
          case 7 : switch(law2)
                        {
                           case 2 : printf("poisson composé log exponentielle\n");
                                         temp=mxGetScalar(prhs[6]);
                                         C=1;
                                         break;
                           case 3 : printf("constante\n");
                                         K=mxGetScalar(prhs[6]);
                                         printf("K=%f\n",K);
                                         C=1;break;
                           default : printf("Erreur\n");exit(0);
                        }
                        break;
          case 8 : switch (law2)
          					{
                           case 1 : printf("poisson composé log normale\n");
                                        mu=mxGetScalar(prhs[6]);
                                        sigma2=mxGetScalar(prhs[7]);
                                        C=1;
                                        break;
                           case 2 : printf("poisson composé log exponentielle\n");
                                        temp=mxGetScalar(prhs[6]);
                                        C=mxGetScalar(prhs[7]);
                                        break;
                           case 3 : printf("constante\n");
                                        K=mxGetScalar(prhs[6]);
                                        C=mxGetScalar(prhs[7]);
                                        break;
                           default : printf("Erreur\n");exit(0);
                        }
                        break;
          case 9 : switch(law2)
                        {
                           case 1 : printf("poisson composé log normale\n");
                                        mu=mxGetScalar(prhs[6]);
                                        sigma2=mxGetScalar(prhs[7]);
                                        C=mxGetScalar(prhs[8]);
                                        break;
                           default : printf("Erreur\n");exit(0);
                        }
                        break;
          default :   printf("usage : nom_executable <Xmax> <Ymax> <Zmax> <resolution> <pas_de_espace> <choix_loi> <param_loi1> (<param_loi2>) (<constante:densifier tirages> sinon constante=1) \n"); exit(0);
      }

/*Calcul du nb de pas en fonction de Xmax,Ymax,Zmax et de deltat*/
      nbpasX=Arrondi(Xmax/deltat);    /* arrondi a l'entier le plus proche */
      nbpasY=Arrondi(Ymax/deltat);
      nbpasZ=Arrondi(Zmax/deltat);

/*Definition umax*/
      umax=1/(res*res*res)-1;

/*Calcul N*/
      N2=((double)C*2*umax*(Xmax+1)*(Ymax+1)*(Zmax+1))/PI;
      printf("N2=%f\n",N2);
      N=Arrondi(N2);
      printf("N=%d\n",N);



   b=(Boule *)malloc(N*sizeof(Boule));
   u=(double *)malloc(N*sizeof(double));
   for (i=0;i<N;i++)
   {
      b[i].xi=-0.5+rand()/(double)RAND_MAX*(Xmax+1);         /* tirage aleatoire des xi dans [-0.5,Xmax+1] */
      b[i].yi=-0.5+rand()/(double)RAND_MAX*(Ymax+1);         /* tirage aleatoire des yi dans [-0.5,Ymax+1] */
      b[i].zi=-0.5+rand()/(double)RAND_MAX*(Zmax+1);         /* tirage aleatoire des zi dans [-0.5,Zmax+1] */
      u[i]=rand()/(double)RAND_MAX*umax;                         /* tirage aleatoire des ri */
      b[i].ri=pow((double)(1/(res*res*res)-u[i]),-(double)1/3);            /* calcul des ri */
   }

   W=(double *)malloc(N*sizeof(double));
   tmp=(double *)malloc(N*sizeof(double));
   switch(law2)
   {
          case 1 :  for (i=0;i<N;i++)
                         {
                             tmp[i]=NormalBoxMuller(mu,sigma2);
                             W[i]=exp(tmp[i]);
                        }
                        H1=C*(1-exp(mu+sigma2/2));
                        break;

          case 2 : for (i=0;i<N;i++)
                        {
                            tmp[i]=rand()/(double)RAND_MAX;
                            W[i]=pow(pow(1+temp,(double)1/temp)*tmp[i],temp);
                        }
                        H1=0;
                        break;

          case 3 : for (i=0;i<N;i++)   W[i]=(double)K;
                        H1=C*(1-K);
                        break;

          default : printf("Erreur lors de l'initialisation de W\n");
   }

   free(tmp);

   Qr=(double ***)malloc(nbpasX*sizeof(double));
   for (i=0;i<nbpasX;i++)
   {
      Qr[i]=(double **)malloc(nbpasY*sizeof(double));
      for (j=0;j<nbpasY;j++)
         Qr[i][j]=(double *)malloc(nbpasZ*sizeof(double));
   }


   plhs[0]=mxCreateDoubleMatrix(nbpasX,1,mxREAL);
   plhs[1]=mxCreateDoubleMatrix(nbpasY,1,mxREAL);
   plhs[2]=mxCreateDoubleMatrix(nbpasZ,1,mxREAL);

   xt=mxGetPr(plhs[0]);
   yt=mxGetPr(plhs[1]);
   zt=mxGetPr(plhs[2]);

   for (i=0;i<nbpasX;i++)
     	xt[i]=i*deltat;
   for (i=0;i<nbpasY;i++)
       yt[i]=i*deltat;
   for (i=0;i<nbpasZ;i++)
       zt[i]=i*deltat;


   for(i=0;i<nbpasX;i++){
     for (j=0;j<nbpasY;j++){
       for (k=0;k<nbpasZ;k++){
            Qr[i][j][k]=1;
       }
     }
   }

  BX=(int **)malloc(N*sizeof(int));
  BY=(int **)malloc(N*sizeof(int));
  BZ=(int **)malloc(N*sizeof(int));
  for(i=0;i<N;i++)
  {
     BX[i]=(int *)malloc(nbpasX*sizeof(int));
     BY[i]=(int *)malloc(nbpasY*sizeof(int));
     BZ[i]=(int *)malloc(nbpasZ*sizeof(int));
  }


   for (l=0;l<N;l++)
   {
     for(i=0;i<nbpasX;i++)
  	 if (fabs(xt[i]-b[l].xi)<=b[l].ri/2)
  	 {
  	BX[l][i]=1;
  	 }
     for(j=0;j<nbpasY;j++)
            if (fabs(yt[j]-b[l].yi)<=b[l].ri/2)
  	 {
  	  BY[l][j]=1;
  	 }
     for(k=0;k<nbpasZ;k++)
            if (fabs(zt[k]-b[l].zi)<=b[l].ri/2)
  	 {
  	  BZ[l][k]=1;
  	 }

   }

     /* calcul des Qr  */
     for (l=0;l<N;l++)
      if (b[l].xi<=(Xmax+b[l].ri/2)&& b[l].yi<=(Ymax+b[l].ri/2)&& b[l].zi<=(Zmax+b[l].ri/2))
       {
        for(i=0;i<nbpasX;i++)
         if (BX[l][i]!=0)
         {
          for (j=0;j<nbpasY;j++)
           if (BY[l][j]!=0)
           {
             for (k=0;k<nbpasZ;k++)
  	           if (BZ[l][k]!=0)
                {
  	               if ((xt[i]-b[l].xi)*(xt[i]-b[l].xi)+ (yt[j]-b[l].yi)*(yt[j]-b[l].yi) + (zt[k]-b[l].zi)*(zt[k]-b[l].zi)<=(b[l].ri/2)*(b[l].ri/2))
                    {
                       Qr[i][j][k]*=W[l];
                    }
                }
           }
  	    }
       }


   /*    for(i=0;i<nbpasX;i++)
          for (j=0;j<nbpasY;j++)
             for (k=0;k<nbpasZ;k++)
                if(Qr[i][j][k]==1)
                {
  	               Qr[i][j][k]=0;
                }
   */
    /* normalisation de Qr   */
    plhs[3]=mxCreateDoubleMatrix(nbpasX*nbpasY*nbpasZ,1,mxREAL);
    qrbis=mxGetPr(plhs[3]);
    l=0;
    for(i=0;i<nbpasX;i++)
     for (j=0;j<nbpasY;j++)
       for (k=0;k<nbpasZ;k++)
       {
           Qr[i][j][k]=Qr[i][j][k]*pow(res,-H1);
           qrbis[l]=Qr[i][j][k];
           l++;
       }

   free(u);

   free(W);

  }
  else
  {
            printf("usage : nom_executable <Xmax> <Ymax> <Zmax> <resolution> <pas_de_espace> <choix_loi> <param_loi1> (<param_loi2>) (<constante:densifier tirages> sinon constante=1) \n"); exit(0);
  }

}
/**********************Fin du main*************************************************/


int Arrondi(double x)
{
	return (((x-(int)x)<=0.5) ? (int)x : (int)x+1);
}


/***********************************
	Tirage aléatoire d'un entier entre 0 et n-1
***********************************/
int TirageAleatoire(int n)
{
	return ((n==0) ? 0 : rand()%n);
}

/***************************************************************************************
 	Simulation d'une loi normale par la méthode de Box Muller : Si
    x1 et x2 suivent une loi uniforme dans]0,1[, y1=sqrt(-2ln(x1))cos(2M_PI.x2)
    et y2=sqrt(-2ln(x1))sin(2M_PI.x2) suivent une loi normale réduite. On en déduit que
    m+sy1 et m+sy2suivent une loi N(m,s)
****************************************************************************************/
double NormalBoxMuller(double mean,double s)
{
	 double 	x1;
    long    pas;

    /* tirage d'un réel uniforme dans [0,1]  */
    pas = 100;
    x1 = (double)TirageAleatoire(pas)/(double)pas;
    x1 = (x1 <1E-6 ? 1E-6 : x1);
    x1 = (x1 ==1 ? 1-1E-6 : x1);

    return (mean + s*sqrt(-2*log(x1))*cos(2*PI*x1));
}



