测绘程序设计实习 CSU

1、实习目的与内容

在这里插入图片描述

2、需求分析

1、 总体需求:控制网平差程序对野外控制网观测数据进行平差数据处理,其目的就是根据最小二乘原理,消除网中的各种几何矛盾,求出全网各待定元素(未知点的平面坐标或 三维坐标)。
2、 功能需求:
(1)优化设计:根据控制网的观测精度与网形,全面评定网的精度。
(2)数据输入:文件数据导入。
(3)概算:未知点近似坐标的推算。
(4)平差计算:对观测数据进行精密平差计算,得到平差后的点位坐标、方向观测值和边
长观测值及相应的改正数,并进行精度评定。
(5)成果输出:概算结果输出,平差结果报表输出,控制网图形输出, 绘制误差椭圆,绘制坐标轴和比例尺。

3、总体设计

1、功能设计
数据输入 坐标概算 平差计算 图形绘制
2、系统架构设计
目总共包括6个类:主对话框类、画图对话框类、控制网平差类以及涉及的角度类、矩阵类和常用功能类。控制网平差类主要包括三个功能模块:数据输入、坐标概算、平差计算和图形绘制。
在这里插入图片描述
在这里插入图片描述
在这里插入图片描述

3、界面设计
在这里插入图片描述

4、主要代码

4.1、 变量名类(ValueName.h)的主要代码:

#pragma once
#include "Angle.h"

/*********************************************************
* 文件名:ValueName.h                                    *
*                                                        *
* 描述:该头文件主要定义一些控制网平差类所用到的结构体   *
		变量或枚举变量。包括:控制点结构体CrtlPoint、角度*
		观测值结构体ObsAngle、距离观测结构体ObaDist和    *
		左右角枚举变量AngleStyle					     *
**
* 历史:** 日期**    理由**  签名**                      *
*		  2021/7/6    无     苑艺博                      *
**
* 外部过程:无                                           *
*********************************************************/

//控制点类型
enum PointStyle
{
	KnownPoint,//已知点
	UnknownPoint//未知点
};

/**********************************************************
* 结构体名:CtrlPoint								      *
*														  *
* 描述:控制点结构体									  *
*														  *
* 使用方法:成员有坐标dx,dy和是否求出的状态bState        *
*														  *
* 历史:**日期** **理由** **签名**						  *
*       2021/7/6    无     苑艺博						  *
*********************************************************/
struct CtrlPoint
{
	CString strID;
	double dx;//坐标x
	double dy;//坐标y;
	bool bState;//是否求出的状态,TRUE求出,FALSE未求出
	int iNum;//编号
	PointStyle ePointStyle;
	double dE, dF;//误差椭圆的长半轴和短半轴
	CAngle dAlfa;//误差椭圆长半轴的方位角
	double dMx, dMy, dMk;//点位误差
};


/********************************************************************
* 结构体名:ObsAngle												*
*																	*
* 描述:角度观测值结构体											*
*																	*
* 使用方法:成员有站点StationPoint和照准点ObjPoint(CtrlPoint指针) *
*			和角度值AngleValue(CAngle类型)						*
*																	*
* 历史:**日期** **理由** **签名**									*
*       2021/7/6    无     苑艺博									*
********************************************************************/
struct ObsAngle
{
	CtrlPoint* StationPoint;//站点
	CtrlPoint* ObjPoint;//照准点
	CAngle AngleValue;//角度值
	double AngleError;//角度改正数
};


/********************************************************************
* 结构体名:ObsDist  												*
*																	*
* 描述:距离观测值结构体											*
*																	*
* 使用方法:成员有起点StartPoint和终点EndPoint(CtrlPoint指针)     *
*			和距离值DistValue(double 类型)						*
*																	*
* 历史:**日期** **理由** **签名**									*
*       2021/7/6    无     苑艺博									*
********************************************************************/
struct ObsDist
{
	CtrlPoint* StartPoint;//起点
	CtrlPoint* EndPoint;//终点
	double DistValue;//距离值
	double DistError;//距离改正数
};

4.2、主对话框类头文件(ControlNetworkDlg.h)的代码:

// ControlNetworkDlg.h: 头文件
//

#pragma once


// CControlNetworkDlg 对话框
class CControlNetworkDlg : public CDialogEx
{
// 构造
public:
	CControlNetworkDlg(CWnd* pParent = nullptr);	// 标准构造函数

// 对话框数据
#ifdef AFX_DESIGN_TIME
	enum { IDD = IDD_CONTROLNETWORK_DIALOG };
#endif

	protected:
	virtual void DoDataExchange(CDataExchange* pDX);	// DDX/DDV 支持


// 实现
private:
	ControlNetworkAdjust CNA_Cal;//定义控制网平差类对象
	ControlNetworkDrawingDlg CND;//定义控制网绘制对话框类
protected:
	HICON m_hIcon;

	// 生成的消息映射函数
	virtual BOOL OnInitDialog();
	afx_msg void OnSysCommand(UINT nID, LPARAM lParam);
	afx_msg void OnPaint();
	afx_msg HCURSOR OnQueryDragIcon();
	DECLARE_MESSAGE_MAP()
public:
	friend void ReadFileContent(CStdioFile& sf,CString& strFileContent);//从文件中读取数据到CString字符串中
	afx_msg void OnBnClickedReadfiledata();
	CString strData;
	CString strOutData;
	afx_msg void OnBnClickedCalapprocoord();
	afx_msg void OnBnClickedAdjustcal();
	afx_msg void OnBnClickedNetdrawing();
};

4.3、主对话框类cpp文件(ControlNetworkDlg.cpp)的主要代码:

//从文件中读取数据到字符串中
void ReadFileContent(CStdioFile& sf,CString& strFileContent)
{
	CString strLine;
	strFileContent.Empty();
	bool bEOF = sf.ReadString(strLine);
	while (bEOF)
	{
		strFileContent += strLine;

		bEOF = sf.ReadString(strLine);
		if (bEOF) strFileContent += _T("\r\n");
	}
}

//读取数据
void CControlNetworkDlg::OnBnClickedReadfiledata()
{
	// TODO: 在此添加控件通知处理程序代码

	CFileDialog dlgFile(TRUE, _T("txt"), NULL,
		OFN_ALLOWMULTISELECT | OFN_EXPLORER,
		_T("(文本文件)|*.dat"));
	if (dlgFile.DoModal() == IDCANCEL) return;
	CString strFileName = dlgFile.GetPathName();
	setlocale(LC_ALL, "chs");
	CStdioFile sf;
	if (!sf.Open(strFileName, CFile::modeRead)) return;


	ReadFileContent(sf, strData);//读取并在文本框中显示文件内容

	sf.SeekToBegin();//文件指针定位到开头
	CNA_Cal.ReadFileData(sf);//从文件读取数据

	UpdateData(false);
}

//坐标概算
void CControlNetworkDlg::OnBnClickedCalapprocoord()
{
	// TODO: 在此添加控件通知处理程序代码

	//计算未知点的近似坐标
	if (CNA_Cal.CalUnknownPoint())
	{
		CNA_Cal.PrintRoughCalResult(strOutData);//输出概算结果到文本框中

		UpdateData(false);
	}
}

//平差计算
void CControlNetworkDlg::OnBnClickedAdjustcal()
{
	// TODO: 在此添加控件通知处理程序代码

	if (CNA_Cal.AdjustCal())
	{

		CNA_Cal.PrintAdjustCalResult(strOutData);//输出平差计算结果到文本框中

		UpdateData(false);
	}
}

//图形绘制
void CControlNetworkDlg::OnBnClickedNetdrawing()
{
	// TODO: 在此添加控件通知处理程序代码
	if (CNA_Cal.JudgeDrawConditon())
	{
		CND.CNA_Draw = &CNA_Cal;//向子对话框中传递控制网平差对象

		CND.DoModal();//显示子对话框
	}
}

4.4、控制网平差类头文件(ControlNetworkAdjust.h)的代码:

#pragma once
#include "ValueName.h"
#include "Matrix.h"

const double EPSILON = 1.0E-12;
const double PI = 4.0 * atan(1.0);


/*******************************************************************
* 类名:ControlNetworkAdjust								     *
*															   *
* 描述:控制网平差类											   *
*															   *
* 使用方法:													   *
*															   * 
* 历史:**日期** **理由** **签名**								   *
*       2021/7/6    无     苑艺博								   *
*															   *
* 外部类:无													   *
*********************************************************************/
class ControlNetworkAdjust
{
private:
	int m_KnownPointCount;//已知点数
	CtrlPoint* m_KnownPoint;//已知点指针
	int m_UnknownPointCount;//未知点数
	CtrlPoint* m_UnknownPoint;//未知点指针
	int m_ObsAngleCount;//观测角度数
	ObsAngle* m_ObsAngle;//观测角度指针
	int m_ObsDistCount;//观测边长数
	ObsDist* m_ObsDist;//观测边长指针

	//以下是绘图用到的变量
	double dXmin;//坐标x的最小值
	double dYmin;//坐标y的最小值
	double dDetX;//坐标x的最大值和最小值的差值
	double dDetY;//坐标y的最大值和最小值的差值
	double dScale;//绘图比例
public:
	ControlNetworkAdjust();//构造函数
	ControlNetworkAdjust(ControlNetworkAdjust& A);//拷贝构造函数
	~ControlNetworkAdjust();//析构函数
private:
	friend void ReadFileContent(CStdioFile& sf, CString& strFileContent);//从文件中读取数据到CString字符串中
	
	void ClearVariables(); //清空变量

	CAngle Azi(const CtrlPoint& P1, const CtrlPoint& P2);//已知两个控制点,求P1->P2的方位角
	double DIST(const CtrlPoint& P1, const CtrlPoint& P2);//已知两个控制点,求P1和P2的距离
	int SplitStringArray(CString str, char split, CStringArray& aStr);// 分割字符串,得到字符串数组


	CtrlPoint* SearchPoint(CString strName);//根据点名寻找符合的已知点或未知点,返回CtrlPoint指针变量
	bool JudgeUnknownPoint();//判断未知点是否全部求出,若全部求出返回真,反之返回假   
	ObsAngle* SearchObsAngle_StationPtIsLookedPt(CString strName);//寻找测站点为所找点且照准点坐标已求出的角度观测值
	ObsAngle* SearchObsAngle(CString strStationName, CString strObjName);//寻找测站点和照准点的点号符合的角度观测值
	ObsDist* SearchObsDist_EndPtIsLookedPt(CString strName);//寻找终点为所找点且起点坐标已求出的边长观测值
	ObsDist* SearchObsDist_StartPtIsLookedPt(CString strName);//寻找起点为所找点且终点坐标已求出的边长观测值
	ObsAngle* SearchZeroDirection(int iNum);//寻找任一角度观测值中测站的零方向所对应的角度观测值


