//airy pattern using the bessel function
//www.timestocome.com

//to compile on OSX
//gcc airy.c -o airy -framework GLUT -framework OpenGL


#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <GLUT/glut.h>



void display (void){

	typedef GLfloat point2[2];			//define a point array

		//table of first ten 1st order bessel numbers

	double bessel[11] = {0.0, 0.44005, 0.57672, 0.33905, -0.06604,
							-0.32757, -0.27668, -0.00468, 0.23463,
							0.24531, 0.04347};

	int i;

	GLfloat size = 3.0;
	GLfloat x, y;
	GLfloat I=0.0, I0=1.0;
	glClear(GL_COLOR_BUFFER_BIT);		//clear window





	//airy pattern

	for(i=0; i<11; i++){

		glColor3f(1.0, 0.0, 0.0);			//draw color

		I = I0*pow(((2*bessel[i])/i),2);

		x = i*25+25;			
		y = I*25+25;						//plot points

		glBegin(GL_POLYGON);				//draw a square pixel

					glVertex2f(x,			y);
					glVertex2f(x + size,	y);
					glVertex2f(x + size,	y + size);
					glVertex2f(x,			y + size);

		glEnd();

	}



	//bessel function	

	for(i=0; i<11; i++){

		glColor3f(0.0, 1.0, 0.0);			//draw color

		x = i*25+25;			
		y = bessel[i]*25+25;						//plot points

		glBegin(GL_POLYGON);				//draw a square pixel

					glVertex2f(x,			y);
					glVertex2f(x + size,	y);
					glVertex2f(x + size,	y + size);
					glVertex2f(x,			y + size);

		glEnd();

	}



	//lagrange estimate of bessel function
	//bessel function	

	for(i=0; i<9; i++){

		glColor3f(0.0, 0.0, 1.0);			//draw color

		x = i*25+25;			

		y = ((bessel[i]-bessel[i+2])/(bessel[i+1]-bessel[i+2])*bessel[i+1] +
			(bessel[i]-bessel[i+1])/(bessel[i+2]-bessel[i+1])*bessel[i+2])*25+25;

		glBegin(GL_POLYGON);				//draw a square pixel
					glVertex2f(x,			y);
					glVertex2f(x + size,	y);
					glVertex2f(x + size,	y + size);
					glVertex2f(x,			y + size);

		glEnd();

	}



	glFlush();							//plot quickly..  only benefit if on a network

}





void myinit (void){

//attributes

	glClearColor(1.0, 1.0, 1.0, 0.0);	//backround
	glColor3f(1.0, 0.0, 0.0);			//draw color



//set up viewing

	glMatrixMode(GL_PROJECTION);
	gluOrtho2D(0.0, 500.0, 0.0, 250.0);	//2d coordinate sys lower left corner 0,0
	glMatrixMode(GL_MODELVIEW);

}





void main(int argc, char **argv)

{

	glutInit(&argc, argv);
    glutInitWindowSize(500, 250);
   	glutInitDisplayMode(GLUT_RGB|GLUT_SINGLE);
    (void)glutCreateWindow("airy pattern");
    glutDisplayFunc(display);
	myinit();
	glutMainLoop();



}