#include "main.h"

void synthese_angle(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	*/
		//printf("generation point ...\n");
		u = (double)randu();
		//point_sphere.theta = 0;
		//point_sphere.theta = T_PI/2.0;//2.0*asin(sqrt(u));
		point_sphere.theta = 2.0*asin(sqrt(u));
		u = (double)randu();
		//point_sphere.phi = 0;
		//point_sphere.phi = T_PI/2.0;;//T_PI/2.0;//u*D_PI;
		point_sphere.phi = u*D_PI;
		if (i==0)
		{
			point_sphere.phi = 0.0;
			point_sphere.theta = 0.0;
		}
		/*else if (i==1)
		{
			point_sphere.phi = 0.0;
			point_sphere.theta = T_PI/4.0;
		}
		else if (i==2)
		{
			point_sphere.phi = T_PI/4.0;
			point_sphere.theta = T_PI/4.0;
		}
		else if (i==3)
		{
			point_sphere.phi = T_PI/2.0;
			point_sphere.theta = T_PI/4.0;
		}
		else if (i==4)
		{
			point_sphere.phi = 0.0;
			point_sphere.theta = T_PI/2.0;
		}
		else if (i==5)
		{
			point_sphere.phi = T_PI/4.0;
			point_sphere.theta = T_PI/2.0;
		}
		else if (i==6)
		{
			point_sphere.phi = T_PI/2.0;
			point_sphere.theta = T_PI/2.0;
		}
		else if (i==7)
		{
			point_sphere.phi = 0.0;
			point_sphere.theta = 3*T_PI/4.0;
		}
		else if (i==8)
		{
			point_sphere.phi = T_PI/4.0;
			point_sphere.theta = 3*T_PI/4.0;
		}
		else if (i==9)
		{
			point_sphere.phi = T_PI/2.0;
			point_sphere.theta = 3*T_PI/4.0;
		}
		else 
		{
			point_sphere.phi = u*D_PI;
			point_sphere.theta = 2.0*asin(sqrt(u));
		}*/

		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	*/
		//printf("transfo coord ...\n");
		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	*/
		//printf("A ...\n");
		coef = Ray/point_sphere.ri;
		//if (ri>=2*T_PI*Ray)
		//	delta = T_PI*Ray;
		//else
			delta = ri/2.0;
		A = 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)
		{
			//printf("subdivision ...\n");
			if ((i>0 && i<10) || (i>nb_points-10 && i<nb_points))	printf("sub ...\n theta_inf = %lf	theta_sup = %lf\n", theta_inf, theta_sup);
			inf_ind = dicho_sub_inf(theta_inf);
			sup_ind = dicho_sub_sup(theta_sup);
		}
		else if (maillage == EQUIANGULAIRE)
		{
			//printf("equiangulaire ...\n");
			inf_ind = dicho_equi_inf(theta_inf);
			sup_ind = dicho_equi_sup(theta_sup);
		}
		else //if (maillage == HEALPIX)
		{
			//printf("healpix ...\n");
			inf_ind = dicho_healpix_inf(theta_inf);
			sup_ind = dicho_healpix_sup(theta_sup);
		}
		
		
			/*	Calcul de la cascade	*/
		if (motif == DISQUE)
		{
			//printf("disque ...\n");
			shaper_disque_angle(pt_car, Wi[i], A, coef, inf_ind, sup_ind);
		}
		else if (motif == CERCLE)
		{
			//shaper_cercle_angle(pt_car, point_sphere, Wi[i], A, ri, coef, inf_ind, sup_ind);
			shaper_etoile_angle(pt_car, point_sphere, coef, Wi[i], A, inf_ind, sup_ind);
		}
		else if (motif == COS)
		{
			//if ((i>0 && i<10) || (i>nb_points-10 && i<nb_points))	printf("cos ...\n inf_ind = %d	sup_ind = %d\n", inf_ind, sup_ind);
			shaper_cos2_angle(pt_car, Wi[i], A, ri, coef, inf_ind, sup_ind);
		}
		else //if (motif == EXP)
		{
			//printf("exp ...\n");
			shaper_exp_angle(pt_car, Wi[i], A, ri, coef, inf_ind, sup_ind);
		}
	}
}

/**************************************************************************************/
/*			Disque																	  */
/**************************************************************************************/
void shaper_disque_angle(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;
	}
}