	void ObsAngleErrorMatrix(CMatrix& B, CMatrix& L, CMatrix& P);//构建观测角度误差方程和权矩阵
	void ObsDistErrorMatrix(CMatrix& B, CMatrix& L, CMatrix& P);//构建观测边长误差方程和权矩阵
	void CalBLP(CMatrix& B, CMatrix& L, CMatrix& P);//计算B、L、P
	void CalAdjustedAngleDist(CMatrix& V);//计算角度和边长改正数以及改正后边长和角度
	void PrecisionEvaluate(CMatrix& Qx,double m0);//精度评定,包括点位误差和误差椭圆参数


	void XYCal();//计算dXmin,dYmin,dDetX,dDetY
	void SetScale(CRect& rect);//设置自适应窗口的绘图比例
	POINT LP2CP(const CtrlPoint& Point, CRect& rect);//从逻辑坐标转换到客户区坐标
	void TrianglePointDrawing(CDC* pDC, const POINT& CentralPoint);//以指定点为中心画小三角形来表示已知点


	void PrintCMatrix(CMatrix& A, CString strFileName);//输出矩阵
public:
	void ReadFileData(CStdioFile& sf);//从文件中读取数据
	bool CalUnknownPoint();//计算未知点近似坐标
	void PrintRoughCalResult(CString& strOut);//向文本框中输出概算坐标
	bool AdjustCal();//平差计算
	void PrintAdjustCalResult(CString& strOut);//向文本框中输出平差坐标,并存入文件中

	bool JudgeDrawConditon();//判断是否达到绘图条件
	void ControlNetworkDrawing(CDC* pDC, CRect& rect);//绘制三角网
};

4.5、控制网平差类cpp文件(ControlNetworkAdjust.cpp)的主要代码:

#include "pch.h"
#include "ControlNetworkAdjust.h"
#include "CommonSurveyFunctions.h"
#include <locale>

//构造函数
ControlNetworkAdjust::ControlNetworkAdjust()
{
	m_KnownPointCount = 0;
	m_KnownPoint = NULL;

	m_UnknownPointCount = 0;
	m_UnknownPoint = NULL;

	m_ObsAngleCount = 0;
	m_ObsAngle = NULL;

	m_ObsDistCount = 0;
	m_ObsDist = NULL;
}

//拷贝构造函数
ControlNetworkAdjust::ControlNetworkAdjust(ControlNetworkAdjust& A)
{
	m_KnownPointCount = A.m_KnownPointCount;
	m_KnownPoint = new CtrlPoint[m_KnownPointCount];
	m_KnownPoint = A.m_KnownPoint;

	m_UnknownPointCount = A.m_UnknownPointCount;
	m_UnknownPoint = new CtrlPoint[m_UnknownPointCount];
	m_UnknownPoint = A.m_UnknownPoint;

	m_ObsAngleCount = A.m_ObsAngleCount;
	m_ObsAngle = new ObsAngle[m_ObsAngleCount];
	m_ObsAngle = A.m_ObsAngle;

	m_ObsDistCount = A.m_ObsDistCount;
	m_ObsDist = new ObsDist[m_ObsDistCount];
	m_ObsDist = A.m_ObsDist;
}

//析构函数
ControlNetworkAdjust::~ControlNetworkAdjust()
{
	if (m_KnownPoint != NULL)
	{
		delete[] m_KnownPoint;
		m_KnownPoint = NULL;
	}

	if (m_UnknownPoint != NULL)
	{
		delete[] m_UnknownPoint;
		m_UnknownPoint = NULL;
	}

	if (m_ObsAngle != NULL)
	{
		delete[] m_ObsAngle;
		m_ObsAngle = NULL;
	}

	if (m_ObsDist != NULL)
	{
		delete[] m_ObsDist;
		m_ObsDist = NULL;
	}
}

//清空变量
void ControlNetworkAdjust::ClearVariables()
{
	if (m_KnownPoint != NULL)
	{
		delete[] m_KnownPoint;
		m_KnownPoint = NULL;
	}

	if (m_UnknownPoint != NULL)
	{
		delete[] m_UnknownPoint;
		m_UnknownPoint = NULL;
	}

	if (m_ObsAngle != NULL)
	{
		delete[] m_ObsAngle;
		m_ObsAngle = NULL;
	}

	if (m_ObsDist != NULL)
	{
		delete[] m_ObsDist;
		m_ObsDist = NULL;
	}
}

//已知两个控制点,求P1->P2的方位角
CAngle ControlNetworkAdjust::Azi(const CtrlPoint& P1, const CtrlPoint& P2)
{
	CAngle angAzi;

	angAzi(RAD) = Azimuth(P1.dx, P1.dy, P2.dx, P2.dy);

	return angAzi;
}

//已知两个控制点,求P1和P2的距离
double ControlNetworkAdjust::DIST(const CtrlPoint& P1, const CtrlPoint& P2)
{
	return Dist(P1.dx, P1.dy, P2.dx, P2.dy);
}

//分割字符串
int ControlNetworkAdjust::SplitStringArray(CString str, char split, CStringArray& aStr)
{
	int startIdx = 0;
	int idx = str.Find(split, startIdx);
	aStr.RemoveAll();//先清空
	while (-1 != idx)
	{
		CString sTmp = str.Mid(startIdx, idx - startIdx);
		aStr.Add(sTmp);
		startIdx = idx + 1;
		idx = str.Find(split, startIdx);
	}
	CString sTmp = str.Right(str.GetLength() - startIdx);
	if (!sTmp.IsEmpty())
		aStr.Add(sTmp);
	return aStr.GetSize();
}

//根据点名寻找符合的已知点或未知点,返回CtrlPoint指针变量
CtrlPoint* ControlNetworkAdjust::SearchPoint(CString strName)
{
	strName.Trim();
	CString strPointID;//修整后的点名

	//遍历已知点
	for (int i = 0; i < m_KnownPointCount; i++)
	{
		strPointID = m_KnownPoint[i].strID.Trim();
		if (strName == strPointID)
		{
			return &m_KnownPoint[i];
		}
	}

	//遍历未知点
	for (int i = 0; i < m_UnknownPointCount; i++)
	{
		strPointID = m_UnknownPoint[i].strID.Trim();
		if (strPointID == strName)
		{
			return &m_UnknownPoint[i];
		}
	}
}


//判断未知点是否全部求出,若全部求出返回真,反之返回假
bool ControlNetworkAdjust::JudgeUnknownPoint()
{
	for (int i = 0; i < m_UnknownPointCount; i++)
	{
		if (m_UnknownPoint[i].bState == FALSE)
		{
			return FALSE;
		}
	}

	return TRUE;
}

//寻找测站点为所找点且照准点坐标已求出的角度观测值
ObsAngle* ControlNetworkAdjust::SearchObsAngle_StationPtIsLookedPt(CString strName)
{
	strName.Trim();
	CString strStationPtID;//修整后角度观测值测站点点名

	//遍历角度观测值
	for (int i = 0; i < m_ObsAngleCount; i++)
	{
		strStationPtID = m_ObsAngle[i].StationPoint->strID.Trim();
		if (strName == strStationPtID && m_ObsAngle[i].ObjPoint->bState == TRUE)
		{
			return &m_ObsAngle[i];
		}
	}

	return NULL;
}

//寻找测站点和照准点的点号符合的角度观测值
ObsAngle* ControlNetworkAdjust::SearchObsAngle(CString strStationName, CString strObjName)
{
	strStationName.Trim();
	strObjName.Trim();

	CString strStationID;//修整后的测站点的点号
	CString strObjID;//修整后的照准点的点号

	//遍历角度观测值
	for (int i = 0; i < m_ObsAngleCount; i++)
	{
		strStationID = m_ObsAngle[i].StationPoint->strID.Trim();
		strObjID = m_ObsAngle[i].ObjPoint->strID.Trim();

		if (strStationID == strStationName && strObjID == strObjName)
		{
			return &m_ObsAngle[i];
		}
	}

	return NULL;
}

//寻找终点为所找点且起点坐标已求出的边长观测值
//输入:strName为所找点点名
//输出:满足条件的边长观测值的指针变量
ObsDist* ControlNetworkAdjust::SearchObsDist_EndPtIsLookedPt(CString strName)
{
	CString strEndID;//修整后边长终点点名

	strName.Trim();

	for (int i = 0; i < m_ObsDistCount; i++)
	{
		strEndID = m_ObsDist[i].EndPoint->strID.Trim();

		if (strName == strEndID && m_ObsDist[i].StartPoint->bState == TRUE)
		{
			return &m_ObsDist[i];
		}
	}

	return NULL;
}

//寻找起点为所找点且终点坐标已求出的边长观测值
//输入:strName为所找点点名
//输出:满足条件的边长观测值的指针变量
ObsDist* ControlNetworkAdjust::SearchObsDist_StartPtIsLookedPt(CString strName)
{
	CString strStartID;//修整后边长终点点名

	strName.Trim();

	for (int i = 0; i < m_ObsDistCount; i++)
	{
		strStartID = m_ObsDist[i].StartPoint->strID.Trim();

		if (strName == strStartID && m_ObsDist[i].EndPoint->bState == TRUE)
		{
			return &m_ObsDist[i];
		}
	}

	return NULL;
}

//寻找任一角度观测值中测站的零方向所对应的角度观测值
//输入:任一角度观测值观测值在数组中的位置iNum
//输出:该角度观测值中测站的零方向所对应的角度观测值,类型为ObsAngle指针变量
ObsAngle* ControlNetworkAdjust::SearchZeroDirection(int iNum)
{
	//向后遍历角度观测值
	for (int i = iNum; i >= 0; i--)
	{
		//如果角度观测值为0,说明该角度观测值所对应的方向即为零方向
		if (m_ObsAngle[i].AngleValue(DMS) == 0)
		{
			return &m_ObsAngle[i];
		}
	}

	return NULL;
}




