/*
 * =====================
 *
 * Author:	Matt Mazzola
 * Date:	12.08.2010
 * Title:	Spectrum Analyzer
 *
 * Because this program uses three tasks,
 * One to collect data, one to process, and
 * one to display, you can swap out different FFT
 * algorithms with minimal modifications to the other
 * tasks, or you could update the display task
 * to a better graph while preserving the other tasks
 * =====================
 */

#include <assert.h>
#include <fcntl.h>
#include <pthread.h>
#include <rtai.h>
#include <rtai_fifos.h>
#include <rtai_lxrt.h>
#include <rtai_sched.h>
#include <rtai_sem.h>
#include <stdio.h>
#include <stdlib.h>
#include <sys/mman.h>
#include <unistd.h>
#include <stdint.h>
#include <math.h>

//#define DISPLAY_DELTA	// by commenting this you can turn the delta display on/off
						// the delta refers to the time difference between the graphs printing
//#define DEBUG

#ifdef DEBUG
# define DEBUG_PRINT(x) printf x
#else
# define DEBUG_PRINT(x) do {} while (0)
#endif


#define TIMESLICE		10000000			// for Round Robin scheduling

#define NUM_OF_BUFFERS	10
#define SAMPLES			128					// this is not used anymore, see mNumOfSamples

#define PRIORITY		0
#define STACK_SIZE		256
#define MAX_MSG_SIZE	512

#define CONTROL_REGISTER	0x10F00000
#define OPTION_REGISTER		0x22400000
#define COMPLETE_REGISTER	0x10800000

#define Y_MAX					20
#define GRAPH_VERTICAL_SCALE	1

#define ASCII_ESC	27

#define TWOPI		6.28318531
#define FOURPI		12.56637061


// -------------- Not Used in last version -------
// getNextBuffer() function previously used status and stage as
// parameters when searching but because of preemption the status and stages
// were not always updated as expected and caused errors.
// i left them in the Buffer structure but don't use them.
// they could be taken out to reduce memory footprint but
// they're only an int and are insignificant and may find use later
typedef enum{
	OPEN	= 0,
	LOCKED
} Buffer_Status;

typedef enum{
	OVERWRITE	= 0,
	RAW_DATA,
	FFT_DATA
} Buffer_Stage;
// ----------------- See Above ----------------------

typedef enum{
	FFT_WIN_TYP_RECTANGLE	= 0x00,	// rectangle (Box car)
	FFT_WIN_TYP_HAMMING,			// hamming
	FFT_WIN_TYP_HANN,				// hann
	FFT_WIN_TYP_TRIANGLE,			// triangle (Bartlett)
	FFT_WIN_TYP_BLACKMAN,			// blackmann
	FFT_WIN_TYP_FLT_TOP,			// flat top
	FFT_WIN_TYP_WELCH,				// welch
} Fft_Window_Type;

typedef enum{
	FFT_REVERSE	= 0x00,
	FFT_FORWARD
} Fft_Direction;

typedef struct{
	uint8_t			id;
	Buffer_Status	status;
	Buffer_Stage	stage;
	double			*vReal;
	double			*vImag;
	double			peak;
} Buffer;

/* global variables */
Buffer 	mRingBuffer[NUM_OF_BUFFERS];		// create array of buffers
int 	mNumOfSamples		= SAMPLES;		// must factor of 2^n where n = 1,2,3...
int 	mSampleFrequency	= 32000;		// Hz ( 44.1kHz is CD quality )
int 	lastWindowType 		= FFT_WIN_TYP_RECTANGLE;

char mTask_1_name[]	= "acquire data";
char mTask_2_name[]	= "compute fft";
char mTask_3_name[]	= "display data";

RTIME mPeriod;
SEM *sem_1_2;
SEM *sem_2_3;
SEM *sem_3_1;

static RTIME previousTime = 0;
static RTIME currentTime = 0;

//============ Prototypes ================
//==================== FFT ===============
void fft_windowing(double *vData, uint16_t samples, uint8_t windowType, uint8_t dir);
void fft_compute(double *vR, double *vI, uint16_t samples, uint8_t dir);
void fft_complexToMagnitude(double *vR, double *vI, uint16_t samples);
double fft_majorPeak(double *vD, uint16_t samples, double samplingFrequency) ;
void swap( double *x, double *y );
uint8_t exponent( uint16_t value );

