#include    <stdio.h>
#include    <stdlib.h>
#include    <math.h>
#include    <unistd.h>
#include    <gsl/gsl_sf_legendre.h>
#include    <plplot/plplot.h>

#ifndef PI
#define PI    3.1415926535897932384
#endif

#define		modulus(a)		(a > 0 ? a : -a)

typedef struct quant_num
{
	int l, m, n;
} quant_num;

const long unsigned int NUM_PTS = 5000;
const PLINT XMAX = 1;
const PLINT YMAX = 1;

void aLegendre(PLFLT *x, PLFLT *y, quant_num *numbers)
{
	long unsigned int i = 0;
	float step = 2.0/NUM_PTS, t = -1.0;
	
	while(t < 1.0 && i < NUM_PTS)
	{
		x[i] = (PLFLT)t;
		y[i] = (PLFLT)gsl_sf_legendre_sphPlm(numbers->l, abs(numbers->m), (const double)t);
		i++;
		t += step;
	}
}

void drawlegend(int n)
{
	PLINT 	nLegend = n+1;
	char 	*text[nLegend];
	int i;
	PLINT	opt_array[nLegend];
	PLINT	text_colors[nLegend];
	PLINT	line_colors[nLegend];
	PLINT	line_styles[nLegend];
	PLINT	line_widths[nLegend];
	PLINT 	symbol_numbers[nLegend];
	PLINT 	symbol_colors[nLegend];
	PLINT 	symbols[nLegend];
	PLFLT	symbol_scales[nLegend];
	PLINT	box_colors[nLegend] = {15};
	
	symbol_numbers[n] = n;
	symbol_colors[n] = n;
	symbols[n] 		 = '.';
	symbol_scales[n] = 1.;
	
	for(i = 0; i <= n; i++)
	{
		text[i] = malloc(15 * sizeof(char));
		if(i>0)
			sprintf(text[i], "l=5, &#124;m&#124;=%d",i);
		else
			sprintf(text[i], "l=5, m=%d",i);
		line_colors[i]	= i+1;
		line_styles[i]	= 1;
		line_widths[i]	= 1;
		text_colors[i] 	= 15;
		opt_array[i] 	= PL_LEGEND_LINE;
	}
	
	pllegend(0, 0.72, .87, .05, 15,
		nLegend, opt_array, 0.5, 0.7, 1.5, 0,
		text_colors, (const char**) text, 
		box_colors, NULL, NULL, 
		line_colors, line_styles, line_widths,
		symbol_colors, symbol_scales, symbol_numbers, symbols);
	for(i = 0; i <= n; i++)
		free(text[i]);
}

int main(void)
{
    PLFLT X[NUM_PTS], Y[NUM_PTS];
    quant_num test_num[6];
    int i;
	
    plscol0(0, 255, 255, 255);
    plinit();
    plscol0(15, 0, 0, 0);
    plcol(15);
    plenv(-XMAX,XMAX, -YMAX,YMAX, 0, 2);

    pllab("x", "P#dl#u#um#d(x)", "Associated Legendre Polynomials");
    plbox( "bcnst", 0, 0, "bcnstv", 0, 0);    plcol(1);
    for(i = 0; i <= 5; i++)
    {
		test_num[i].l = 5; 	test_num[i].m = i; 	test_num[i].n = 1; 
		orbital(X, Y, &test_num[i]); plcol(i+1);
		plline(NUM_PTS, X, Y);
    }

    drawlegend(test_num[0].l);
    
    plend();

    return 0;
}