#include "main.h"

void synthese_eclairement(int nb_points, double umax, double *Wi)
{
	// Donnees calculees 
	int i;
	double V, W, g;
	
	// Indices des points correspondants aux bornes de la zone d'influence
	int inf_ind, sup_ind;
	double theta_inf, theta_sup, angle_theta;
	
	// Points uniformement repartis sur la sphere
	point4d_t pt_car;
	pointSphere_t point_sphere;
	double u, ri, delta, coef, A;
	
	printf("synthese ...\n");
	
	for (i=0; i<nb_points; i++)
	{
		/*	Generation de points uniformement repartis sur la sphere	*/
		u = (double)randu();
		point_sphere.theta = 2.0*asin(sqrt(u));
		u = (double)randu();
		point_sphere.phi = u*D_PI;
		u = (double)randu();
        V = 1.0/sqrt(umax-(umax-1.0)*u);
        u = (double)randu();
        W = u*2.0*(4.0/(T_PI*pow(V,3)));
        g = 1/(2.0*T_PI*(1-cos(V/2.0))*CARRE(Ray+V)*V);
        while(W>g)
        {
			u = (double)randu();
            V = 1.0/sqrt(umax-(umax-1.0)*u);
            u = (double)randu();
            W = u*2.0*(4.0/(T_PI*pow(V,3)));
            g = 1.0/(2.0*T_PI*(1-cos(V/2.0))*CARRE(Ray+V)*V);
		}
        ri = L*V;
        //printf("V = %f, W = %f, g = %f et ri = %f\n", V, W, g, ri);
        
		// ri = 1.0/sqrt(umax-(umax-1.0)*u);
		//ri = L/sqrt(umax-(umax-1.0)*u);	//Calcul de la gamme d'echelle de la texture : ri_max pour u->1 et ri_min pour u->0
		point_sphere.ri = Ray+ri;

		/*	Transformation des coordonnees spheriques en coordonnees cartesiennes	*/
		pt_car.x = point_sphere.ri*sin(point_sphere.theta)*cos(point_sphere.phi);
		pt_car.y = point_sphere.ri*sin(point_sphere.theta)*sin(point_sphere.phi);
		pt_car.z = point_sphere.ri*cos(point_sphere.theta);
		
		/*	Calcul de la zone d'influence de chaque point en fonction de la version utilisee	*/
		coef = Ray/point_sphere.ri;
		if (ri>=2*Ray)
			delta = Ray;
		else
			delta = ri/2.0;
		A = asin(delta/Ray);
		
		if (pt_car.z == 0)
			angle_theta = T_PI/2.0;
		else
			angle_theta = atan2(sqrt(CARRE(pt_car.x)+CARRE(pt_car.y)),pt_car.z);
		theta_inf = angle_theta-A;		theta_sup = angle_theta+A;

			// Recherche dichotomique des theta les plus proches de theta_inf et theta_sup
		if (maillage == SUBDIVISION)
		{
			inf_ind = dicho_sub_inf(theta_inf);
			sup_ind = dicho_sub_sup(theta_sup);
		}
		else if (maillage == EQUIANGULAIRE)
		{
			inf_ind = dicho_equi_inf(theta_inf);
			sup_ind = dicho_equi_sup(theta_sup);
		}
		else if (maillage == HEALPIX)
		{
			inf_ind = dicho_healpix_inf(theta_inf);
			sup_ind = dicho_healpix_sup(theta_sup);
		}
		
		
			/*	Calcul de la cascade	*/
		if (motif == DISQUE)
			shaper_disque_eclairement(pt_car, Wi[i], A, coef, inf_ind, sup_ind);
		else if (motif == COS)
			shaper_cos2_eclairement(pt_car, Wi[i], A, ri, coef, inf_ind, sup_ind);
		if (motif == EXP)
			shaper_exp_eclairement(pt_car, Wi[i], A, ri, coef, inf_ind, sup_ind);
		
	}
}

