//儀器系數自動標定算法。
//標定的結果:
//系數: COEFF[3][3]={{4.287,-2.185,-0.1945},
//                {-0.3181,10.56,-0.2134},
//                {-0.3802,-0.1864,0.3893}} 
//
//本底: cps0[3]={0.2782,0.1101,0.6737}

float former[5][3]={{137.4,9.61,0.59},
         {5.05,284.3,0.37},
         {2.43,8.42,5.71},
         {118.9,270.9,3.05},
         {2.89,7.3,0.9}}; //五個模型的含量數據（已知）

float cps[5][3]={{35.396,2.828,37.787},
	{17.154,28.164,31.537},
	{2.236,1.313,17.828},
	{45.424,28.446,66.166},
	{1.543,0.919,4.608}}; //五個模型的測量數據（已知）

float coef[3][3];  //9個小寫字母系數（中間數據）
float COEFF[3][3];  //9個大寫字母系數（待標定）
float cps0[3];    //3個本底計數率（待標定）

float A[3][3] ;	//存放三元一次方程組的系數矩陣
float B[3];	//存放三元一次方程組的常數項
float X[3];	//存放三元一次方程組的解:X1=1.537067,X2=-0.896159,X3=3.592657

float DETA ( ) //行列式計算
{
	return (A[0][0]*A[1][1]*A[2][2]+A[1][0]*A[2][1]*A[0][2]+A[2][0]*A[0][1]*A[1][2]
	 -A[2][0]*A[1][1]*A[0][2]-A[0][0]*A[2][1]*A[1][2]-A[1][0]*A[0][1]*A[2][2]);	
}

void swap (int k)//常數項與k列系數交換
{
	int i;
	float t;
	for (i=0;i<3;i++) {
		t=B[i];
		B[i]=A[i][k];
		A[i][k]=t;
		}
}	

void XYZ ( )	//求解三元一次方程組
{
	float d;
	d=DETA( );	//求系數行列式的值
	swap (0) ;	//將x1的系數與常數項交換
	X[0]=DETA()/d;	//解出x1
	swap (0) ;	//恢復原來的系數與常數項
	swap (1) ;	//將x2的系數與常數項交換
	X[1]=DETA()/d;	//解出x2
	swap (1) ;	//恢復原來的系數與常數項
	swap (2) ;	//將x3的系數與常數項交換
	X[2]=DETA()/d;	//解出x3
	swap (2) ;	//恢復原來的系數與常數項
}

void RTK (int k)	//求解第k行的3個小寫字母系數
{
	int i,j;
	for (i=0;i<3;i++)	//準備方程組的系數
		for (j=0;j<3;j++) A[i][j]=former[i][j];
	for (i=0;i<3;i++) B[i]=cps[i][k]-cps0[k];//準備方程組的常數項
	XYZ ();	//解出3個小寫字母系數
	for (i=0;i<3;i++) coef[k][i]=X[i];//保存3個小寫字母系數
}

void conv ( )	//將9個小寫字母系數轉換成9個大寫字母系數
{
	float d;
	int i,j;
	for (i=0;i<3;i++) //取9個小寫字母系數	
		for (j=0;j<3;j++) A[i][j]=coef[i][j];
	d=DETA( );	//求系數行列式的值
	COEFF[0][0]=(A[1][1]*A[2][2]-A[2][1]*A[1][2])/d;//計算9個大寫字母系數
	COEFF[0][1]=(A[2][1]*A[0][2]-A[0][1]*A[2][2])/d;
	COEFF[0][2]=(A[0][1]*A[1][2]-A[0][2]*A[1][1])/d;
	COEFF[1][0]=(A[1][2]*A[2][0]-A[1][0]*A[2][2])/d;
	COEFF[1][1]=(A[0][0]*A[2][2]-A[2][0]*A[0][2])/d;
	COEFF[1][2]=(A[1][0]*A[0][2]-A[0][0]*A[1][2])/d;
	COEFF[2][0]=(A[1][0]*A[2][1]-A[1][1]*A[2][0])/d;
	COEFF[2][1]=(A[0][1]*A[2][0]-A[0][0]*A[2][1])/d;
	COEFF[2][2]=(A[0][0]*A[1][1]-A[0][1]*A[1][0])/d;
}
	
 main ( )
{
	int i,j,k=1;
	float c;
	for (i=0;i<3;i++) cps0[i]=0; //3個本底計數率初始化為0。
	for (i=0;i<5;i++) {  //進行5次迭代運算。
		for (j=0;j<3;j++) RTK (j);//計算9個小寫字母系數
		for (j=0;j<3;j++) 	//計算零值模型各道的凈計數率
			B[j]=coef[j][0]*former[4][0]+
			coef[j][1]*former[4][1]+coef[j][2]*former[4][2];
		for (j=0;j<3;j++) cps0[j]=cps[4][j]-B[j];//計算各道的本底計數率
		}
	conv ( );	//將9個小寫字母系數轉換成9個大寫字母系數。
	for (i=0;i<3;i++) 	//計算混合模型各元素的含量
		B[i]=COEFF[i][0]*(cps[3][0]-cps0[0])+
		COEFF[i][1]*(cps[3][1]-cps0[1])+COEFF[i][2]*(cps[3][2]-cps0[2]);
	c=former[3][0]/B[0];
	if ( c<0.96 || c>1.04 ) k=0;//鈾含量誤差超過百分之四,標定失敗
	c=former[3][1]/B[1];
	if ( c<0.94 || c>1.06 ) k=0;//釷含量誤差超過百分之六,標定失敗
	c=former[3][2]/B[2];
	if ( c<0.88 || c>1.12 ) k=0;//鉀含量誤差超過百分之十二,標定失敗
	while (1) ; //在這一行設置斷點,中止程序運行,以便觀察程序運行的結果 
}