/*
 * Program to calculate points on an ellipse
 * 
 * Extract from:
 * 'Programmers Guide to PC&PS/2 Video Systems", Richard wilton
 * Microsoft Press, ISBN 1-55615-103-9
 *
 * This program calculates the points of an ellipse by calculating the
 * points in the first quadrant and reflecting them into the other 3
 * quadrants.  There are two different algorithms used, one for when
 * dy/dx < -1 and one for when dy/dx >-1.  Refer to the book for
 * a derivation of the algorithm; it is very similar to the derivation
 * of bresenham's line algorithm.  There are problems related with drawing
 * very small ellipses; the recommends using a more intense graphics
 * mode.  Also, the iteration will not terminate if a0 or b0 is 0,
 * You should check for this.
 *
 */

#define TESTING

#ifdef TESTING
void SetPixel(x,y,n)
{
  printf("%4d %4d  ",x,y);
}
#endif TESTING
void Set4Pixels(int x,int y,int xc,int yc,int n)
{
	SetPixel(xc + x, yc + y, n);
	SetPixel(xc - x, yc + y, n);
	SetPixel(xc + x, yc - y, n);
	SetPixel(xc - x, yc - y, n);
#ifdef TESTING
	printf("\n");
#endif TESTING

}

/*
 * xc,yc = center 
 * a0,b0 = major and minor axis ( radius) 
 */
void Ellipse(int xc,int yc,int a0,int b0, int PixelValue)
{
	int x = 0;
	int y = b0;

	long a= a0;
	long b = b0;


	long Asquared = a*a;
	long TwoAsquared = 2 * Asquared;
	long Bsquared = b*b;
 	long TwoBsquared = 2 * Bsquared;

	long d;
	long dx,dy;

	d = Bsquared - Asquared*b + Asquared/4L;
	dx = 0;
	dy = TwoAsquared *b;

  	while (dx < dy) {
  		Set4Pixels(x,y,xc,yc,PixelValue);
	
		if (d > 0L) {
			--y;
			dy -= TwoAsquared;
			d-=dy;
		}

		++x;
		dx+=TwoBsquared;
		d += Bsquared + dx;
	}

	d += (3L *(Asquared-Bsquared)/2L -(dx+dy)) / 2L;

 	while (y >=0) {
		Set4Pixels(x,y,xc,yc,PixelValue);

		if (d < 0L) {
			++x;
			dx+= TwoBsquared;
			d+= dx;
		}

		--y;
		dy -= TwoAsquared;
		d+= Asquared - dy;
	}
}
#ifdef TESTING
main()
{
	/* Draw a circle of radius 40 at center 100,100 */
	Ellipse( 100,100,40,40,0);
}
#endif TESTING