#include "main.h"


void Ylm(int order, int level)
{
	int i, B, j, ind1, ind2, inf_ind, sup_ind, ind_m, k, kk, N;
	double theta, theta_1, val, val_test, coeff, x, y, z, norm, angle_theta, angle_phi;
	
	double *workspace, *storeplm;
	
	FILE *file;
	
	if(level<order)
	{
		fprintf(stderr, "Erreur : level >= order\n");
		exit(1);
	}
	
	B = (int)sqrt(nb_sommets)/2;
	workspace = (double *)malloc(16*B * sizeof(double));
	storeplm = (double *)malloc(2*B*(B-order) * sizeof(double));
	
	PmlTableGen(B, order, storeplm, workspace);

	
	printf("niveau et ordre de l'harmonique spherique : %d, %d\n", level, order);
	
	coeff = sqrt((2*level+1)/(4*T_PI));
	angle_theta = 0;//T_PI/2;
	angle_phi = 0;//T_PI/2;
	N = 25;
	x = sin(angle_theta)*cos(angle_phi);
	y = sin(angle_theta)*sin(angle_phi);
	z = cos(angle_theta);
	
	norm = sqrt(CARRE(x)+CARRE(y)+CARRE(z));
	
	file = fopen("ind.txt","w");
	for (j=0; j<(B); j++)
	{
		theta = (2*j+1)*T_PI/(N*B)-angle_theta;
				
		ind1 = 0;
		ind2 = nb_sommets-1;
		inf_ind = -10;
		sup_ind = -10;
		while(ind1<ind2)
		{
			ind_m = (int)((ind1+ind2)/2.0);
			val_test = liste_triee_y[ind_m];
			if (val_test==theta)
				inf_ind = ind_m;
			else if (val_test>theta)
				ind2 = ind_m-1;
			else if (val_test<theta)
				ind1 = ind_m+1;
		}
		while ((ind2-1)>0 && ind2<nb_sommets && liste_triee_y[ind2]==liste_triee_y[ind2-1])
			ind2 = ind2-1;
		inf_ind = ind2;
		
		theta_1 = (2*(j+1)+1)*T_PI/(N*B)+angle_theta;
		ind1 = 0;
		ind2 = nb_sommets-1;
		while(ind1<ind2)
		{
			ind_m = (int)((ind1+ind2)/2.0);
			val_test = liste_triee_y[ind_m];
			if (val_test==theta_1)
				inf_ind = ind_m;
			else if (val_test>theta_1)
				ind2 = ind_m-1;
			else if (val_test<theta_1)
				ind1 = ind_m+1;
		}
		while ((ind2-1)>0 && ind2<nb_sommets && liste_triee_y[ind2]==liste_triee_y[ind2-1])
			ind2 = ind2-1;
		sup_ind = ind2-1;
		
		val = storeplm[(level-order)*(2*B)+j];
		/*if(val<0)
			val=-1;
		else if(val>0)
			val=1;*/

		fprintf(file,"%f %d, %d %f\n", theta, inf_ind, sup_ind, val);
		
		k=0;
		kk=inf_ind;
		if(order!=0)
			ind_m=inf_ind+(int)((2*B)/(2)/order);
		
		if(order==0)
		{
			for (i=inf_ind; i<=sup_ind; i++)
			{
				val_test = acos((liste_sommets[i].x*x+liste_sommets[i].y*y+liste_sommets[i].z*z)/(sqrt(CARRE(liste_sommets[i].x)+CARRE(liste_sommets[i].y)+CARRE(liste_sommets[i].z))*norm));
				if(val_test<=((2*(j+1)+1)*T_PI/(N*B)) && val_test>=((2*j+1)*T_PI/(N*B)))
					liste_sommets[i].t = val*coeff;
			}
		}
		else
		{
			//i=inf_ind;
			ind_m=inf_ind+(int)((2*B)/(2)/order);
			kk=inf_ind;
			k=0;
			while (ind_m<nb_sommets && i<=sup_ind)
			{
				for (i=kk; i<=ind_m; i++)
				{
					if(k%2==0)
						liste_sommets[i].t = val*coeff;
					else
						liste_sommets[i].t = -val*coeff;
				}
				kk = ind_m;				
				ind_m = ind_m+(int)((2*B)/(2)/order);
				k++;
			}
		}
	}
	
	fclose(file);
	
	free(workspace);
	free(storeplm);

}



