#include "main.h"


/****************************************/
/*	Lecture et ecriture des fichiers	*/
/****************************************/

void lecture_equi(FILE *file1)
{
	int i;
	double x, y, z;
	
	printf("Maillage equiangulaire\n");
	fscanf(file1, "%d\n", &nb_sommets);
	nb_sommets = CARRE(nb_sommets);
	printf("%d sommets sont utilises\n", nb_sommets);
	liste_sommets = (point4d_t *)malloc(nb_sommets*sizeof(point4d_t));
	for (i=0; i<nb_sommets; i++)
	{
		fscanf(file1, "%lf %lf %lf\n", &x, &y, &z);
		liste_sommets[i].x = x*Ray;
		liste_sommets[i].y = y*Ray;
		liste_sommets[i].z = z*Ray;
		liste_sommets[i].t = 1.0;
	}
	liste_triee_y = (double *)malloc(nb_sommets*sizeof(double));
	for (i=0; i<nb_sommets; i++)
		fscanf(file1, "%lf\n", &liste_triee_y[i]);
}

void ecriture_equi(FILE *file2)
{
	int i;

	fprintf(file2, "%d\n", (int)sqrt(nb_sommets));
	//Coordonnees des points de chaque anneau
	for (i=0; i<nb_sommets; i++)
	{
		fprintf(file2, "%f %f %f %f\n", liste_sommets[i].x/Ray, liste_sommets[i].y/Ray, liste_sommets[i].z/Ray, liste_sommets[i].t);
	}

}


/********************************/
/*	Recherche dichotomique		*/
/********************************/
int dicho_equi_inf(double theta_inf)
{
	int ind1 = 0, ind2 = nb_sommets-1, inf_ind = -10, ind_m;
	double valeur_test;
	
	while (ind1<ind2)
	{
		ind_m = (int)((ind1+ind2)/2.0);
		valeur_test = liste_triee_y[ind_m];			
		if (valeur_test==theta_inf)
			inf_ind = ind_m;
		else if (valeur_test>theta_inf)
			ind2 = ind_m-1;
		else if (valeur_test<theta_inf)
			ind1 = ind_m+1;
	}	
	while (ind2>1 && ind2<nb_sommets && liste_triee_y[ind2]==liste_triee_y[ind2-1])
		ind2 = ind2-1;
	if (inf_ind == -10)	inf_ind = ind2;
	
	return (int) inf_ind;
}

int dicho_equi_sup(double theta_sup)
{
	int ind1 = 0, ind2 = nb_sommets-1, sup_ind = -10, ind_m;	
	double valeur_test;
		
	while (ind1<ind2)
	{
		ind_m = (int)((ind1+ind2)/2.0);
		valeur_test = liste_triee_y[ind_m];				
		if (valeur_test==theta_sup)
			sup_ind = ind_m;
		else if (valeur_test>theta_sup)
			ind2 = ind_m-1;
		else if (valeur_test<theta_sup)
			ind1 = ind_m+1;
	}
	while (ind1>0 && ind1<(nb_sommets-1) && liste_triee_y[ind1]==liste_triee_y[ind1+1])
		ind1 = ind1+1;
	if (sup_ind == -10)	sup_ind = ind1;
	
	return (int) sup_ind;
}


/******************************************************************/
/*	Calcul de la valeur du multiplieur en fonction de sa position */
/******************************************************************/
double W_values(double theta, double phi, int i, int ny)
{
	double *val_Wi, Wi;
	
	val_Wi = (double*)malloc(4*sizeof(double));
	val_Wi[0] = 0.001;
	val_Wi[1] = 0.03;
	val_Wi[2] = 0.01;
	val_Wi[3] = 0.005;
	
	if (theta<3.1415/2.0 && phi<3.1415/2.0)
		Wi = 0.001;
	else if (theta<3.1415/2.0 && phi>3.1416/2.0)
		Wi = 0.01;
	else if (theta>3.1416/2.0 && phi<3.1415/2.0)
		Wi = 0.03;
	else if (theta>3.1416/2.0 && phi>3.1416/2.0)
		Wi = 0.005;
	else if (theta>3.1415/2.0 && theta<3.1416/2.0 && phi>3.1415/2.0 && phi<3.1416/2.0)
		Wi = 0.007;
	else if (theta<3.1415/2.0 && phi>3.141/2.0 && phi<3.142/2.0)
	{
		Wi = val_Wi[k1];
		//printf("k1 %d\n", k1);
		k1 += 2;
		if (k1==4)
		{	k1 = 0; }
		//printf("Wi = %f\n", Wi);
	}
	else if (theta>3.1416/2.0 && phi>3.1415/2.0 && phi<3.1416/2.0)
	{
		Wi = val_Wi[k2];
		k2 += 2;
		if (k2==5)
			k2 = 1;
	}
	else if (theta>3.1415/2.0 && theta<3.1416/2.0 && phi<3.1415/2.0)
	{
		Wi = val_Wi[k3];
		k3 ++;
		if (k3==2)
			k3 = 0;
	}
	else if (theta>3.1415/2.0 && theta<3.1416/2.0 && phi>3.1416/2.0)
	{
		Wi = val_Wi[k4];
		k4 ++;
		if (k4==4)
			k4 = 2;
	}
	//Wi = 0.5;
	
	if (i==0 || i==ny-1)
		Wi = Wi;
	else if (i%2==0)
		Wi = Wi;
	else
		Wi = -Wi;
	
	return Wi;
	
}

