
#include "tirages.h"


/******************************************************************************/
/* VALEUR TIREE SELON UNE LOI NORMALE CENTREE REDUITE                         */
/******************************************************************************/

double randn(void)
{
	//	variables locales
	double x1,x2,w;
	static unsigned char valeur_disponible = 0;
	static double autre_valeur;

	if (!valeur_disponible)
	{
		// tirage de deux valeurs suivant une loi normale centree reduite et retour d'une d'elle
		do
		{
			x1 = 2.0*randu()-1.0;
			x2 = 2.0*randu()-1.0;
			w = CARRE(x1)+CARRE(x2);
		} while (w>=1.0);

		w = sqrt((-2.0*log(w))/w);
		autre_valeur = x1*w;
		valeur_disponible = 1;
		return x2*w;
	}
	else
	{
		//	retour d'une valeur precedemment calculee
		valeur_disponible = 0;
		return autre_valeur;
	}
}


/******************************************************************************/
/* ZONES D'INFLUENCE                                                          */
/*                                                                            */
/* Determination des points de la matrice Qr contenus dans un carre de centre */
/* (xi,yi,zi) et de cote ri (l'angle du carre par rapport a Qr etant nul).    */
/******************************************************************************/

long zone_influence(double xi, double yi, double zi, double ri, long nx, long ny, long nz,
                    double dx, positions_t *positions)
{

	//	variables locales
	long dix,fix,diy,fiy,diz,fiz,retour=0;

	//	limites de la zone d'influence
	dix = MAX(0,(long)ceil ((xi-ri/2.0)/dx));
	fix = MIN(nx-1,(long)floor((xi+ri/2.0)/dx));

	//	si le point a une zone d'inflence selon x, on continue
	if (dix<=fix)
	{
		diy = MAX(0,(long)ceil ((yi-ri/2.0)/dx));
		fiy = MIN(ny-1,(long)floor((yi+ri/2.0)/dx));
		if (diy<=fiy)
		{
			diz = MAX(0,(long)ceil ((zi-ri/2.0)/dx));
			fiz = MIN(nz-1,(long)floor((zi+ri/2.0)/dx));

			//	si le point a une zone d'inflence selon y, on continue
			if (diz<=fiz)
			{
				positions->tableau[positions->n].x = xi;
				positions->tableau[positions->n].y = yi;
				positions->tableau[positions->n].z = zi;
				positions->tableau[positions->n].r = ri;
				positions->tableau[positions->n].dix = dix;
				positions->tableau[positions->n].fix = fix;
				positions->tableau[positions->n].diy = diy;
				positions->tableau[positions->n].fiy = fiy;
				positions->tableau[positions->n].diz = diz;
				positions->tableau[positions->n].fiz = fiz;
				positions->n++;
				retour=1;
			}
		}
	}
	return retour;
}


/******************************************************************************/
/* TIRAGE DES POSITIONS                                                       */
/*                                                                            */
/* Generation d'un point (xi,yi,zi,ri) et determination de sa zone d'influence.  */
/******************************************************************************/

void tirer_positions(long N, double x_min, double x_max, double y_min, double y_max, double z_min, double z_max,
                     double u_max, double r_min, double bord, double dx,
                     positions_t *positions, long nx, long ny, long nz)
{

	// variables locales
	long i;
	double xi,yi,zi,ui,ri,delta_x=x_max-x_min,delta_y=y_max-y_min,delta_z=z_max-z_min;

	//mexPrintf("N = %d\n", N);
	// tirage des positions
	for (i=0; i<N; i++)
	{
		xi = x_min-0.5*bord+(delta_x+bord)*randu();
		yi = y_min-0.5*bord+(delta_y+bord)*randu();
		zi = z_min-0.5*bord+(delta_z+bord)*randu();
		ui = u_max*randu();
		ri = pow(1.0/CUBE(r_min)-ui,-(double)1.0/3);

		zone_influence(xi, yi, zi, ri, nx, ny, nz, dx, positions);
	}
}


/******************************************************************************/
/* TIRAGE DES MULTIPLIEURS                                                    */
/*                                                                            */
/* Generation d'un multiplieur pour chaque point (xi,yi,ri) en fonction de la */
/* loi, ainsi que d'un angle de rotation pour le motif. On construit la liste */
/* des variables aleatoires retournees sous Matlab.                           */
/******************************************************************************/

void tirer_multiplieurs(positions_t positions, tirage_t *tirages, loi_t loi,
                        motif_t type_motif, angle_t angle, numero_fonction_t numero_fonction,
                        long numero_motif,long numero_zone)
{

	//	variables locales
	long i;
	double cste;
	parametres_t p=loi.parametres;

	// recopie des nouvelles positions, tirage des angles de rotation et affectation des caracteristiques de la zone
	for (i=0; i<positions.n; i++)
	{
		memcpy(&(tirages[i]),&(positions.tableau[i]),4*sizeof(double));
		tirages[i].theta = 2.*T_PI*randu();
		tirages[i].type_motif = type_motif;
		tirages[i].angle = angle;
		tirages[i].numero_fonction = numero_fonction;
		tirages[i].numero_motif   =numero_motif;
		tirages[i].numero_zone    =numero_zone;
	}

	// tirage des multiplieurs
	switch (loi.type)
	{
		case log_normal:
			cste=sqrt(p.ln.sigma2);
			for (i=0; i<positions.n; i++)
				tirages[i].w=exp(cste*randn()+p.ln.mu);
				break;
		case log_exponential:
			cste = pow(1.+p.le.Temp,1./p.le.Temp);
			for (i=0; i<positions.n; i++)
				tirages[i].w = pow(cste*randu(),p.le.Temp);
				break;
		case log_poisson:
			for (i=0; i<positions.n; i++)
				tirages[i].w = p.lp.W0;
			break;
	}
}