/**************************************************************************************/
/*			Disque																	  */
/**************************************************************************************/
void shaper_disque_eclairement(point4d_t pt_car, double Wi, double A, double coef, int inf_ind, int sup_ind)
{
	int j, indic;
	double x, y, z, norm, den, num, angle;

	x = coef*pt_car.x;
	y = coef*pt_car.y;
	z = coef*pt_car.z;
	norm = sqrt(CARRE(x)+CARRE(y)+CARRE(z));

	for(j=inf_ind; j<=sup_ind; j++)
	{
		if (maillage == 0)
			indic = liste_triee_x[j];
		else if (maillage == 1 || maillage == 2)
			indic = j;
		den = sqrt(CARRE(liste_sommets[indic].x)+CARRE(liste_sommets[indic].y)+CARRE(liste_sommets[indic].z))*norm;
		num = liste_sommets[indic].x*x+liste_sommets[indic].y*y+liste_sommets[indic].z*z;
		angle = acos(num/den);
		if (angle<=A)
			liste_sommets[indic].t *= Wi;
	}
}

/**************************************************************************************/
/*			Shaper cos^2															  */
/**************************************************************************************/
void shaper_cos2_eclairement(point4d_t pt_car, double Wi, double A, double ri, double coef, int inf_ind, int sup_ind)
{
	int j, indic;
	double coord, OH, x, y, z, norm, den, num, angle;

	x = coef*pt_car.x;
	y = coef*pt_car.y;
	z = coef*pt_car.z;
	norm = sqrt(CARRE(x)+CARRE(y)+CARRE(z));

	for(j=inf_ind; j<=sup_ind; j++)
	{
		if (maillage == SUBDIVISION)
			indic = liste_triee_x[j];
		else if (maillage == EQUIANGULAIRE || maillage == HEALPIX)
			indic = j;
		den = sqrt(CARRE(liste_sommets[indic].x)+CARRE(liste_sommets[indic].y)+CARRE(liste_sommets[indic].z))*norm;
		num = liste_sommets[indic].x*x+liste_sommets[indic].y*y+liste_sommets[indic].z*z;
		angle = acos(num/den);
		if (angle<=A)
		{
			OH = cos(A)*Ray;
			coord = tan(angle)*OH;
			liste_sommets[indic].t *= exp(log(Wi?Wi:EPSILON)*CARRE(cos(T_PI*coord/ri)));
		}
	}
}


/**************************************************************************************/
/*			Shaper exp																  */
/**************************************************************************************/
void shaper_exp_eclairement(point4d_t pt_car, double Wi, double A, double ri, double coef, int inf_ind, int sup_ind)
{
	int j, indic;
	double coord, ri2_4, OH, x, y, z, norm, den, num, angle;


	ri2_4 = CARRE(ri)*0.25;

	x = coef*pt_car.x;
	y = coef*pt_car.y;
	z = coef*pt_car.z;
	norm = sqrt(CARRE(x)+CARRE(y)+CARRE(z));

	for(j=inf_ind; j<=sup_ind; j++)
	{
		if (maillage == 0)
			indic = liste_triee_x[j];
		else if (maillage == 1 || maillage == 2)
			indic = j;
		den = sqrt(CARRE(liste_sommets[indic].x)+CARRE(liste_sommets[indic].y)+CARRE(liste_sommets[indic].z))*norm;
		num = liste_sommets[indic].x*x+liste_sommets[indic].y*y+liste_sommets[indic].z*z;
		angle = acos(num/den);
		if (angle<=A)
		{
			OH = cos(A)*Ray;
			coord = CARRE(tan(angle)*OH);
			liste_sommets[indic].t *= exp(log(Wi?Wi:EPSILON)*exp(1.0-1.0/(1.0+EPSILON-CARRE(coord/ri2_4))));
		}
	}
}