//构建观测角度误差方程和权矩阵
//输入:无
//输出:B、L均为误差方程V=BX+L中角度和边长合矩阵,P为角度和边长合权阵,且在B、L、和P中角度在前,边长在后
void ControlNetworkAdjust::ObsAngleErrorMatrix(CMatrix& B, CMatrix& L, CMatrix& P)
{
	double p = 180 / PI * 3600;//使弧度制转为秒需乘以的常数

	//遍历角度观测值,构建观测角度误差方程和权矩阵
	for (int i = 0; i < m_ObsAngleCount; i++)
	{
		//计算照准方向的方位角近似值
		CAngle SightAzi = Azi(*m_ObsAngle[i].StationPoint, *m_ObsAngle[i].ObjPoint);

		//计算零方向的方位角近似值
		ObsAngle* ZeroDirection = SearchZeroDirection(i);
		CAngle ZeroAzi;
		if (ZeroDirection != NULL)
		{
			ZeroAzi = Azi(*ZeroDirection->StationPoint, *ZeroDirection->ObjPoint);
		}

		//计算照准方向的边长
		double dDist = DIST(*m_ObsAngle[i].StationPoint, *m_ObsAngle[i].ObjPoint);

		//计算误差方程B的系数
		double dA = sin(SightAzi(RAD)) / dDist * p / 1000;
		double dB = cos(SightAzi(RAD)) / dDist * (-p) / 1000;

		//计算常数项
		CAngle Angle = SightAzi - ZeroAzi;//计算夹角
		//如果计算结果小于0
		if (Angle(DEG) < 0)
		{
			Angle(DEG) = Angle(DEG) + 360;
		}
		CAngle AngleError = Angle - m_ObsAngle[i].AngleValue;
		double dAngleErrorSec = AngleError(DEG) * 3600;//转为秒

		//构建误差方程和权阵
		{
			//构建B
			//如果测站点是未知点
            if (m_ObsAngle[i].StationPoint->ePointStyle == UnknownPoint)
			{
				//计算各参数在误差方程中B中的列数
				int iNumx = 2 * m_ObsAngle[i].StationPoint->iNum;//坐标x在误差方程中B中的列数
				int iNumy = iNumx + 1;//坐标y在误差方程中B中的列数

				B(i, iNumx) = dA;
				B(i, iNumy) = dB;
			}
			//如果照准点是未知点
			if (m_ObsAngle[i].ObjPoint->ePointStyle == UnknownPoint)
			{
				//计算各参数在误差方程中B中的列数
				int iNumx = 2 * m_ObsAngle[i].ObjPoint->iNum;//坐标x在误差方程中B中的列数
				int iNumy = iNumx + 1;//坐标y在误差方程中B中的列数

				B(i, iNumx) = -dA;
				B(i, iNumy) = -dB;
			}

			//构建L
			L(i, 0) = dAngleErrorSec;

			//构建P
			P(i, i) = 1;
		}
	}
}

//构建观测边长误差方程和权矩阵
//输入:无
//输出:B、L均为误差方程V=BX+L中角度和边长合矩阵,P为角度和边长合权阵,且在B、L、和P中角度在前,边长在后
void ControlNetworkAdjust::ObsDistErrorMatrix(CMatrix& B, CMatrix& L, CMatrix& P)
{
	//遍历角度观测值,构建观测角度误差方程和权矩阵
	for (int i = 0; i < m_ObsDistCount; i++)
	{
		//计算边长的方位角近似值
		CAngle DistAzi = Azi(*m_ObsDist[i].StartPoint, *m_ObsDist[i].EndPoint);

		//计算边长近似值
		double dDist = DIST(*m_ObsDist[i].StartPoint, *m_ObsDist[i].EndPoint);

		//计算常数项
		double dDistErrorMm = (dDist - m_ObsDist[i].DistValue) * 1000;//转为毫米

		//构建观测边长误差方程和权矩阵
		{
			//构建B
			//如果起点为未知点
			if (m_ObsDist[i].StartPoint->ePointStyle == UnknownPoint)
			{
				//计算各参数在误差方程中B中的列数
				int iNumx = 2 * m_ObsDist[i].StartPoint->iNum;//坐标x在误差方程中B中的列数
				int iNumy = iNumx + 1;//坐标y在误差方程中B中的列数

				B(m_ObsAngleCount + i, iNumx) = -cos(DistAzi(RAD));
				B(m_ObsAngleCount + i, iNumy) = -sin(DistAzi(RAD));
			}

			//如果终点为未知点
			if (m_ObsDist[i].EndPoint->ePointStyle == UnknownPoint)
			{
				//计算各参数在误差方程中B中的列数
				int iNumx = 2 * m_ObsDist[i].EndPoint->iNum;//坐标x在误差方程中B中的列数
				int iNumy = iNumx + 1;//坐标y在误差方程中B中的列数

				B(m_ObsAngleCount + i, iNumx) = cos(DistAzi(RAD));
				B(m_ObsAngleCount + i, iNumy) = sin(DistAzi(RAD));
			}

			//构建L
			L(m_ObsAngleCount + i, 0) = dDistErrorMm;

			//构建P
			P(m_ObsAngleCount + i, m_ObsAngleCount + i) = 1000 / m_ObsDist[i].DistValue;
		}
	}
}

//计算B、L、P
void ControlNetworkAdjust::CalBLP(CMatrix& B, CMatrix& L, CMatrix& P)
{
	//计算B中角度观测值零方向的系统误差Z的系数
	int iNumz = 2 * m_UnknownPointCount;//Z在B中的列数

	for (int i = 0; i < m_ObsAngleCount; i++)
	{
		//如果是零方向
		if (m_ObsAngle[i].AngleValue(DEG) == 0)
		{
			B(i, iNumz) = -1;

			//把该测站上其他观测值的Z的系数赋值
			for (int j = i + 1; j < m_ObsAngleCount; j++)
			{
				if (m_ObsAngle[j].AngleValue(DEG) != 0)
				{
					B(j, iNumz) = -1;
				}
				else
				{
					break;
				}
			}

			iNumz++;
		}
	}

	//计算误差方程和权阵
	ObsAngleErrorMatrix(B, L, P);
	ObsDistErrorMatrix(B, L, P);
}

//计算角度和边长改正数以及改正后的值
void ControlNetworkAdjust::CalAdjustedAngleDist(CMatrix& V)
{
	//计算改正后的角度观测值
	for (int i = 0; i < m_ObsAngleCount; i++)
	{
		//计算角度改正数,单位秒,保留两位小数
		double dAngleError = V(i, 0);
		CString strAngleError;
		strAngleError.Format(_T("%.2f"), dAngleError);
		m_ObsAngle[i].AngleError = _tstof(strAngleError);

		//计算改正后的角度值
		//由于角度改正数一般不大于60秒,如果大于,则说明观测值存在问题
		if (m_ObsAngle[i].AngleError < 60)
		{
			m_ObsAngle[i].AngleValue(DMS) = m_ObsAngle[i].AngleValue(DMS) + m_ObsAngle[i].AngleError / 10000;//计算改正后角度值
		}
		else
		{
			MessageBox(NULL, _T("角度改正数大于60秒,角度观测值存在问题,请检查!"), _T("警告"), MB_OK);
		}
	}

	//计算改正后的边长观测值
	for (int i = 0; i < m_ObsDistCount; i++)
	{
		//计算边长改正数,单位m,保留四位小数
		double dDistError= V(m_ObsAngleCount + i, 0) / 1000;
		CString strDistError;
		strDistError.Format(_T("%.4f"), dDistError);
		m_ObsDist[i].DistError = _tstof(strDistError);

		//计算改正后的边长
		m_ObsDist[i].DistValue = m_ObsDist[i].DistValue + m_ObsDist[i].DistError;
	}
}

//精度评定,包括点位误差和误差椭圆参数
void ControlNetworkAdjust::PrecisionEvaluate(CMatrix& Qx, double m0)
{
	for (int i = 0; i < m_UnknownPointCount; i++)
	{
		//计算未知点坐标x,y在矩阵Qx中的行数
		int iNumx = 2 * m_UnknownPoint[i].iNum;
		int iNumy = iNumx + 1;

		double Qxx = Qx(iNumx, iNumx);
		double Qyy = Qx(iNumy, iNumy);
		double Qxy = Qx(iNumx, iNumy);

		//计算点位误差
		m_UnknownPoint[i].dMx = m0 * sqrt(Qxx);
		m_UnknownPoint[i].dMy = m0 * sqrt(Qyy);
		m_UnknownPoint[i].dMk = sqrt(m_UnknownPoint[i].dMx * m_UnknownPoint[i].dMx + m_UnknownPoint[i].dMy * m_UnknownPoint[i].dMy);

		//计算误差椭圆参数
		double dAlfa;//长轴方位角
		dAlfa = 0.5 * atan(2 * Qxy / (Qxx - Qyy));
		if (dAlfa < 0)
		{
			dAlfa += PI;
		}
		m_UnknownPoint[i].dAlfa(RAD) = dAlfa;
		m_UnknownPoint[i].dE = m0 * sqrt(Qxx + Qxy * tan(m_UnknownPoint[i].dAlfa(RAD)));
		m_UnknownPoint[i].dF = m0 * sqrt(Qxx + Qxy * tan(m_UnknownPoint[i].dAlfa(RAD) + PI / 2));
	}
}





//计算dXmin,dYmin,dDetX,dDetY
void ControlNetworkAdjust::XYCal()
{
	double dXmax, dYmax;

	dXmin = m_KnownPoint[0].dx;
	dXmax = m_KnownPoint[0].dx;
	dYmin = m_KnownPoint[0].dy;
	dYmax = m_KnownPoint[0].dy;

	//遍历已知点
	for (int i = 0; i < m_KnownPointCount; i++)
	{
		dXmin = m_KnownPoint[i].dx < dXmin ? m_KnownPoint[i].dx : dXmin;
		dXmax = m_KnownPoint[i].dx > dXmax ? m_KnownPoint[i].dx : dXmax;
		dYmin = m_KnownPoint[i].dy < dYmin ? m_KnownPoint[i].dy : dYmin;
		dYmax = m_KnownPoint[i].dy > dYmax ? m_KnownPoint[i].dy : dYmax;
	}

	//遍历未知点
	for (int i = 0; i < m_UnknownPointCount; i++)
	{
		dXmin = m_UnknownPoint[i].dx < dXmin ? m_UnknownPoint[i].dx : dXmin;
		dXmax = m_UnknownPoint[i].dx > dXmax ? m_UnknownPoint[i].dx : dXmax;
		dYmin = m_UnknownPoint[i].dy < dYmin ? m_UnknownPoint[i].dy : dYmin;
		dYmax = m_UnknownPoint[i].dy > dYmax ? m_UnknownPoint[i].dy : dYmax;
	}

	dDetX = dXmax - dXmin;
	dDetY = dYmax - dYmin;
}

