#include "main.h"

void synthese_map(int maillage)
{
	int i, j, count, test = 0;
	double Wi, *phi_int;
	
	pointSphere_t point_sphere;
	double ri;
	
	double x_j, y_i = -1, max_x = 2, max_y = 1, min_x = -2;
	
	int nx, ny;
	double delta_mire, tmp;
	

	for (i=0; i<nb_sommets; i++)
		liste_sommets[i].t = 0.0;
	
	/*******************	Synthese	*****************************/
	
	ny = 21;
	delta_mire = 2.0/(ny-1);
	//printf("delta_mire = %f\n", delta_mire);
	
	ri = delta_mire;
	phi_int = (double *)malloc(2*ny*sizeof(double));
	for (i=0; i<2*ny; i++)
		phi_int[i] = 0.0;
		
	k1 = 0;
	k2 = 1;
	k3 = 0;
	k4 = 2;
	
	printf("synthese ...\n");
	
	for (i=0; i<ny; i++)
	{
		//printf("y_i : %f,	i = %d\n",y_i, i);
		point_sphere.theta = acos(y_i);
		
		//printf("theta = %f\n", point_sphere.theta);
		
		if (i==0 || i==ny-1)
			x_j = 0.0;
		else if (i%2==0)
			x_j = 0.0;
		else
			x_j = delta_mire;
		//x_j = -2;
		
		count = 0;
		while (x_j<max_x)
		{		
			if (sin(point_sphere.theta) <= 0.0001)	
				point_sphere.phi = 0.0;
			else
			{
				tmp = x_j/sin(point_sphere.theta);
				/*if (tmp>1)
					tmp = 1;
				else if (tmp<-1)
					tmp = -1;
				point_sphere.phi = acos(tmp);*/
				//printf("	tmp = %f\n", tmp);
				if (tmp>1 && tmp<2)
					point_sphere.phi = acos(tmp-1)+3.1415;
				else if (tmp>=2)
					point_sphere.phi = 3*3.1415/2.0;
				else
					point_sphere.phi = acos(tmp);
			}
			//printf("	phi = %f", point_sphere.phi);
			
			Wi = W_values(point_sphere.theta, point_sphere.phi, i, ny);
			
			phi_int[count] = point_sphere.phi;
			//printf("	phi_int = %f, count = %d", phi_int[count], count);
			if (i!=0 || i!=ny-1)
			{for (j=1; j<count+1; j++)
			{
				if ((phi_int[count]>=phi_int[j-1]-0.001 && phi_int[count]<=phi_int[j-1]+0.001))
				//(phi_int[count] == phi_int[j-1])
					test = 1;
			}}
			count ++;
			//printf("	test = %d", test);
				
			if (test == 1)
			{
				//printf("	break");
				if (i==0 || i==ny-1)
					x_j = max_x;
				else
					x_j += 2*delta_mire;
			}
			else
			{
				last_step(ri, Wi, point_sphere.theta, point_sphere.phi);

				/*if (i==0 || i==ny-1)
					x_j = 1;
				else
					x_j += 2*delta_mire;*/
				if (i==0 || i==ny-1)
					x_j = max_x;
				else
					x_j += 2*delta_mire;
			}
			//printf("\n");
			test = 0;
		}
		
		if (i==0 || i==ny-1)
			x_j = min_x;
		else if (i%2==0)
			x_j = -2*delta_mire;
		else
			x_j = -delta_mire;
		//x_j = -2;
		
		/*for (j=0; j<2*ny; j++)
			phi_int[j] = 0.0;*/
		
		//count = 0;
		test = 0;
		while (x_j>min_x)
		{		
			if (sin(point_sphere.theta) <= 0.0001)	
				point_sphere.phi = 0.0;
			else
			{
				tmp = x_j/sin(point_sphere.theta);
				/*if (tmp>1)
					tmp = 1;
				else if (tmp<-1)
					tmp = -1;
				point_sphere.phi = acos(tmp);*/
				//printf("	tmp = %f\n", tmp);
				if (tmp<-1 && tmp>-2)
					point_sphere.phi = acos(tmp+1)-3.1415;
				else if (tmp<=-2)
					point_sphere.phi = -T_PI/2.0;
				else
					point_sphere.phi = acos(tmp);
			}

			//printf("	phi = %f", point_sphere.phi);
			
			Wi = W_values(point_sphere.theta, point_sphere.phi, i, ny);
			
			phi_int[count] = point_sphere.phi;
			//printf("	phi_int = %f, count = %d", phi_int[count], count);
			if (i!=0 || i!=ny-1)
			{for (j=1; j<count+1; j++)
			{
				if ((phi_int[count]>=phi_int[j-1]-0.001 && phi_int[count]<=phi_int[j-1]+0.001) || (phi_int[count]>=phi_int[j-1]-0.001 && phi_int[count]<=phi_int[j-1]+0.001))
				//(phi_int[count] == phi_int[j-1])
					test = 1;
			}}
			count ++;
				
			if (test == 1)
			{
				//printf("	break");
				if (i==0 || i==ny-1)
					x_j = min_x;
				else
					x_j -= 2*delta_mire;
			}
			else
			{
				last_step(ri, Wi, point_sphere.theta, point_sphere.phi);

			
				/*if (i==0 || i==ny-1)
					x_j = 1;
				else
					x_j += 2*delta_mire;*/
				if (i==0 || i==ny-1)
					x_j = min_x;
				else
					x_j -= 2*delta_mire;
			}
			//printf("\n");
			test = 0;
		}
		
		y_i += delta_mire;
	}
		
}


void last_step(double ri, double Wi, double theta, double phi)
{
	double delta, coef, A, ri_int;
	point4d_t pt_car;
	// Indices des points correspondants aux bornes de la zone d'influence
	int inf_ind, sup_ind;
	double theta_inf, theta_sup, angle_theta;

	ri_int = Ray+ri;
	pt_car.x = ri*sin(theta)*cos(phi);
	pt_car.y = ri*sin(theta)*sin(phi);
	pt_car.z = ri*cos(theta);

	coef = Ray/ri;
	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;

	inf_ind = dicho_equi_inf(theta_inf);
	sup_ind = dicho_equi_sup(theta_sup);

	if (motif == DISQUE)
		shaper_disque(pt_car, Wi, A, coef, inf_ind, sup_ind);
	else if (motif == COS)
		shaper_cos2(pt_car, Wi, A, ri, coef, inf_ind, sup_ind);
	else //if (motif == EXP)
		shaper_exp(pt_car, Wi, A, ri, coef, inf_ind, sup_ind);
}


/**************************************************************************************/
/*			Disque																	  */
/**************************************************************************************/
void shaper_disque(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++)
	{
		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(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, val;

	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++)
	{
		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;
			if(Wi<0)
			{
				val = -Wi;
				liste_sommets[indic].t = val*(pow(val,CARRE(cos(T_PI*coord/ri)))-1);//pow((Wi?Wi:EPSILON),CARRE(cos(T_PI*coord/ri)));//exp(log(Wi?Wi:EPSILON)*CARRE(1-cos(T_PI*coord/ri)));
			}
			else
			{
				liste_sommets[indic].t += Wi*(1-pow(Wi,CARRE(cos(T_PI*coord/ri))));//pow((Wi?Wi:EPSILON),CARRE(cos(T_PI*coord/ri)));//exp(log(Wi?Wi:EPSILON)*CARRE(1-cos(T_PI*coord/ri)));
			}	
		}
	}
}


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


