//-------------------------------------------------------------------
//	fourdemo.cpp
//
//	This program illustrates the approximation of a function on
//	the interval from 0 to 2*PI by partial sums of terms in the 
//	the interval from 0 to 2*PI by a succession of partial sums 
//	of terms from the corresponding Fourier Series:
//
//	 A0 + A1*cos(x) + B1*sin(x) + A2*cos(2x) + B2*sin(2x) + ... 
//
//	where  A0 = mean-value of f(x) over interval [0,2*PI]
//	
//	Ak = 2 * mean-value of f(x)*cos(kx) over interval [0,2*PI] 
//
//	Bk = 2 * mean-value of f(x)*sin(kx) over interval [0,2*PI] 
//
//	et cetra.  The 'myfunction()' routine can be edited to let
//	us watch alternative curves being dynamically approximated 
//	by the succession of trigonometric polynomials.  
//
//		to compile:  $ g++ fourdemo.cpp -o fourdemo
//		to execute:  $ ./fourdemo
//
//	NOTE: This program is designed for execution in our Harney
//	classroom using the instructor's Gateway NV53A laptop (and
//	the installed overhead display at 1024-by-768 resolution),
//	hence requiring services of our 'gateway.c' kernel module.
//
//	programmer: ALLAN CRUSE
//	date begun: 04 APR 2011
//	revised on: 11 APR 2011 -- for Gateway laptop's radeon GPU
//-------------------------------------------------------------------

#include <stdio.h>	// for printf(), perror() 
#include <fcntl.h>	// for open() 
#include <stdlib.h>	// for exit() 
#include <unistd.h>	// for lseek(), usleep() 
#include <string.h>	// for memcpy(), memset()
#include <math.h>	// for sin(), cos(), exp()
#include <sys/mman.h>	// for mmap()
#include <sys/ioctl.h>	// for ioctl()

#define VRAM_BASE	0xA0000000	// suitable address for mmap
#define GET_CRTC_START_ADDRESS	0	// number for ioctl service
#define PI		3.14159265	// value to be used for PI 
#define N		640		// number of subintervals

const unsigned int	hpitch = 1024;	// native display's horiz pitch
const unsigned int	maxcol = 1024;	// native horizontal resolution  
const unsigned int	maxrow = 768; 	// native vertical resolution
unsigned int		*fb;		// pointer to native frame-buffer

unsigned int	Screensave[ hpitch * maxrow ];	// to preserve screen
unsigned int	Screenback[ hpitch * maxrow ];	// to hold background

unsigned int	window_wide = 640;
unsigned int	window_high = 480;
unsigned int	window_upper = 40;
unsigned int	window_lower = window_upper + window_high;
unsigned int	window_edgeL = (maxcol - window_wide) / 2;
unsigned int	window_edgeR = window_edgeL + window_wide;

double		scaleX = 32.0;
double		scaleY = -32.0;
double		transX = window_edgeL + window_wide / 2.0;
double		transY = window_upper + window_high / 2.0;

double		domain[ N + 1 ];	// collection of end-points
double		thesum[ N + 1 ];	// array of function-heights

char devname[] = "/dev/vram";		// name of the device-file

const unsigned int	white 	= 0x00FFFFFF;
const unsigned int	gray	= 0x003F3F3F;
const unsigned int	yellow	= 0x00FFFF00;
const unsigned int	red	= 0x00FF0000;
const unsigned int	green	= 0x0000FF00;
const unsigned int	cyan	= 0x0000FFFF;

double myfunction( double x )
{
	return	6.0 * exp( - x*x / 32 );	// <--- Nice bell-curve

	// several alternative curves 
//	return	1.0 +  x * (3.0 - x)*(5.0 - x) / 8.0;	
//	return	x;
//	return	cos( 3*x );
//	return	sin( 5*x );
//	return	x * x / 16.0;
//	return	fabs( x - PI/2.0 );
//	return	( ( x - PI ) * ( x - PI ) ) / 8.0;
}

void draw_thesum( int color );
void draw_curve( double function( double ), int freq, int color );
double mean_value( double function( double ), int k, double trig( double ) );


