#include <string.h>
#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <math.h>

#define T_PI 3.14159265358979				//pi
#define D_PI 6.28318530717958647692			//2*pi

int main(int argc, char **argv)
{
	int i, ncap, Nside, npix;
	double x, y, z, *theta, *phi, fodd;
	long iring, iphi, ip, nl2;
	
	FILE *f_coord;
	
	if (argc != 3)
	{
		printf("Usage : ./healpix Nside fichier\n");
		exit(1);
	}
	
	Nside = atoi(argv[1]);
	f_coord = fopen(argv[2],"w");
	npix = 12*Nside*Nside;
	theta = (double *)calloc(npix, sizeof(double));
	phi = (double *)calloc(npix, sizeof(double));
	
	fprintf(f_coord, "%d\n", 2);
	fprintf(f_coord, "%d\n", npix);
	
	ncap = 2*(Nside*Nside-Nside);
	
	for (i=0; i<npix; i++)
	{
		if(i<ncap)
		{
			iring = (long)(0.5*(1.0+sqrt(1.0+2.0*(double)i)));
			iphi = (long)(i+1)-2.0*iring*(iring-1.0);
			theta[i] = acos(1.0-(double)(iring*iring)/(double)(3*Nside*Nside));
			phi[i] = ((double)iphi-0.5)*T_PI/(2.0*(double)iring);
		}
		else if(i<(12*Nside*Nside-ncap))
		{
			ip = (long)i-(long)ncap;
			iring = (long)(ip/(4.0*(long)Nside))+(long)Nside;
			iphi = ip % (4*(long)Nside)+1;
			fodd = ((iring+(long)Nside)&1) ? 1.0 : 0.5;
			nl2 = 2.0*(long)Nside;
			theta[i] = acos((double)(nl2-iring)*2.0/(double)(3.0*Nside));
			phi[i] = ((double)iphi-fodd)*T_PI/(double)nl2;
		}
		else
		{
			ip = (long)(12*Nside*Nside-i);
			iring = (long)(0.5*(1.0+sqrt((double)(2.0*ip-1.0))));
			iphi = 4.0*iring+1.0-(ip-2.0*iring*(iring-1));
			theta[i] = acos(-1.0+(double)(iring*iring)/(double)(3.0*Nside*Nside));
			phi[i] = ((double)iphi-0.5)*T_PI/(double)(2.0*iring);
		}
		x = sin(theta[i])*cos(phi[i]);
		y = sin(theta[i])*sin(phi[i]);
		z = cos(theta[i]);
		
		fprintf(f_coord, "%f %f %f\n", x, y, z);
	}
	for (i=0; i<npix; i++)
		fprintf(f_coord, "%f\n", theta[i]);
	for (i=0; i<npix; i++)
		fprintf(f_coord, "%f\n", phi[i]);
	
	fclose(f_coord);
	free(theta);

	return 0;
}

