/*****************************************************************************
	CHUA.C - A Chua chaos osscillator simulator
	Written by Bertie Genade 1995
*****************************************************************************/

#include<conio.h>
#include<dos.h>
#include<errno.h>
#include<stdio.h>
#include<stdlib.h>
#include<string.h>
#include<math.h>
#include<graphics.h>


double C1,C2,L,G1;			/* Two terminal devices parameters */
double Bp1,Bp2,m0,m1,m2;	/* Non-linear resitor parameters */
double Vc1,Vc2,Il;			/* State variables */
double dVc1,dVc2,dIl;		/* Derivatives of state variables */
double Ls,C1s,C2s;			/* Time scaled values of L,C1 and C2 */
double tscale;				/* Used for time scaling */
double frequency;			/* Frequency of fundamental harmonic */
double tstep;				/* Time increment of each differiantial step */
double size;				/* Size of graphics window, use 3 for X-Y plot */
int samplerate;				/* Samplerate at which .WAV file is compiled */
long int datalength;		/* Total length of sound data written to .WAV file */
long int datidx;			/* Current write pointer position in sound data */
int waitc;					/* Number of calculations per sample */

FILE *fp;

char buffer[16000];			/* Memorybuffer for making .WAV file */
char str[100];
char fname[10];

enum PLOTMODES {PLOTXY,PLOTVC1,PLOTVC2,PLOTIL,PLOTALL};
enum WAVMODES {WAVVC1,WAVVC2,WAVIL};
int plotmode, wavmode;
char *wmodes[]={"Vc1","Vc2","Il"};

struct waveformat{
	char		rID[4];
	long int	rLen;
	char		wID[4];
	char		fID[4];
	long int	fLen;
	int			WformatTag;
	int			nChannels;
	long		nSamplesPerSec;
	long		nAvgBytesPerSec;
	int			nBlockAlign;
	int			FormatSpecific;
	char		dID[4];
	long int	dLen;
	} wavef;

void savewaveheader(void){
	strcpy(wavef.rID,"RIFF");
	strcpy(wavef.wID,"WAVE");
	strcpy(wavef.fID,"fmt ");
	strcpy(wavef.dID,"data");
	wavef.rLen = datalength + 0x18;
	wavef.fLen = 12;
	wavef.WformatTag = 1;
	wavef.nChannels = 1;
	wavef.nSamplesPerSec = samplerate;
	wavef.nAvgBytesPerSec = samplerate;
	wavef.nBlockAlign = 1;
	wavef.FormatSpecific = 1;
	wavef.dLen = datalength;
	if((fp=fopen("CHUA.WAV","wb"))==NULL) exit(1);
	fwrite(wavef.rID,sizeof(struct waveformat),1,fp);
	}

void plot(double x, double y, int mode){
	int ax,ay;

	switch(mode){
	case PLOTXY:ax=320+x/size*240;
				ay=240-y/size*240;
				putpixel(ax,ay,15);
				break;
	case PLOTVC1:ax=(double)datidx/datalength*640;
				ay=240-Vc1/size*240.0;
                putpixel(ax,ay,15);
				break;
	case PLOTVC2:ax=(double)datidx/datalength*640;
				ay=240-Vc2/size*240.0;
                putpixel(ax,ay,15);
				break;
	case PLOTIL:ax=(double)datidx/datalength*640;
				ay=240-Il/size*240.0;
                putpixel(ax,ay,15);
				break;
	case PLOTALL:ax=(double)datidx/datalength*640;
				ay=80-Vc1/size*80.0;
                putpixel(ax,ay,15);
				ay=240-Vc2/size*80.0;
                putpixel(ax,ay,15);
				ay=400-Il/size*80.0;
				putpixel(ax,ay,15);
				break;
				}
	}