//设置自适应窗口的绘图比例
void ControlNetworkAdjust::SetScale(CRect& rect)
{
	double a = rect.Width();
	double b = rect.Height();
	double ry = double(rect.Width()) * 3 / 4 / dDetY;//宽度扩大到2/3,方便在右边画坐标轴和比例尺
	double rx = double(rect.Height()) * 3 / 4 / dDetX;
	dScale = rx < ry ? rx : ry;
}

//从逻辑坐标转换到客户区坐标
POINT ControlNetworkAdjust::LP2CP(const CtrlPoint& Point, CRect& rect)
{
	POINT P;
	P.y = rect.Height() - (Point.dx - dXmin) * dScale;
	P.x = (Point.dy - dYmin) * dScale;
	return P;
}

//以指定点为中心画小三角形来表示已知点
void ControlNetworkAdjust::TrianglePointDrawing(CDC* pDC, const POINT& CentralPoint)
{
	POINT Vertice[3];//三角形三个顶点

	Vertice[0].x = CentralPoint.x;
	Vertice[0].y = CentralPoint.y - 10;

	Vertice[1].x = CentralPoint.x + 8;
	Vertice[1].y = CentralPoint.y + 6;

	Vertice[2].x = CentralPoint.x - 8;
	Vertice[2].y = CentralPoint.y + 6;

	for (int i = 0; i < 2; i++)
	{
		pDC->MoveTo(Vertice[i].x, Vertice[i].y);
		pDC->LineTo(Vertice[i + 1].x, Vertice[i + 1].y);
	}

	pDC->MoveTo(Vertice[0].x, Vertice[0].y);
	pDC->LineTo(Vertice[2].x, Vertice[2].y);
}

//绘制控制网
void ControlNetworkAdjust::ControlNetworkDrawing(CDC* pDC, CRect& rect)
{
	XYCal();//计算dXmin,dYmin,dDetX,dDetY
	SetScale(rect);//设置自适应窗口的绘图比例

	//三角网的原点,可以控制三角网的整体平移
	double dOrgNetX = rect.Width() / 4;
	double dOrgNetY = -rect.Height() / 5;

	CPen pen(PS_SOLID, 1, RGB(0, 0, 0));
	CPen* pOldPen = pDC->SelectObject(&pen);

	POINT StartPoint;//起点在图上的坐标
	POINT EndPoint;//终点在图上的坐标

	//在已知点上画三角形
	StartPoint = LP2CP(m_KnownPoint[0], rect);
	EndPoint = LP2CP(m_KnownPoint[1], rect);
	StartPoint.x = StartPoint.x + dOrgNetX;
	StartPoint.y = StartPoint.y + dOrgNetY;
	EndPoint.x = EndPoint.x + dOrgNetX;
	EndPoint.y = EndPoint.y + dOrgNetY;
	TrianglePointDrawing(pDC, StartPoint);
	TrianglePointDrawing(pDC, EndPoint);

	//根据已知点之间的斜率来画另一条平行线
	int dx = StartPoint.x - EndPoint.x;
	int dy = StartPoint.y - EndPoint.y;
	double dAngle;//垂直于已知点之间直线的直线与X轴的夹角
	if (dx == 0)
	{
		dAngle = PI / 2;
		StartPoint.x = StartPoint.x + 3;
		StartPoint.y = StartPoint.y;
		EndPoint.x = EndPoint.x + 3;
		EndPoint.y = EndPoint.y;
	}
	else
	{
		dAngle = PI / 2 - atan(dy / dx);
		StartPoint.x = StartPoint.x + 3 * cos(dAngle);
		StartPoint.y = StartPoint.y - 3 * sin(dAngle);
		EndPoint.x = EndPoint.x + 3 * cos(dAngle);
		EndPoint.y = EndPoint.y - 3 * sin(dAngle);
	}

	pDC->MoveTo(StartPoint);
	pDC->LineTo(EndPoint);

	//遍历角度观测值,画出照准边
	for (int i = 0; i < m_ObsAngleCount; i++)
	{
		StartPoint = LP2CP(*m_ObsAngle[i].StationPoint, rect);
		EndPoint = LP2CP(*m_ObsAngle[i].ObjPoint, rect);
		StartPoint.x = StartPoint.x + dOrgNetX;
		StartPoint.y = StartPoint.y + dOrgNetY;
		EndPoint.x = EndPoint.x + dOrgNetX;
		EndPoint.y = EndPoint.y + dOrgNetY;

		pDC->MoveTo(StartPoint);
		pDC->LineTo(EndPoint);
	}

	//画误差椭圆
	
	//误差椭圆的中心
	double dOrgX;
	double dOrgY;

	//遍历未知点,画误差椭圆
	for (int i = 0; i < m_UnknownPointCount; i++)
	{
		double dStartX, dStartY, dEndX, dEndY;

		//误差椭圆的中心坐标
		POINT Org;
		Org = LP2CP(m_UnknownPoint[i], rect);//坐标转换
		dOrgX = Org.x + dOrgNetX;
		dOrgY = Org.y + dOrgNetY;


		double dF = m_UnknownPoint[i].dF * 5000;
		double dE = m_UnknownPoint[i].dE * 5000;
		double dAlfa = m_UnknownPoint[i].dAlfa(RAD);

		//绘制短半轴
		dStartX = (dF * sin(dAlfa)) * dScale + dOrgX;
		dStartY = (-dF * cos(dAlfa)) * dScale + dOrgY;
		dEndX = (-dF * sin(dAlfa)) * dScale + dOrgX;
		dEndY = (dF * cos(dAlfa)) * dScale + dOrgY;
		pDC->MoveTo(dStartX, dStartY);
		pDC->LineTo(dEndX, dEndY);

		//绘制长半轴
		dStartX = (-dE * cos(dAlfa)) * dScale + dOrgX;
		dStartY = (-dE * sin(dAlfa)) * dScale + dOrgY;
		dEndX = (dE * cos(dAlfa)) * dScale + dOrgX;
		dEndY = (dE * sin(dAlfa)) * dScale + dOrgY;
		pDC->MoveTo(dStartX, dStartY);
		pDC->LineTo(dEndX, dEndY);

		double ex, fy;
		ex = dE;
		fy = 0;
		//转换到长半轴方向上
		dStartX = (ex * cos(dAlfa)
			- fy * sin(dAlfa)
			) * dScale + dOrgX;
		dStartY = (fy * cos(dAlfa)
			+ ex * sin(dAlfa)
			) * dScale + dOrgY;
		pDC->MoveTo(dStartX, dStartY);

		for (int i = 6; i <= 360; i += 6)
		{
			//在坐标轴方向的坐标
			ex = dE * cos((i / 180.0) * PI);
			fy = dF * sin((i / 180.0) * PI);

			//转换到长半轴方向上
			dEndX = (ex * cos(dAlfa)
				- fy * sin(dAlfa)
				) * dScale + dOrgX;
			dEndY = (fy * cos(dAlfa)
				+ ex * sin(dAlfa)
				) * dScale + dOrgY;
			pDC->LineTo(dEndX, dEndY);
		}
	}
	



	//画坐标轴
	double dLenX = double(rect.Height()) / 10;//坐标轴X轴的长度
	double dLenY = double(rect.Width())  / 10;//坐标轴Y轴的长度

	//坐标轴的原点坐标
	double dOrgCoordX =  rect.left + double(rect.Width()) / 10;
	double dOrgCoordY = rect.bottom - double(rect.Height()) / 10;

	//画Y轴
	pDC->MoveTo(dOrgCoordX, dOrgCoordY);
	pDC->LineTo(dOrgCoordX + dLenY, dOrgCoordY);
	pDC->LineTo(dOrgCoordX + dLenY - 5, dOrgCoordY - 5);
	pDC->MoveTo(dOrgCoordX + dLenY, dOrgCoordY);
	pDC->LineTo(dOrgCoordX + dLenY - 5, dOrgCoordY + 5);

	//画X轴
	pDC->MoveTo(dOrgCoordX, dOrgCoordY);
	pDC->LineTo(dOrgCoordX, dOrgCoordY - dLenX);
	pDC->LineTo(dOrgCoordX - 5, dOrgCoordY - dLenX + 5);
	pDC->MoveTo(dOrgCoordX, dOrgCoordY - dLenX);
	pDC->LineTo(dOrgCoordX + 5, dOrgCoordY - dLenX + 5);



	//画比例尺
	//比例尺的起点
	double dOrgScaleX = dOrgCoordX + dLenX + double(rect.Width()) / 10;
	double dOrgScaleY = dOrgCoordY;
	pDC->MoveTo(dOrgScaleX, dOrgScaleY);
	pDC->LineTo(dOrgScaleX + 200 * dScale, dOrgScaleY);
	pDC->LineTo(dOrgScaleX + 200 * dScale, dOrgScaleY - 2);
	pDC->MoveTo(dOrgScaleX, dOrgScaleY);
	pDC->LineTo(dOrgScaleX, dOrgScaleY - 2);


	pDC->SelectObject(pOldPen);
	pen.DeleteObject();




	CFont Font;//创建字体
	VERIFY(Font.CreateFont(
		10,                        // nHeight
		0,                         // nWidth
		0,                         // nEscapement
		0,                         // nOrientation*
		FW_NORMAL,                 // nWeight
		FALSE,                     // bItalic
		FALSE,                     // bUnderline
		0,                         // cStrikeOut
		ANSI_CHARSET,              // nCharSet
		OUT_DEFAULT_PRECIS,        // nOutPrecision
		CLIP_DEFAULT_PRECIS,       // nClipPrecision
		DEFAULT_QUALITY,           // nQuality
		DEFAULT_PITCH | FF_SWISS,  // nPitchAndFamily
		_T("宋体")));                 // lpszFacename

	//pDC->SetTextColor(RGB(0, 0, 0));///改变字体颜色

	CFont* pOldFont = pDC->SelectObject(&Font);

	//写点名
	POINT Point;
	for (int i = 0; i < m_KnownPointCount; i++)
	{
		Point = LP2CP(m_KnownPoint[i], rect);
		Point.x = Point.x + dOrgNetX;
		Point.y = Point.y + dOrgNetY;

		pDC->TextOutW(Point.x, Point.y + 10, m_KnownPoint[i].strID);
	}

	for (int i = 0; i < m_UnknownPointCount; i++)
	{
		Point = LP2CP(m_UnknownPoint[i], rect);
		Point.x = Point.x + dOrgNetX;
		Point.y = Point.y + dOrgNetY;

		pDC->TextOutW(Point.x, Point.y + 10, m_UnknownPoint[i].strID);
	}

	//写坐标轴名
	pDC->TextOutW(dOrgCoordX - 12, dOrgCoordY - dLenX - 2, _T("X"));
	pDC->TextOutW(dOrgCoordX + dLenY - 3, dOrgCoordY + 5, _T("Y"));

	pDC->SelectObject(pOldPen);
	Font.DeleteObject();


	//写比例尺名
	pDC->TextOutW(dOrgScaleX + 70 * dScale, dOrgScaleY + 6, _T("200m"));
}




