/***************************************************************/
/*Copyright (c) by Continocean Tech Inc. 	       	       */
/*12 Mountain Ave., Montville, New Jersey, USA 07045           */
/*Tel: 201-257-1912 Fax: 201-257-9634			       */
/*Email: continocean@webexpert.net			       */
/*All rights reserved. No part of this work covered by the     */
/*copyright hereon may be reproduced or used in any form or by */
/*any means----graphic, electronic, or mechanical, including   */
/*photocopying, recording, taping, or information storage and  */
/*retrieval systems----without written permission of 	       */
/*Continocean Tech inc.					       */
/***************************************************************/

                     ENGSIM for MS VC++

                  Engineering Simulation Library
                    for Microsoft Visual C++

                   Usage of ENGSIM Instructions


IMPORTANT MESSAGE:
    Please use the Microsoft Development Studio's File-OpenWorkspace to
open \engsim\sim project space. Then use File-Open to open file main.cpp
for editing. There are several examples under \engsim\examples which are 
included in main.cpp. You can edit main.cpp and the examples for exercises.


PART I. COMPLEX, MATRIX, COMPLEX MATRIX, INTEGER MATRIX AND STRING MATRIX

Complex                declare a complex number
e.g.                   Complex a;

matrix	  		declare a matrix
			matrix a(2,2), h(10,20), p(200);
			//a(1,1)...... a(2,2)
			//p(1) ...... p(200)

matrixn			declare a matrix with negative indexing
			matrixn a(2,2); //a(-2,-2) ...... a(2,2)
			matrixn p(100); //p(-100) ......p(100)

cmatrix	  	declare a complex matrix
			cmatrix a(2,2),b(2,2),p(100);

cmatrixn	declare a complex matrix with negative indexing
			cmatrixn a(2,2),b(2,2),p(100);
			//a(-2,-2)......a(2,2), p(-100)......p(100)

imatrix	   	declare an integer matrix
			imatrix a(2,2), b(2,2), p(100);
			//a(1,1) ...... a(2,2)
			//p(1) ...... p(100)

imatrixn	declare an integer matrix with negative indexing
			imatrixn a(2,2),b(2,2),p(100);
			//a(-2,-2) ... a(2,2)
			//p(-100) ... p(100)

smatrix	  	declare a string matrix
			smatrix a(2,2),p(100);
			//a(1,1) ... a(2,2)
			//p(1) ... p(100)

J           pure unit imaginary for Complex
            a= 1.0+2.0*J;

identity	identity matrix (diagonal elements are all 1's, others all 0's)
		    (for matrix only)
		    a=identity(10);	   //10x10 identity matrix

