//本程序是針對(duì)G(s)=(8*e^(-4s))/(1+2s)設(shè)計(jì)的大林算法，可以進(jìn)行無振鈴和有振鈴實(shí)驗(yàn)，以及有振鈴時(shí)消除振鈴的實(shí)驗(yàn)
//L=20,用勞斯判據(jù)判斷特征根是否在z平面的左半平面的簡單方法，特征方程為1-e^(-T/Tm_)*z^(-1)-(1-e^(-T/Tm_))*z^(-21)=0
//即對(duì)特征方程為z^(21)-e^(-T/Tm_)*z^(20)-(1-e^(-T/Tm_))=0直接列勞斯陣列表
//判斷有幾次符號(hào)改變，即可知道右半平面極點(diǎn)數(shù)，反過來也就知道左半平面極點(diǎn)數(shù)，然后用消除振鈴現(xiàn)象的方法，進(jìn)行消除振鈴現(xiàn)象的實(shí)驗(yàn)
#include	<stdio.h>
#include	<dos.h>
#include	<math.h>
#include	<bios.h>
#include	<conio.h>
#include	<graphics.h>
#define	BASE_LINE1	190
#define	BASE_LINE2	320
#define	BASE_LINE3	420

#define M 10

double	a,b,c,d,p,q,t,y1_t,y1_kT,y1_kT_[20],y1_t_[200];
double	y_kT,e_kT,e_kT_T,u_kT,u_kT_T,u_kT_21T,u_kT_[22],Tm,Tm_,T,l,K;

int	i,j,m,gdrive=VGA,gmode=VGAHI;
unsigned char   exit_f;
int	scan_keyboard_result;

void	ji_suan(void){
	p=1.0-exp(-0.5*t);	q=exp(-0.5*t);
}

void	work0(void)
{	i=1;	j=1;
	do{
		if(j==M)	t=T;
		else	t=T*((double)j/(double)M);
		ji_suan();
		y1_t=p*8.0*u_kT+q*y1_kT;
		for(m=198;m>=0;m--)	y1_t_[m+1]=y1_t_[m];
		y1_t_[0]=y1_t;
		if(j==M){
			y1_kT=y1_t;
			for(m=18;m>=0;m--)	y1_kT_[m+1]=y1_kT_[m];
			y1_kT_[0]=y1_kT;
			y_kT=y1_kT_[19];
			e_kT=1.0-y_kT;
			for(m=20;m>=0;m--)	u_kT_[m+1]=u_kT_[m];
			u_kT_21T=u_kT_[21];
			u_kT_T=u_kT_[1];
			u_kT=a*e_kT+b*e_kT_T+c*u_kT_T+d*u_kT_21T;
			u_kT_[0]=u_kT;
			e_kT_T=e_kT;
		}
		setcolor(GREEN);	moveto(i,250-(int)(150.0*y1_t_[199]));lineto(i,250-(int)(150.0*y1_t_[199]+1));	//對(duì)象連續(xù)輸出響應(yīng)
		setcolor(BLUE);		moveto(i,250-(int)(150.0*e_kT));lineto(i,250-(int)(150.0*e_kT+1.0));	//離散誤差輸出響應(yīng)
		setcolor(YELLOW);	moveto(i,400-(int)(150.0*u_kT));lineto(i,400-(int)(150.0*u_kT+1.0));
		i=i+1;	j+=1;
		if(j==M+1)	j=1;
		if(i>639)	{	exit_f=1;	}
		delay(20);
	}while(exit_f == 0);
	printf("Please press ESC key & Exit!\n");
	exit_f=0;
	do{	delay(20);
		if(bioskey(1)!=0){
			scan_keyboard_result=bioskey(0);
			if(scan_keyboard_result==0x011b)	exit_f=1;
		}
	}while(exit_f == 0);
	closegraph();
}

void	main(void)
{       int n;
	i=0;
	j=0;
	exit_f=0;
//	Tm=2.0;Tm_=0.5;T=0.2;l=20.0;K=8.0;	//Tm_<Tm時(shí)，會(huì)產(chǎn)生振玲現(xiàn)象，本行是有振玲條件，同學(xué)還通過改變Tm_觀察波形變化
	Tm=2.0;Tm_=3.0;T=0.2;l=20.0;K=8.0;	//Tm_>Tm時(shí)，不會(huì)產(chǎn)生振玲現(xiàn)象，本行是無振玲條件，同學(xué)還通過改變Tm_觀察波形變化
	a=(1-exp(-T/Tm_))/(K*(1-exp(-T/Tm)));
	b=-a*exp(-T/Tm);
	c=exp(-T/Tm_);
	d=(1-exp(-T/Tm_));
	for(n=0;n<200;n++)	y1_t_[n]=0.0;
	for(n=0;n<20;n++)	y1_kT_[n]=0.0;
	for(n=0;n<22;n++)	u_kT_[n]=0.0;
	initgraph(&gdrive,&gmode,"e:\\bc\\bgi");
	u_kT_21T=u_kT_[21];
	u_kT_T=u_kT_[1];
	e_kT_T=0.0;
	y_kT=y1_kT_[19];
	e_kT=1.0-y_kT;
	u_kT=a*e_kT+b*e_kT_T+c*u_kT_T+d*u_kT_21T;
	e_kT_T=e_kT;
	y1_t=0.0;
	y1_kT=0.0;
	y1_t_[0]=y1_t;
	y1_kT_[0]=y1_kT;
	u_kT_[0]=u_kT;
	setcolor(GREEN);	moveto(i,250-(int)(150.0*y1_t));lineto(i,250-(int)(150.0*y1_t+1));	//對(duì)象連續(xù)輸出響應(yīng)
	setcolor(BLUE);		moveto(i,250-(int)(150.0*e_kT));lineto(i,250-(int)(150.0*e_kT+1.0));	//離散誤差輸出響應(yīng)
	setcolor(YELLOW);	moveto(i,400-(int)(150.0*u_kT));lineto(i,400-(int)(150.0*u_kT+1.0));
	work0();
}