经典平差程序

// 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