//从文件中读取数据
void ControlNetworkAdjust::ReadFileData(CStdioFile& sf)
{
	//清空变量
	ClearVariables();
	
	CString strLine;
	CStringArray aStrTmp;

	//读取已知点
	sf.ReadString(strLine);
	m_KnownPointCount = _ttoi(strLine);
	m_KnownPoint = new CtrlPoint[m_KnownPointCount];
	for (int i = 0; i < m_KnownPointCount; i++)
	{
		sf.ReadString(strLine);
		SplitStringArray(strLine, ',', aStrTmp);
		m_KnownPoint[i].strID = aStrTmp[0];
		m_KnownPoint[i].dx = _tstof(aStrTmp[1]);
		m_KnownPoint[i].dy = _tstof(aStrTmp[2]);
		m_KnownPoint[i].bState = TRUE;
		m_KnownPoint[i].ePointStyle = KnownPoint;
	}

	//读取未知点
	sf.ReadString(strLine);
	m_UnknownPointCount = _ttoi(strLine);
	m_UnknownPoint = new CtrlPoint[m_UnknownPointCount];
	sf.ReadString(strLine);
	SplitStringArray(strLine, ',', aStrTmp);
	for (int i = 0; i < m_UnknownPointCount; i++)
	{
		m_UnknownPoint[i].strID = aStrTmp[i];
		m_UnknownPoint[i].dx = -1;
		m_UnknownPoint[i].dy = -1;//坐标x、y赋值为-1,方便后面判断是否已经进行概算
		m_UnknownPoint[i].bState = FALSE;
		m_UnknownPoint[i].iNum = i;
		m_UnknownPoint[i].ePointStyle = UnknownPoint;
		m_UnknownPoint[i].dE = -1;//长轴赋值为-1,方便后面判断是否已经平差计算
	}

	//读取边长观测值
	sf.ReadString(strLine);
	m_ObsDistCount = _ttoi(strLine);
	m_ObsDist = new ObsDist[m_ObsDistCount];
	for (int i = 0; i < m_ObsDistCount; i++)
	{
		sf.ReadString(strLine);
		SplitStringArray(strLine, ',', aStrTmp);
		m_ObsDist[i].StartPoint = SearchPoint(aStrTmp[0]);
		m_ObsDist[i].EndPoint = SearchPoint(aStrTmp[1]);
		m_ObsDist[i].DistValue = _tstof(aStrTmp[2]);
	}

	//读取角度观测值
	sf.ReadString(strLine);
	m_ObsAngleCount = _ttoi(strLine);
	m_ObsAngle = new ObsAngle[m_ObsAngleCount];
	for (int i = 0; i < m_ObsAngleCount; i++) 
	{
		sf.ReadString(strLine);
		SplitStringArray(strLine, ',', aStrTmp);
		m_ObsAngle[i].StationPoint = SearchPoint(aStrTmp[0]);
		m_ObsAngle[i].ObjPoint = SearchPoint(aStrTmp[1]);
		m_ObsAngle[i].AngleValue(DMS)= _tstof(aStrTmp[2]);
	}
}

//计算未知点近似坐标
bool ControlNetworkAdjust::CalUnknownPoint()
{
	if (m_UnknownPointCount == 0)
	{
		MessageBox(NULL, _T("未读取数据!"), _T("警告"), MB_OK);
		return FALSE;
	}
	else if (m_UnknownPoint[0].dx >= 0)
	{
		MessageBox(NULL, _T("已经完成坐标概算!"), _T("警告"), MB_OK);
		return FALSE;
	}
	else
	{
		do
		{
			//遍历未知点,计算未知点的近似坐标
			for (int i = 0; i < m_UnknownPointCount; i++)
			{
				//如果该未知点的坐标未求出
				if (m_UnknownPoint[i].bState == FALSE)
				{
					ObsAngle* LookedAngle2 = NULL;//测站点和照准点点均已求出的角度观测值
					ObsAngle* LookedAngle1 = NULL;//测站点坐标已求出且照准点为所求未知点的角度观测值

					//寻找终点为未知点且起点坐标已求出的边长观测值
					ObsDist* LookedDist = SearchObsDist_EndPtIsLookedPt(m_UnknownPoint[i].strID);
					
					//如果边长观测值存在
					if (LookedDist != NULL)
					{
						//寻找测站点为Dist起点且照准点坐标已求出的角度观测值
						LookedAngle2 = SearchObsAngle_StationPtIsLookedPt(LookedDist->StartPoint->strID);

						//寻找照准边为Dist所在边的角度观测值
						LookedAngle1 = SearchObsAngle(LookedDist->StartPoint->strID, LookedDist->EndPoint->strID);
					}
					else
					{
						//寻找起点为未知点且终点坐标已求出的边长观测值
						LookedDist = SearchObsDist_StartPtIsLookedPt(m_UnknownPoint[i].strID);

						if (LookedDist != NULL)
						{
							//寻找测站点为Dist起点且照准点坐标已求出的角度观测值
							LookedAngle2 = SearchObsAngle_StationPtIsLookedPt(LookedDist->EndPoint->strID);

							//寻找照准边为Dist所在边的角度观测值
							LookedAngle1 = SearchObsAngle(LookedDist->EndPoint->strID, LookedDist->StartPoint->strID);
						}
					}
						
					//如果两个角度观测值存在,证明该未知点坐标可求出
					if (LookedAngle1 != NULL && LookedAngle2 != NULL)
					{
						//计算LookedAngle1和LookedAngle2之间的夹角
						CAngle Angle = LookedAngle1->AngleValue - LookedAngle2->AngleValue;

						//计算LookedAngle2所在照准边的方位角
						CAngle LookedAngle2Azi = Azi(*LookedAngle2->StationPoint, *LookedAngle2->ObjPoint);

						//计算LookedAngle1所在照准边的方位角
						CAngle LookedAngle1Azi = LookedAngle2Azi + Angle;

						//如果求得LookedAngle1所在照准边的方位角大于等于360,减去360
						if (LookedAngle1Azi(DEG) > 360)
						{
							LookedAngle1Azi(DEG) = LookedAngle1Azi(DEG) - 360;
						}

						double dx, dy;//坐标增量

						//计算坐标增量
						dx = LookedDist->DistValue * cos(LookedAngle1Azi(RAD));
						dy = LookedDist->DistValue * sin(LookedAngle1Azi(RAD));

						//计算照准点的近似坐标
						LookedAngle1->ObjPoint->dx = LookedAngle1->StationPoint->dx + dx;
						LookedAngle1->ObjPoint->dy = LookedAngle1->StationPoint->dy + dy;
						LookedAngle1->ObjPoint->bState = TRUE;
					}
				}
			}
		} while (!JudgeUnknownPoint());

		return TRUE;
	}
}

//输出矩阵到文件
void ControlNetworkAdjust::PrintCMatrix(CMatrix& A, CString strFileName)
{
	CStdioFile sf;
	if (!sf.Open(strFileName, CFile::modeWrite | CFile::modeCreate)) return;

	CString strResult;
	CString strContent;

	for (int i = 0; i < A.Row(); i++)
	{
		for (int j = 0; j < A.Col(); j++)
		{
			strResult.Format(_T("%-15.4f"), A(i, j));
			strContent += strResult;
		}
		strContent += _T("\n");
	}

	sf.WriteString(strContent);

	sf.Close();
}

//向文本框中输出概算坐标
void ControlNetworkAdjust::PrintRoughCalResult(CString& strOut)
{
	strOut.Empty();

	CString strResult;

	strOut += (_T("*****************坐标概算*****************\r\n\r\n"));

	strResult.Format(_T("%-20s%-22s%-22s\r\n\r\n"), _T("点号"), _T("X坐标"), _T("Y坐标"));
	strOut += strResult;

	//输出已知点坐标
	strOut += (_T("*****************已知点*******************\r\n"));
	for (int i = 0; i < m_KnownPointCount; i++)
	{
		strResult.Format(_T("%-15s%-20.6f%-20.6f\r\n"), m_KnownPoint[i].strID, m_KnownPoint[i].dx, m_KnownPoint[i].dy);
		strOut += strResult;
	}

	//输出未知点坐标
	strOut += (_T("\r\n*****************未知点*******************\r\n"));
	for (int i = 0; i < m_UnknownPointCount; i++)
	{
		strResult.Format(_T("%-15s%-20.6f%-20.6f\r\n"), m_UnknownPoint[i].strID, m_UnknownPoint[i].dx, m_UnknownPoint[i].dy);
		strOut += strResult;
	}

	return;
}