/**************************************************************************************/
/*			Cercle																	  */
/**************************************************************************************/
void shaper_cercle_angle(point4d_t pt_car, pointSphere_t point_sphere, double Wi, double A, double ri, double coef, int inf_ind, int sup_ind)
{
	int j, indic, i=0;
	double x, y, z, norm, den, num, angle, coord, r;
	double temp, phi_P;
	double phi, theta, dist;
	double x1, y1;
	//double x1, y1, z1, norm1, den1, num1, angle1, u;
	//double x2, y2, z2, norm2, den2, num2, angle2;

	x = sin(point_sphere.theta)*cos(point_sphere.phi);
	y = sin(point_sphere.theta)*sin(point_sphere.phi);
	z = cos(point_sphere.theta);
	//printf("x = %f, y : %f, z = %f\n", x, y, z);
	/*x = coef*pt_car.x;
	y = coef*pt_car.y;
	z = coef*pt_car.z;*/
	norm = sqrt(CARRE(x)+CARRE(y)+CARRE(z));
	//printf("%f, %f, %f, %f, %d, %d\n", x, y, z, A, inf_ind, sup_ind);
	/*u = (double)randu();
	x1 = coef*point_sphere.ri*sin(point_sphere.theta+A)*cos(point_sphere.phi);
	y1 = coef*point_sphere.ri*sin(point_sphere.theta+A)*sin(point_sphere.phi);
	z1 = coef*point_sphere.ri*cos(point_sphere.theta+A);
	norm1 = sqrt(CARRE(x1)+CARRE(y1)+CARRE(z1));
	x2 = coef*point_sphere.ri*sin(point_sphere.theta)*cos(point_sphere.phi+A);
	y2 = coef*point_sphere.ri*sin(point_sphere.theta)*sin(point_sphere.phi+A);
	z2 = coef*point_sphere.ri*cos(point_sphere.theta);
	norm2 = sqrt(CARRE(x2)+CARRE(y2)+CARRE(z2));*/

	for(j=inf_ind; j<=sup_ind; j++)
	{
		if (maillage == 0)
			indic = liste_triee_x[j];
		else //if (maillage == 1 || maillage == 2)
			indic = j;
		r = sqrt(CARRE(liste_sommets[indic].x)+CARRE(liste_sommets[indic].y)+CARRE(liste_sommets[indic].z));
		den = r*norm;
		num = liste_sommets[indic].x*x+liste_sommets[indic].y*y+liste_sommets[indic].z*z;
		angle = acos(num/den);
		/*den1 = r*norm1;
		num1 = liste_sommets[indic].x*x1+liste_sommets[indic].y*y1+liste_sommets[indic].z*z1;
		angle1 = acos(num1/den1);
		den2 = r*norm2;
		num2 = liste_sommets[indic].x*x2+liste_sommets[indic].y*y2+liste_sommets[indic].z*z2;
		angle2 = acos(num2/den2);*/
		/*if ((angle<=A && angle>=(A/1.5))||(angle<=(A/4) && angle>=(A/6)))
		{
			coord = angle*Ray;
			liste_sommets[indic].t *= Wi;//*CARRE(cos(T_PI*sqrt(coord)/A));
		}*/
		/*if (angle<=A)
		{
			coord = angle*Ray;
			liste_sommets[indic].t *= Wi*CARRE(cos(T_PI*sqrt(coord)/A));
		}*/
		/*if (angle<=A && angle1<=A)
		{
			//coord = angle*Ray;
			liste_sommets[indic].t *= Wi;//*CARRE(cos(T_PI*sqrt(coord)/A));
		}*/
		/*if (angle<=A && angle1<=A && angle2<=A)
		{
			coord = angle*Ray;
			liste_sommets[indic].t *= Wi*CARRE(cos(T_PI*sqrt(coord)/A));
		}*/
		
		x1 = (liste_sommets[indic].x-x);//cos(point_sphere.theta)*(liste_sommets[indic].x-x) - sin(point_sphere.theta)*(liste_sommets[indic].y-y);
		y1 = cos(point_sphere.theta)*(liste_sommets[indic].y-y) - sin(point_sphere.theta)*(liste_sommets[indic].z-z);//sin(point_sphere.theta)*(liste_sommets[indic].x-x) + cos(point_sphere.theta)*(liste_sommets[indic].y-y);
		
		phi_P = acos((x1)/sqrt(CARRE(x1)+CARRE(y1)));
		//temp = A;
		//dist = sqrt(1-CARRE(cos(angle)));
		dist = sqrt(CARRE(liste_sommets[indic].x-x)+CARRE(liste_sommets[indic].y-y)+CARRE(liste_sommets[indic].z-z));
		temp = A*(0.1*cos(8*phi_P)+0.5);
		if (dist<=temp)
		{
			i ++;
			//if(i<=5)
			//	printf("dist = %f, temp = %f\n", dist, temp);
			//coord = angle*Ray;
			liste_sommets[indic].t *= Wi;//*CARRE(cos(T_PI*sqrt(coord)/A));
		}
		
		/*phi = acos(liste_sommets[indic].x/(sqrt(CARRE(liste_sommets[indic].x)+CARRE(liste_sommets[indic].y))));
		theta = acos(liste_sommets[indic].z/r);
		dist = phi/cos(theta);
		if (dist<=(point_sphere.theta+A) && dist>=(point_sphere.theta-A))
		{
			liste_sommets[indic].t *= Wi;
		}*/
	}
}