void fillbuffer(int length){
	double g2, wc;
	int c;
	for(c=0 ; c<length ; ++c){
		for(wc=0 ; wc<waitc ; ++wc){
			if(Vc1<Bp1) g2 = m0*(Vc1-Bp1) + m1*Bp1;
			if((Vc1>=Bp1)&&(Vc1<=Bp2)) g2 = m1*Vc1;
			if(Vc1>Bp2) g2 = m2*(Vc1-Bp2) + m1*Bp2;
			dVc1 = (G1 * (Vc2 - Vc1) - g2) / C1s;
			dVc2 = (G1 * (Vc1 - Vc2) + Il) / C2s;
			dIl = -Vc2 / Ls;
			Vc1 = Vc1 + dVc1 * tstep;
			Vc2 = Vc2 + dVc2 * tstep;
			Il = Il + dIl * tstep;
			}
		switch(wavmode){
		case WAVVC1:buffer[c] = (char)(Vc1 * 42.0);
					break;
		case WAVVC2:buffer[c] = (char)(Vc2 * 42.0);
					break;
		case WAVIL:	buffer[c] = (char)(Il * 42.0);
					break;
					}
		plot(Vc1,Vc2,plotmode);
		++datidx;
		}
	}


void fillblock(int length, double freq){
	int bc;

	tscale = 0.325/freq;
	tstep = 1.0/((double)waitc) * 1/16000.0;

	Vc1 = 0;						/* Initial conditions */
	Vc2 = .1;
	Il = 0;

	Ls = L * tscale;				/* Derive time scaled component values */
	C1s = C1 * tscale;
	C2s = C2 * tscale;

	for(bc=0 ; bc<=(length/16000) ; ++bc){
		fillbuffer(length%16000);
		fwrite(buffer,length%16000,1,fp);
		}
	}
	


void calculate(void){
	int gdriver=VGA, gmode=VGAHI;
	initgraph(&gdriver,&gmode,"");
	setcolor(14);
	switch(plotmode){
	case PLOTXY:line(0,240,640,240);
				line(320,0,320,4809);
				break;
	case PLOTVC1:
	case PLOTVC2:
	case PLOTIL:line(0,240,640,240);
				break;
	case PLOTALL:line(0,80,640,80);
				line(0,240,640,240);
				line(0,400,640,400);
				break;
				}
	setcolor(15);

	savewaveheader();
	datidx = 0;
	fillblock(15999,frequency);
	sprintf(str,"L=%g C1=%g C2=%g G=%g Bp1=%g Bp2=%g m0=%g m1=%g m2=%g F=%g",
			L,C1,C2,G1,Bp1,Bp2,m0,m1,m2,frequency);
	outtextxy(0,0,str);
	switch(plotmode){
	case PLOTXY:outtextxy(0,470,"Vc2 versus Vc1");
				break;
	case PLOTVC1:outtextxy(0,470,"Vc1");
				break;
	case PLOTVC2:outtextxy(0,470,"Vc2");
				break;
	case PLOTIL:outtextxy(0,470,"Il");
				break;
	case PLOTALL:outtextxy(0,470,"Vc1, Vc2 and Il");
				break;
				}
	sprintf(str,"Y-axis +%g to -%g per graph",size,size);
	outtextxy(130,470,str);
	if(plotmode==PLOTXY) outtextxy(470,470,"Aspect 1:1");
			else outtextxy(370,470,"X-axis totals 1 second.");
	getch();
	fclose(fp);
	closegraph();
	}

void enter_string(char *str){
	char rp, l=0;
	*str=0;
	while((rp=getch())!=13){
		if((rp>31)&&(rp<128)&&(l<30)){
			*str=rp;
			++str, ++l;
			putch(rp);
			}
		if((rp==8)&&(l>0)){
			--str, --l;
			cprintf("\b \b");
			}
		}
	*str=0;
	}

float enter_float(void){
	enter_string(str);
	return(atof(str));
	}

