#include "tirages.h"
#include "calculQr.h"

void mexFunction(int nb_sorties, mxArray *sorties[], int nb_entrees,const mxArray *entrees[]) {

	// variables locales
	long nx, ny, nz, L, l, N, indice_debut_non_shaper, i, j, k, tmp;
	double x_max, y_max, z_max, r_min, u_max, dx, lambda, H1, Ttheo2, *Qr_out, ***Qr_tmp, quotient;
	unsigned long  graineinit;
	double *grainefin;
	tirage_t *tirages_zoom;
	positions_t positions_zoom;
	motif_t motif;
	angle_t angle;
	numero_fonction_t numero_fonction;
	loi_t loi;
	positions_t positions;
	tirage_t *tirages;

		//	//	//	Recuperation des donnees d'entree

	// verification des entrees et sorties
	if (nb_entrees!=10)
		mexErrMsgTxt("Il faut dix arguments en entree !");
	if (nb_sorties!=2)
		mexErrMsgTxt("Il faut deux arguments en sortie !");

	// taille de la texture resultante
	if (mxGetM(entrees[0])!=1 || mxGetN(entrees[0])!=1 || (x_max=*mxGetPr(entrees[0]))<=0.0)
		mexErrMsgTxt("x_max doit etre un scalaire strictement positif !");
	if (mxGetM(entrees[1])!=1 || mxGetN(entrees[1])!=1 || (y_max=*mxGetPr(entrees[1]))<=0.0)
		mexErrMsgTxt("y_max doit etre un scalaire strictement positif !");
	if (mxGetM(entrees[2])!=1 || mxGetN(entrees[2])!=1 || (z_max=*mxGetPr(entrees[2]))<=0.0)
		mexErrMsgTxt("z_max doit etre un scalaire strictement positif !");
	if (mxGetM(entrees[3])!=1 || mxGetN(entrees[3])!=1 || (r_min=*mxGetPr(entrees[3]))<=0.0)
		mexErrMsgTxt("r_min doit etre un scalaire strictement positif !");
	if (r_min>=1.0)
		mexErrMsgTxt("r_min doit etre strictement inferieur a 1 !");
	if (mxGetM(entrees[4])!=1 || mxGetN(entrees[4])!=1 || (dx=*mxGetPr(entrees[4]))<=0.0)
		mexErrMsgTxt("dx doit etre un scalaire strictement positif !");


	//	loi de la texture
	if (mxGetM(entrees[5])!=1 || mxGetN(entrees[5])!=1)
		mexErrMsgTxt("type_loi doit etre un scalaire !");
	loi.type = (type_loi_t)*mxGetPr(entrees[5]);
	if (loi.type<1 || loi.type>3)
		mexErrMsgTxt("Loi inconnue !");
	L = MAX(mxGetM(entrees[6]),mxGetN(entrees[6]));
	l = MIN(mxGetM(entrees[6]),mxGetN(entrees[6]));
	if (l==1 && L==sizeof(parametres_t)/sizeof(double)+1)
		memcpy(&(loi.parametres),mxGetPr(entrees[6]),L*sizeof(double));
	else
		mexErrMsgTxt("Mauvais nombre de parametres de loi !");
	if (loi.C<=0.)
		mexErrMsgTxt("La parametre de densite de la loi doit etre strictement positif !");

	switch (loi.type)
	{
		case log_normal:
			if (loi.parametres.ln.sigma2<0.0)
				mexErrMsgTxt("Le parametre sigma2 doit etre positif !");
			H1 = loi.C*(1.0-exp(loi.parametres.ln.mu+loi.parametres.ln.sigma2/2.0));
			Ttheo2 = -loi.C*(exp(2.0*(loi.parametres.ln.sigma2+loi.parametres.ln.mu))-1.0);
			break;
		case log_exponential:
			if (loi.parametres.le.Temp<=0.0)
				mexErrMsgTxt("Le parametre Temp doit etre strictement positif !");
			H1 = 0.0;
			Ttheo2 = -loi.C*(CARRE(1.0+loi.parametres.le.Temp)/(1.0+2.0*loi.parametres.le.Temp)-1.0);
			break;
		case log_poisson:
			if (loi.parametres.lp.W0<0.0 || loi.parametres.lp.W0>=1.0)
				mexErrMsgTxt("Le parametre W0 doit etre positif et strictement inferieur a 1 !");
			H1 = loi.C*(1.0-loi.parametres.lp.W0);
			Ttheo2 = loi.C*(-2.0*(1.0-loi.parametres.lp.W0)+1.0-CARRE(loi.parametres.lp.W0));
			break;
	}
	if (Ttheo2<-1.)
		mexWarnMsgTxt("La convergence n'est pas guarantie !");

	//	motif de la texture
	if (mxGetM(entrees[7])!=1 || mxGetN(entrees[7])!=1)
		mexErrMsgTxt("motif doit etre un scalaire !");
	motif = (motif_t)*mxGetPr(entrees[7]);
	if (motif<1 || motif>4)
		mexErrMsgTxt("type de motif inconnu !");
	if (mxGetM(entrees[8])!=1 || mxGetN(entrees[8])!=1)
		mexErrMsgTxt("angle doit etre un scalaire !");
	angle = (angle_t)*mxGetPr(entrees[8]);
	if (angle!=fixe && angle!=variable)
		mexErrMsgTxt("angle doit etre egal a 0 ou 1 !");
	
	//graineinit=*(unsigned long *)mxGetData(entrees[9]);
    graineinit=*(double *)mxGetData(entrees[9]);
	//mexPrintf("graineinit = %u \n",graineinit);

	switch (motif)
	{
		case m_carre:
			if (angle==fixe)
				numero_fonction = n_carre_fixe;
			else
				numero_fonction = n_carre_variable;
			break;
		case m_cercle:
			numero_fonction = n_cercle;
			break;
		case m_shaper_cos_carre:
			numero_fonction = n_shaper_cos_carre;
			break;
		case m_shaper_exp:
			numero_fonction = n_shaper_exp;
			break;
	}

		//	//	//	Calcul de la cascade

	// tirage des positions des multiplieurs
	u_max = 1.0/CUBE(r_min)-1.0;
	nx = (long)(x_max/dx)+1;
	ny = (long)(y_max/dx)+1;
	nz = (long)(z_max/dx)+1;
	lambda = loi.C*2.0/T_PI*u_max*(x_max+1.0)*(y_max+1.0)*(z_max+1.0);
	N = (long)(lambda+0.5);
	positions.tableau = (position_t *)malloc(N*sizeof(position_t));
	positions.n = 0;
	//mexPrintf("Allocation OK\n");

	if (graineinit == 0)
		rand_init(time(NULL));
	else
		rand_init(graineinit);

	tirer_positions(N, 0.0, x_max, 0.0, y_max, 0.0, z_max, u_max, r_min, 1.0, dx, &positions, nx, ny, nz);
	if (est_shaper(motif))
		indice_debut_non_shaper = positions.n;
	else
		indice_debut_non_shaper = 0;
	//mexPrintf("Tirage positions OK\n");

	// allocation memoire pour les sorties
	Qr_tmp = (double ***)malloc(nx*sizeof(double**));
	for (i=0; i<nx; i++)
	{
        Qr_tmp[i]=(double **)malloc(ny*sizeof(double*));
		for(j=0; j<ny; j++)
		{
			Qr_tmp[i][j]=(double *)malloc(nz*sizeof(double));
		}
	}
	allocation_sorties(sorties, &Qr_out, nx, ny, nz, grainefin);
	sorties[1]=mxCreateDoubleMatrix(1,1,mxREAL); 
	grainefin=mxGetPr(sorties[1]);
	//mexPrintf("Allocation sorties OK\n");
	
	tirages = (tirage_t*)malloc(positions.n*sizeof(tirage_t));

	// tirage des multiplieurs (des wi ainsi que des thetai)
	tirer_multiplieurs(positions, tirages, loi, motif, angle, numero_fonction,0,0);
	//mexPrintf("Tirage Wi OK\n");

	// calcul de Qr global (de 1 a r_min)
	calcul_Qr(dx, nx, ny, nz, Qr_tmp, indice_debut_non_shaper, positions, tirages);
	//mexPrintf("Calcul Qr OK\n");

	// normalisation de Qr
	quotient = pow(r_min,H1);

	tmp = 0.0;
	for (i=0; i<nx; i++)
    {
        for (j=0; j<ny; j++)
        {
			for (k=0; k<nz; k++)
			{	
				Qr_out[tmp] = Qr_tmp[i][j][k]/quotient;
				tmp++;
			}
        }
    }
	
	*grainefin = (double) get_mt();
	//printf("graine sortie = %f\n", grainefin);
	
	
    // liberation de la memoire temporaire
	free(positions.tableau);
	free(tirages);
	free(grainefin);
	for (i=0; i<nx; i++)
	{
        for(j=0; j<ny; j++)
		{
			free(Qr_tmp[i][j]);
		}
		free(Qr_tmp[i]);
	}
	free(Qr_tmp);
	
	free(Qr_out);
}