//向文本框中输出平差坐标,并存入文件中
void ControlNetworkAdjust::PrintAdjustCalResult(CString& strOut)
{
	CString strFileData;
	CString strResult;

	strFileData += (_T("******************************平差计算************************************\n"));

	//输出方向观测成果
	strFileData += (_T("\n******************************方向观测成果********************************\n"));
	strResult.Format(_T("%-13s%-13s%-13s%-10s%-17s%-10s\n"), _T("测站"), _T("照准"), _T("方向值(dms)"), _T("改正数(s)"),
		_T("平差后值(dms)"), _T("备注"));
	strFileData += strResult;
	for (int i = 0; i < m_ObsAngleCount; i++)
	{
		//如果是零方向
		if (m_ObsAngle[i].AngleValue(DEG) == 0)
		{
			strResult.Format(_T("%-15s%-15s%-18.6f\n"), m_ObsAngle[i].StationPoint->strID,
				m_ObsAngle[i].ObjPoint->strID, m_ObsAngle[i].AngleValue(DMS));
			strFileData += strResult;
		}
		else//如果不是零方向
		{
			strResult.Format(_T("%-15s%-15s%-18.6f%-15.2f%-18.6f\n"), m_ObsAngle[i].StationPoint->strID, m_ObsAngle[i].ObjPoint->strID,
				m_ObsAngle[i].AngleValue(DMS) - m_ObsAngle[i].AngleError / 10000, m_ObsAngle[i].AngleError, m_ObsAngle[i].AngleValue(DMS));
			strFileData += strResult;
		}
	}

	//输出距离观测成果
	strFileData += (_T("\n******************************距离观测成果********************************\n"));
	strResult.Format(_T("%-13s%-15s%-14s%-10s%-13s%-20s\n"), _T("测站"), _T("照准"), _T("距离(m)"), _T("改正数(m)"), _T("平差后值(m)"),
		_T("方位角(dms)"));
	strFileData += strResult;
	for (int i = 0; i < m_ObsDistCount; i++)
	{
		CAngle DistAzi = Azi(*m_ObsDist[i].StartPoint, *m_ObsDist[i].EndPoint);//计算方位角
		strResult.Format(_T("%-15s%-15s%-18.4f%-15.4f%-18.4f%-18.6f\n"), m_ObsDist[i].StartPoint->strID, m_ObsDist[i].EndPoint->strID,
			m_ObsDist[i].DistValue - m_ObsDist[i].DistError, m_ObsDist[i].DistError, m_ObsDist[i].DistValue, DistAzi(DMS));
		strFileData += strResult;
	}

	//输出平面点位误差
	strFileData += (_T("\n******************************平面点位误差********************************\n"));
	strResult.Format(_T("%-13s%-12s%-10s%-15s%-16s%-20s\n"), _T("点名"), _T("长轴(m)"), _T("短轴(m)"), _T("长轴方位角(dms)"),
		_T("点位中误差(m)"), _T("备注"));
	strFileData += strResult;
	for (int i = 0; i < m_UnknownPointCount; i++)
	{
		strResult.Format(_T("%-15s%-15.4f%-15.4f%-23.6f%-18.4f\n"), m_UnknownPoint[i].strID, m_UnknownPoint[i].dE,
			m_UnknownPoint[i].dF, m_UnknownPoint[i].dAlfa(DMS), m_UnknownPoint[i].dMk);
		strFileData += strResult;
	}

	//输出控制点成果
	strFileData += (_T("\n******************************控制点成果*********************************\n"));
	strResult.Format(_T("%-19s%-20s%-20s%-13s%-20s\n"), _T("点名"), _T("X(m)"), _T("Y(m)"), _T("H(m)"), _T("备注"));
	strFileData += strResult;
	for (int i = 0; i < m_KnownPointCount; i++)
	{
		strResult.Format(_T("%-15s%-18.4f%-18.4f%-17s%-20s\n"), m_KnownPoint[i].strID, m_KnownPoint[i].dx, m_KnownPoint[i].dy, _T(" "), _T("未知点"));
		strFileData += strResult;
	}
	for (int i = 0; i < m_UnknownPointCount; i++)
	{
		strResult.Format(_T("%-15s%-18.4f%-18.4f%-17s%-20s\n"), m_UnknownPoint[i].strID, m_UnknownPoint[i].dx, m_UnknownPoint[i].dy, _T(" "), _T(" "));
		strFileData += strResult;
	}

	//输出到文件中
	CString strFileName = _T("AdjustReult.dat");
	CStdioFile sf;
	if (!sf.Open(strFileName, CFile::modeReadWrite | CFile::modeCreate)) return;

	sf.WriteString(strFileData);
	MessageBox(NULL, _T("平差结果已保存到“AdjustReult.dat”文件中!"), _T("提示"), MB_OK);

	sf.SeekToBegin();//文件指针定位到开头

	ReadFileContent(sf,strOut);//从文件读取数据

	sf.Close();
}


//判断是否达到绘图条件
bool ControlNetworkAdjust::JudgeDrawConditon()
{
	if (m_UnknownPointCount == 0)
	{
		MessageBox(NULL, _T("未读取数据!"), _T("警告"), MB_OK);
		return FALSE;
	}
	else if (m_UnknownPoint[0].dx < 0)
	{
		MessageBox(NULL, _T("未坐标概算!"), _T("警告"), MB_OK);
		return FALSE;
	}
	else if (m_UnknownPoint[0].dE < 0)
	{
		MessageBox(NULL, _T("未平差计算!"), _T("警告"), MB_OK);
		return FALSE;
	}
	else
	{
		return TRUE;
	}
}

//平差计算
bool ControlNetworkAdjust::AdjustCal()
{
	if (m_UnknownPointCount == 0)
	{
		MessageBox(NULL, _T("未读取数据!"), _T("警告"), MB_OK);
		return FALSE;
	}
	else if (m_UnknownPoint[0].dx < 0)
	{
		MessageBox(NULL, _T("未坐标概算!"), _T("警告"), MB_OK);
		return FALSE;
	}
	else if (m_UnknownPoint[0].dE >= 0)
	{
		MessageBox(NULL, _T("已经完成平差计算!"), _T("警告"), MB_OK);
		return FALSE;
	}
	else
	{
		//定义误差方程和权阵
		CMatrix X(m_UnknownPointCount * 3 + m_KnownPointCount, 1);
		CMatrix B(m_ObsAngleCount + m_ObsDistCount, m_UnknownPointCount * 3 + m_KnownPointCount);
		CMatrix L(m_ObsAngleCount + m_ObsDistCount, 1);
		CMatrix P(m_ObsAngleCount + m_ObsDistCount, m_ObsAngleCount + m_ObsDistCount);
		CMatrix V(m_ObsAngleCount + m_ObsDistCount, 1);


		CMatrix VT;
		CMatrix BT;//B的转置矩阵
		CMatrix NBB;
		CMatrix NBB1;//NBB的逆矩阵

		//平差计算直到满足精度要求
		double dXmax;//X中坐标x或y增量的最大值
		do
		{
			//计算B、L、P
			CalBLP(B, L, P);

			BT = ~B;//求B的转置矩阵
			NBB = BT * P * B;//求NBB
			NBB1 = NBB.Inv();//求NBB的逆矩阵
			X = -1 * NBB1 * BT * P * L;//求X
			V = B * X + L;//计算改正数V

			//求未知点坐标的近似平差值
			for (int i = 0; i < m_UnknownPointCount; i++)
			{
				m_UnknownPoint[i].dx += X(2 * i, 0) / 1000;
				m_UnknownPoint[i].dy += X(2 * i + 1, 0) / 1000;
			}

			//求dXmax
			dXmax = fabs(X(0, 0));
			for (int i = 0; i < 2 * m_UnknownPointCount; i++)
			{
				dXmax = fabs(X(i, 0)) > dXmax ? fabs(X(i, 0)) : dXmax;
			}
		} while (dXmax > 0.1);

		//将零方向的角度改正归算到该测站其他照准边上
		double V0;//记录零方向的角度改正数
		for (int i = 0; i < m_ObsAngleCount; i++)
		{
			//如果是零方向
			if (m_ObsAngle[i].AngleValue(DEG) == 0)
			{
				V0 = V(i, 0);
				V(i, 0) = 0;//零方向改正数化为0

				//如果是该测站上其他照准方向
				for (int j = i + 1; j < m_ObsAngleCount; j++)
				{
					if (m_ObsAngle[j].AngleValue(DEG) != 0)
					{
						V(j, 0) = V(j, 0) - V0;
					}
					else
					{
						break;
					}
				}
			}
		}

		//计算边长和角度观测值的改正数以及改正后的值
		CalAdjustedAngleDist(V);

		//计算单位权中误差
		VT = ~V;
		CMatrix VV = VT * P * V;
		double m0 = sqrt(VV(0, 0) / double(V.Row() - X.Row())) / 1000;

		//精度评定,计算点位误差和误差椭圆参数
		PrecisionEvaluate(NBB1, m0);

		return TRUE;
	}
}

4.6、画图对话框类头文件(ControlNetworkDrawingDlg.h)的代码:

#pragma once


// ControlNetworkDrawingDlg 对话框

class ControlNetworkDrawingDlg : public CDialogEx
{
	DECLARE_DYNAMIC(ControlNetworkDrawingDlg)

public:
	ControlNetworkDrawingDlg(CWnd* pParent = nullptr);   // 标准构造函数
	virtual ~ControlNetworkDrawingDlg();

// 对话框数据
#ifdef AFX_DESIGN_TIME
	enum { IDD = IDD_CONTROLNETWORKDRAWING_DIALOG };
#endif

protected:
	virtual void DoDataExchange(CDataExchange* pDX);    // DDX/DDV 支持

	DECLARE_MESSAGE_MAP()
public:
	afx_msg void OnPaint();

//实现
public:
	ControlNetworkAdjust* CNA_Draw;//控制网平差类对象
};

4.7、画图对话框类cpp文件(ControlNetworkDrawingDlg.cpp)的主要代码:

void ControlNetworkDrawingDlg::OnPaint()
{
	CPaintDC dc(this); // device context for painting
					   // TODO: 在此处添加消息处理程序代码
					   // 不为绘图消息调用 CDialogEx::OnPaint()

	CWnd* pWin = GetDlgItem(IDC_CTRLNET);//获取Picture控件的指针
	CRect rect;
	pWin->GetClientRect(rect);//把控件的长宽、坐标等信息保存在rect里
	CDC* pDC = pWin->GetDC();//获取该控件的画布

	CNA_Draw->ControlNetworkDrawing(pDC, rect);//画图
}

4.8、角度类头文件(Angle.h)的主要代码:

#pragma once

//枚举数据类型,用于代表角度形式
enum AngleStyle
{
	DEG,//十进制度
	DMS,//度分秒
	RAD//弧度制
};

/*******************************************************************
* 类名:CAngle													   *
*																   *
* 描述:角度类													   *
*																   *
* 使用方法:1、初始化:CAngle(double value=0,AngleStyle style=DMS) *
*			即默认值为0,角度为度分秒,如果value为a,AngleStyle为  *
*			Deg或RAD,则类对象的值便为a,类型为度或弧度。          *
*			2、获取指定类型的角度值的函数有两个					   *
*			第一个为double& operator() (AngleStyle style),使用该  *
*			函数,返回引用类对象,类对象的值和类型改变(指定类型) *
*			第二个为double operator() (AngleStyle style) const,   *
*			const CAngle类型变量调用,该函数返回的类对象不能引用, *
*			且类对象的值和类型不变(原来类型)					   * 
*			3、角度加减,即CAngle operator + (const CAngle& m1,	   *
			const CAngle& m2),CAngle operator - (const CAngle& m1,*
			const CAngle& m2),该函数返回的类对象类型为RAD。       *
*																   *
* 历史:**日期** **理由** **签名**								   *
*       2021/7/6    无     苑艺博								   *
*																   *
* 外部类:无													   *
*******************************************************************/
class CAngle
{
public:
	CAngle(double value=0,AngleStyle style=DMS);
	~CAngle(void);
private:
	double dValue;//角度值
	AngleStyle  nCurStyle;//当前角度值类型
private:
	//设置常成员函数的作用:1.类成员不会被改变
	//2.可以被常类对象调用
	double Deg(double dDms) const;
	double Dms(double dDeg) const;
	
public:
	//获取指定的类型获取角度值,
    //由于返回的是dValue的引用,所以该值大小可以改变,即可以进行赋值
	double& operator() (AngleStyle style);   //CAanle ang1,   ang1.operator(DEG); angl1(DEG)