//==================== Thread Routines ===
void* get_samples( void *aPtr );
void* compute_fft( void *aPtr );
void* display_data( void *aPtr );

// =================== Helping Functions =
void print_values( double *aVector, int size );
void print_graph( double *aVector, int size, int height );

// notice this takes parameters aStatus/aStage, these are disregarded
void getNextBuffer(
	Buffer 			**aCurrentBufferPtr,
	Buffer 	        *aRingBuffer,
	Buffer_Status 	aStatusCondition,
	Buffer_Stage 	aStageCondition,
	int 			aTotalNumOfBuffers );

// used for testing when not connected to signal generator
void generateSignal(
		double* aVector,
		int aAmplitude,
		int aSignalFrequency,
		int aSampleFrequency,
		int aSamples
		);

int main( int argc, char *argv[] )
{

	if( argc != 2 )
	{
		printf("Error: Enter Sample Number As Argument\n");
		return EXIT_FAILURE;
	}
	else
	{
		mNumOfSamples = (int)atoi( argv[1] );
		printf("%d\n", mNumOfSamples );
	}

	// initialize buffer id's
	int i;
	for( i = 0; i < NUM_OF_BUFFERS; i++ )
	{
		mRingBuffer[i].id		= i+1;
		mRingBuffer[i].status	= OPEN;
		mRingBuffer[i].stage	= OVERWRITE;
		mRingBuffer[i].vReal	= malloc( sizeof( double) * mNumOfSamples );
		mRingBuffer[i].vImag	= malloc( sizeof( double) * mNumOfSamples );
		memset( mRingBuffer[i].vReal, 0, sizeof( double) * mNumOfSamples );
		memset( mRingBuffer[i].vImag, 0, sizeof( double) * mNumOfSamples );
		mRingBuffer[i].peak		= 0;
	}
	DEBUG_PRINT(("%d Buffers Initialized\n", i));


	int mTamplePeriodInNanoSeconds = 1000000000/mSampleFrequency;

    sem_1_2 = rt_typed_sem_init( nam2num("sem1"), 0, CNT_SEM );
    sem_2_3 = rt_typed_sem_init( nam2num("sem2"), 0, CNT_SEM );
    sem_3_1 = rt_typed_sem_init( nam2num("sem3"), NUM_OF_BUFFERS, CNT_SEM );

    DEBUG_PRINT(("Period: %d\n", mTamplePeriodInNanoSeconds));
	mPeriod = start_rt_timer( nano2count( mTamplePeriodInNanoSeconds ) );

	pthread_t mAcquireData;
	pthread_t mComputeFFT;
	pthread_t mDisplayData;
	
	pthread_create( &mAcquireData, NULL, get_samples, NULL );
	pthread_create( &mComputeFFT, NULL, compute_fft, NULL );
	pthread_create( &mDisplayData, NULL, display_data, NULL );


	while(1)
	{
		// Keep Parent Alive
		sleep(5);
	}
	
	return EXIT_SUCCESS;
}

