#include "calculQr.h"

/******************************************************************************/
/* VARIABLES GLOBALES                                                         */
/******************************************************************************/

double    dx;
long      ny;
double    ***Qr_tmp;

/******************************************************************************/
/* FONCTION DE MOTIF (CARRE, ANGLE FIXE)                                      */
/******************************************************************************/

void f_carre_fixe(position_t *position, tirage_t *tirage)
{
	//	variables locales
	double xi,yi,zi,ri,wi;
	long i,j,k;

	//mexPrintf("carre \n");
	xi = tirage->x;
	yi = tirage->y;
	zi = tirage->z;
	ri = tirage->r;
	wi = tirage->w;
	for (i=position->dix; i<=position->fix; i++)
	{
		for (j=position->diy; j<=position->fiy; j++)
		{
			for (k=position->diz; k<=position->fiz; k++)
			{
				Qr_tmp[i][j][k] *= wi;
			}
		}
	}
}


/******************************************************************************/
/* FONCTION DE MOTIF (CARRE, ANGLE VARIABLE)                                  */
/******************************************************************************/

void f_carre_variable(position_t *position, tirage_t *tirage)
{
	//	variables locales
	double xi,yi,zi,ri_2rc2,wi,cos_thetai,sin_thetai;
	long i,j,k;

	xi = tirage->x;
	yi = tirage->y;
	zi = tirage->z;
	ri_2rc2 = tirage->r*C_1_2RC2;
	wi = tirage->w;
	cos_thetai = cos(tirage->theta);
	sin_thetai = sin(tirage->theta);
	for (i=position->dix; i<=position->fix; i++)
	{
		for (j=position->diy; j<=position->fiy; j++)
		{
			for (k=position->diz; k<=position->fiz; k++)
			{
				if (fabs(cos_thetai*(i*dx-xi)+sin_thetai*(j*dx-yi))<=ri_2rc2 &&
					fabs(sin_thetai*(i*dx-xi)-cos_thetai*(j*dx-yi))<=ri_2rc2)
					Qr_tmp[i][j][k] *= wi;
			}
		}
	}
}


/******************************************************************************/
/* FONCTION DE MOTIF (CERCLE)                                                 */
/******************************************************************************/

void f_cercle(position_t *position, tirage_t *tirage)
{
	//	variables locales
	double xi,yi,zi,ri2_4,wi,r2;
	long i,j,k;
	
	//mexPrintf("cercle \n");
	
	xi = tirage->x;
	yi = tirage->y;
	zi = tirage->z;
	ri2_4 = 0.25*CARRE(tirage->r);
	wi = tirage->w;
	for (i=position->dix; i<=position->fix; i++)
	{
		for (j=position->diy; j<=position->fiy; j++)
		{
			for (k=position->diz; k<=position->fiz; k++)
			{
				r2 = CARRE(i*dx-xi)+CARRE(j*dx-yi)+CARRE(k*dx-zi);
				if (r2<ri2_4)
				{
					Qr_tmp[i][j][k] *= wi;
				}
			}
		}
	}
}


/******************************************************************************/
/* FONCTION DE MOTIF (SHAPER cos^2(pi.r) )                                    */
/******************************************************************************/

void f_shaper_cos_carre(position_t *position, tirage_t *tirage)
{
	//	variables locales
	double xi,yi,zi,ri,ri2_4,wi,r2;
	long i,j,k;

	xi = tirage->x;
	yi = tirage->y;
	zi = tirage->z;
	ri = tirage->r;
	ri2_4 = 0.25*CARRE(ri);
	wi = tirage->w;
	for (i=position->dix; i<=position->fix; i++)
	{
		for (j=position->diy; j<=position->fiy; j++)
		{
			for (k=position->diz; k<=position->fiz; k++)
			{
				r2 = CARRE(i*dx-xi)+CARRE(j*dx-yi)+CARRE(k*dx-zi);
				if (r2<ri2_4)
					Qr_tmp[i][j][k] += log((wi?wi:EPSILON))*CARRE(cos(T_PI*sqrt(r2)/ri));
			}
		}
	}
}


/******************************************************************************/
/* FONCTION DE MOTIF (SHAPER exp[ 1 - 1/(1-(2r)^4) ])                         */
/******************************************************************************/

void f_shaper_exp(position_t *position,tirage_t *tirage)
{
	//	variables locales
	double xi,yi,zi,ri,ri2_4,wi,r2;
	long i,j,k;

	xi = tirage->x;
	yi = tirage->y;
	zi = tirage->z;
	ri = tirage->r;
	ri2_4 = 0.25*CARRE(ri);
	wi = tirage->w;
	for (i=position->dix; i<=position->fix; i++)
	{
		for (j=position->diy; j<=position->fiy; j++)
		{
			for (k=position->diz; k<=position->fiz; k++)
			{
				r2 = CARRE(i*dx-xi)+CARRE(j*dx-yi)+CARRE(k*dx-zi);
				if (r2<ri2_4)
					Qr_tmp[i][j][k] += log((wi?wi:EPSILON))*exp(1.-1./(1.+EPSILON-CARRE(r2/ri2_4)));
			}
		}
	}
}


/******************************************************************************/
/* CALCUL DE QR                                                               */
/*                                                                            */
/* Construction de la matrice Qr en effectuant pour chacun de ses points le   */
/* produit des multiplieurs des points generes aleatoirement dont le motif    */
/* 'eclaire' le point. On commence par parcourir les motifs de type fonction  */
/* lissante.                                                                  */
/******************************************************************************/

void calcul_Qr(double _dx, long nx, long _ny, long nz, double ***_Qr_tmp, long debut_non_shaper,
               positions_t positions, tirage_t *tirages)
{
	//	variables locales
	long i,j,k;
	void (*f_motif[9])(position_t *, tirage_t *) = {
		f_carre_fixe,f_carre_variable,f_cercle,f_shaper_cos_carre,f_shaper_exp
	};

	//	initialisation des variables globales
	dx = _dx;
	ny = _ny;
	Qr_tmp = _Qr_tmp;
	//mexPrintf("initialisation locale OK\n");

	//	fonctions de motif de type fonction lissante
	if (debut_non_shaper>0)
	{
		for (i=0; i<nx; i++)
		{
			for (j=0;j<ny; j++)
			{
				for (k=0; k<nz; k++)
					Qr_tmp[i][j][k] = 0.0;
			}
		}
		//mexPrintf("initialisation OK\n");
		for (i=0; i<debut_non_shaper; i++)
			f_motif[tirages[i].numero_fonction](positions.tableau+i, tirages+i);
		for (i=0; i<nx; i++)
		{
			for (j=0;j<ny; j++)
			{
				for (k=0; k<nz; k++)
					Qr_tmp[i][j][k] = exp(Qr_tmp[i][j][k]);
			}
		}
		//mexPrintf("calcul OK\n");
	}
	else
	{	
		for (i=0; i<nx; i++)
		{
			for (j=0;j<ny; j++)
			{
				for (k=0; k<nz; k++)
					Qr_tmp[i][j][k] = 1.0;
			}
		}
	}
	//mexPrintf("initialisation OK\n");
	// autres fonctions de motif
	for (i=debut_non_shaper; i<positions.n; i++)
		f_motif[tirages[i].numero_fonction](positions.tableau+i, tirages+i);
	//mexPrintf("calcul OK\n");
}