int main( int argc, char **argv )
{
	// map video display-memory into application address-space
	int	fd = open( devname, O_RDWR );
	if ( fd < 0 ) { perror( devname ); exit(1); }
	int	size = lseek( fd, 0, SEEK_END );
	int	prot = PROT_READ | PROT_WRITE;
	int	flag = MAP_FIXED | MAP_SHARED;
	void	*mm = (void *)VRAM_BASE;
	if ( mmap( mm, size, prot, flag, fd, 0 ) == MAP_FAILED )
		{ perror( "mmap" ); exit(1); }
	unsigned int	crtc_start_address = 0;
	if ( ioctl( fd, GET_CRTC_START_ADDRESS, &crtc_start_address ) < 0 )
		{ perror( "ioctl" ); exit(1); }
	unsigned int	fb_offset = crtc_start_address & 0x3FFFFFFF;
	fb = (unsigned int *)( VRAM_BASE + fb_offset );

	// save the screen-image to be restored later
	memcpy( Screensave, &fb[0], sizeof( Screensave ) );

	// erase the desktop 
	memset( &fb[0], 0x00, 4 * hpitch * 560 );
	
	// draw the window-border
	int	rowmin = window_upper, rowmax = window_lower;
	int	colmin = window_edgeL, colmax = window_edgeR;
	int	r = rowmin, c = colmin;
	do { fb[ r * hpitch + c ] = white; ++c; } while ( c < colmax );
	do { fb[ r * hpitch + c ] = white; ++r; } while ( r < rowmax );
	do { fb[ r * hpitch + c ] = white; --c; } while ( c > colmin );
	do { fb[ r * hpitch + c ] = white; --r; } while ( r > rowmin );
	getchar();
	
	// paint the interval from zero to 2*PI
	r = rowmin + 5;
	do 	{
		int	cL = (int)(    0 * scaleX + transX );
		int	cR = (int)( 2*PI * scaleX + transX );
		c = cL;
		do { fb[ r * hpitch + c ] = gray; ++c; } while ( c < cR );
		++r;
		}
	while ( r < rowmax - 5 );  
	getchar();	

	// draw the window's vertical and horizontal axes
	r = (rowmin + rowmax)/2, c = colmin + 5;
	do { fb[ r * hpitch + c ] = white; ++c; } while ( c < colmax-5 );
	c = (colmin + colmax)/2, r = rowmin + 5;
	do { fb[ r * hpitch + c ] = white; ++r; } while ( r < rowmax-5 );
	getchar();

	// initialize the 'domain[]' array
	double	xx = 0.0;
	double	xinc = 10.0 * PI / N;
	for (int i = 0; i <= N/2; i++)
		{
		domain[ N/2 + i ] = xx;
		domain[ N/2 - i ] = -xx;
		xx += xinc;
		}

	// draw the graph of the function we want to approximate
	draw_curve( myfunction, 1, yellow );
	getchar();

	// save this screen for reuse as the standard background
	memcpy( Screenback, &fb[0], sizeof( Screenback ) );

		//===========================//

	// compute the fourier series constant-term
	double	coeff = mean_value( myfunction, 0, cos ); 
	for (int i = 0; i <=N; i++) thesum[ i ] = coeff;

	// loop to compute and display successive partial sums
	for (int k = 1; k < 180; k++)
		{
		memcpy( &fb[0], Screenback, 4 * 550 * hpitch );
		draw_thesum( cyan );
		
		// add the fourier-series freq=k terms to 'thesum[]'
		coeff = 2.0 * mean_value( myfunction, k, cos );
		for (int i = 0; i <= N; i++)
			{
			double	x = domain[ i ];
			double	y = cos( k*x ) * coeff;
			thesum[ i ] += y;
			}
		coeff = 2.0 * mean_value( myfunction, k, sin );
		for (int i = 0; i <= N; i++)
			{
			double	x = domain[ i ];
			double	y = sin( k*x ) * coeff;
			thesum[ i ] += y;
			}

		double	delta = 2*PI / N;
		double	maxdiff = 0.0;
		for (int i = 0; i <= N; i++)
			{
			double	x = domain[ i ];
			if (( x <= delta ) ||( x >= 2*PI - delta )) continue;
			double	y = myfunction( x );
			double	diff = fabs( thesum[ i ] - y );
			if ( diff > maxdiff ) maxdiff = diff;
			}
		// printf( " k=%d  maxdiff = %1.6f \n", k, maxdiff );
		if ( maxdiff < 0.05 ) break;
		usleep( 500000 );
		}
	printf( "\n" );
	getchar();

	// redraw the graph of the function we want to approximate
	memcpy( &fb[0], Screenback, sizeof( Screenback ) );
	draw_curve( myfunction, 1, yellow );
	getchar();
	
	// recover the original saved desktop screen-image
	memcpy( &fb[0], Screensave, sizeof( Screensave ) );
}