// gets voltage reading from ACD chip
// adds value to vector in buffer
void* get_samples( void *aPtr )
{
	// create pointer to first buffer in the ring
	Buffer* mBufferPtr 		= mRingBuffer;
	int		mSamplesCounter	= 0;

	// open memory mappings to ACD registers
	static unsigned char  *control, *option, *complete;
	static short *reader;

	int devmem = open("/dev/mem", O_RDWR|O_SYNC);
	assert(devmem != -1);

	control = ( unsigned char *)mmap(0, getpagesize(), PROT_READ|PROT_WRITE, MAP_SHARED, devmem, CONTROL_REGISTER);
	assert(&control != MAP_FAILED);

	reader = (short *)control;

	option = (unsigned char*)mmap(0, getpagesize(), PROT_READ|PROT_WRITE, MAP_SHARED, devmem, OPTION_REGISTER);
	assert(&option != MAP_FAILED);

	complete = (unsigned char*)mmap(0, getpagesize(), PROT_READ|PROT_WRITE, MAP_SHARED, devmem, COMPLETE_REGISTER);
	assert(&complete != MAP_FAILED);

	RT_TASK* get_samples_task = rt_task_init( nam2num( mTask_1_name ), PRIORITY, STACK_SIZE, MAX_MSG_SIZE );
	rt_task_make_periodic( get_samples_task, rt_get_time(), mPeriod );
	rt_set_sched_policy( get_samples_task, RT_SCHED_RR, TIMESLICE );

	while( 1 )
	{
		// initialized count to number of available buffers since
		// is incremented every time the display module finishes
		rt_sem_wait( sem_3_1 );
		DEBUG_PRINT(("Get Samples: Grabbed Semaphore\n"));

		// change buffer status to LOCKED
		mBufferPtr->status	= LOCKED;
		rt_task_make_periodic( get_samples_task, rt_get_time(), mPeriod );

		// if we haven't collected 512 samples continue sampling
		while( mSamplesCounter < mNumOfSamples )
		{
			// collect sample
				// initialize ACD chip
				// Bit (0-2)	= Channel Select ( Binary value 0 to 7 )
				// Bit (3)		= 0 Unipolar/ 1 Bipolar
				// Bit (4)		= Range Select 0 5V, 1 10V
				// Bit (5-7)	= Set to ( 010 )

				// since I wanted setup for channel 0 with 5v Bipolar:
				// 765 4   3 210
				// 010 0   1 000
				//   4       8
				*control = 0x48;

				// wait for conversion to occur
				if(rt_task_wait_period() < 0){
					printf("RTAI: Period trouble\n");
				}

				// a period over 12us should be sufficient for conversion to occur
				// but we check if it is complete just to be sure
				while( (*complete & 0x80) != 0x00 )
				{
					// No Body
				}

				// store result as integer into vector
				// **************************************************
				// ** suppose to be  = ( 5.0/2048.0 )*(double)*(reader);
				// ** but I was never able to test this correctly
				// ** because I don't know how to dynamically
				// ** adjust the output of the FFT to optimally
				// ** display on the graph
				// ** in general it seams the higher the data fed to the fft
				// ** the higher the output but the graph could be
				// ** optimized if you find the right value to convert
				// ** i chose 60 from trial and error.

				// ** other notes about general conversion techniques
				// ** remember i'm using bipolar so 12th bit is sign
				// ** so eq is (range)/(2^11) = (5.0/2048.0)
				// ** if unipolar eq = (range)/(2^12)
				// **************************************************
//				mBufferPtr->vReal[ mSamplesCounter ] = ( 5.0/2048.0 )*((double)(*reader));
				mBufferPtr->vReal[ mSamplesCounter ] = (60)*( 5.0/2048.0 )*((double)(*reader));

				mSamplesCounter++;
		}

		DEBUG_PRINT(("Generate Samples on Buffer: "));
		DEBUG_PRINT(("BUFFER: %2d, %2d, %2d\n", mBufferPtr->id, mBufferPtr->status, mBufferPtr->stage ));

//		generateSignal( mBufferPtr->vReal, 5, 1000, 8000, mNumOfSamples );

		DEBUG_PRINT(("Get Samples: Collected %d Samples\n", mNumOfSamples ));
//		print_values( mBufferPtr->vReal, mNumOfSamples ); // for testing./

		// reset counter
		mSamplesCounter = 0;

		// when 512 samples update status to OPEN // update stage to RAW_DATA
		mBufferPtr->stage	= RAW_DATA;
		mBufferPtr->status	= OPEN;

		// post semaphore for FFT computing thread
		// that through should sequentially continue through the available buffers
		// until it encounters the buffer we just filled with data
		DEBUG_PRINT(("Get Samples: Signaling FFT Semaphore, Moving Buffer\n\n"));
		rt_sem_signal( sem_1_2 );\

		// sleep for small amount so lower priority task can grab semaphore immediately
		rt_sleep( nano2count( 1000 ) );

		// move buffer pointer to next buffer that has OPEN status and OVERWRITE stage
		getNextBuffer( &mBufferPtr, mRingBuffer, OPEN, OVERWRITE, NUM_OF_BUFFERS );
	}
}