void menu(void){
	window(1,1,80,25);
	clrscr();
	textattr(4*16+15);
	cprintf("Chua circuit simulator v1.00        Written by Bertie Genade     Copyright 1995");
	textbackground(0);
	window(5,2,80,25);
	textcolor(14);
	cprintf("A  Change L\n\rB  Change C1\n\rC  Change C2\n\rD  Change G\n\r");
	cprintf("E  Change Bp1\n\rF  Change Bp2\n\rG  Change m0\n\r");
	cprintf("H  Change m1\n\rI  Change m2\n\r");
	cprintf("J  Change frequency\n\rK  Change graph size\n\r");
	textcolor(10);
	cprintf("L  Display Vc2 versus Vc1\n\rM  Display Vc1\n\r");
	cprintf("N  Display Vc2\n\rO  Display Il\n\r");
	cprintf("P  Display Vc1,Vc2,Il together\n\r");
	textcolor(13);
	cprintf("Q  Select .WAV file output\n\r");
	cprintf("R  Play .WAV file\n\rS  Save .WAV file\n\r");
	textcolor(12);
	cprintf("T  Exit to DOS\n\n\r");
	textcolor(15);
	cprintf("Press A - T to select item");
	window(40,2,80,25);
	textcolor(11);
	cprintf("L = %g\n\rC1 = %g\n\rC2 = %g\n\rG = %g\n\r",L,C1,C2,G1);
	cprintf("Bp1 = %g\n\rBp2 = %g\n\rm0 = %g\n\rm1 = %g\n\rm2 = %g\n\r",Bp1,Bp2,m0,m1,m2);
	cprintf("Frequency = %g\n\rGraph size = %g\n\r",frequency,size);
	cprintf("\n\n\n\n\nWave output = %s\n\n\r",wmodes[wavmode]);
	textcolor(15);
	}

void main(){
	int rp;

	samplerate=16000;
	datalength=16000;
	frequency=100;

	waitc = 20;
	size = 3;

	C1 = 1/9.0;						/* 1/9.0 = double scroll */
	L =  1/7.0;						/* 1/7.0 = double scroll */
	C2 = 1;			                /* 1 = double scroll */
	G1 = 0.7;						/* 0.7 = double scroll */
	Bp1 = -1;						/* Non-linear resistor parameters */
	Bp2 = 1;
	m0 = -0.5;
	m1 = -0.8;
	m2 = -0.5;

	while(1){
		menu();
		rp=getch();
		if(rp>96) rp-=32;
		switch(rp){
		case 65:cprintf("Enter new value for L : "); L=enter_float();
				break;
		case 66:cprintf("Enter new value for C1 : "); C1=enter_float();
				break;
		case 67:cprintf("Enter new value for C2 : "); C2=enter_float();
				break;
		case 68:cprintf("Enter new value for G : "); G1=enter_float();
				break;
		case 69:cprintf("Enter new value for Bp1 : "); Bp1=enter_float();
				break;
		case 70:cprintf("Enter new value for Bp2 : "); Bp2=enter_float();
				break;
		case 71:cprintf("Enter new value for m0 : "); m0=enter_float();
				break;
		case 72:cprintf("Enter new value for m1 : "); m1=enter_float();
				break;
		case 73:cprintf("Enter new value for m2 : "); m2=enter_float();
				break;
		case 74:cprintf("Enter new value for frequency : ");
				frequency=enter_float();
				break;
		case 75:cprintf("Enter new value for graph size : "); size=enter_float();
				break;
		case 76:plotmode=PLOTXY;
				calculate();
				break;
		case 77:plotmode=PLOTVC1;
    			calculate();
				break;
		case 78:plotmode=PLOTVC2;
    			calculate();
				break;
		case 79:plotmode=PLOTIL;
    			calculate();
				break;
		case 80:plotmode=PLOTALL;
    			calculate();
				break;
		case 81:cprintf("1  Vc1\n\r2  Vc2\n\r3  Il\n\rSelect 1,2 or 3.");
				switch(getch()){
				case '1':wavmode=WAVVC1;
						break;
				case '2':wavmode=WAVVC2;
						break;
				case '3':wavmode=WAVIL;
						break;
						}
				break;
		case 82:system("wave chua.wav");
				break;
		case 83:cprintf("Enter name of .WAV file :"); enter_string(fname);
				sprintf(str,"copy chua.wav %s.wav",fname);
				system(str);
				break;
		case 27:
		case 84:exit(0);
				}
		}
	}