void draw_pixel( int x, int y, int color )
{
	// perform clipping 
	if (( x < window_edgeL+5 )||( x > window_edgeR-5 )) return; 
	if (( y < window_upper+5 )||( y > window_lower-5 )) return; 
	fb[ y * hpitch + x ] = color;
}

void exchange_ends( int &x1, int &y1, int &x2, int &y2 )
{ 
	int	xx, yy;

	xx = x2; x2 = x1; x1 = xx; yy = y2; y2 = y1; y1 = yy; 
}

void octant_odd( int x1, int y1, int x2, int y2, int color )
{
	// make sure  x1 <= x2 
	if ( x1 > x2 ) exchange_ends( x1, y1, x2, y2 );

	int	xinc = ( x2 > x1 ) ? 1 : -1;
	int	yinc = ( y2 > y1 ) ? 1 : -1;

	int	dx = ( x2 - x1 ) * xinc;
	int	dy = ( y2 - y1 ) * yinc;
	int	errorterm = -dx;
	for (int x = x1, y = y1; x <= x2; x++)
		{
		draw_pixel( x, y, color );
		errorterm += ( dy + dy );
		if ( errorterm >= 0 ) 
			{
			y += yinc;
			errorterm -= ( dx + dx );
			}
		}
}

void octant_even( int x1, int y1, int x2, int y2, int color )
{
	// make sure  y1 <= y2 
	if ( y1 > y2 ) exchange_ends( x1, y1, x2, y2 );

	int	xinc = ( x2 > x1 ) ? 1 : -1;
	int	yinc = ( y2 > y1 ) ? 1 : -1;

	int	dx = ( x2 - x1 ) * xinc;
	int	dy = ( y2 - y1 ) * yinc;
	int	errorterm = -dy;
	for (int x = x1, y = y1; y <= y2; y++)
		{
		draw_pixel( x, y, color );
		errorterm += ( dx + dx );
		if ( errorterm >= 0 ) 
			{
			x += xinc;
			errorterm -= ( dy + dy );
			}
		}
}

void draw_line( int x1, int y1, int x2, int y2, int color )
{
	// use the Bresenham line-drawing algorithm to draw 
	// the straight line connecting (x1,y1) to (x2,y2)
	int	deltaX = ( x1 < x2 ) ? x2-x1 : x1-x2;
	int	deltaY = ( y1 < y2 ) ? y2-y1 : y1-y2;

	if ( deltaX > deltaY )	// slant is less than 45-degrees
		octant_odd( x1, y1, x2, y2, color );
	else			// slant is more than 45-degrees
		octant_even( x1, y1, x2, y2, color );
}

void draw_curve( double function( double ), int freq, int color )
{
	for (int i = 0; i < N; i++)
		{
		double	x1 = domain[ i ];
		double	y1 = function( freq * x1 );
		double	x2 = domain[ i+1 ];
		double	y2 = function( freq * x2 );

		int	h1 = (int)( x1 * scaleX + transX );
		int	v1 = (int)( y1 * scaleY + transY );
		int	h2 = (int)( x2 * scaleX + transX );
		int	v2 = (int)( y2 * scaleY + transY );

		draw_line( h1, v1, h2, v2, color );
		}
}

void draw_thesum( int color )
{
	for (int i = 0; i < N; i++)
		{
		double	x1 = domain[ i ];
		double	y1 = thesum[ i ];
		double	x2 = domain[ i+1 ];
		double	y2 = thesum[ i+1 ];

		int	h1 = (int)( x1 * scaleX + transX );
		int	v1 = (int)( y1 * scaleY + transY );
		int	h2 = (int)( x2 * scaleX + transX );
		int	v2 = (int)( y2 * scaleY + transY );

		draw_line( h1, v1, h2, v2, color );
		}
}

double mean_value( double function( double ), int k, double trig( double ) )
{
	// This function returns the mean-value of the product-function
	//  y = f(x) * trig( k*x )  over the interval from  0  to  2*PI  
	// computed using the Trapezoid Rule for numerical integration. 
	
	double	twoPI = 2.0 * PI;
	double	delta = twoPI / N;

	double	sum = 0.0;
	double	xL = 0.0, xR = delta;
	do	{
		double	edgeL = function( xL ) * trig( k * xL ); 		
		double	edgeR = function( xR ) * trig( k * xR ); 		
		double	area = ((edgeL + edgeR) / 2.0 ) * delta;
		sum += area;
		xL += delta;
		xR += delta;
		}
	while ( xL < twoPI );
	return	sum / twoPI;	
}