void* compute_fft( void *aPtr )
{
	// create pointer to first buffer in the ring
	Buffer* mBufferPtr 		= mRingBuffer;

	RT_TASK* compute_fft_task = rt_task_init( nam2num( mTask_2_name ), PRIORITY, STACK_SIZE, MAX_MSG_SIZE );
	rt_task_make_periodic( compute_fft_task, rt_get_time(), mPeriod );
	rt_set_sched_policy( compute_fft_task, RT_SCHED_RR, TIMESLICE );

	while( 1 )
	{
		// wait for the get_samples() thread to post semaphore
		rt_sem_wait( sem_1_2 );
		DEBUG_PRINT(("Compute FFT: Grabbed Semaphore\n"));

		// remember this is initialized to the first buffer ( mRingBuffer )
		// lock current buffer;
		mBufferPtr->status	= LOCKED;

		DEBUG_PRINT(("Compute FFT on: "));
		DEBUG_PRINT(("BUFFER: %2d, %2d, %2d\n", mBufferPtr->id, mBufferPtr->status, mBufferPtr->stage ));

		//compute FFT on data in place
		fft_windowing(mBufferPtr->vReal, mNumOfSamples, FFT_WIN_TYP_HAMMING, FFT_FORWARD);		// Weigh data
		fft_compute(mBufferPtr->vReal, mBufferPtr->vImag, mNumOfSamples, FFT_FORWARD);		// Compute FFT
		fft_complexToMagnitude(mBufferPtr->vReal, mBufferPtr->vImag, mNumOfSamples);			// Compute magnitudes

		mBufferPtr->peak = fft_majorPeak( mBufferPtr->vReal, mNumOfSamples, mSampleFrequency );

		//update stage to FFT_DATA
		mBufferPtr->stage	= FFT_DATA;
		mBufferPtr->status	= OPEN;

		// post semaphore for display_data() thread
		// that through should sequentially continue through the available buffers
		// until it encounters the buffer we just processed
		DEBUG_PRINT(("Compute FFT: Signaling Display Semaphore\n\n"));
		rt_sem_signal( sem_2_3 );

		// sleep for small amount so lower priority task can grab semaphore immediately
		rt_sleep( nano2count( 1000 ) );

		// move buffer pointer to next buffer that has OPEN status and RAW_DATA stage
		getNextBuffer( &mBufferPtr, mRingBuffer, OPEN, RAW_DATA, NUM_OF_BUFFERS );
	}
}

void* display_data( void *aPtr )
{
	Buffer* mBufferPtr 		= mRingBuffer;

	RT_TASK* display_data_task = rt_task_init( nam2num( mTask_3_name ), PRIORITY, STACK_SIZE, MAX_MSG_SIZE );
	rt_task_make_periodic( display_data_task, rt_get_time(), mPeriod );
	rt_set_sched_policy( display_data_task, RT_SCHED_RR, TIMESLICE );

	while( 1 )
	{
		// wait for the get_samples() thread to post semaphore
		rt_sem_wait( sem_2_3 );
		DEBUG_PRINT(("Display Data: Grabbed Semaphore\n"));

		// change buffer status to LOCKED
		mBufferPtr->status = LOCKED;

		DEBUG_PRINT(("Display Data on: "));
		DEBUG_PRINT(("BUFFER: %2d, %2d, %2d\n", mBufferPtr->id, mBufferPtr->status, mBufferPtr->stage ));

#ifdef DISPLAY_DELTA
		currentTime = rt_get_time_ns();
		printf("%20.2f\n", (float)( currentTime - previousTime ) );
		previousTime = currentTime;
#endif

		printf("Peak: %8.2f (Hz)\n", mBufferPtr->peak );

//		print_values( mBufferPtr->vReal, mNumOfSamples );

		print_graph( mBufferPtr->vReal, mNumOfSamples/2.0, Y_MAX );

		printf("%c[23A", ASCII_ESC );

		// reset vImag to 0 to reset for next processing
		memset( mBufferPtr->vImag, 0, sizeof( double) * mNumOfSamples );

		// update buffer stage and status
		mBufferPtr->stage	= OVERWRITE;
		mBufferPtr->status	= OPEN;

		// post semaphore for get_samples() thread
		DEBUG_PRINT(("Display Data: Signaling 'Get Samples' Semaphore\n\n"));
		rt_sem_signal( sem_3_1 );

		rt_sleep( nano2count( 1000 ) );

		// move buffer pointer to next buffer that has OPEN status and FFT_DATA stage
		getNextBuffer( &mBufferPtr, mRingBuffer, OPEN, FFT_DATA, NUM_OF_BUFFERS );
	}
}