	//重载,获取指定的类型获取角度值,该值不可改变,const CAngle类型变量调用
	double operator() (AngleStyle style) const;
	//重载运算符+/-
    friend CAngle operator + (const CAngle& m1,const CAngle& m2);
	friend CAngle operator - (const CAngle& m1,const CAngle& m2);
};

4.9、角度类cpp文件(Angle.cpp)的主要代码:

#include "pch.h"
#include "Angle.h"
#include "math.h"

//重载构造函数,有缺省值
CAngle::CAngle(double value,AngleStyle style)
{
	dValue=value;
	nCurStyle=style;
}

CAngle::~CAngle(void)
{
}
//重载()函数
double& CAngle::operator() (AngleStyle style) //指定的类型获取角度值
{
	//double dAngleValue;
    if(style==DMS)
	{
		if(nCurStyle==DEG)
		{
            dValue=Dms(dValue);
		}
		else if(nCurStyle==RAD)
		{
			dValue=Dms(dValue*180.0/PI);
		}
		nCurStyle=DMS;
        		
	}
	else if(style==DEG)
	{
	    if(nCurStyle==DMS)
		{
            dValue=Deg(dValue);
		}
		else if(nCurStyle==RAD)
		{
			dValue=dValue*180.0/PI;
		}
       nCurStyle=DEG;
	}
	else
	{
       if(nCurStyle==DMS)
		{
            dValue=Deg(dValue)*PI/180;
		}
		else if(nCurStyle==DEG)
		{
			dValue=dValue*PI/180;
		}
        nCurStyle=RAD;
	}
	return dValue;
}
//重载()函数,该函数是常函数,只能被常CAngle对象使用
double CAngle::operator() (AngleStyle style) const //指定的类型获取角度值
{
	double dAngleValue;
    if(style==DMS)
	{
		if(nCurStyle==DEG)
		{
            dAngleValue=Dms(dValue);
		}
		else if(nCurStyle==RAD)
		{
			dAngleValue=Dms(dValue*180.0/PI);
		}
		else
		{
			dAngleValue=dValue;
		}
        		
	}
	else if(style==DEG)
	{
	    if(nCurStyle==DMS)
		{
            dAngleValue=Deg(dValue);
		}
		else if(nCurStyle==RAD)
		{
			dAngleValue=dValue*180.0/PI;
		}
       else
		{
			dAngleValue=dValue;
		}
	}
	else
	{
       if(nCurStyle==DMS)
		{
            dAngleValue=Deg(dValue)*PI/180;
		}
		else if(nCurStyle==DEG)
		{
			dAngleValue=dValue*PI/180;
		}
        else
		{
			dAngleValue=dValue;
		}
	}
	return dAngleValue;
}


//私有成员,度分秒向十进制度转换
double CAngle::Deg(double dDms) const
{
	int iDeg,iMin;
	double dSec;

	iDeg = int(dDms + EPSILON);//度//加一个很小的数,以防止取整时的出错
	iMin = int((dDms - iDeg) * 100+ EPSILON);//分
	dSec = ((dDms - iDeg) * 100 - iMin) * 100 ;//秒

	return iDeg + (double)iMin / 60 + dSec / 3600;
}

//私有成员,十进制度向度分秒转换
double CAngle::Dms(double dDeg) const
{
	int iDeg,iMin;
	double dSec;
	double dTmp;

	iDeg = int(dDeg + EPSILON);//整数部分度
	dTmp = (dDeg - iDeg) * 60;//小数部分转换成分
	iMin = int(dTmp+ EPSILON);//取分的整数部分
	dSec = (dTmp - iMin) * 60;//截取秒

	return iDeg + (double)iMin / 100 + dSec / 10000;
}

//友元重载+函数
CAngle operator + (const CAngle& m1,const CAngle& m2)
{
       CAngle addAngle(0,RAD);
	   addAngle(RAD)=m1(RAD)+m2(RAD);
	   return addAngle;
}
//友元重载-函数
CAngle operator - (const CAngle& m1,const CAngle& m2)
{
       CAngle subAngle(0,RAD);
	   subAngle(RAD)=m1(RAD)-m2(RAD);
	   return subAngle;
}

4.10、矩阵类头文件(Matrix.h)的主要代码:

#pragma once

class CMatrix
{
public:
	CMatrix(int row=3,int col=3);
	// copy constructor

    CMatrix (const CMatrix& m);

	~CMatrix(void);
private:
	double **dMatData;//保存矩阵元素数据的二维数组
	int iRow;//矩阵的行
	int iCol;//矩阵的列
    
public:
	int Row() const {return iRow;}//返回行
	int Col() const {return iCol;}//返回列
    void SetSize (int row, int col);//调整数组的大小,原有数据不变(未测试)
	void setzeros(CMatrix& B);//矩阵置零

	double& operator () (int row, int col);//获取矩阵元素
    double  operator () (int row, int col) const;//重载获取矩阵元素函数,只有const对象能访问
	CMatrix& operator = (const CMatrix& m) ;
    
    //注意:友元函数并不是类自己的成员函数
    friend CMatrix operator + (const CMatrix& m1,const CMatrix& m2);
	friend CMatrix operator - (const CMatrix& m1,const CMatrix& m2);
	friend CMatrix operator * (const CMatrix& m1,const CMatrix& m2);
	friend CMatrix operator * (const double& num, const CMatrix& m1);
	friend CMatrix operator * (const CMatrix& m1,const double& num);

	friend CMatrix operator ~ (const CMatrix& m);//矩阵转置
	CMatrix Inv();//矩阵求逆
	void Unit();//生成单位矩阵
};

4.11、矩阵类cpp文件(Matrix.cpp)的主要代码:

#include "pch.h"
#include "pch.h"
#include "Matrix.h"
#include "math.h"

/*******************************************************************
* 类名:CMatrix												       *
*																   *
* 描述:矩阵类													   *
*																   *
* 使用方法:													   *
*																   *
* 历史:**日期** **理由** **签名**								   *
*       2021/7/8    无     苑艺博								   *
*																   *
* 外部类:无													   *
*******************************************************************/
CMatrix::CMatrix(int row,int col)
{
	iRow=row;
	iCol=col;
    dMatData = new double*[row];

	for (int i=0; i < row; i++)
	{
		dMatData[i]= new double[col];
		for(int j=0;j<col;j++)
		{
				dMatData[i][j]=0;
		}
	 }
}

// copy constructor,
//拷贝构造函数的作用:
//(1)以类对象作为函数参数传值调用时;
//(2)函数返回值为类对象;
//(3)用一个已定义的对象去初始化一个新对象时;

CMatrix::CMatrix (const CMatrix& m)
{
   	iRow=m.Row();
	iCol=m.Col();
    dMatData = new double*[iRow];

	for (int i=0; i < iRow; i++)
	{
		dMatData[i]= new double[iCol];
	//	for(int j=0;j<iCol;j++)
		{
			memcpy(dMatData[i],m.dMatData[i],sizeof(double)*iCol);
		}
	 }
   
}

CMatrix::~CMatrix(void)
{
    for (int i=0; i < iRow; i++)
	{
		delete[] dMatData[i];
	 }
	delete[] dMatData;
}

//返回数组元素(引用返回)
double& CMatrix::operator () (int row, int col)
{
    if (row >= iRow || col >= iCol)
	{
      throw( "CMatrix::operator(): Index out of range!");
	}
	
    return dMatData[row][col]; 
}

返回数组元素(重载)
double CMatrix::operator () (int row, int col) const
{
    if (row >= iRow || col >= iCol)
	{
      throw( "CMatrix::operator(): Index out of range!");
	}
	
    return dMatData[row][col]; 
}


//重载预算符+
CMatrix operator + (const CMatrix& m1,const CMatrix& m2)
{

   if((m1.Col()!=m2.Col()) ||(m1.Row()!=m2.Row()) )
   {
       throw( "CMatrix::operator+: The two matrix have different size!");
   }

   CMatrix matTmp(m1.Row(),m1.Col());
   for(int i=0;i<m1.Row();i++)
   {
	   for(int j=0;j<m1.Col();j++)
	   {
             matTmp(i,j)=m1(i,j)+m2(i,j);     
	   }
   }
   return matTmp;
}

//重载赋值运算符=,当左右两边矩阵的大小不相等时,
//以右边的大小为基准,调整左边矩阵的大小

CMatrix &CMatrix::operator = (const CMatrix& m) 
{
	//revised in 2011-4-1, by Daiwujiao
 //   if(iRow!=m.Row()||iCol!=m.Col())
	//{
 //       throw( "CMatrix::operator=: The two matrix have different size!");
	//}

	if(iRow!=m.Row()||iCol!=m.Col())
	{
		SetSize(m.Row(),m.Col());
	}
	for (int i=0; i < iRow; i++)
	{
		
		for(int j=0;j<iCol;j++)
		{
				dMatData[i][j]=m(i,j);
		}
	 }
    return *this;
}


//调整矩阵大小,原有值不变
void CMatrix::SetSize (int row, int col)
{
   if (row == iRow && col == iCol)
   {
      return;
   }

   double **rsData = new double*[row];
   	for (int i=0; i < row; i++)
	{
		rsData[i]= new double[col];
		for(int j=0;j<col;j++)
		{
				rsData[i][j]=0;
		}
	 }

	int minRow=(iRow>row)?row:iRow;
    int minCol= (iCol>col)?col:iCol;
    int  colSize = minCol * sizeof(double);
    

   for (int i=0; i < minRow; i++)
   {
      memcpy( rsData[i], dMatData[i], colSize);
   }

    for (int i=0; i < minRow; i++)
	{
         delete[] dMatData[i];
	}
	delete[] dMatData;
	dMatData=rsData;
    iRow=row;
	iCol=col;
    return;
}
//重载预算符-
CMatrix operator - (const CMatrix& m1,const CMatrix& m2)
{

   if((m1.Col()!=m2.Col()) ||(m1.Row()!=m2.Row()) )
   {
       throw( "CMatrix::operator-: The two matrix have different size!");
   }

   CMatrix matTmp(m1.Row(),m1.Col());
   for(int i=0;i<m1.Row();i++)
   {
	   for(int j=0;j<m1.Col();j++)
	   {
             matTmp(i,j)=m1(i,j)-m2(i,j);     
	   }
   }
   return matTmp;
}