unitary		unitary matrix (all elements are all 1's)
		    (for matrix only)
		    a=unitary(10);	   //10x10 unitary matrix

identityn	identity matrix with negative indexing
		    (for matrix with negative indexing only)
		    a=identityn(10);	   //21x21 identity matrixn

unitaryn	unitary matrix with negative indexing
		    (for matrix with negative indexing only)
		    a=unitaryn(10);	   //21x21 unitary matrixn

real()      real of a complex number
		    double da;
            da=real(a);
			real of Complex matrix
			matrix mc(2,2);
			mc=real(a);
			matrixn mc(2,2);
			mc=real(a);

imag()      imaginary of a complex number
            double da;
            da=imag(a);
		    imaginary of Complex matrix
			matrix mc(2,2);
			mc=imag(a);
			matrixn mc(2,2);
			mc=imag(a);

=           equal operator for complex, matrix, matrixn, 
		    cmatrix, cmatrixn, imatrix, imatrixn
            b= 1.0+2.0*J;
		    b=a;

+           plus operator for complex, matrix, matrixn, 
		    cmatrix, cmatrixn, imatrix, imatrixn
            c= a+b;

-           minus operator for complex, matrix, matrixn, 
		    cmatrix, cmatrixn, imatrix, imatrixn
            c= a-b;
            c=-b;

*           multiplication operator for complex, matrix, matrixn, 
		    cmatrix, cmatrixn, imatrix, imatrixn
            c= a*b;

/		    division operator for complex, matrix, matrixn, 
		    cmatrix, cmatrixn
		    c=a/b;

~           conjugate operator for complex, cmatrix, cmatrixn
            c= (~a);

==		    logic equal operator for complex, matrix, matrixn, 
		    cmatrix, cmatrixn, imatrix, imatrixn, smatrix
		    if(a==b){...};

!=		    logic not equal operator for complex, matrix, matrixn, 
		    cmatrix, cmatrixn, imatrix, imatrixn, smatrix
		    while(a!=b){...};

exp()		exponential for complex
		    exp(J*2.0*PI);

abs()		absolute for complex
		    sqrt(sum a(i,j)*a(i,j)) for matrix, matrixn
		    sqrt(sum a(i,j)*(~a(i,j))) for cmatrix, cmatrixn
			double dc;
			dc=abs(a);

normsq()		real*real+imag*imag for complex	
		    sum a(i,j)*a(i,j) for matrix, matrixn, imatrix, imatrixn
		    sum a(i,j)*(~a(i,j)) for cmatrix, cmatrin
			double dc;
			dc=norm(a);

phase()		phase of a complex number  between (0,2*PI)
		    phase(c);

arg()		same as phase()

pow( , )	power for complex number
		    Complex c;
		    pow(c,2.0);
		    pow(c,3);

printw()	print to screen for complex, matrix, matrixn, 
		    cmatrix, cmatrixn, imatrix, imatrixn, smatrix,
		    double, integer, character
		    printw(a);
		    printw(2.0+3.0*J);
		    printw(1.0);
		    printw(1);
		    printw('c');
		    printw("test=%g, dog=%d\n",0.5, 100);

fprintf(,)	print to a disk file for complex, matrix, matrixn, 
		    cmatrix, cmatrixn, imatrix, imatrixn, smatrix,
		    double, integer, character
		    fprintf(fp, c);

fscanf( , )	read from a disk file for complex, matrix, matrixn, 
		    cmatrix, cmatrixn, imatrix, imatrixn, smatrix,
		    double, integer, character
		    fscanf(fp, a);

( , )		an element of a matrix, matrixn, cmatrix, cmatrixn
			imatrix, imatrixn, smatrix
			a(2,-5)=10.1+2.1*J;
			p(-5)=3.0+9.9*J;  

init(, ...)	initialize a matrix, matrixn, cmatrix, cmatrixn
			imatrix, imatrixn, smatrix
			matrix a(2,2);
			init(a, 1.0,2.0,
				4.0,5.0
				);
			matrixn a(2,2);
			a=init(a,1.0,2.0,3.0,4.0,5.0,
				 2.0,2.0,3.0,4.0,5.0,
				 3.0,2.0,3.0,4.0,5.0,
				 4.0,2.0,3.0,4.0,5.0,
				 5.0,2.0,3.0,4.0,5.0
				 );
			cmatrix a(2,2);
			a=init(a,
				1.0+2.0*J, 2.0+3.0*J,
				1.1+2.1*J, 2.1+3.1*J
				);
			cmatrixn a(2,2);
			a=init(a,
				1.0+2.0*J, 2.0+3.0*J, 1.0+2.0*J, 2.0+3.0*J, 1.0+2.0*J,
				1.1+2.1*J, 2.1+3.1*J, 1.0+2.0*J, 2.0+3.0*J, 1.0+2.0*J,
				1.0+2.0*J, 2.0+3.0*J, 1.0+2.0*J, 2.0+3.0*J, 1.0+2.0*J,
				1.1+2.1*J, 2.1+3.1*J, 1.0+2.0*J, 2.0+3.0*J, 1.0+2.0*J,
				1.0+2.0*J, 2.0+3.0*J, 1.0+2.0*J, 2.0+3.0*J, 1.0+2.0*J
				);
			imatrix a(2,2);
			init(a,
		   		1,2,
				4,5
				);
			imatrixn a(2,2);
			init(a,
		   		1,2,3,2,3,
				4,5,6,2,3,
				7,8,9,2,3,
				4,5,6,2,3,
				7,8,9,2,3
				);
			smatrix a(2,2);
			a=init(a,
				'H','i',
				'O','K');


joinr(a, b);
	right join matrices a,b as (a b) for matrix, cmatrix, imatrix
	
joind(a, b);
	down join matrices a,b as 
	   a(1,1)......
	   ............ 
	(  b(1,1)......  )
	   ............
	for matrix, cmatrix, imatrix

shiftl( )	shift left by 1 for matrix, matrixn, 
		    cmatrix, cmatrixn, imatrix, imatrixn

shiftr( )	shift right by 1 for matrix, matrixn, 
		    cmatrix, cmatrixn, imatrix, imatrixn

shiftu()	shift up by 1 for matrix, matrixn, 
		    cmatrix, cmatrixn, imatrix, imatrixn

shiftd()	shift down by 1 for matrix, matrixn, 
		    cmatrix, cmatrixn, imatrix, imatrixn

sub(a, m1, n1, m2, n2)	
			Take submatrix beginning at a(m1,n1) and ending at a(m2,n2)
			for matrix, matrixn, cmatrix, imatrix		
			(matrixn becomes a matrix)
			sub(a,m1,n1,m2,n2);

TT()		transpose for matrix, matrixn, 
		    cmatrix, cmatrixn, imatrix, imatrixn, smatrix
			c=TT(a);

HH			hermitian for cmatrix, cmatrixn
			c=HH(a);


Max()		find maximum of all elements for matrix
			double a_max=Max(a);

Min()		find minimum of all elements for matrix
			double a_min=Min(a);

Max( , , ,)	find maximum of all elements for matrix
			Max(a, y_max, i_max, j_max); //y_max=a(i_max,j_max)

Min( , , ,)	find minimum of all elements for matrix
			Min(a, y_min, i_min, j_min); //y_min=a(i_min,j_min) 
 
corrm(a, b, n)	
			n-by-n correlation matrix of a and b for cmatrix, cmatrixn
			c=corrm(a,b,2);

det()		determinant of a matrix, matrixn, cmatrix, cmatrixn
			imatrix, imatrixn
			double dc;
			dc=det(a);

|			or for imatrix, imatrixn
			c=a|b;

&			and for imatrix, imatrixn
			c=a&b;

^			X-or for imatrix, imatrixn
			c=a^b;

int2string(ia, sa, len)	
			convert integer ia to string sa(1) ... sa(len)

float2string(fa, sa, len) 
			convert float number fa to string sa(1) ... sa(len)


PART II. GRAPHICS

plot(y)		plot array y 
			matrix y(100);
			plot(y);

plot(x, y)	plot array, x=X-axis, y=Y-axis
			matrix x(100), y(110);
			plot(x, y);

plot(x, y1, y2)		
			plot array, x=X-axis, y1,y2,...=Y-axis
			matrix x(100), y1(110), y2(100);
			plot(x, y1, y2);

plot(x, y1, y2, y3)
plot(x, y1, y2, y3, y4);
plot(x, y1, y2, y3, y4, y5);

plot(y, "legend");
plot(x, y, "legend");
plot(x, y1, "legend1", y2, "legend2");
plot(x, y1, "legend1", y2, "legend2", y3, "legend3");
plot(x, y1, "legend1", y2, "legend2", y3, "legend3", y4, "legend4");
plot(x, y1, "legend1", y2, "legend2", y3, "legend3", y4, "legend4", y5, "legend5");

set_picture_size(x);
	scale picture
	x=1.0 default size
	0.0<x<1.0  shrink by x
	1.0<x<2.0  expand by x
	set_picture_size(1.0); //default
	set_picture_size(1.5);
	set_picture_size(0.5);

	
set_picture_position(x, y);
	relocate picture position
	-10000<x,y<10000
	set_picture_position(0, 0);	//default
	set_picture_position(2000, 1000); //move right by 2000, up by 1000
	set_picture_position(-2000, -1000); //move left by 2000, down by 1000

set_axis_title("x_title", "y_title");

set_text_color("color");
	color=black (default), blue, green, cyan, red, magenta, yellow, white 
	      darkblue, darkgreen, darkcyan, darkred, darkmagenta, darkyellow
	       darkgrey, lightgrey

set_text_position(x, y);	
	x=row, y=column
	x,y>0
	set_text_position(100, 20);
	

set_text_size(x);
	set text font size
	set_text_size(10); //default
	set_text_size(15);




PART III. COMMUNICATIONS ENGINEERING AND SIGNAL PROCESSING SIMULATION



convol(x, y)
	convolution of arrays x and y for matrix, matrixn, cmatrix and cmatrixn
	matrix x(100),y(100),z(200);
	z=convol(x,y);
	matrixn xn(100),yn(100),zn(200);
	zn=convol(xn,yn);
	cmatrix cx(100),cy(100),cz(200);
	cz=convol(cx,cy);
	cmatrixn cxn(100),cyn(100),czn(200);
	czn=convol(cxn,cyn);

corr(x, y)
	correlation of arrays x and y for matrix (array) and cmatrix (array)
	matrix x(100), y(100);
	matrixn z(100);	//z(-100) ... z(100)
	z=corr(x,y);
	cmatrix x(100), y(100);
	cmatrixn z(100);	//z(-100) ... z(100)
	z=corr(x,y);

corrnorm(x, y)
	normalized correlation of arrays x and y for matrix (array) and cmatrix (array)
	matrix x(100), y(100);
	matrixn z(100);	//z(-100) ... z(100)
	z=corrnorm(x,y);
	cmatrix x(100), y(100);
	cmatrixn z(100);	//z(-100) ... z(100)
	z=corrnorm(x,y);


crandom();
	generate uniform-distributed integer in (0, 2147483647)
	int a;
	a=crandom();

fft(x, sp_freq, sp_amp, sp_phase, n, Ts);
	compute spectrum of an array x using FFT
	n=pow(2,r)=number of points in FFT, 
	dimension of sp_freq (frequency), sp_amp (amplitude) and sp_phase (phase) is n/2
	Ts=sampling-period
	int n=4098;	//2^12
	matrix x(n), freq(n/2), amp(n/2), phase(n/2);
	float Ts=0.01;
	fft(x, freq, amp, phase, n, Ts);
	plot(freq, amp, "Amplitude", phase, "Phase");

integral(func, a, b, step)	    
        integrate function func from a to b with step-size step
	double function(double x){
	       double y;
	       ......
	       return y;
	}
	double a=-1.0,b=2.0,step=0.01;
	integral(func,a,b,step);


lpf(fmT, n, h);
	generate low pass filter coefficient h(1) ...... h(n) with cut-off frequency fmT
	int n=20;
	float fmT=1000;
	matrix h(n);
	lpf(fmT, n, h);

norm(m, v)
	generate Gaussian (Normal) distributed random variable with mean m and variance v*v
	float a;
	a=norm(1.0, 2.0); //Gaussian distributed random variable with mean 1.0 and variance 4.0
 
phasecontinue(p)
	make the phase array p to be continuous by subtracting or adding 2PI.
	matrix p(100);
	phasecontinue(p);
	plot(p);	//the plot  will be continuous

phasefollow(p1, p2)
	make phase array p1 follow phase array p2 by subtracting or adding 2PI.
	matrix p1(100),p2(100),x(100);
	phasefollow(p1,p2);
	......
	plot(x,p1,p2);	//

qpiqpskmod(a2, a1)
	PI/4-QPSK modulation. 
	binary sequence (a2 a1)to QPSK mapping: 00=PI/4, 01=3*PI/4, 11=5PI/4, 10=7PI/4.
	cmatrix x(100);
	imatrix a(200);
	a(2)=0;
	a(1)=1;
	x(1)=qpiqpskmod(a(2), a(1));

qpiqpskdem(x, a2, a1);
	PI/4-QPSK demodulation. 
	binary sequence (a2 a1)to QPSK mapping: 00=PI/4, 01=3*PI/4, 11=5PI/4, 10=7PI/4.
	cmatrix x(100);
	imatrix a(200);
	x(1)=exp(J*PI/4.0);
	qpiqpskdem(x(1), a(2), a(1));

qpidqpskmod(a2, a1, x1);
	PI/4-DQPSK modulation. 
	binary sequence (a2 a1)to QPSK mapping: 00=PI/4, 01=3*PI/4, 11=5PI/4, 10=7PI/4.
	cmatrixn x(100);
	imatrixn a(200);
	a(0)=0;
	a(-1)=1;
	x(-1)=exp(J*PI/4.0);
	x(0)=qpidqpskmod(a(0), a(-1),x(-1));

qpidqpskdem(y0, y1, a2, a1);
	PI/4-DQPSK demodulation. 
	binary sequence (a2 a1)to QPSK mapping: 00=PI/4, 01=3*PI/4, 11=5PI/4, 10=7PI/4.
	cmatrix y(100);
	imatrix a(200);
	y(-1)=exp(J*PI/4.0);
	y(0)=exp(J*3.0*PI/4.0);
	qpidqpskdem(y(-1), y(0), a(0), a(-1));


ray(m, v)
	generate Rayleigh (m=0) or Rician (m!=0) distributed variable
	float a;
	a=ray(1.0, 2.0);

ray_fade(fm, T, k);
	generate Jakes Rayleigh fading channel coefficient
	Doppler frequency is fm Hz, T is symbol period, 
	k is any integer number greater than 0 
	(setting k=1 and k=1000 respectively gives two-ray channel or diversity channel)
	cmatrixn c1(100), c2(100);
	float fm=75.0, T=41.0e-6;
	c1(0)=ray_fade(fm, T, 1);
	c2(0)=ray_fade(fm, T, 1000); //c1 and c2 are two independent rays.

sinc(x);
	sinc function sin(x)/x
	double y;
	y=sinc(1.5*PI);

spectrum(func, sp_freq, sp_amp, sp_phase, n, T);
	compute spectrum of function func
	T=function period, from 0 to T
	n=pow(2,r)=number of points in FFT, 
	dimension of sp_freq (frequency), sp_amp (amplitude) and sp_phase (phase) is n/2
	sampling-period=T/n
	float function(float x)
	{...... }
	int n=4098;	//2^12
	matrix freq(n/2), amp(n/2), phase(n/2);
	float T=100.0;
	spectrum(function, freq, amp, phase, n, T);
	plot(freq, amp, "Amplitude", phase, "Phase");

tail(func, a)
	integrate the tail of the function func from a to +infinity
	double function(double x){
	       double y;
	       ......
	       return y;
	}
	double a=1.0;
	tail(func,a);