void getNextBuffer(
	Buffer 			**aCurrentBufferPtr,
	Buffer 			*aRingBuffer,
	Buffer_Status 	aStatusCondition,
	Buffer_Stage 	aStageCondition,
	int 			aTotalNumOfBuffers )
{
	// move buffer pointer to next buffer
//	do{
		if( (*aCurrentBufferPtr)->id == aTotalNumOfBuffers )
		{
			// set pointer back to first to simulate ring
			(*aCurrentBufferPtr) = aRingBuffer;
		}
		else
		{
			(*aCurrentBufferPtr)++;
		}

//	} while ( (*aCurrentBufferPtr)->status != aStatusCondition && (*aCurrentBufferPtr)->stage != aStageCondition );
}


void fft_windowing(double *vData, uint16_t samples, uint8_t windowType, uint8_t dir)
{
	// Weighing factors are computed once before multiple use of FFT
	// The weighing function is symmetric; half the weighs are recorded
	double samplesMinusOne = ( (double)(samples) - 1.0 );
	uint16_t i;
	for (i = 0; i < (samples >> 1); i++) {
		double indexMinusOne = (double)(i);
		double ratio = (indexMinusOne / samplesMinusOne);
		double weighingFactor = 1.0;
		// compute and record weighting factor
		switch (windowType) {
		case FFT_WIN_TYP_RECTANGLE: 	// rectangle (box car)
			weighingFactor = 1.0;
			break;
		case FFT_WIN_TYP_HAMMING: 		// hamming
			weighingFactor = 0.54 - (0.46 * cos(TWOPI * ratio));
			break;
		case FFT_WIN_TYP_HANN: 			// hann
			weighingFactor = 0.54 * (1.0 - cos(TWOPI * ratio));
			break;
		case FFT_WIN_TYP_TRIANGLE: 		// triangle (Bartlett)
			weighingFactor = 1.0 - ((2.0 * abs((int)(indexMinusOne - (samplesMinusOne / 2.0)))) / samplesMinusOne);
			break;
		case FFT_WIN_TYP_BLACKMAN: 		// blackmann
			weighingFactor = 0.42323 - (0.49755 * (cos(TWOPI * ratio))) + (0.07922 * (cos(FOURPI * ratio)));
			break;
		case FFT_WIN_TYP_FLT_TOP: 		// flat top
			weighingFactor = 0.2810639 - (0.5208972 * cos(TWOPI * ratio)) + (0.1980399 * cos(FOURPI * ratio));
			break;
		case FFT_WIN_TYP_WELCH: 		// welch
			weighingFactor = 1.0 - pow( ((indexMinusOne - samplesMinusOne / 2.0) / (samplesMinusOne / 2.0)) ,2 );
			break;
		}
		if (dir == FFT_FORWARD) {
			vData[i] *= weighingFactor;
			vData[samples - (i + 1)] *= weighingFactor;
		}
		else {
			vData[i] /= weighingFactor;		
			vData[samples - (i + 1)] /= weighingFactor;
		}
	}
}

void fft_compute(double *vR, double *vI, uint16_t samples, uint8_t dir)
{
	// Computes in-place complex-to-complex FFT
	// Reverse bits
	uint16_t j = 0;
	uint16_t i;
	for ( i = 0; i < (samples - 1); i++)
	{
		if (i < j) {
			 swap(&vR[i], &vR[j]);
			 swap(&vI[i], &vI[j]);
		}
		uint16_t k = (samples >> 1);
		while (k <= j)
		{
			 j -= k;
			 k >>= 1;
		}
		j += k;
	}

	// Compute the FFT
	double c1 = -1.0;
	double c2 = 0.0;
	uint8_t l2 = 1;
	uint8_t l;
	for ( l = 0; l < exponent(samples); l++) {
		uint8_t l1 = l2;
		l2 <<= 1;
		double u1 = 1.0;
		double u2 = 0.0;
		for (j = 0; j < l1; j++) {
			uint16_t g;
			for ( g = j; g < samples; g += l2) {
				uint16_t i1 = g + l1;
				double t1 = u1 * vR[i1] - u2 * vI[i1];
				double t2 = u1 * vI[i1] + u2 * vR[i1];
				vR[i1] = vR[g] - t1;
				vI[i1] = vI[g] - t2;
				vR[g] += t1;
				vI[g] += t2;
			}
			double z = (u1 * c1) - (u2 * c2);
			u2 = (u1 * c2) + (u2 * c1);
			u1 = z;
		}
		c2 = sqrt((1.0 - c1) / 2.0);
		if (dir == FFT_FORWARD) c2 = -c2;
		c1 = sqrt((1.0 + c1) / 2.0);
	}

	// Scaling for forward transform
	if (dir == FFT_FORWARD) {
		uint16_t i;
		for ( i = 0; i < samples; i++) {
			 vR[i] /= samples;
			 vI[i] /= samples;
		}
	}
}