//重载预算符*,两个矩阵相乘,m1的列要等于m2的行
CMatrix operator * (const CMatrix& m1,const CMatrix& m2)
{

   if((m1.Col()!=m2.Row()))
   {
       throw( "CMatrix::operator*: The col of matrix m1 doesn't equ to row of m2 !");
   }

   CMatrix matTmp(m1.Row(),m2.Col());
   for(int i=0;i<m1.Row();i++)
   {
	   for(int j=0;j<m2.Col();j++)
	   {
		   for(int k=0;k<m2.Row();k++)
		   {
             matTmp(i,j)+=m1(i,k)*m2(k,j);     
		   }
	   }
   }
   return matTmp;
}

//重载预算符*,矩阵右乘一个数
CMatrix operator * (const CMatrix& m1,const double& num)
{
   CMatrix matTmp(m1.Row(),m1.Col());
   for(int i=0;i<m1.Row();i++)
   {
	   for(int j=0;j<m1.Col();j++)
	   {
             matTmp(i,j)=m1(i,j)*num;     

	   }
   }
   return matTmp;
}

//重载预算符*,矩阵左乘一个数
CMatrix operator * (const double& num, const CMatrix& m1)
{
   CMatrix matTmp(m1.Row(),m1.Col());
   for(int i=0;i<m1.Row();i++)
   {
	   for(int j=0;j<m1.Col();j++)
	   {
             matTmp(i,j)=m1(i,j)*num;     
	   }
   }
   return matTmp;
}

//矩阵转置
CMatrix operator ~ (const CMatrix& m)
{
  CMatrix matTmp(m.Col(),m.Row());

   for (int i=0; i < m.Row(); i++)
      for (int j=0; j < m.Col(); j++)
      {
         matTmp(j,i) = m(i,j);
      }
   return matTmp;
}


//矩阵求逆
//采用选全主元法
CMatrix CMatrix::Inv()
{
    if (iRow!=iCol)
	{
        throw("待求逆的矩阵行列不相等!");
	}
    
    int i, j, k, vv;
 
    CMatrix InvMat(iRow,iRow);

    //复制矩阵
	InvMat=*this;
   
    int* MainRow=new int[iRow];
	int* MainCol=new int[iRow];//用于记录主元素的行和列

    double dMainCell;//主元元素的值
    double dTemp;//临时变量

    for(k = 0;k<iRow;k++)
	{
        dMainCell = 0;
        //选全主元
        for( i = k;i<iRow ;i++)
		{
            for( j = k;j<iRow;j++)
			{
                dTemp = fabs(InvMat(i, j));
                if(dTemp > dMainCell)
				{
                    dMainCell = dTemp;
                    MainRow[k] = i;
                    MainCol[k] = j;
				}
			}
		}

		if( fabs(dMainCell) < 0.0000000000001)//矩阵秩亏,不能求逆
		{
            throw("矩阵秩亏");
		}

        if(MainRow[k] != k)//交换行
		{
            for( j = 0 ;j<iRow;j++)
			{
                vv = MainRow[k];
                dTemp = InvMat(k, j);
                InvMat(k, j) = InvMat(vv, j);
                InvMat(vv, j) = dTemp;
			}
		}

        if(MainCol[k] != k)//交换列
		{
            for(i = 0;i<iRow;i++)
			{
                vv = MainCol[k];
                dTemp = InvMat(i, k);
                InvMat(i, k) = InvMat(i, vv);
                InvMat(i, vv) = dTemp;
			}
		}
        InvMat(k, k) = 1.0 / InvMat(k, k);//计算乘数

        for( j = 0;j< iRow;j++) //计算主行
		{
            if(j != k)
			{
                InvMat(k, j) = InvMat(k, j) * InvMat(k, k);
			}
		}
        for(i = 0;i<iRow;i++)//消元
		{
            if( i !=k)
			{
                for(j = 0;j<iRow;j++)
				{
					if(j != k)
					{
                        InvMat(i, j) -= InvMat(i, k) * InvMat(k, j);
					}
				}
			}
		}
        for( i = 0;i< iRow;i++ )//计算主列
		{
            if( i != k)
			{
				InvMat(i, k) = -InvMat(i, k) * InvMat(k, k);
			}
		}
	}

    for( k = iRow - 1;k>=0;k--)
	{
        if(MainCol[k] != k)// 交换行
		{
            for( j = 0;j<iRow;j++)
			{
                vv = MainCol[k];
                dTemp = InvMat(k, j);
                InvMat(k, j) = InvMat(vv, j);
                InvMat(vv, j) = dTemp;
			}
		}

        if(MainRow[k] != k)//交换列
		{
            for( i = 0;i<iRow;i++)
			{
				vv = MainRow[k];
                dTemp = InvMat(i, k);
                InvMat(i, k) = InvMat(i, vv);
                InvMat(i, vv) = dTemp;
			}
		}
	}
	delete[] MainRow;
	delete[] MainCol;
    return InvMat;
}
//单位化矩阵
void CMatrix::Unit()
{
     for(int i=0;i<iRow;i++)
	 {
		 for(int j=0;j<iCol;j++)
		 {
			 dMatData[i][j]=(i==j)?1:0;
		 }
	 }
}
//矩阵置零
void CMatrix::setzeros(CMatrix&B) 
{
	for(int i=0;i<B.iCol;i++)
		for (int j = 0; j < B.iCol; j++)
		{
			B(i,j)= 0;
		}
}

4.12、常用功能类头文件(CommonSurveyFunctions.h)的主要代码:

//常用测量计算函数
#pragma once

// 求平面上两点间距离
double Dist(double X1, double Y1,double X2, double Y2);

//重载,求空间上两点间距离
double Dist(double X1, double Y1, double Z1, 
			double X2, double Y2, double Z2);

//求两点的方位角
double  Azimuth(const double &X1, const double &Y1,
				const double &X2,const double &Y2);
//符号函数
int sgn(double x);

4.13、常用功能类cpp文件(CommonSurveyFunctions.cpp)的主要代码:

#include "pch.h"//预编译头
#include "CommonSurveyFunctions.h"
#include "math.h"
#include "Angle.h"


// 求平面上两点间距离
double Dist(double X1, double Y1,double X2, double Y2)
{
	double d;
	d=sqrt((X2-X1)*(X2-X1)+(Y2-Y1)*(Y2-Y1));
	return d;
}

//重载,求空间上两点间距离
double Dist(double X1, double Y1, double Z1, 
			double X2, double Y2, double Z2)
{
	double d;
	d=sqrt((X2-X1)*(X2-X1)+(Y2-Y1)*(Y2-Y1)+(Z2-Z1)*(Z2-Z1));
	return d;
}

//求两点的方位角
double  Azimuth(const double &X1, const double &Y1,
				const double &X2,const double &Y2)
{
		double dx, dy;
        dx = X2 - X1;
        dy = Y2 - Y1 + EPSILON;
        return PI - sgn(dy) * PI / 2 - atan(dx / dy);
}
//符号函数
int sgn(double x)
{
    if (x>=0) 
		return(1);//非负则返回1
	else
		return(-1);
}

5、 总结

(1)能力考查:本次课程设计主要考查了我们对控制网平差知识原理的掌握程度以及编程能力,而控制网平差所要求的高精度也需要我们具有缜密的心思和通过调试程序来解决问题的能力,而且,这次课程设计任务量较大、用时较长,更需要我们对软件和系统有整体设计和详细安排。总之,本次课程设计综合提升了我们测绘程序设计的能力。
(2)问题总结:下面对于程序不同模块的相关问题进行总结。
1、 数据读取:这部分较为简单,但是需要解决一个问题,即对一个文件既要读取数据到事先设置好的数组中,还要输出到对话框中。对此,我先用文件读取函数将文本内容读取到CString变量中,这样便实现了输出到对话框,接着将文本指针重新定位到开头,再逐行读取数据到数组中。这样只需通过文件对话框选择一次文件便可同时满足两种需求。
2、 坐标概算:这部分可采用多种方法,但是每种方法得到的坐标精度可能不同。我尝试了遍历角度、遍历边长和遍历未知点等多种方法,最终得到了一种算法较为简单且精度较高的坐标概算算法。而且,不同的坐标概算算法对数据的规范程度要求不同,如果角度和距离观测值非必要部分缺少,设计的算法也要能够成功进行坐标概算。
3、 平差计算:这部分的难点在于B、L、P矩阵的计算,这需要对相关公式的熟练掌握。同时,这部分也需要解决一个问题,即如何在把平差结果输出到对话框的同时保存到文件中,因为保存到文件时换行符只需\n,而输出到对话框的CString变量时换行符需要\r\n。为了解决这个矛盾,我先将平差结果格式化后保存到文件中,然后将文件指针定位到开头,再通过文件读取函数把存入文件中的平差结果文本输出到对话框中。这样便同时实现了两种功能。
4、 图形绘制:这部分包括三角网形绘制、误差椭圆绘制以及坐标轴和比例尺等元素的绘制,同时还要求各部分的比例大小和相对位置等适中,也要尽可能地满足处理不同数据均能够绘制合适的图形,对画图的程序设计和编程能力要求较高。
5、 用户体验:首先对程序的健壮性和容错性等要求较高,需要应对用户的多种不规范使用行为,比如不按照顺序点击按钮、重复点击按钮以及重复读入数据等等。同时,要求程序界面设计相对简洁美观。而且,要求程序使用方法简单,在不丧失程序必要功能的前提下,尽可能降低程序使用的复杂度。这些要求我自认为处理得不错,基本满足要求。
(3)存在问题:在定义X矩阵的大小时,行数定义为(3 * 未知点数 + 已知点数),其实应为(2 * 未知点数 + 角度测站数),忽略了数据不规范的情况下(未知点数 + 已知点数 != 角度测站数)的情况。
(4)实习收获:本次课程设计实习,让我对控制网平差相关原理和知识有了更深的掌握,锻炼了我的算法设计、处理bug和代码编写等编程能力,让我认识到测绘程序软件设计不同于其他一般软件设计的特点,即测绘程序设计关键在于对数据的处理,具有较高的精度要求,也让我对测绘程序设计这门学科有了更深的理解和热爱。

  • 16
    点赞
  • 67
    收藏
    觉得还不错? 一键收藏
  • 6
    评论
评论 6
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值