/**************************************************************************************/
/*			Shaper cos^2															  */
/**************************************************************************************/
void shaper_cos2_angle(point4d_t pt_car, double Wi, double A, double ri, double coef, int inf_ind, int sup_ind)
{
	int j, indic;
	double coord, 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)
		{
			coord = angle*Ray;
			liste_sommets[indic].t *= exp(log(Wi?Wi:EPSILON)*CARRE(cos(T_PI*coord/ri)));
		}
	}
}


/**************************************************************************************/
/*			Shaper exp																  */
/**************************************************************************************/
void shaper_exp_angle(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, 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)
		{
			coord = CARRE(angle*Ray);
			liste_sommets[indic].t *= exp(log(Wi?Wi:EPSILON)*exp(1.0-1.0/(1.0+EPSILON-CARRE(coord/ri2_4))));
		}
	}
}


/**************************************************************************************/
/*			Shaper etoile															  */
/**************************************************************************************/
void shaper_etoile_angle(point4d_t pt_car, pointSphere_t point_sphere, double coef, double Wi, double A, int inf_ind, int sup_ind)
{
	int j, indic;
	double x, y, z, x1, y1, z1, x2, y2;
	double temp, phi_P, dist;
	double x_P, y_P, z_P;
	double sin1, cos1;
	double angle;

	sin1 = sin(point_sphere.theta);
	cos1 = cos(point_sphere.theta);

	x = coef*pt_car.x;
	y = coef*pt_car.y;
	z = coef*pt_car.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;
		
		x_P = liste_sommets[indic].x-x;
		y_P = liste_sommets[indic].y-y;
		z_P = liste_sommets[indic].z-z;
		
		dist = sqrt(CARRE(x_P)+CARRE(y_P)+CARRE(z_P));
		
		x1 = cos(-point_sphere.phi)*liste_sommets[indic].x - sin(-point_sphere.phi)*liste_sommets[indic].y;
		y1 = sin(-point_sphere.phi)*liste_sommets[indic].x + cos(-point_sphere.phi)*liste_sommets[indic].y;
		z1 = liste_sommets[indic].z;
		
		x2 = cos(-point_sphere.theta)*x1 + sin(-point_sphere.theta)*z1;
		y2 = y1;

		phi_P = acos((x2)/sqrt(CARRE(x2)+CARRE(y2)));
		
		// Fonction de l'angle
			//Fleur
		temp = A*(0.3*cos(8*(phi_P))+0.35);
			
			//Coeur
		//temp = A*((0.3*cos(1*(phi_P))+0.5)-(0.1*sin(3*(phi_P))+0.3));
			
		if (dist<=temp)
			liste_sommets[indic].t *= Wi;
			
			//Coeur qui bave ...
		/*angle = D_PI*(double)randu();
		temp = A*((0.3*cos(1*(phi_P))+0.5)-(0.1*sin(3*(phi_P))+0.3));
		if(fabs(cos(angle)*dist)<=temp && fabs(sin(angle)*dist)<=temp)
			liste_sommets[indic].t *= Wi;*/
			
			//Disque qui bave
		/*angle = D_PI*(double)randu();
		if(fabs(cos(angle)*dist)<=0.5*A && fabs(sin(angle)*dist)<=0.5*A)
			liste_sommets[indic].t *= Wi;*/
		
		/*angle = D_PI*(double)randu();
		temp = cos(angle)*A*((0.3*cos(1*(phi_P))+0.5)-(0.1*sin(3*(phi_P))+0.3));
		if(dist<=temp)
			liste_sommets[indic].t *= Wi;*/
		
	}
}