void fft_complexToMagnitude(double *vR, double *vI, uint16_t samples) {
// vM is half the size of vR and vI
	uint8_t i;
	for ( i = 0; i < samples; i++) {
		vR[i] = sqrt( pow( vR[i], 2 ) + pow( vI[i], 2 ) );
	}
}


double fft_majorPeak(double *vD, uint16_t samples, double samplingFrequency) {
	double maxY = 0;
	uint16_t IndexOfMaxY = 0;
	for (uint16_t i = 1; i < ((samples >> 1) - 1); i++) {
		if ((vD[i-1] < vD[i]) && (vD[i] > vD[i+1])) {
			if (vD[i] > maxY) {
				maxY = vD[i];
				IndexOfMaxY = i;
			}
		}
	}
	double delta = 0.5 * ((vD[IndexOfMaxY-1] - vD[IndexOfMaxY+1]) / (vD[IndexOfMaxY-1] - (2.0 * vD[IndexOfMaxY]) + vD[IndexOfMaxY+1]));
	double interpolatedX = ((IndexOfMaxY + delta)  * samplingFrequency) / (samples-1);
	// retuned value: interpolated frequency peak apex
	return(interpolatedX);
}


void swap(double *x, double *y) {
	double temp = *x;
	*x = *y;
	*y = temp;
}

uint8_t exponent(uint16_t value) {
	// computes the exponent of a powered 2 value
	uint8_t result = 0;
	while (((value >> result) & 1) != 1) {
		result++;
	}
	return(result);
}

void generateSignal(
		double* aVector,
		int aAmplitude,
		int aSignalFrequency,
		int aSampleFrequency,
		int aSamples
		)
{
	int i;
	for( i = 0; i < aSamples; i++ )
	{
		aVector[i] = ( aAmplitude * ( sin( 2.0 * 3.14159 * ( (double)aSignalFrequency / (double)aSampleFrequency ) * i ) / 2.0 + 0.5 ) );
	}
	return;
}

void print_values( double *aVector, int size )
{
	int i;
	for( i = 0; i < size; i++ )
	{
		printf("v[%2d]: %5.4f\n", i, aVector[i] );
	}
	return;
}

void print_graph( double *aVector, int width, int height )
{
	int mYmax	= height;
	int mXmax	= width;

	// setup pointers for sprintf
	// sprintf is much more  efficient than printf for this case
	// or printing many small characters and concatenating together
	char graph[ (mYmax + 2) * ( mXmax + 4 ) + 1 ];
	char* writeLocation = graph;

	// for loop indecies
	int i;
	int j;

	writeLocation += 1;
	// compute string to print to screen
	// count down from yMax to 0 to preserve + up coordinate system
	for( i = mYmax; i > 0; i-- )
	{
		if( i%5 == 0 )
		{
			writeLocation += sprintf( ( writeLocation -1 ), "%2d|", i );
		}
		else
		{
			writeLocation += sprintf( ( writeLocation -1 ), "  |" );
		}
		// count up from 0 to preserve + right coordinate system
		for( j = 0; j < mXmax; j++ )
		{
			// if the current value for this column is higher than the position
			// fill this column in with # sign to create pseudo bar graph
			if( i <= aVector[j] )
			{
				writeLocation += sprintf( ( writeLocation -1 ), "#" );
			}
			else
			{
				writeLocation += sprintf( ( writeLocation -1 ), " " );
			}
		}
		// attach newline character
		writeLocation += sprintf( ( writeLocation -1 ), "\n" );
	}

	writeLocation += sprintf( ( writeLocation -1 ), "  +" );
	for( j = 0; j < mXmax; j++ )
	{
		writeLocation += sprintf( ( writeLocation -1 ), "-" );
	}
	writeLocation += sprintf( ( writeLocation -1 ), "\n   " );

	for( j = 0; j < mXmax; j += 5 )
	{
		writeLocation += sprintf( ( writeLocation -1 ), "%-5d", j );
	}
	writeLocation += sprintf( ( writeLocation -1 ), "\n" );

	printf("%s", graph );
}