//?Classical.cpp:?implementation?of?the?CClassical?class.
//
//////////////////////////////////////////////////////////////////////

#include?“stdafx.h“
#include?“Classical.h“
#include?“public.h“

//////////////////////////////////////////////////////////////////////
//?Construction/Destruction
//////////////////////////////////////////////////////////////////////



//????＝＝＝?《測量平差程序設計》－第二章?經典平差?源程序?＝＝＝
//????＝＝＝?《測量平差程序設計》－第二章?經典平差?源程序?＝＝＝


//////////////////////////////////////////////////////////////////////
//??參數平差程序（觀測值獨立）
//?參數：
//?int?n?-?觀測值個數
//?int?t?-?參數個數
//?double?A[]?-?誤差方程系數矩陣，數組長度為n×t
//?double?L[]?-?觀測值向量數組長度為n
//?double?P[]?-?觀測值權矩陣的對角線元素數組長度為n
//?double?X[]?-?參數向量函數調用前為參數的先驗值調用后為參數平差值，
//??????????????數組長度為t
//?double?Qx[]?-?函數調用前為參數的先驗權逆陣，調用后為參數平差值的權逆陣，
//???????????????僅存下三角矩陣，數組長度為t×（t+1）/2?
//?double?V[]?-?殘差向量，數組長度為n
//?返回值：驗后單位權中誤差
//////////////////////////////////////////////////////////////////////////
double?Adjust（int?nint?t?double?A[]double?L[]double?P[]
			???double?X[]double?N[]double?V[]）
{
	double?*U=new?double?[t];?//臨時數組，存法方程自由項
	
	//??法方程系數矩陣與自由項清零?
	for（int?i=0;i	for（int?j=0;j	
	//??組成法方程式?
	for（int?k=0;k	{
		for（int?i=0;i		{
			U[i]+=A[k*t+i]*P[k]*L[k];
			for（int?j=0;?j<=i;?j++）
			{
				N[ij（ij）]+=A[k*t+i]*P[k]*A[k*t+j];
			}
		}
	}

	//??計?算?X?	
	if（!inverse（Nt））
	{
		delete[]U;?
		return?-1.0;
	}
	for（i=0;i	{
		X[i]=0.0;
		for（int?j=0;j	}
	delete?[]U;
	
	//???計?算?V?
	double?pvv=0;
	for（i=0;i	{
		double?vi=-L[i];
		for（int?j=0;j		V[i]=vi;
		pvv+=vi*vi*P[i];
	}
	
	if（n==t）return?0.0;?//n=t無多余觀測
	return?sqrt（pvv/（n-t））;

}


//////////////////////////////////////////////////////////////////////////
//?????參數平差例題（觀測值獨立）
void?Adjust_Example（）?
{
	char?*file=“例題\\參數平差\\data.txt“;
	FILE?*fp=fopen（file“r“）;
	if（fp==NULL）
	{
		MyBreak（“數據文件打不開“）;
		return;
	}

	int?nt;
	fscanf（fp“%d%d“&n&t）;

????double?*A=new?double[n*t];
????double??*L=new?double[n];
	double??*P=new?double[n];
	double??*X=new?double[t];
	double??*N=new?double[t*（t+1）/2];
	double??*V=new?double[n];
	
	for（int?i=0;i	{
		for（int?j=0;j		{
			fscanf（fp“%lf“A+i*t+j）;
		}
		fscanf（fp“%lf“L+i）;?//?讀取誤差方程自由項
	}

	for（i=0;i	{
		fscanf（fp“%lf“P+i）;
	}
	fclose（fp）;
	
	fp=fopen（“例題\\參數平差\\result.txt““w“）;
	
	fprintf（fp“===?參數平差（觀測值獨立）計算?===\n“）;
	fprintf（fp“\n誤差方程:\n“）;
	
	for（i=0;i	{
		for（int?j=0;j			fprintf（fp“%7.3lf?“A[i*t+j]）;
		fprintf（fp“%7.3lf?\n“L[i]）;
	}

	double?m=Adjust（?n?t?A?L??P??X?N?V）;
	
	fprintf（fp“參數平差值及其中誤差：\n“）;
	for（i=0;i	{
		fprintf（fp“%4d?%10.4lf?%10.4lf\n“i+1X[i]sqrt（N[ij（ii）]）*m）;
	}
	
	fprintf（fp“\n最小二乘殘差（V）：\n“）;
	for（i=0;i