水准网平差c++代码

实验一 matlab完成水准网平差实验数据：水准网有2个已知点，3个未知点，7个测段。

已知点高程H1=5.016M H2=6.016h1=1.359; h2=2.009; h3=0.363; h4=1.012;h5=0.657; h6=0.238; h7=-0.595;S1=1.1 S2=1.7 S3=2.3 S4=2.7 S5=2.4 S6=1.4 S7=2.6求解（1）求个待定点高程,H5的高差中误差；3、4号点的高程中误差。

课程设计内容1、平差程序设计思路：使用间接平差法求解（1） 由题意知必要观测数t=3，选取3、4、5号点高程X1、X2、X3为参数。

（2） 误差方程：V1=x1v2=x2v3=x1v4=x2v5=x2-x1+h2-h1-h5v6=x3-x1v7=-x3(3)取1M 的观测高程为单位权观测，即 p=1/s；(4)求法方程：Nbbx-W=0 Nbb=b’pbW=b’pl(5)求的平差值x=Nbb^-1*W L=l+V V=bx-l（6）高差权函数式： k=-x1+x2(6)求中误差： 单位权中误差δ0，协因数阵Nbb^-1.求得中误差δ2、平差程序流程代码说明 ：h1=1.359;h2=2.009;h3=0.363;h4=1.012;h5=0.657;h6=0.238;h7=-0.595;H1=5.016H2=6.016h=[h1 h2 h3 h4 h5 h6 h7]'s=[1.1 1.7 2.3 2.7 2.4 1.4 2.6]' B=[1 0 0 ;0 1 0;1 0 0;0 1 0 ; -1 1 0 ;-1 0 1 ;0 0 -1 ]p=diag(1./s)l=[0;0;4;3;7;2;0]W=B'*p*lNbb=B'*p*Bx=inv(Nbb)*WV=(B*x-l)H=h+V/1000Q=inv(Nbb)n=7;t=3;j=V'*p*Vd= sqrt(j/4)f=[-1 1 0]'q=f'*Q*fD=d*sqrt(q)D1=d*sqrt(Q)(3) 平差程序流程代码说明:clccleardisp(‘观测高差，单位m’)h1=1.359;h2=2.009;h3=0.363;h4=1.012;h5=0.657;h6=0.238;h7=-0.595;H1=5.016 % 已知点高程，单位mH2=6.016 % 已知点高程，单位mh=[h1 h2 h3 h4 h5 h6 h7]'s=[1.1 1.7 2.3 2.7 2.4 1.4 2.6]' %S是线路长度disp(‘系数矩阵B、l’)B=[1 0 0 ;0 1 0;1 0 0;0 1 0 ; -1 1 0 ;-1 0 1 ;0 0 -1 ]p=diag(1./s) %定义权阵l=[0;0;4;3;7;2;0]W=B'*p*lNbb=B'*p*Bdisp(‘参数的解’)x=inv(Nbb)*WV=(B*x-l) % 误差方程(mm)H=h+V/1000 % 观测值的平差值Q=inv(Nbb) %观测值协因数阵n=7 %观测值数t=3 % 必要观测数j=V'*p*V %计算单位权中误差的参数d= sqrt(j/4) %单位权中误差f=[-1 1 0]' %权函数系数阵q=f'*Q*f %权函数协因数阵D=d*sqrt(q) %高差中误差D1=d*sqrt(Q) %高程中误差（4）计算结果： p =0.9091 0 0 0 0 0 00 0.5882 0 0 0 0 00 0 0.4348 0 0 0 00 0 0 0.3704 0 0 00 0 0 0 0.4167 00 0 0 0 0 0.7143 00 0 0 0 0 00.3846W =-2.60614.02781.4286Nbb =2.4748 -0.4167 -0.7143-0.4167 1.3753 0-0.7143 0 1.0989x =-0.24272.85521.1423V =-0.24272.8552-4.2427-0.1448-3.9021-0.6151-1.1423H =1.35882.01190.35881.01190.65310.2374-0.5961成果检核：H1+H5-H2=0；H3+H5-H4=0;H6+H7+H3=0即1.3588+0.6531-2.0119=0;0.3588+0.6531-1.0119=0；0.2374-0.5961+0.3588=0等式成立 检核通过0.5307 0.1608 0.3450 0.1608 0.7758 0.1045 0.3450 0.1045 1.1342 j = 19.7994d = 2.2248f = -11q =0.9850D = 2.2080D1 =1.6208 0.8921 1.30670.8921 1.9597 0.71921.3067 0.71922.3694。

常用测量程序设计代码（1）用全站仪在A点观测了B点斜边和垂直角，求A到B的高差。

（提示：22sin(1)cos2ABDh D a k a i vR=+-+-，D--斜边，a--垂直角，i--仪器高，v--反光镜高，k--大气折光系数）using System;using System.Collections.Generic;using System.Text;namespace ConsoleApplication1{class Application{staticvoid Main(string[] args){Console.Write("请输入斜边=");double D = double.Parse(Console.ReadLine());Console.Write("请输入垂直角[ddd.mmss]=");double a = DEG(double.Parse(Console.ReadLine()));Console.Write("请输入仪器高=");double i = double.Parse(Console.ReadLine());Console.Write("请输入反光镜高=");double v = double.Parse(Console.ReadLine());double h = D * Math.Sin(a) + (1 - 0.13) * D / 6371000.0 * D / 6371000.0 * Math.Cos(a) * Math.Cos(a) / 2.0 + i - v;Console.WriteLine("高差为{0}",h);}//将ddd.mmss转为弧度staticpublicdouble DEG(double ang){int fuhao = (int)(ang / Math.Abs(ang));ang = Math.Abs(ang);int d = (int)ang;int m = ((int)(ang * 100)) - d * 100;double s = ang * 10000 - m * 100 - d * 10000;return ((d + m / 60.0 + s / 3600.0) * fuhao) / 180.0 * Math.PI;}}}（2）如图所示，已知A点的坐标与A点到B点的边长与方位角，计算B点的坐标。

误差理论与测量平差础）课程设计报告系（部）：土木工程系实习单位：山东交通学院班级：测绘084学生姓名：田忠星学号080712420带队教师：夏小裕﹑周宝兴时间：10 年12 月13 日到10 年12 月19 日山东交通学院目录：1.摘要P32.概述P33.水准网间接平差程序设计思路P3—P44.平差程序流程图P4—P65.程序源代码及说明P7—P236.计算结果P23—P267.总结P26—P27一：摘要在测量工作中，为了能及时发现错误和提高测量成果的精度，常作多余观测，这就产生了平差问题。

在一个平差问题中，当所选的独立参数X?的个数等于必要观测数t 时，可将每个观测值表达成这t 个参数的函数，组成观测方程，这种以观测方程为函数模型的平差方法，就是间接平差。

二：概述：该课程设计的主要目是对水准网进行间接平差，在输入数据后依次计算高程近似值﹑误差方程和平差计算。

三：水准网间接平差程序设计思路1．根据平差问题的性质，选择t 个独立量（既未知点的高程）作为参数X?2.将每一个观测量的平差值（既观测的高程差值）分别表达成L L V3．由误差方程系数 B 和自由项组成法方程，法方程个数等于参数的个数t ；4. 解算法方程，求出参数X?,计算参数（高程）的平差值X?=X0 +x? ；5．由误差方程计算V，求出观测量（高差）平差值L L V 6.评定精度单位权中误差平差值函数的中误差四：平差程序流程图1． 已知数据的输入 需要输入的数据包括水准网中已知点数﹑未知点数以及这些点 的点号， 已知高程和高差观测值﹑距离观测值。

程序采用文件方 式进行输入，约定文件输入的格式如下： 第一行：已知点数﹑未知点数﹑观测值个数 第二行：点号（已知点在前，未知点在后） 第三行：已知高程（顺序与上一行的点号对应） 第四行：高差观测值，按“起点点号，终点点号。

高差观测值， 距离观测值”的顺序输入。

本节中使用的算例的数据格式如下2,3,7 1,2,3,4,5 5.016,6.016 1,3,1.359,1.11,4,2.009,1.7 2,3,0.363,2.3 2,4,1.012,2.7 3,4,0.657,2.4 3,5,0.238,1.4 5,2,-0.595,2.6 2.平差计算过程V TPV rV TPVnus（1）近似高程的计算。

《测绘程序设计（）》上机实验报告（Visual C++.Net）班级：学号：姓名：序号：二零一一年五月实验8 平差程序设计基础1.实验目的：1.1 巩固过程的定义与调用；1.2 巩固类的创建于使用；1.3 巩固间接平差模型与平差计算；1.4 掌握平差程序设计的基本技巧与步骤。

2.实验内容：水准网平差程序设计。

设计一个水准网平差的程序，要求数据从文件中读取。

计算部分也界面无关。

3.设计思路：在本次的实验中，我着重想表现的是一种面向对象的编程思想。

于是，在程序中我设计了4个类：CPoint、CObserve、CMatrix、Leveling,分别定义点的属性、观测数据属性、矩阵和水准网平差计算的属性与方法。

水准网平差计算一般步骤为：(1)读取观测数据和已知数据；(2)计算未知点高程近似值；(3)列高差观测值误差方程；(4)根据水准路线长度计算高差观测值的权；(5)组成法方程；(6)解法方程，求得未知点高程改正数及平差后高程值；(7)求高差观测值残差及平差后高差观测值；(8)精度评定；(9)输出平差结果。

水准网高程近似值计算算法4.界面设计：仅添加了一个button按钮，单击后读取数据，并进行水准网平差计算，计算结果保存在记事本中5.主要代码：文件一： CPoint.h代码：class ControlPoint{public:ControlPoint(void){};~ControlPoint(void){};public:CString pointID; //点号double H; //高程};class CObserve{public:CObserve(void){};~CObserve(void){};public:ControlPoint *pStartObs; //后视点ControlPoint *pEndObs; //前视点double h; //路线长度double dDist; //高差};文件二：Leveling.h代码：#pragma once#include"CPoint.h"#include"Matrix.h"class Leveling{public:Leveling(void);~Leveling(void);private:ControlPoint *m_pKnownPoint; //已知点数组int m_iKnownPointCount; //已知点个数ControlPoint *m_pUnknownPoint; //待测点数组int m_iUnknownPointCount; //待测点个数CObserve *m_pObsData; //观测数据数组int m_iObsDataCount; //观测数据个数public:bool LoadData(const CString &FileName); //从文件中导入数据void OutMatrixToFile(const CMatrix& mat,CStdioFile& SF); //把矩阵输出到文件中void SetKownPointSize(int n); //设置已知点数据大小void SetUnknownPointSize(int n);void SetObsDataSize(int n);CString* SplitString(CString str , char split, int iSubStrs); //分割字符串public://根据点号从已知点数组中找到控制点，并返回该点的指针ControlPoint* SearchKnownPointUsingID(CString pointID);//根据点号从未知点数组中找到控制点，并返回该点的指针ControlPoint* SearchUnknownPointUsingID(CString pointID);//根据点号从未知点和已知点数组中找到控制点，并返回该点的指针ControlPoint* SearchPointUsingID(CString pointID);void ComputeApproximateH(void); //求待测点高差近似值public://组成误差方程,B 为系数矩阵，f为常数项向量void FormErrorEquations(CMatrix& B, CMatrix& f);void Weight(CMatrix& P); //求高差观测值的权阵void IndirectlyAdjust(const CString& strFileName);//间接平差计算主函数};文件三：Leveling.cpp代码：#include"StdAfx.h"#include"Leveling.h"#include"math.h"#include<locale.h>Leveling::Leveling(void){m_pKnownPoint=NULL;m_iKnownPointCount=0;m_pUnknownPoint=NULL;m_iUnknownPointCount=0;m_pObsData=NULL;m_iObsDataCount=0;}Leveling::~Leveling(void){//释放动态数组内存if(m_pUnknownPoint!=NULL){delete[] m_pUnknownPoint;m_pUnknownPoint=NULL;}if(m_pKnownPoint!=NULL){delete[] m_pKnownPoint;m_pKnownPoint=NULL;}if(m_pObsData!=NULL){delete[] m_pObsData;m_pObsData=NULL;}}bool Leveling::LoadData(const CString &FileName){CStdioFile sf; //创建文件对象//以读的形式打开文件，如果打开失败则返回if(!sf.Open(FileName, CFile::modeRead)) return false;CString strLine;BOOL bEOF=sf.ReadString(strLine);//读取第一行,已知点数m_iKnownPointCount= _ttoi((strLine)); //把读取的第一行字符串转换为数值型SetKownPointSize(m_iKnownPointCount);//设置已知点数组大小int n=0;//读取并保存已知点数据for(int i=0;i<m_iKnownPointCount;i++){sf.ReadString(strLine);CString *pstrData=SplitString(strLine,',',n);m_pKnownPoint[i].pointID=pstrData[0];m_pKnownPoint[i].H=_tstof(pstrData[1]);delete[] pstrData;pstrData=NULL;}//开始读取未知知点数据sf.ReadString(strLine);//未知点个数m_iUnknownPointCount= _ttoi((strLine));SetUnknownPointSize(m_iUnknownPointCount); //设置未知点数组大小sf.ReadString(strLine);//未知点点号//读取并保存未知点数据CString *pstrData=SplitString(strLine,',',n);for(int i=0;i<m_iUnknownPointCount;i++){m_pUnknownPoint[i].pointID=pstrData[i];m_pUnknownPoint[i].H=0;}delete[] pstrData;pstrData=NULL;//开始读取观测数据个数sf.ReadString(strLine);//观测数个数m_iObsDataCount= _ttoi((strLine));SetObsDataSize(m_iObsDataCount);//读取并保存观测数据for(int i=0;i<m_iObsDataCount;i++){sf.ReadString(strLine);CString *pstrData=SplitString(strLine,',',n);m_pObsData[i].pStartObs=SearchPointUsingID(pstrData[0]);m_pObsData[i].pEndObs=SearchPointUsingID(pstrData[1]);m_pObsData[i].h=_tstof(pstrData[2]);m_pObsData[i].dDist=_tstof(pstrData[3]);delete[] pstrData;pstrData=NULL;}sf.Close();return true;}void Leveling::SetKownPointSize(int n){m_iKnownPointCount=n;m_pKnownPoint=new ControlPoint[n];}void Leveling::SetUnknownPointSize(int n){m_iUnknownPointCount=n;m_pUnknownPoint=new ControlPoint[n];}void Leveling::SetObsDataSize(int n){m_iObsDataCount=n;m_pObsData=new CObserve[n];}//字符串分割函数CString* Leveling::SplitString(CString str , char split, int iSubStrs) {int iPos = 0; //分割符位置int iNums = 0; //分割符的总数CString strTemp = str;CString strRight;//先计算子字符串的数量while (iPos != -1){iPos = strTemp.Find(split);if (iPos == -1){break;}strRight = strTemp.Mid(iPos + 1, str.GetLength());strTemp = strRight;iNums++;}if (iNums == 0) //没有找到分割符{//子字符串数就是字符串本身iSubStrs = 1;return NULL;}//子字符串数组iSubStrs = iNums + 1; //子串的数量= 分割符数量+ 1CString* pStrSplit;pStrSplit = new CString[iSubStrs];strTemp = str;CString strLeft;for (int i = 0; i < iNums; i++){iPos = strTemp.Find(split);//左子串strLeft = strTemp.Left(iPos);//右子串strRight = strTemp.Mid(iPos + 1, strTemp.GetLength()); strTemp = strRight;pStrSplit[i] = strLeft;}pStrSplit[iNums] = strTemp;return pStrSplit;}//根据点号从已知点数组中找到控制点，并返回该点的指针ControlPoint* Leveling::SearchKnownPointUsingID(CString pointID) {for(int i=0;i<m_iKnownPointCount;i++){if(pointID==m_pKnownPoint[i].pointID){return &m_pKnownPoint[i];}}return NULL;}//根据点号从未知点数组中找到控制点，并返回该点的指针ControlPoint* Leveling::SearchUnknownPointUsingID(CString pointID) {for(int i=0;i<m_iUnknownPointCount;i++){if(pointID==m_pUnknownPoint[i].pointID){return &m_pUnknownPoint[i];}}return NULL;}//根据点号从未知点和已知点数组中找到控制点，并返回该点的指针ControlPoint* Leveling::SearchPointUsingID(CString pointID){ControlPoint* pCP=NULL;pCP=SearchKnownPointUsingID(pointID);if(pCP==NULL){pCP=SearchUnknownPointUsingID(pointID);}return pCP;}void Leveling::ComputeApproximateH(void) //求待测点近似高程{for(int i=0;i<m_iUnknownPointCount;i++){if(m_pUnknownPoint[i].H!=0)continue;for(int j=0;j<m_iObsDataCount;j++){if(m_pUnknownPoint[i].pointID==m_pObsData[j].pStartObs->pointID&&m_pObsData[j].pEndObs->H>0)m_pUnknownPoint[i].H=m_pObsData[j].pEndObs->H-m_pObsData[j].h;elseif(m_pUnknownPoint[i].pointID==m_pObsData[j].pEndObs->pointID&&m_pObsData[j].pStartObs->H>0)m_pUnknownPoint[i].H=m_pObsData[j].pStartObs->H+m_pObsData[j].h;}}}void Leveling::FormErrorEquations(CMatrix& B, CMatrix& f){B.SetSize(m_iObsDataCount, m_iUnknownPointCount);f.SetSize(m_iObsDataCount, 1);//计算系数阵for(int i=0;i<m_iObsDataCount;i++)for(int j=0;j<m_iUnknownPointCount;j++){if(m_pObsData[i].pStartObs->pointID==m_pUnknownPoint[j].pointID) B(i,j)=-1;if(m_pObsData[i].pEndObs->pointID==m_pUnknownPoint[j].pointID)B(i,j)=1;}//计算常数阵for(int i=0;i<m_iObsDataCount;i++)f(i,0)=m_pObsData[i].h-(m_pObsData[i].pEndObs->H-m_pObsData[i].pStartObs->H);}void Leveling::Weight(CMatrix& P) //计算权阵{P.SetSize(m_iObsDataCount,m_iObsDataCount);for(int i=0;i<m_iObsDataCount;i++)for(int j=0;j<m_iObsDataCount;j++){if(i==j)P(i,j)=1.0/m_pObsData[i].dDist;}}void Leveling::IndirectlyAdjust(const CString& strFileName){//使用最小二乘法计算平差值//NBB=BT*P*BCMatrix B; //系数矩阵CMatrix BT; //B的转置矩阵CMatrix NBB; //NBBCMatrix invNbb; //NBB逆阵CMatrix W; //W=BT*P*fCMatrix f; //常数项矩阵CMatrix P; //权阵CMatrix x; //近似值改正项矩阵CMatrix V; //改正项矩阵CMatrix Omiga; //方差CMatrix X; //平差值ComputeApproximateH();FormErrorEquations(B,f); //系数阵Weight(P); //权阵BT=~B; //B的转置矩阵NBB=BT*P*B;invNbb=NBB.Inv();; //NBB逆阵W=BT*P*f;x=invNbb*W;X.SetSize(m_iUnknownPointCount,1);for(int i=0;i<m_iUnknownPointCount;i++){X(i,0)=m_pUnknownPoint[i].H+x(i,0);}//精度评定V=B*x-f;//开始输出间接平差的结果CStdioFile SF;CString strLine;setlocale(LC_ALL,"");if(!SF.Open(strFileName, CFile::modeCreate|CFile::modeWrite)) return;//开始写数据SF.WriteString(_T("－－－－－－－－－－水准网间接平差结果－－－－－－－－－－\n"));//写已知点数据strLine.Format(_T("已知点个数：%d\n"),m_iKnownPointCount);SF.WriteString(strLine);for(int i=0;i<m_iKnownPointCount;i++){strLine.Format(_T("%s,%4f\n"),m_pKnownPoint[i].pointID,m_pKnownPoint[i].H);SF.WriteString(strLine);}//待测点平差值strLine.Format(_T("待测点个数：%d\n"),m_iUnknownPointCount);SF.WriteString(strLine);strLine.Format(_T("平差结果：\n"));SF.WriteString(strLine);for(int i=0;i<m_iUnknownPointCount;i++){strLine.Format(_T("%s,%4f\n"),m_pUnknownPoint[i].pointID,X(i,0));SF.WriteString(strLine);}//输出系数阵和常数项SF.WriteString(_T("B矩阵：\r\n"));OutMatrixToFile(B,SF);//输出权矩阵SF.WriteString(_T("\r\nP矩阵(对角阵)：\r\n"));OutMatrixToFile(P,SF);SF.WriteString(_T("高差改正数x(mm)：\r\n"));OutMatrixToFile(x,SF);SF.WriteString(_T("NBB矩阵：\r\n"));OutMatrixToFile(NBB,SF);SF.WriteString(_T("NBB矩阵的逆矩阵：\r\n"));OutMatrixToFile(invNbb,SF);SF.WriteString(_T("W矩阵：\r\n"));OutMatrixToFile(W,SF);//输出观测值残差SF.WriteString(_T("观测值残差(mm)：\r\n"));OutMatrixToFile(V,SF);//计算单位权中误差，并输出Omiga=~V*P*V;double Sigma;Sigma = sqrt(Omiga(0, 0) / (m_iObsDataCount - m_iUnknownPointCount));strLine.Format(_T("单位权中误差(mm)：%.4f\r\n"),Sigma);SF.WriteString(strLine);double Qx;SF.WriteString(_T("点位误差(mm)：\r\n"));for(int i=0;i<invNbb.Row();i++){Qx=sqrt(invNbb(i,i))*Sigma;strLine.Format(_T("%.4f "),Qx);SF.WriteString(strLine);}SF.WriteString(_T("\r\n"));SF.Close();}//把矩阵输出到文件中void Leveling::OutMatrixToFile(const CMatrix& mat,CStdioFile& SF){CString strLine,strTmp;for(int i=0;i<mat.Row();i++){strLine.Empty();for(int j=0;j<mat.Col();j++){strTmp.Format(_T("%.4f "),mat(i,j));strLine=strLine+strTmp;}SF.WriteString(strLine+_T("\r\n"));}}文件四：Matrix.h代码：#pragma onceclass CMatrix{public:CMatrix(int row=3,int col=3);// copy constructorCMatrix (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);//调整数组的大小，原有数据不变（未测试）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();//生成单位矩阵};文件五：LevelingAdjustDlg.cpp代码：void CLevelingAdjustDlg::OnBnClickedButton2(){// TODO: 在此添加控件通知处理程序代码CFileDialog dlgFile(TRUE,_T("txt"),NULL,OFN_EXPLORER,_T("(文本文件)|*.txt")); //创建打开文件对话框if(dlgFile.DoModal()==IDCANCEL) return;//如果选择取消按钮，则退出CString strFileName=dlgFile.GetPathName();//获取选择的文件的名称Leveling ST;ST.LoadData(strFileName);ST.IndirectlyAdjust(_T("IndirectlyAdjustResult.txt"));MessageBox(_T("计算完毕！\r\n结果保存在Result.txt文件中！"));ShellExecute(this->m_hWnd,_T("open"),_T("notepad.exe"),_T("IndirectlyAdjustResult.txt"),_T(""),SW_SHOW );}6.运行结果：实验的运行结果如下图所示：8.实验总结这次的实验让我对于面对对象的编程思想有了进一步的理解。

B=NaN(45,30);P=NaN(45,45);Si=NaN(19,1);X0=NaN(18,1);Y0=NaN(18,1);XY=NaN(30,1);S0=NaN(19,1);L=NaN(45,1);r=15;rou=206265;Li=NaN(26,1);L0=NaN(26,1);a=NaN(26,1);b=NaN(26,1);c=NaN(26,1);d=NaN(26,1);j1=NaN(26,1);j2=NaN(26,1);sgm0=6;X1=NaN(15,1);Y1=NaN(15,1);S1=NaN(19,1);L1=NaN(26,1);L1dms=NaN(26,3);Sdev=NaN(19,1);Qx=NaN(15,1);Qy=NaN(15,1);Qxy=NaN(15,1);K=NaN(15,1);E=NaN(15,1);F=NaN(15,1);Qe=NaN(15,1);Fe=NaN(15,1);Fai=NaN(15,3);Li(1,1)=(134+35/60+52.7/3600)/180*pi;Li(2,1)=(202+19/60+10.7/3600)/180*pi;Li(3,1)=(129+30/60+14.9/3600)/180*pi;Li(4,1)=(162+9/60+14.2/3600)/180*pi;Li(5,1)=(138+43/60+48.6/3600)/180*pi;Li(6,1)=(175+52/60+30.1/3600)/180*pi;Li(7,1)=(144+48/60+40.2/3600)/180*pi;Li(8,1)=(134+33/60+3.1/3600)/180*pi;Li(9,1)=(189+13/60+5.9/3600)/180*pi;Li(10,1)=(170+33/60+51.6/3600)/180*pi;Li(11,1)=(112+44/60+18.3/3600)/180*pi;Li(12,1)=(171+13/60+58/3600)/180*pi;Li(13,1)=(146+38/60+42.3/3600)/180*pi;Li(14,1)=(360-(212+55/60+25.8/3600))/180*pi;Li(15,1)=(50+26/60+14.2/3600)/180*pi;Li(16,1)=(360-(297+41/60+39.7/3600))/180*pi;Li(17,1)=(167+26/60+28.8/3600)/180*pi;Li(18,1)=(360-(245+3/60+31.1/3600))/180*pi;Li(19,1)=(91+23/60+26.5/3600)/180*pi;Li(20,1)=(360-(206+19/60+41.2/3600))/180*pi;Li(21,1)=(360-(189+32/60+52.5/3600))/180*pi;Li(22,1)=(66+59/60+34.1/3600)/180*pi;Li(23,1)=(360-(279+55/60+6.8/3600))/180*pi;Li(24,1)=(155+58/60+21.4/3600)/180*pi;Li(25,1)=(96+35/60+25.7/3600)/180*pi;Li(26,1)=(360-(294+25/60+59.4/3600))/180*pi;Si(1,1)=45.878;Si(2,1)=25.083;Si(3,1)=58.859;Si(4,1)=33.166;Si(5,1)=43.223;Si(6,1)=62.642;Si(7,1)=37.399;Si(8,1)=46.764;Si(9,1)=39.761;Si(10,1)=35.981;Si(11,1)=40.878;Si(12,1)=53.546;Si(13,1)=45.165;Si(14,1)=64.024;Si(15,1)=44.644;Si(16,1)=56.887;Si(17,1)=42.365;Si(18,1)=52.129;Si(19,1)=56.626;X0(1,1)=256.832;Y0(1,1)=288.001;X0(2,1)=267.182;Y0(2,1)=265.153;X0(3,1)=310.011;Y0(3,1)=224.779;X0(4,1)=307.803;Y0(4,1)=191.667;X0(5,1)=291.856;Y0(5,1)=151.494;X0(6,1)=236.081;Y0(6,1)=122.978;X0(7,1)=201.643;Y0(7,1)=108.393;X0(8,1)=155.941;Y0(8,1)=118.303;X0(9,1)=134.685;Y0(9,1)=151.906;X0(10,1)=110.827;Y0(10,1)=178.840;X0(11,1)=89.129;Y0(11,1)=213.509;X0(12,1)=119.980;Y0(12,1)=257.274; X0(13,1)=151.325;Y0(13,1)=289.791;X0(14,1)=251.244;Y0(14,1)=188.915; X0(15,1)=190.059;Y0(15,1)=212.716;X0(16,1)=133.395;Y0(16,1)=207.687; X0(17,1)=204.453;Y0(17,1)=252.576;X0(18,1)=213.783;Y0(18,1)=303.863; %{X0(1,1)=256.829;Y0(1,1)=287.987;X0(2,1)=267.176;Y0(2,1)=265.120;X0(3,1)=310.005;Y0(3,1)=224.746;X0(4,1)=307.804;Y0(4,1)=191.653;X0(5,1)=291.857;Y0(5,1)=151.480;X0(6,1)=236.082;Y0(6,1)=122.964;X0(7,1)=201.667;Y0(7,1)=108.433;X0(8,1)=155.965;Y0(8,1)=118.343;X0(9,1)=134.709;Y0(9,1)=151.946;X0(10,1)=110.851;Y0(10,1)=178.880;X0(11,1)=89.129;Y0(11,1)=213.509;X0(12,1)=119.980;Y0(12,1)=257.274; X0(13,1)=151.325;Y0(13,1)=289.791;X0(14,1)=251.244;Y0(14,1)=188.915; X0(15,1)=190.059;Y0(15,1)=212.716;X0(16,1)=133.395;Y0(16,1)=207.687; X0(17,1)=204.453;Y0(17,1)=252.576;X0(18,1)=213.780;Y0(18,1)=303.849; %}%X0(18,1) is the coordinate of point 0 ;S0(1,1)=sqrt((X0(18,1)-X0(1,1))^2+(Y0(18,1)-Y0(1,1))^2);S0(14,1)=sqrt((X0(13,1)-X0(18,1))^2+(Y0(13,1)-Y0(18,1))^2);S0(15,1)=sqrt((X0(11,1)-X0(16,1))^2+(Y0(11,1)-Y0(16,1))^2);S0(16,1)=sqrt((X0(15,1)-X0(16,1))^2+(Y0(15,1)-Y0(16,1))^2);S0(17,1)=sqrt((X0(15,1)-X0(17,1))^2+(Y0(15,1)-Y0(17,1))^2);S0(18,1)=sqrt((X0(18,1)-X0(17,1))^2+(Y0(18,1)-Y0(17,1))^2);S0(19,1)=sqrt((X0(4,1)-X0(14,1))^2+(Y0(4,1)-Y0(14,1))^2);for i=1:12S0(i+1,1)=sqrt((X0(i,1)-X0(i+1,1))^2+(Y0(i,1)-Y0(i+1,1))^2);end%the angle operationa(1,1)=Y0(2,1)-Y0(1,1); b(1,1)=X0(2,1)-X0(1,1);c(1,1)=Y0(18,1)-Y0(1,1); d(1,1)=X0(18,1)-X0(1,1);for i=2:12a(i,1)=Y0(i+1,1)-Y0(i,1); b(i,1)=X0(i+1,1)-X0(i,1);c(i,1)=Y0(i-1,1)-Y0(i,1); d(i,1)=X0(i-1,1)-X0(i,1);enda(13,1)=Y0(18,1)-Y0(13,1); b(13,1)=X0(18,1)-X0(13,1);c(13,1)=Y0(12,1)-Y0(13,1); d(13,1)=X0(12,1)-X0(13,1);% the next is the 14th angle ;a(14,1)=Y0(1,1)-Y0(18,1); b(14,1)=X0(1,1)-X0(18,1);c(14,1)=Y0(13,1)-Y0(18,1); d(14,1)=X0(13,1)-X0(18,1);a(15,1)=Y0(16,1)-Y0(11,1); b(15,1)=X0(16,1)-X0(11,1);c(15,1)=Y0(10,1)-Y0(11,1); d(15,1)=X0(10,1)-X0(11,1);a(16,1)=Y0(12,1)-Y0(11,1); b(16,1)=X0(12,1)-X0(11,1);c(16,1)=Y0(16,1)-Y0(11,1); d(16,1)=X0(16,1)-X0(11,1);a(17,1)=Y0(11,1)-Y0(16,1); b(17,1)=X0(11,1)-X0(16,1);c(17,1)=Y0(15,1)-Y0(16,1); d(17,1)=X0(15,1)-X0(16,1);a(18,1)=Y0(16,1)-Y0(15,1); b(18,1)=X0(16,1)-X0(15,1);c(18,1)=Y0(17,1)-Y0(15,1); d(18,1)=X0(17,1)-X0(15,1);a(19,1)=Y0(17,1)-Y0(15,1); b(19,1)=X0(17,1)-X0(15,1);c(19,1)=Y0(14,1)-Y0(15,1); d(19,1)=X0(14,1)-X0(15,1);a(20,1)=Y0(14,1)-Y0(15,1); b(20,1)=X0(14,1)-X0(15,1);c(20,1)=Y0(16,1)-Y0(15,1); d(20,1)=X0(16,1)-X0(15,1);a(21,1)=Y0(15,1)-Y0(17,1); b(21,1)=X0(15,1)-X0(17,1);c(21,1)=Y0(18,1)-Y0(17,1); d(21,1)=X0(18,1)-X0(17,1);a(22,1)=Y0(17,1)-Y0(18,1); b(22,1)=X0(17,1)-X0(18,1);c(22,1)=Y0(13,1)-Y0(18,1); d(22,1)=X0(13,1)-X0(18,1);a(23,1)=Y0(1,1)-Y0(18,1); b(23,1)=X0(1,1)-X0(18,1);c(23,1)=Y0(17,1)-Y0(18,1); d(23,1)=X0(17,1)-X0(18,1);a(24,1)=Y0(15,1)-Y0(14,1); b(24,1)=X0(15,1)-X0(14,1);c(24,1)=Y0(4,1)-Y0(14,1); d(24,1)=X0(4,1)-X0(14,1);a(25,1)=Y0(14,1)-Y0(4,1); b(25,1)=X0(14,1)-X0(4,1);c(25,1)=Y0(3,1)-Y0(4,1); d(25,1)=X0(3,1)-X0(4,1);a(26,1)=Y0(5,1)-Y0(4,1); b(26,1)=X0(5,1)-X0(4,1);c(26,1)=Y0(14,1)-Y0(4,1); d(26,1)=X0(14,1)-X0(4,1);for i=1:26if(a(i,1)>0&&b(i,1)>0)j1(i,1)=atan(a(i,1)/b(i,1));else if(a(i,1)>0&&b(i,1)<0)j1(i,1)=atan(a(i,1)/b(i,1))+pi;else if(a(i,1)<0&&b(i,1)<0)j1(i,1)=atan(a(i,1)/b(i,1))+pi;elsej1(i,1)=atan(a(i,1)/b(i,1))+2*pi;endendendif(c(i,1)>0&&d(i,1)>0)j2(i,1)=atan(c(i,1)/d(i,1));else if(c(i,1)>0&&d(i,1)<0)j2(i,1)=atan(c(i,1)/d(i,1))+pi;else if(c(i,1)<0&&d(i,1)<0)j2(i,1)=atan(c(i,1)/d(i,1))+pi;elsej2(i,1)=atan(c(i,1)/d(i,1))+2*pi;endendendL0(i,1)=j1(i,1)-j2(i,1);if(L0(i,1)<0)L0(i,1)= L0(i,1)+2*pi;endendfor i=20:45L(i,1)=Li(i-19,1)-L0(i-19,1);L(i,1)=L(i,1)*rou;endB=zeros(45,30);B(1,21)=-(X0(1,1)-X0(18,1))/S0(1,1);B(1,22)=-(Y0(1,1)-Y0(18,1))/S0(1,1);B(1,1)=(X0(1,1)-X0(18,1))/S0(1,1);B(1,2)=(Y0(1,1)-Y0(18,1))/S0(1,1);for i=2:10B(i,2*i-3)=-(X0(i,1)-X0(i-1,1))/S0(i,1);B(i,2*i-2)=-(Y0(i,1)-Y0(i-1,1))/S0(i,1);B(i,2*i-1)=(X0(i,1)-X0(i-1,1))/S0(i,1);B(i,2*i)=(Y0(i,1)-Y0(i-1,1))/S0(i,1);endfor i=1:19L(i,1)=(Si(i,1)-S0(i,1))*1000;endB(11,19)=-(X0(11,1)-X0(10,1))/S0(11,1);B(11,20)=-(Y0(11,1)-Y0(10,1))/S0(11,1);B(12,23)=(X0(12,1)-X0(11,1))/S0(12,1);B(12,24)=(Y0(12,1)-Y0(11,1))/S0(12,1);B(13,23)=-(X0(13,1)-X0(12,1))/S0(13,1);B(13,24)=-(Y0(13,1)-Y0(12,1))/S0(13,1);B(13,25)=(X0(13,1)-X0(12,1))/S0(13,1);B(13,26)=(Y0(13,1)-Y0(12,1))/S0(13,1);B(14,25)=-(X0(18,1)-X0(13,1))/S0(14,1);B(14,26)=-(Y0(18,1)-Y0(13,1))/S0(14,1);B(14,21)=(X0(18,1)-X0(13,1))/S0(14,1);B(14,22)=(Y0(18,1)-Y0(13,1))/S0(14,1);B(15,27)=(X0(16,1)-X0(11,1))/S0(15,1);B(15,28)=(Y0(16,1)-Y0(11,1))/S0(15,1);B(16,27)=-(X0(15,1)-X0(16,1))/S0(16,1);B(16,28)=-(Y0(15,1)-Y0(16,1))/S0(16,1);B(17,29)=(X0(17,1)-X0(15,1))/S0(17,1);B(17,30)=(Y0(17,1)-Y0(15,1))/S0(17,1);B(18,29)=-(X0(18,1)-X0(17,1))/S0(18,1);B(18,30)=-(Y0(18,1)-Y0(17,1))/S0(18,1);B(18,21)=(X0(18,1)-X0(13,1))/S0(18,1);B(18,22)=(Y0(18,1)-Y0(13,1))/S0(18,1);B(19,7)=(X0(4,1)-X0(14,1))/S0(19,1);B(19,8)=(Y0(4,1)-Y0(14,1))/S0(19,1);%the rim is over;B(20,1)=rou*((Y0(2,1)-Y0(1,1))/(S0(2,1)^2)-(Y0(18,1)-Y0(1,1))/(S0(1,1)^ 2));B(20,2)=-rou*((X0(2,1)-X0(1,1))/(S0(2,1)^2)-(X0(18,1)-X0(1,1))/(S0(1,1) ^2));B(20,3)=-rou*((Y0(2,1)-Y0(1,1))/(S0(2,1)^2));B(20,4)=rou*((X0(2,1)-X0(1,1))/(S0(2,1)^2));B(20,21)=rou*(Y0(18,1)-Y0(1,1))/(S0(1,1)^2);B(20,22)=-rou*((X0(18,1)-X0(1,1))/(S0(1,1)^2));for i=2:9B(i+19,2*i-1)=rou*((Y0(i+1,1)-Y0(i,1))/(S0(i+1,1)^2)-(Y0(i-1,1)-Y0(i,1) )/(S0(i,1)^2));B(i+19,2*i)=-rou*((X0(i+1,1)-X0(i,1))/(S0(i+1,1)^2)-(X0(i-1,1)-X0(i,1)) /(S0(i,1)^2));B(i+19,2*i-3)=rou*((Y0(i-1,1)-Y0(i,1))/(S0(i,1)^2));B(i+19,2*i-2)=-rou*((X0(i-1,1)-X0(i,1))/(S0(i,1)^2));B(i+19,2*i+1)=-rou*((Y0(i+1,1)-Y0(i,1))/(S0(i+1,1)^2));B(i+19,2*i+2)=rou*((X0(i+1,1)-X0(i,1))/(S0(i+1,1)^2));endB(29,19)=rou*((Y0(11,1)-Y0(10,1))/(S0(11,1)^2)-(Y0(9,1)-Y0(10,1))/(S0(1 0,1)^2));B(29,20)=-rou*((X0(11,1)-X0(10,1))/(S0(11,1)^2)-(X0(9,1)-X0(10,1))/(S0( 10,1)^2));B(29,17)=rou*((Y0(9,1)-Y0(10,1))/(S0(10,1)^2));B(29,18)=-rou*((X0(9,1)-X0(10,1))/(S0(10,1)^2));B(30,19)=rou*((Y0(10,1)-Y0(11,1))/(S0(11,1)^2));B(30,20)=-rou*((X0(10,1)-X0(11,1))/(S0(11,1)^2));B(30,23)=-rou*((Y0(12,1)-Y0(11,1))/(S0(12,1)^2));B(30,24)=rou*((X0(12,1)-X0(11,1))/(S0(12,1)^2));B(31,23)=rou*((Y0(13,1)-Y0(12,1))/(S0(13,1)^2)-(Y0(11,1)-Y0(12,1))/(S0( 12,1)^2));B(31,24)=-rou*((X0(13,1)-X0(12,1))/(S0(13,1)^2)-(X0(11,1)-X0(12,1))/(S0 (12,1)^2));B(31,25)=-rou*((Y0(13,1)-Y0(12,1))/(S0(13,1)^2));B(31,26)=rou*((X0(13,1)-X0(12,1))/(S0(13,1)^2));B(32,25)=rou*((Y0(18,1)-Y0(13,1))/(S0(14,1)^2)-(Y0(12,1)-Y0(13,1))/(S0( 13,1)^2));B(32,26)=-rou*((X0(18,1)-X0(13,1))/(S0(14,1)^2)-(X0(12,1)-X0(13,1))/(S0 (13,1)^2));B(32,23)=rou*((Y0(12,1)-Y0(13,1))/(S0(13,1)^2));B(32,24)=-rou*((X0(12,1)-X0(13,1))/(S0(13,1)^2));B(32,21)=-rou*(Y0(18,1)-Y0(13,1))/(S0(14,1)^2);B(32,22)=rou*((X0(18,1)-X0(13,1))/(S0(14,1)^2));B(33,21)=rou*((Y0(1,1)-Y0(18,1))/(S0(1,1)^2)-(Y0(13,1)-Y0(18,1))/(S0(14 ,1)^2));B(33,22)=-rou*((X0(1,1)-X0(18,1))/(S0(1,1)^2)-(X0(13,1)-X0(18,1))/(S0(1 4,1)^2));B(33,25)=rou*((Y0(13,1)-Y0(18,1))/(S0(14,1)^2));B(33,26)=-rou*((X0(13,1)-X0(18,1))/(S0(14,1)^2));B(33,1)=-rou*(Y0(1,1)-Y0(18,1))/(S0(1,1)^2);B(33,2)=rou*((X0(1,1)-X0(18,1))/(S0(1,1)^2));B(34,19)=rou*((Y0(10,1)-Y0(11,1))/(S0(11,1)^2));B(34,20)=-rou*((X0(10,1)-X0(11,1))/(S0(11,1)^2));B(34,27)=-rou*((Y0(16,1)-Y0(11,1))/(S0(15,1)^2));B(34,28)=rou*((X0(16,1)-X0(11,1))/(S0(15,1)^2));B(35,27)=rou*((Y0(16,1)-Y0(11,1))/(S0(15,1)^2));B(35,28)=-rou*((X0(16,1)-X0(11,1))/(S0(15,1)^2));B(35,23)=-rou*((Y0(12,1)-Y0(11,1))/(S0(12,1)^2));B(35,24)=rou*((X0(12,1)-X0(11,1))/(S0(12,1)^2));B(36,27)=rou*((Y0(11,1)-Y0(16,1))/(S0(15,1)^2)-(Y0(15,1)-Y0(16,1))/(S0( 16,1)^2));B(36,28)=-rou*((X0(11,1)-X0(16,1))/(S0(15,1)^2)-(X0(15,1)-X0(16,1))/(S0 (16,1)^2));B(37,29)=rou*((Y0(17,1)-Y0(15,1))/(S0(17,1)^2));B(37,30)=-rou*((X0(17,1)-X0(15,1))/(S0(17,1)^2));B(37,27)=-rou*((Y0(16,1)-Y0(15,1))/(S0(16,1)^2));B(37,28)=rou*((X0(16,1)-X0(15,1))/(S0(16,1)^2));B(38,29)=-rou*((Y0(17,1)-Y0(15,1))/(S0(17,1)^2));B(38,30)=rou*((X0(17,1)-X0(15,1))/(S0(17,1)^2));B(39,27)=rou*((Y0(16,1)-Y0(15,1))/(S0(16,1)^2));B(39,28)=-rou*((X0(16,1)-X0(15,1))/(S0(16,1)^2));B(40,29)=rou*((Y0(15,1)-Y0(17,1))/(S0(17,1)^2)-(Y0(18,1)-Y0(17,1))/(S0(18,1)^2));B(40,30)=-rou*((X0(15,1)-X0(17,1))/(S0(17,1)^2)-(X0(18,1)-X0(17,1))/(S0 (18,1)^2));B(40,21)=rou*((Y0(18,1)-Y0(17,1))/(S0(18,1)^2));B(40,22)=-rou*((X0(18,1)-X0(17,1))/(S0(18,1)^2));B(41,21)=rou*((Y0(17,1)-Y0(18,1))/(S0(18,1)^2)-(Y0(13,1)-Y0(18,1))/(S0( 14,1)^2));B(41,22)=-rou*((X0(17,1)-X0(18,1))/(S0(18,1)^2)-(X0(13,1)-X0(18,1))/(S0 (14,1)^2));B(41,25)=rou*((Y0(13,1)-Y0(18,1))/(S0(14,1)^2));B(41,26)=-rou*((X0(13,1)-X0(18,1))/(S0(14,1)^2));B(41,29)=-rou*(Y0(17,1)-Y0(18,1))/(S0(18,1)^2);B(41,30)=rou*((X0(17,1)-X0(18,1))/(S0(18,1)^2));B(42,21)=rou*((Y0(1,1)-Y0(18,1))/(S0(1,1)^2)-(Y0(17,1)-Y0(18,1))/(S0(18 ,1)^2));B(42,22)=-rou*((X0(1,1)-X0(18,1))/(S0(1,1)^2)-(X0(17,1)-X0(18,1))/(S0(1 8,1)^2));B(42,29)=rou*((Y0(17,1)-Y0(18,1))/(S0(18,1)^2));B(42,30)=-rou*((X0(17,1)-X0(18,1))/(S0(18,1)^2));B(42,1)=-rou*(Y0(1,1)-Y0(18,1))/(S0(1,1)^2);B(42,2)=rou*((X0(1,1)-X0(18,1))/(S0(1,1)^2));B(43,7)=rou*(Y0(4,1)-Y0(14,1))/(S0(19,1)^2);B(43,8)=-rou*((X0(4,1)-X0(14,1))/(S0(19,1)^2));B(44,7)=rou*((Y0(14,1)-Y0(4,1))/(S0(19,1)^2)-(Y0(3,1)-Y0(4,1))/(S0(4,1) ^2));B(44,8)=-rou*((X0(14,1)-X0(4,1))/(S0(19,1)^2)-(X0(3,1)-X0(4,1))/(S0(4,1 )^2));B(44,5)=rou*((Y0(3,1)-Y0(4,1))/(S0(4,1)^2));B(44,6)=-rou*((X0(3,1)-X0(4,1))/(S0(4,1)^2));B(45,7)=rou*((Y0(5,1)-Y0(4,1))/(S0(5,1)^2)-(Y0(14,1)-Y0(4,1))/(S0(19,1) ^2));B(45,8)=-rou*((X0(5,1)-X0(4,1))/(S0(5,1)^2)-(X0(14,1)-X0(4,1))/(S0(19,1 )^2));B(45,9)=-rou*((Y0(5,1)-Y0(4,1))/(S0(5,1)^2));B(45,10)=rou*((X0(5,1)-X0(4,1))/(S0(5,1)^2));%µ¥Î»»¯³ÉͳһµÄºÁÃ×for i=20:45for j=1:30B(i,j)=B(i,j)/1000;endend%P matrixP=zeros(45,45);for i=1:19P(i,i)=sgm0^2/(S0(i,1)/2000*1000)^2;endfor i=20:45P(i,i)=sgm0^2/36;endXY=inv(B'*P*B)*B'*P*L;V=B*XY-L;sgm1=sqrt(V'*P*V/r);Q=inv(B'*P*B);Dx=sgm1^2*Q;for i=1:15X1(i,1)=XY(2*i-1,1)/1000;Y1(i,1)=XY(2*i,1)/1000;endfor i=1:10X1(i,1)= X1(i,1)+ X0(i,1);Y1(i,1)= Y1(i,1)+ Y0(i,1);endX1(11,1)= X1(11,1)+ X0(18,1);Y1(11,1)= Y1(11,1)+ Y0(18,1);for i=12:13X1(i,1)= X1(i,1)+ X0(i,1);Y1(i,1)= Y1(i,1)+ Y0(i,1);endX1(14,1)= X1(14,1)+ X0(16,1);Y1(14,1)= Y1(14,1)+ Y0(16,1);X1(15,1)= X1(15,1)+ X0(17,1);Y1(15,1)= Y1(15,1)+ Y0(17,1);% question 3S1(1,1)=sqrt((X1(11,1)-X1(1,1))^2+(Y1(11,1)-Y1(1,1))^2);S1(11,1)=sqrt((X1(10,1)-X0(11,1))^2+(Y1(10,1)-Y0(11,1))^2);S1(12,1)=sqrt((X0(11,1)-X1(12,1))^2+(Y0(11,1)-Y1(12,1))^2);S1(13,1)=sqrt((X1(12,1)-X1(13,1))^2+(Y1(12,1)-Y1(13,1))^2);S1(14,1)=sqrt((X1(13,1)-X1(11,1))^2+(Y1(13,1)-Y1(11,1))^2);S1(15,1)=sqrt((X1(14,1)-X0(11,1))^2+(Y1(14,1)-Y0(11,1))^2);S1(16,1)=sqrt((X1(14,1)-X0(15,1))^2+(Y1(14,1)-Y0(15,1))^2);S1(17,1)=sqrt((X1(15,1)-X0(15,1))^2+(Y1(15,1)-Y0(15,1))^2);S1(18,1)=sqrt((X1(15,1)-X1(11,1))^2+(Y1(15,1)-Y1(11,1))^2);S1(19,1)=sqrt((X1(4,1)-X0(14,1))^2+(Y1(4,1)-Y0(14,1))^2);for i=1:9S1(i+1,1)=sqrt((X1(i,1)-X1(i+1,1))^2+(Y1(i,1)-Y1(i+1,1))^2); end%¼ÆËã¸÷½Ç¶ÈµÄ¹Û²âÖµµÄƽ²îÖµ%******¿ÉÄÜ»áÓÐ360¶ÈµÄ²îÒì*****a(1,1)=Y1(2,1)-Y1(1,1); b(1,1)=X1(2,1)-X1(1,1);c(1,1)=Y1(11,1)-Y1(1,1); d(1,1)=X1(11,1)-X1(1,1);for i=2:9a(i,1)=Y1(i+1,1)-Y1(i,1); b(i,1)=X1(i+1,1)-X1(i,1);c(i,1)=Y1(i-1,1)-Y1(i,1); d(i,1)=X1(i-1,1)-X1(i,1);enda(10,1)=Y0(11,1)-Y1(10,1); b(10,1)=X0(11,1)-X1(10,1);c(10,1)=Y1(9,1)-Y1(10,1); d(10,1)=X1(9,1)-X1(10,1);a(11,1)=Y1(12,1)-Y0(11,1); b(11,1)=X1(12,1)-X0(11,1);c(11,1)=Y1(10,1)-Y0(11,1); d(11,1)=X1(10,1)-X0(11,1);a(12,1)=Y1(13,1)-Y1(12,1); b(12,1)=X1(13,1)-X1(12,1);c(12,1)=Y0(11,1)-Y1(12,1); d(12,1)=X0(11,1)-X1(12,1);a(13,1)=Y1(11,1)-Y1(13,1); b(13,1)=X1(11,1)-X1(13,1);c(13,1)=Y1(12,1)-Y1(13,1); d(13,1)=X1(12,1)-X1(13,1);a(14,1)=Y1(1,1)-Y1(11,1); b(14,1)=X1(1,1)-X1(11,1);c(14,1)=Y1(13,1)-Y1(11,1); d(14,1)=X1(13,1)-X1(11,1);a(15,1)=Y1(14,1)-Y0(11,1); b(15,1)=X1(14,1)-X0(11,1);c(15,1)=Y1(10,1)-Y0(11,1); d(15,1)=X1(10,1)-X0(11,1);a(16,1)=Y1(12,1)-Y0(11,1); b(16,1)=X1(12,1)-X0(11,1);c(16,1)=Y1(14,1)-Y0(11,1); d(16,1)=X1(14,1)-X0(11,1);a(17,1)=Y0(11,1)-Y1(14,1); b(17,1)=X0(11,1)-X1(14,1); c(17,1)=Y0(15,1)-Y1(14,1); d(17,1)=X0(15,1)-X1(14,1); a(18,1)=Y1(14,1)-Y0(15,1); b(18,1)=X1(14,1)-X0(15,1); c(18,1)=Y1(15,1)-Y0(15,1); d(18,1)=X1(15,1)-X0(15,1); a(19,1)=Y1(15,1)-Y0(15,1); b(19,1)=X1(15,1)-X0(15,1); c(19,1)=Y0(14,1)-Y0(15,1); d(19,1)=X0(14,1)-X0(15,1); a(20,1)=Y0(14,1)-Y0(15,1); b(20,1)=X0(14,1)-X0(15,1); c(20,1)=Y1(14,1)-Y0(15,1); d(20,1)=X1(14,1)-X0(15,1); a(21,1)=Y0(15,1)-Y1(15,1); b(21,1)=X0(15,1)-X1(15,1); c(21,1)=Y1(11,1)-Y1(15,1); d(21,1)=X1(11,1)-X1(15,1); a(22,1)=Y1(15,1)-Y1(11,1); b(22,1)=X1(15,1)-X1(11,1); c(22,1)=Y1(13,1)-Y1(11,1); d(22,1)=X1(13,1)-X1(11,1); a(23,1)=Y1(1,1)-Y1(11,1); b(23,1)=X1(1,1)-X1(11,1);c(23,1)=Y1(15,1)-Y1(11,1); d(23,1)=X1(15,1)-X1(11,1); a(24,1)=Y0(15,1)-Y0(14,1); b(24,1)=X0(15,1)-X0(14,1); c(24,1)=Y1(4,1)-Y0(14,1); d(24,1)=X1(4,1)-X0(14,1);a(25,1)=Y0(14,1)-Y1(4,1); b(25,1)=X0(14,1)-X1(4,1);c(25,1)=Y1(3,1)-Y1(4,1); d(25,1)=X1(3,1)-X1(4,1);a(26,1)=Y1(5,1)-Y1(4,1); b(26,1)=X1(5,1)-X1(4,1);c(26,1)=Y0(14,1)-Y1(4,1); d(26,1)=X0(14,1)-X1(4,1); for i=1:26if(a(i,1)>0&&b(i,1)>0)j1(i,1)=atan(a(i,1)/b(i,1));else if(a(i,1)>0&&b(i,1)<0)j1(i,1)=atan(a(i,1)/b(i,1))+pi;else if(a(i,1)<0&&b(i,1)<0)j1(i,1)=atan(a(i,1)/b(i,1))+pi;elsej1(i,1)=atan(a(i,1)/b(i,1))+2*pi;endendendif(c(i,1)>0&&d(i,1)>0)j2(i,1)=atan(c(i,1)/d(i,1));else if(c(i,1)>0&&d(i,1)<0)j2(i,1)=atan(c(i,1)/d(i,1))+pi;else if(c(i,1)<0&&d(i,1)<0)j2(i,1)=atan(c(i,1)/d(i,1))+pi;elsej2(i,1)=atan(c(i,1)/d(i,1))+2*pi;endendendL1(i,1)=j1(i,1)-j2(i,1);if(L1(i,1)<0)L1(i,1)= L1(i,1)+2*pi;endend%***»¡¶Èת»¯Îª½Ç¶È***for i=1:26L1dms(i,:)=degrees2dms(rad2deg(L1(i,1)));end%******±ß³¤Ïà¶ÔÖÐÎó²î*****for i=1:19Sdev(i,1)=sqrt(Dx(i,i))/S1(i,1)/1000;end%*******Îó²îÍÖÔ²****for i=1:15Qx(i,1)=Q(2*i-1,2*i-1);Qy(i,1)=Q(2*i,2*i);Qxy(i,1)=Q(2*i-1,2*i);K(i,1)=sqrt((Qx(i,1)-Qy(i,1))^2+4*Qxy(i,1)^2);E(i,1)=sqrt(0.5*sgm1^2*(Qx(i,1)+Qy(i,1)+K(i,1))); F(i,1)=sqrt(0.5*sgm1^2*(Qx(i,1)+Qy(i,1)-K(i,1))); Qe(i,1)=0.5*(Qx(i,1)+Qy(i,1)+K(i,1));Fe(i,1)=atan((Qe(i,1)-Qx(i,1))/Qxy(i,1));if(Fe(i,1)<0)Fe(i,1)= Fe(i,1)+2*pi;endFe(i,1)=rad2deg(Fe(i,1));Fai(i,:)=degrees2dms( Fe(i,1));end% the next is used for iteration%{for i= 1:10X0(i,1)=X1(i,1);Y0(i,1)=Y1(i,1);endX0(18,1)=X1(11,1);Y0(18,1)=Y1(11,1);for i=12:13X0(i,1)=X1(i,1);Y0(i,1)=Y1(i,1);endX0(16,1)=X1(14,1);Y0(16,1)=Y1(14,1);X0(17,1)=X1(15,1);Y0(17,1)=Y1(15,1);%}。

教材《误差理论与测量平差基础》第二版武汉大学出版社P108页的例7—1的运行结果：源程序：＃define N 5 /*N是观测值个数*/#define T 3 /*T是必要观测数*/＃include<stdio.h>#include<math。

h>float Nbb［T］［T］,Nb[T］[T]，W[T］［1],x［T]［1］；main（）｛float D（float a［T]［N］，float b［N][N］，float c[N][T]）;float K（float a[T][N］,float b［N］[N］，float c［N］[1])；float G(float a[T］[T］）；float F（float ca［T-1］[T—1］);float DM(float a[1][N］，float b［N］[N] ，float c［N］[1]）；int i,j,m,n;float B［N]［T］,BT［T］［N],V［N]［1］,VT［1］[N］，P［N］［N],C[N］[1]，Bx[N][1],f，g，h，x1; printf（"请输入V的系数B［N]［T]：\n”);for（i=0；i<N;i++）for(j=0;j〈T;j++)scanf（”％8f”，＆B［i］[j］)；printf（"请输入观测值的权阵P［N］[N］:\n");for（i=0；i<N；i++)for(j=0;j<N；j++)scanf（"%8f”,＆P[i]［j])；printf（"请输入常数C［N][1］:\n”）；for(i=0;i〈N；i++）for（j=0;j<1；j++)scanf(”％8f”，＆C[i][j］)；for(i=0；i〈N；i++）for（j=0;j<T；j++）BT[j]［i]=B［i］［j];g=D(BT, P， B);h=K（BT， P， C）；f=G(Nbb);for（i=0；i〈T;i++）for（j=0；j〈1;j++）{x[i］［j］=Nb［i］［0]*W［0］[j]；for(m=1；m〈T；m++)x［i］［j]+=（Nb[i］［m］＊W［m］［j］);｝for（i=0;i〈T;i++）x［i］［0]=x［i］[0]/f;for（i=0;i〈N;i++)for（j=0;j〈1；j++）{Bx［i][j]=B[i］[0］＊x［0］[j]；for(m=1;m〈T；m++)Bx［i］[j]+=(B［i］［m]＊x［m][j］);｝for(i=0;i〈N；i++)V[i］[0]=（Bx[i][0]-C[i][0]）;for(i=0；i<N；i++)for(j=0；j〈1;j++)VT[j][i]=V[i]［j]；x1=DM（VT，P,V)；x1=x1/（N-T）；printf("参数x［T］［1]=\n");for(i=0；i〈T;i++）printf("%15f"，x［i］[0]);printf（"\n"）；printf("改正数V［N][1］=\n”)；for（i=0；i<N；i++)printf（"%15f”,V[i］［0］）；printf（"\n单位权中误差x1=%15f”,sqrt(x1）); printf("\n协因数阵Qxx[T］[T]:\n"）；for(i=0;i<T;i++）｛for（j=0；j〈T；j++)printf（”%15f”，Nb［i］[j］/f);printf（”\n")；｝}float G（float a[T］［T］){int i，j，m，n;float c[T—1］[T—1］,y=0；for(i=0;i<T;i++）for（j=0;j<T;j++）{for（m=0；m<T;m++)for（n=0；n〈T；n++）{if(m〈i＆＆n<j）c[m][n］=a[m］［n］;if（m>i&&n〈j)c[m-1］[n]=a［m］[n］;if（m<i＆&n>j）c［m］[n—1]=a[m］［n]；if（m>i&&n>j）c[m-1]［n-1]=a［m]［n］；｝if（（i+j)%2==0)Nb[j］[i]=F（c);elseNb[j］［i］=(-1)＊F(c）;｝for(m=0；m<T；m++)y+=(a［0][m］＊Nb[m][0］）；return (y);}float F(float ca[T—1］[T—1]){int i,j,m,n，s，t，k=1;float f=1，c，x，sn;for （i=0，j=0；i<T—1&＆j〈T—1;i++，j++）｛if （ca［i］[j］==0){for （m=i；ca[m］［j]==0；m++）;if (m==T—1）｛sn=0;return （sn）;｝elsefor （n=j;n<T—1;n++)｛c=ca[i］[n］；ca［i］[n］=ca［m］［n];ca[m］［n］=c；｝k＊=（-1)；｝for (s=T—2；s>i；s--）｛x=ca[s][j];for (t=j;t〈T-1;t++）ca［s][t］—=ca[i][t]*（x/ca［i][j])；}｝for （i=0;i<T-1;i++)f*=ca［i]［i];sn=k＊f;return (sn);}float D（float a［T］[N］,float b[N][N] ,float c［N］［T])｛int i,j，m；float d［T］[N];for（i=0;i<T；i++)for（j=0;j〈N;j++）｛d[i][j］=a［i］[0]*b［0][j］;for(m=1；m〈N;m++）d[i］[j]+=(a［i][m］*b[m］[j]）；｝for（i=0;i〈T；i++）for（j=0;j〈T；j++)｛Nbb［i]［j］=d[i][0］*c［0］［j］；for(m=1；m〈N;m++)Nbb［i][j］+=（d[i］[m]＊c［m][j])；}return （Nbb［0］［0］）；｝float K(float a[T］[N］,float b[N]［N],float c［N]［1］）｛int i，j,m;float d［T]［N］;for（i=0；i<T;i++）for(j=0;j<N;j++）{d[i］［j］=a［i]［0］＊b［0]［j]；for(m=1;m<N；m++)d[i]［j]+=（a[i］[m］＊b[m］[j])；｝for（i=0;i<T;i++）for(j=0；j〈1；j++){W[i]［j]=d[i]［0]*c［0］[j]；for(m=1;m〈N；m++)W［i］[j］+=(d［i]［m］*c[m］[j］)；}return （W［0］[0］)；｝float DM（float a[1]［N],float b[N]［N] ，float c[N][1]）{int i，j,m;float d[1][N］，x;for（i=0;i〈1;i++）for(j=0;j<N；j++）｛d[i]［j］=a［i］[0]＊b[0］[j]；for（m=1；m<N;m++）d［i][j］+=（a［i］[m］*b［m］［j]);}for(i=0；i<1;i++)for（j=0；j<1;j++)｛x=d[i］[0］＊c［0］[j]；for（m=1；m<N;m++）x+=(d［i][m］*c[m][j]);｝return (x）；}程序说明：1) 用该程序前，根据具体情况输入N和T；2）该程序中的（N-T）自由度(即多余观测数）必须大于等于2，不然程序运行时会出错，原因在于求行列式的逆时，有语句for (s=R-2；s>i；s—-)，R=1时s=-1。

《测绘程序设计》上机实验报告（Visual C++.Net）班级：测绘0901班学号： 04姓名：代娅琴2012年4月29日实验八平差程序设计基础一、实验目的巩固过程的定义与调用巩固类的创建与使用巩固间接平差模型及平差计算掌握平差程序设计的基本技巧与步骤二、实验内容水准网平差程序设计。

设计一个水准网平差的程序，要求数据从文件中读取，计算部分与界面无关。

1.水准网间接平差模型：2.计算示例：近似高程计算：3.水准网平差计算一般步骤(1)读取观测数据和已知数据；(2)计算未知点高程近似值；(3)列高差观测值误差方程；(4)根据水准路线长度计算高差观测值的权；(5)组成法方程；(6)解法方程，求得未知点高程改正数及平差后高程值；(7)求高差观测值残差及平差后高差观测值；(8)精度评定；(9)输出平差结果。

4.水准网高程近似值计算算法5.输入数据格式示例实验代码：#pragma onceclass LevelControlPoint{public:LevelControlPoint(void);~LevelControlPoint(void);public:CString strName;trName=pstrData[0];m_pKnownPoint[i].strID=pstrData[0];m_pKnownPoint[i].H=_tstof(pstrData[1]);m_pKnownPoint[i].flag=1;trName=pstrData[i];m_pUnknownPoint[i].strID=pstrData[i];m_pUnknownPoint[i].H=0;lag=0;pBackObj=SearchPointUsingID(pstrData[0]);pFrontObj=SearchPointUsingI D(pstrData[1]);ObsValue=_tstof(pstrData[2]);ist=_tstof(pstrData[3]);trID==ID){return &m_pKnownPoint[i];}}return NULL;}trID==ID){return &m_pUnknownPoint[i];}}return NULL;}LevelControlPoint* AdjustLevel::SearchPointUsingID(CString ID){LevelControlPoint* cp;cp=SearchKnownPointUsingID(ID);if(cp==NULL)cp=SearchUnknownPointUsingID(ID);return cp;}void AdjustLevel::ApproHeignt(void)lag!=1){pFrontObj->strID==m_pUnknownPoint[i].strID)&& m_pDhObs[j].cpBackObj->flag==1 ){ =m_pDhObs[i].cpBackObj->H - m_pDhObs[i].ObsValue;*/m_pUnknownPoint[i].H=m_pDhObs[j].cpBackObj->H + m_pDhObs[j].HObsValue;m_pUnknownPoint[i].flag=1;break;}}if(m_pUnknownPoint[i].flag!=1)pBackObj->strID==m_pUnknownPoint[i].strID)&& m_pDhObs[j].cpFrontObj->flag==1 ){ =m_pDhObs[j].cpFrontObj->H-m_pDhObs[j].HObsValue;/* m_pUnknownPoint[i].H=m_pDhObs[i].cpFrontObj->H+m_pDhObs[i].ObsValue;*/ m_pUnknownPoint[i].flag=1;break;}}}}if(i==m_iUnknownPointCount-1)lag!=1)ist);p(i,i)=value;}return p;}void AdjustLevel::FormErrorEquation(CMatrix &B, CMatrix &L){(m_iDhObsCount,m_iUnknownPointCount);(m_iDhObsCount,1);for(int i=0;i<m_iDhObsCount;i++)pBackObj->strID);tmpFront=SearchPointUsingID(m_pDhObs[i].cpFrontObj->strID);trID==tmpBack->strID)trID==tmpFront->strID)bsValue-(m_pDhObs[i].cpBackObj->H-m_pDhObs[i].cpFrontO bj->H);*/L(i,0)=m_pDhObs[i].HObsValue-(m_pDhObs[i].cpFrontObj->H - m_pDhObs[i].cpBackObj->H);(_T("%.3f"),L(i,0));L(i,0)=_tstof(tmp);L(i,0)=L(i,0)*1000;+=x(i,0);xt"));xt"));if()==IDCANCEL) return;CString strFileName=();setlocale(LC_ALL,"");CStdioFile sf;if(!(strFileName, CFile::modeCreate|CFile::modeWrite)) return;(LevleContent);();UpdateData(FALSE);}void CIndircLelveDlg::OnBnClickedComputelevel(){f\r\n"), [i].strID,[i].H);LevleContent+=Temp;}(_T("单位权中误差：%.1f mm\r\n"),r0*1000);LevleContent+=Temp;LevleContent+=_T("未知点高程中误差(mm):\r\n");for(int i=0;i< ;i++){();(_T("%s,%.1f\r\n"),[i].strName,Qx[i]*1000);LevleContent+=Temp;}UpdateData(false);}void CIndircLelveDlg::OnBnClickedSavelevleresult(){xt"));if()==IDCANCEL) return;CString strFileName=();setlocale(LC_ALL,"");CStdioFile sf;if(!(strFileName, CFile::modeCreate|CFile::modeWrite)) return;(LevleContent);();UpdateData(FALSE);}三、实验结果打开文件数据：平差结果：四、实验心得这从实验是我们测绘程序设计的最后一次实验，虽然这个学期我们做了好几次相关的实验，但是我却发现自己学的东西也越来越模糊，感觉很多内容都不理解。

程序源代码C语言程序：#include <windows.h>#include <stdarg.h>#include <stdio.h>#include <stdlib.h>#include<math.h>#define PI 3.1415926535898#define p 206264.806247#define MAX 50//矩阵的乘法运算MatrixMutiply(double Matrix1[MAX][MAX],double Matrix2[MAX][MAX],double MatrixResult[MAX][MAX],int m1,int m2,int m3){int i,j,k;double Sum;for(i=0;i<m1;i++){for(j=0;j<m3;j++){/*按照矩阵乘法的规则计算结果矩阵的i*j元素*/Sum=0;for(k=0;k<m2;k++)Sum+=Matrix1[i][k]*Matrix2[k][j];MatrixResult[i][j]=Sum;}}//return MatrixResult;}//矩阵的转置运算MatrixT(double Matrix1[MAX][MAX],double T[MAX][MAX],int m1,int m2) {//m1,m2分别是矩阵的行列数int i,j;//double T[50][50];for(i=0;i<m1;i++){for(j=0;j<m2;j++){T[j][i]=Matrix1[i][j];}}//return T;}//矩阵的逆运算void swap(double a,double b){double c=a;a=b;b=c;}double DinV(double A[50][50],int n) //n代表阶数{int i,j,k;double d;int JS[50],IS[50];for (k=0;k<n;k++){d=0;for (i=k;i<n;i++)for (j=k;j<n;j++){if (fabs(A[i][j])>d){d=fabs(A[i][j]);IS[k]=i;JS[k]=j;}}if (d+1.0==1.0)return 0;if (IS[k]!=k)for (j=0;j<n;j++)swap(A[k][j],A[IS[k]][j]);if (JS[k]!=k)for (i=0;i<n;i++)swap(A[i][k],A[i][JS[k]]);A[k][k]=1/A[k][k];for (j=0;j<n;j++)if (j!=k) A[k][j]=A[k][j]*A[k][k];for (i=0;i<n;i++)if (i!=k)for (j=0;j<n;j++)if (j!=k) A[i][j]=A[i][j]-A[i][k]*A[k][j];for (i=0;i<n;i++)if (i!=k) A[i][k]=-A[i][k]*A[k][k]; }for (k=n-1;k>=0;k--){for (j=0;j<n;j++)if (JS[k]!=k) swap(A[k][j],A[JS[k]][j]);for (i=0;i<n;i++)if (IS[k]!=k)swap(A[i][k],A[i][IS[k]]);}}//将度分秒连写的角度化为弧度值double jd_hd(double D,double F,double M){//int dd=(int)((int)B/10000); //提取度值//int ff=(int)(((int)B-dd*10000)/100); //提取分值//double mm=B-dd*10000-ff*100;//提取秒值double B;B=(D*3600.0+F*60.0+M)/p;//角度化弧度return B;}main(){double D,F,M,sigma_beta=2.0,sigma_s;//scanf("%f%f%f",&D,&F,&M);//printf("%f",jd_hd(D,F,M));doublebeta[50],alf[50],alfo[50],s[50],so[50],Xo[50],Yo[50],B[50][50]={0.0},L[50][50],P[50][50 ]={0.0},c[50][50]={0.0},wx=0.0;//此处为Xo Yo，B矩阵赋初值为零beta代表夹角，alf方位角，alfo方位角近似值，s距离观测值，so距离近似值//Xo Yo坐标近似值，B[50][50]误差矩阵，L[50]为L矩阵，P[50][50]为权阵,c[1][30]代表限制条件的系数阵,w代表限制条件常数项doubleNbb[50][50],Ncc[50][50],W[50][50],Ks,xgu[50][50],Xgu[30][1],Ygu[30][1],V[50][50],sigma_gu,Q[50][50],sigma_xy[50][50];alf[0]=PI;alfo[0]=PI;int i,j,m1,m2,m3;//将测量数据导入，并存入相应数组FILE *fp3;char strline[100]; //读取文件每行的bufferint du[100],fen[100],miao[100];double bian[100];i=0,j=1;if((fp3=fopen("D:\\111\\测量数据.txt","r"))==NULL) //文件位置和文件名{printf("文件不存在！");return 0;}while(!feof(fp3)) //判断文件是否已到末尾{fgets(strline,100,fp3); //读取一行sscanf(strline,"%d %d %d %lf",&du[i],&fen[i],&miao[i],&bian[i]); //从文件读取到的一行数据分别存放在两个数组中i++;}fclose(fp3);while(1){//printf("\n%d\t%d\t%d\t%lf",du[j],fen[j],miao[j],bian[j]);beta[j-1]=jd_hd(du[j],fen[j],miao[j]);s[j-1]=bian[j];j++;if(j>=i){break;}}Xo[0]=5000.0;Yo[0]=5000.0;Xo[15]=5000.0;Yo[15]=5000.0;so[0]=s[0];//用来求未知点坐标近似值for(i=1;i<15;i++){alf[i]=alf[i-1]+beta[i]-PI;if(alf[i]>=(2*PI)){alf[i]=alf[i]-2*PI;}Xo[i]=Xo[i-1]+s[i-1]*cos(alf[i-1]);Yo[i]=Yo[i-1]+s[i-1]*sin(alf[i-1]);//printf("X=%f\t",Xo[i]);}for(i=1;i<15;i++){//求近似距离so[i]=sqrt((Yo[i+1]-Yo[i])*(Yo[i+1]-Yo[i])+(Xo[i+1]-Xo[i])*(Xo[i+1]-Xo[i]));//求近似方位角，分象限进行讨论if((Yo[i+1]-Yo[i])>0&&(Xo[i+1]-Xo[i])>0)//第一象限{alfo[i]=atan((Yo[i+1]-Yo[i])/(Xo[i+1]-Xo[i]));}else if((Yo[i+1]-Yo[i])>0&&(Xo[i+1]-Xo[i])<0)//第二象限{alfo[i]=atan((Yo[i+1]-Yo[i])/(Xo[i+1]-Xo[i]))+PI;}else if((Yo[i+1]-Yo[i])<0&&(Xo[i+1]-Xo[i])>0)//第三象限{alfo[i]=atan((Yo[i+1]-Yo[i])/(Xo[i+1]-Xo[i]))+2*PI;}else //((Yo[i]-Yo[i-1])<0&&(Xo[i]-Xo[i-1])<0)//第四象限{alfo[i]=atan((Yo[i+1]-Yo[i])/(Xo[i+1]-Xo[i]))+PI;}//printf("alf=%f\t",alfo[i]);//printf("so=%f\t",so[i]);}//求B矩阵//将第一个角度的系数单独算出j=0;//B[0][j]=0.0;((Yo[14]-Yo[0])/(so[14]*so[14]))*p/1000.0;//B[0][j+1]=0.0;-1*((Xo[14]-Xo[0])/(so[14]*so[14]))*p/1000.0;B[0][j]=-1*(Yo[1]-Yo[0])/(so[0]*so[0])*p/1000.0;//B[0][j+1]=(Xo[1]-Xo[0])/(so[0]*so[0])*p/1000.0;B[0][26]=((Yo[14]-Yo[0])/(so[14]*so[14]))*p/1000.0;B[0][27]=-1*(Xo[14]-Xo[0])/(so[14]*so[14])*p/1000.0;//将第二个角度的系数单独算出B[1][j]=((Yo[2]-Yo[1])/(so[1]*so[1])-(Yo[0]-Yo[1])/(so[0]*so[0]))*p/1000.0;//B[1][j+1]=-1*((Xo[2]-Xo[1])/(so[1]*so[1])-(Xo[0]-Xo[1])/(so[0]*so[0]))*p/1000.0 ;B[1][j+2]=-1*((Yo[2]-Yo[1])/(so[1]*so[1]))*p/1000.0;B[1][j+3]=((Xo[2]-Xo[1])/(so[1]*so[1]))*p/1000.0;//求其他角度改正的系数for(i=2;i<15;i++){if(i<14){B[i][j]=((Yo[i-1]-Yo[i])/(so[i-1]*so[i-1]))*p/1000.0;B[i][j+1]=-1*((Xo[i-1]-Xo[i])/(so[i-1]*so[i-1]))*p/1000.0;B[i][j+2]=((Yo[i+1]-Yo[i])/(so[i]*so[i])-(Yo[i-1]-Yo[i])/(so[i-1]*so[i-1]))*p/1000.0;B[i][j+3]=-1*((Xo[i+1]-Xo[i])/(so[i]*so[i])-(Xo[i-1]-Xo[i])/(so[i-1]*so[i-1]))*p/1000. 0;B[i][j+4]=-1*((Yo[i+1]-Yo[i])/(so[i]*so[i]))*p/1000.0;B[i][j+5]=(Xo[i+1]-Xo[i])/(so[i]*so[i])*p/1000.0;}else{B[i][j]=(Yo[i-1]-Yo[i])/(so[i-1]*so[i-1])*p/1000.0;B[i][j+1]=-1*((Xo[i-1]-Xo[i])/(so[i-1]*so[i-1]))*p/1000.0;B[i][j+2]=((Yo[i+1]-Yo[i])/(so[i]*so[i])-(Yo[i-1]-Yo[i])/(so[i-1]*so[i-1]))*p/1000.0;B[i][j+3]=-1*((Xo[i+1]-Xo[i])/(so[i]*so[i])-(Xo[i-1]-Xo[i])/(so[i-1]*so[i-1]))*p/1000. 0;}j=j+2;}B[2][1]=0.0;//求边长改正的系数j=0;//将第一个边长的系数单独算出B[i][j]=(Xo[1]-Xo[0])/so[0];//B[i][j+1]=(Yo[1]-Yo[0])/so[0]; //(Yo[1]-Yo[0])/(so[0]*so[0]); i=i+1;for(i;i<30;i++){if(i<29){B[i][j]=-1*(Xo[i-14]-Xo[i-15])/so[i-15];B[i][j+1]=-1*(Yo[i-14]-Yo[i-15])/so[i-15];B[i][j+2]=-1*B[i][j];B[i][j+3]=-1*B[i][j+1];}else{B[i][j]=-1*(Xo[i-14]-Xo[i-15])/so[i-15];B[i][j+1]=-1*(Yo[i-14]-Yo[i-15])/so[i-15];}j=j+2;}B[16][1]=0.0;for(j=1;j<27;j++){for(i=0;i<30;i++){B[i][j]=B[i][j+1];}}//求L矩阵L[0][0]=(beta[0]-(alfo[0]-alfo[14]+PI))*p;for(i=1;i<30;i++){if(i<15){L[i][0]=(beta[i]-(alfo[i]-alfo[i-1]+PI))*p;if(L[i][0]>PI*p){L[i][0]=(L[i][0]-2*PI*p);}L[i][0]=L[i][0];}else{L[i][0]=(s[i-15]-so[i-15])*1000.0;}}//求权阵Pfor(i=0;i<30;i++){if(i<15){P[i][i]=1;}else{sigma_s=5+10*0.000001*s[i-15]*1000; //单位为（''/mm）的平方P[i][i]=sigma_beta*sigma_beta/(sigma_s*sigma_s);}}double Tb[50][50],Tc[50][50],MatrixResult[50][50],TV[50][50];doubletemp1[MAX][MAX],temp2[MAX][MAX],temp3[MAX][MAX],temp4[MAX][MAX],temp5[MAX][MAX];//计算Nbb矩阵,W矩阵m1=30;m2=28;MatrixT(B,Tb,m1,m2);m3=30;MatrixMutiply(Tb,P,temp1,m2,m1,m3);m3=m2;m1=m2;m2=30;MatrixMutiply(temp1,B,Nbb,m1,m2,m3);m1=27;m2=30;m3=1;MatrixMutiply(temp1,L,W,m1,m2,m3);//矩阵输出D:\\111\\导线网输出数据3027.txt文本FILE *fp;fp=fopen("D:\\111\\导线网输出数据3027.txt","w");if(fp!=NULL){fprintf(fp,"距离近似值so(单位：m):\t");fprintf(fp,"方位角近似值alfo(单位：弧度):\n");for(i=0;i<15;i++){fprintf(fp,"%4.12lf\t",so[i]);fprintf(fp,"%.12lf\n",alfo[i]);}fprintf(fp,"Xo(单位：mm)：\t");fprintf(fp,"Yo(单位：mm)：");fprintf(fp,"\n");for(i=0;i<15;i++){fprintf(fp,"%.6lf\t",Xo[i]);fprintf(fp,"%.6lf",Yo[i]);fprintf(fp,"\n");}fprintf(fp,"B矩阵：");fprintf(fp,"\n");for(i=0;i<30;i++){for(j=0;j<27;j++){fprintf(fp,"%.6f ",B[i][j]);}fprintf(fp,"\n");}fprintf(fp,"L矩阵(单位：秒和mm)："); fprintf(fp,"\n");for(i=0;i<30;i++){fprintf(fp,"%.6lf ",L[i][0]);fprintf(fp,"\n");}fprintf(fp,"P矩阵：");fprintf(fp,"\n");for(i=0;i<30;i++){for(j=0;j<30;j++){fprintf(fp,"%.6f ",P[i][j]);}fprintf(fp,"\n");}fprintf(fp,"Nbb矩阵：");fprintf(fp,"\n");for(i=0;i<27;i++){for(j=0;j<27;j++){fprintf(fp,"%.12f ",Nbb[i][j]);}fprintf(fp,"\n");}fclose(fp); //写入完毕，关闭文件}DinV(Nbb,27); //MatrixResult=c * Nbb的逆,此时Nbb已经变成Nbb的逆//计算x^m1=27;m2=27;m3=1;MatrixMutiply(Nbb,W,xgu,m1,m2,m3);double xgu28[50][50],sigma_xy28[50][50];xgu28[0][0]=xgu[0][0];for(i=1;i<27;i++){xgu28[i+1][0]=xgu[i][0];}xgu28[1][0]=0.0;//计算X^(即Xgu估值)Xgu[0][0]=5000.0;Ygu[0][0]=5000.0;for(i=0;i<14;i++){Xgu[i+1][0]=Xo[i+1]+xgu28[2*i][0]/1000.0;Ygu[i+1][0]=Yo[i+1]+xgu28[2*i+1][0]/1000.0; }//精度评定m1=30;m2=27;m3=1;MatrixMutiply(B,xgu,temp5,m1,m2,m3);for(i=0;i<30;i++){if(i<15){V[i][0]=(temp5[i][0]-L[i][0]);}else{V[i][0]=(temp5[i][0]-L[i][0]);}}m1=30;m2=1;MatrixT(V,TV,m1,m2);m1=1;m2=30;m3=30;MatrixMutiply(TV,P,temp4,m1,m2,m3);m1=1;m2=30;m3=1;MatrixMutiply(temp4,V,temp4,m1,m2,m3);sigma_gu=sqrt(temp4[0][0]/3); //单位权中误差double vv=0.0;for (i=0;i<15;i++){for (j=0;j<1;j++){vv=vv+V[i][j];}//puts("");}printf("%lf\t",vv);for(i=0;i<27;i++){for(j=0;j<27;j++){Q[i][j]=Nbb[i][j];}}for(i=0;i<27;i++){sigma_xy[i][0]=sqrt(Q[i][i])*sigma_gu; //坐标平差值中误差//printf("%lf\n",sigma_xy[i][0]);}sigma_xy28[0][0]=sigma_xy[0][0];for(i=1;i<27;i++){sigma_xy28[i+1][0]=sigma_xy[i][0];}sigma_xy28[1][0]=0.0;//printf("%.10lf\n",Ncc[0][0]);FILE *fp1;fp1=fopen("D:\\111\\导线网输出数据3027.txt","a"); fprintf(fp1,"Nbb的逆：");fprintf(fp1,"\n");for(i=0;i<27;i++){for(j=0;j<27;j++){fprintf(fp1,"%lf ",Nbb[i][j]);}fprintf(fp1,"\n");}fprintf(fp1,"W：");fprintf(fp1,"\n");for(i=0;i<27;i++){for(j=0;j<1;j++){fprintf(fp1,"%.12lf ",W[i][j]);}fprintf(fp1,"\n");}fprintf(fp1,"x^(单位：mm)：");fprintf(fp1,"\n");for(i=0;i<28;i++){for(j=0;j<1;j++){fprintf(fp1,"%.10lf ",xgu28[i][j]);}fprintf(fp1,"\n");}fprintf(fp1,"V(单位：秒和mm)：");fprintf(fp1,"\n");for(i=0;i<30;i++){for(j=0;j<1;j++){fprintf(fp1,"%.10lf ",V[i][j]);}fprintf(fp1,"\n");}fprintf(fp1,"X^(单位：m)：\t");fprintf(fp1,"Y^(单位：m)：");fprintf(fp1,"\n");for(i=0;i<15;i++){for(j=0;j<1;j++){fprintf(fp1,"%.6lf\t",Xgu[i][j]);fprintf(fp1,"%.6lf",Ygu[i][j]);}fprintf(fp1,"\n");}fprintf(fp1,"单位权中误差(单位：mm)：\n"); fprintf(fp1,"%.10lf\n",sigma_gu);fprintf(fp1,"X坐标平差值中误差(单位：mm):\t");fprintf(fp1,"Y坐标平差值中误差(单位：mm)：");fprintf(fp1,"\n");for(i=0;i<13;i++){fprintf(fp1,"%lf\t\t",sigma_xy28[2*i][0]);fprintf(fp1,"%lf",sigma_xy28[2*i+1][0]);fprintf(fp1,"\n");}}C语言画图程序//导线网概略图#include <stdio.h>#include <windows.h>#include <math.h>#define NUM 30LRESULT CALLBACK WndProc (HWND, UINT, WPARAM, LPARAM) ;int WINAPI WinMain (HINSTANCE hInstance, HINSTANCE hPrevInstance, PSTR szCmdLine, int iCmdShow){static char szAppName[] = "SineWave" ;HWND hwnd ;MSG msg ;WNDCLASSEX wndclass ;wndclass.cbSize = sizeof (wndclass) ;wndclass.style = CS_HREDRAW | CS_VREDRAW ;wndclass.lpfnWndProc = WndProc ;wndclass.cbClsExtra = 0 ;wndclass.cbWndExtra = 0 ;wndclass.hInstance = hInstance ;wndclass.hIcon = LoadIcon (NULL, IDI_APPLICATION) ;wndclass.hCursor = LoadCursor (NULL, IDC_ARROW) ;wndclass.hbrBackground = (HBRUSH) GetStockObject (WHITE_BRUSH) ; wndclass.lpszMenuName = NULL ;wndclass.lpszClassName = szAppName ;wndclass.hIconSm = LoadIcon (NULL, IDI_APPLICATION) ;RegisterClassEx (&wndclass) ;hwnd = CreateWindow (szAppName, "控制网图",WS_OVERLAPPEDWINDOW,CW_USEDEFAULT, CW_USEDEFAULT,CW_USEDEFAULT, CW_USEDEFAULT,NULL, NULL, hInstance, NULL) ;ShowWindow (hwnd, iCmdShow) ;UpdateWindow (hwnd) ;while (GetMessage (&msg, NULL, 0, 0)){TranslateMessage (&msg) ;DispatchMessage (&msg) ;}return msg.wParam ;}LRESULT CALLBACK WndProc (HWND hwnd, UINT iMsg, WPARAM wParam, LPARAM lParam){static int cxClient, cyClient ;HDC hdc ;int i ;PAINTSTRUCT ps ;POINT pt [NUM] ;switch (iMsg){case WM_SIZE:cxClient = LOWORD (lParam) ;cyClient = HIWORD (lParam) ;return 0 ;case WM_PAINT:hdc = BeginPaint (hwnd, &ps) ;doubley[16]={500.0,363.2211,250.8730,245.2811,343.9816,399.8183,506.0272,596.9804,6 90.1846,686.1801,688.9648,550.4837,552.8629,499.4146,497.8342,500.0}; doublex[16]={500.0,500.0,506.038,269.917,264.999,268.990,346.426,388.125,405.382,466. 024,553.062,555.201,651.049,652.317,591.707,500.0};for (i = 0 ; i < 16 ; i++){pt[i].x = x[i];pt[i].y = y[i];}Polyline (hdc, pt, 16);return 0 ;case WM_DESTROY:PostQuitMessage (0) ;return 0 ;}return DefWindowProc (hwnd, iMsg, wParam, lParam) ;}C#程序：Form1:using System;using System.Collections.Generic;using ponentModel;using System.Data;using System.Drawing;using System.Linq;using System.Text;using System.Threading.Tasks;using System.Windows.Forms;using System.IO;using System.Collections.Generic;namespace导线控制网{public partial class Form1 : Form{public Form1(){InitializeComponent();}public string jisuanjieguo;public const double p = 180*3600/Math.PI;//定义角度转弧度的函数public double jiao_hu(double du,double fen,double miao){return (du*3600.0+fen*60.0+miao)/p;}//定义矩阵运算的类（包括矩阵的加、减、乘、转置和求逆）public class matrix_yusuan{//矩阵相加public static double[,] matrix_jia(double[,] Arry, double[,] Arry1) {int m = Arry.GetLength(0);int n = Arry.GetLength(1);int s = Arry1.GetLength(0);int t = Arry1.GetLength(1);double[,] temp = new double[m, n];double[,] tem = {{0}};if (m == s && n == t){for (int i = 0; i < m; i++){for (int j = 0; j < n; j++){temp[i, j] = Arry[i, j] + Arry1[i, j];}}return temp;}else{Console.WriteLine("两个矩阵大小不同");return tem ;}}//矩阵相减public static double[,] matrix_jian(double[,] Arry, double[,] Arry1) {int m = Arry.GetLength(0);int n = Arry.GetLength(1);int s = Arry1.GetLength(0);int t = Arry1.GetLength(1);double[,] temp = new double[m, n];double[,] tem = { { 0 } };if (m == s && n == t){for (int i = 0; i < m; i++){for (int j = 0; j < n; j++){temp[i, j] = Arry[i, j] - Arry1[i, j];}}return temp;}else{Console.WriteLine("两个矩阵大小不同");return tem;}}//矩阵转置public static double[,] matrix_t(double[,] Arry){int m = Arry.GetLength(0);int n = Arry.GetLength(1);double[,] temp = new double[n, m];for (int i = 0; i < n; i++){for (int j = 0; j < m; j++){temp[i,j] = Arry[j,i];}}return temp;}//矩阵相乘public static double[,] matrix_cheng(double[,] Arry, double[,] Arry1) {int m = Arry.GetLength(0);//矩阵Arry的行数int r = Arry.GetLength(1);//矩阵Arry的列数int k = Arry1.GetLength(0);//矩阵Arry的行数int n = Arry1.GetLength(1);//矩阵Arry1的列数double[,] temp = new double[m, n];double[,] tem = { { 0 } };if (r == k){for (int i = 0; i < m; i++){for (int j = 0; j < n; j++){for (int t = 0; t < r; t++){temp[i, j] += Arry[i, t] * Arry1[t, j];}}}return temp;}else{Console.WriteLine("两个矩阵无法相乘");return tem;}}//求矩阵Arry的逆矩阵public static double[,] matrix_ni(double[,] Arryni){int Level = Arryni.GetLength(1);double[,] NArry = new double[Level, Level];double HLS = matrix_yusuan.matrix_hls(Arryni);double[,] BArry = matrix_bansui(Arryni);for (int i = 0; i < Level; i++){for (int j = 0; j < Level; j++){NArry[i, j] = BArry[i, j] / HLS;}}return NArry;}//求矩阵Arry的伴随矩阵public static double[,] matrix_bansui(double[,] Arryni){int Level = Arryni.GetLength(1);double[,] BArry = new double[Level,Level];for (int m = 0; m < Level; m++){for (int n = 0; n < Level; n++){BArry[m, n] = matrix_yusuan.matrix_yuzi(Arryni, n, m);}}return BArry;}//求矩阵Arry的元素Arry[i,j]的余子式public static double matrix_yuzi(double[,] Arryni, int i, int j)//第i行，第j列，起始值为0{int Level = Arryni.GetLength(1);double[,] Arry1 = new double[Level-1,Level-1];for (int m = 0; m < Level - 1; m++){for (int n = 0; n < Level - 1; n++){if (m < i && n < j){Arry1[m, n] = Arryni[m, n];}else if(m < i && n >= j){Arry1[m, n] = Arryni[m, n + 1];}else if (m >= i && n < j){Arry1[m, n] = Arryni[m + 1, n];}else if (m >= i && n >= j){Arry1[m, n] = Arryni[m + 1, n + 1];}}}//根据矩阵元素的不同位置返回不同的值if ((i + j) % 2 != 0){return (-1)*matrix_yusuan.matrix_hls( Arry1);}else{return matrix_yusuan.matrix_hls(Arry1);}}//求行列式public static double matrix_hls(double[,] Arryni){int Level = Arryni.GetLength(1);//简单来说，对于常用的二维数组，取0或者1，分别获取列和行的长度(行数和列数)double[,] dArry = new double[Level, Level];for (int i = 0; i < Level; i++){for (int j = 0; j < Level; j++){dArry[i, j] = Arryni[i, j];}}int sign = 1;for (int i = 0, j = 0; i < Level && j < Level; i++, j++){//判断该行的正对角元素dArry[i,j]是否为0，若为0，则寻找i行以下的行m（m>i,且dArry[m,j]!=0）进行交换if (dArry[i, j] == 0){//判断是否达到了矩阵的最大行if (i == Level - 1){return 0;}int m = i + 1;//获取一个dArry[m,j]不为零的行for (; dArry[m, j] == 0; m++){if (m == Level - 1){return 0;}}//把i行和m行互换double temp;for (int n = j; n < Level; n++){temp = dArry[i, n];dArry[i, n] = dArry[m, n];dArry[m, n] = temp;}sign *= (-1);}//把当前行以下的行所对应的列变成0double tmp;//for (int s = Level - 1; s > i; s--)for (int s = i + 1; s < Level; s++){tmp = dArry[s, j];for (int t = j; t < Level; t++){dArry[s, t] -= dArry[i, t] * (tmp / dArry[i, j]); }}}double result = 1;for (int i = 0; i < Level; i++){if (dArry[i, i] != 0){result *= dArry[i, i];}else{return 0;}}return sign * result;}}public class aa//用于判断文件导入和计算是否完成{public static int panduan1 = 0;//panduan1用来判断是否进行了计算public static int panduan2 = 0;//panduan2用来判断是否导入文件}int i;int j;//定义数组double[,] hudu = new double[30, 1];double[,] bian = new double[30, 1];double[] alf = new double[30];double[] alfo = new double[30];double[] so = new double[30];double[] Xo = new double[30];double[] Yo = new double[30];double[,] B = new double[30, 27];double[,] L = new double[30, 1];double[,] P = new double[30, 30];double[,] W = new double[27, 1];double[,] temp = new double[30, 30];//此处为Xo Yo，B矩阵赋初值为零 hudu代表夹角，alf方位角，alfo方位角近似值，s 距离观测值，so距离近似值//doubleNbb[50][50],xgu[50][50],Xgu[30][1],Ygu[30][1],V[50][50],sigma_gu,Q[50][50],sigma_xy[50] [50];double[,] Nbb = new double[27, 27];double[,] xgu = new double[27, 1];double[,] xgu28 = new double[28, 1];double[,] Xgu = new double[30, 1];double[,] Ygu = new double[30, 1];double[,] V = new double[30, 1];double[,] Q = new double[27, 27];double[,] sigma_xy = new double[30, 30];double[,] sigma_xy28 = new double[30, 30];double sigma_gu;double D, F, M, sigma_hudu = 2.0, sigma_s;//用来求未知点坐标近似值private void toolStripMenuItem7_Click(object sender, EventArgs e) {if (aa.panduan1 == 1){textBox2.Text = "";textBox2.Text = "近似方位角alfo" + "\t\t" + "\r\n";for (i = 0; i < 15; i++){textBox2.Text += Math.Round(alfo[i], 6);textBox2.Text += "\r\n\r\n";}}else{MessageBox.Show("请先进行计算！", "系统提示");}}private void toolStripMenuItem10_Click(object sender, EventArgs e) {if (aa.panduan1 == 1){textBox2.Text = "";textBox2.Text = "P矩阵" + "\t\t" + "\r\n";for (i = 0; i < P.GetLength(0); i++){for (j = 0; j < P.GetLength(1); j++){textBox2.Text += Math.Round(P[i, j], 6) + "\t";}textBox2.Text += "\r\n\r\n";}}else{MessageBox.Show("请先进行计算！", "系统提示");}}private void toolStripMenuItem9_Click(object sender, EventArgs e) {if (aa.panduan1 == 1){textBox2.Text = "B矩阵" + "\t\t" + "\r\n";for (i = 0; i < B.GetLength(0); i++){for (j = 0; j < B.GetLength(1); j++){textBox2.Text += Math.Round(B[i, j], 6) + "\t";}textBox2.Text += "\r\n\r\n";}}else{MessageBox.Show("请先进行计算！", "系统提示");}}private void button15_Click(object sender, EventArgs e){if (aa.panduan2 == 1){alf[0] = Math.PI;alfo[0] = Math.PI;Xo[0] = 5000.0;Yo[0] = 5000.0;Xo[15] = 5000.0;Yo[15] = 5000.0;so[0] = bian[0, 0];for (i = 1; i < 15; i++){alf[i] = alf[i - 1] + hudu[i, 0] - Math.PI;if (alf[i] >= (2 * Math.PI)){alf[i] = alf[i] - 2 * Math.PI;}Xo[i] = Xo[i - 1] + bian[i - 1, 0] * Math.Cos(alf[i - 1]);Yo[i] = Yo[i - 1] + bian[i - 1, 0] * Math.Sin(alf[i - 1]);}jisuanjieguo += "近似距离(单位：m)：" + "\r" + "近似方位角(单位：弧度)：" + "\r\n";jisuanjieguo += Math.Round(so[0], 3) + "\t\t" + Math.Round(alfo[0], 6) + "\r\n";for (i = 1; i < 15; i++){//求近似距离so[i] = Math.Sqrt((Yo[i + 1] - Yo[i]) * (Yo[i + 1] - Yo[i]) + (Xo[i + 1] - Xo[i]) * (Xo[i + 1] - Xo[i]));//求近似方位角，分象限进行讨论if ((Yo[i + 1] - Yo[i]) > 0 && (Xo[i + 1] - Xo[i]) > 0)//第一象限{alfo[i] = Math.Atan((Yo[i + 1] - Yo[i]) / (Xo[i + 1] - Xo[i])); }else if ((Yo[i + 1] - Yo[i]) > 0 && (Xo[i + 1] - Xo[i]) < 0)//第二象限{alfo[i] = Math.Atan((Yo[i + 1] - Yo[i]) / (Xo[i + 1] - Xo[i])) + Math.PI;}else if ((Yo[i + 1] - Yo[i]) < 0 && (Xo[i + 1] - Xo[i]) > 0)//第三象限{alfo[i] = Math.Atan((Yo[i + 1] - Yo[i]) / (Xo[i + 1] - Xo[i])) + 2 * Math.PI;}else//((Yo[i]-Yo[i-1])<0&&(Xo[i]-Xo[i-1])<0)//第四象限{alfo[i] = Math.Atan((Yo[i + 1] - Yo[i]) / (Xo[i + 1] - Xo[i])) + Math.PI;}jisuanjieguo += Math.Round(so[i], 3) + "\t\t"+ Math.Round(alfo[i], 6) + "\r\n";}//将第一个角度的系数单独算出j = 0;//B[0][j]=0.0;((Yo[14]-Yo[0])/(so[14]*so[14]))*p/1000.0;//B[0][j+1]=0.0;-1*((Xo[14]-Xo[0])/(so[14]*so[14]))*p/1000.0;B[0, j] = -1 * (Yo[1] - Yo[0]) / (so[0] * so[0]) * p / 1000.0;//B[0][j+1]=(Xo[1]-Xo[0])/(so[0]*so[0])*p/1000.0;B[0, 25] = ((Yo[14] - Yo[0]) / (so[14] * so[14])) * p / 1000.0;B[0, 26] = -1 * (Xo[14] - Xo[0]) / (so[14] * so[14]) * p / 1000.0;//将第二个角度的系数单独算出B[1, j] = ((Yo[2] - Yo[1]) / (so[1] * so[1]) - (Yo[0] - Yo[1]) / (so[0] * so[0])) * p / 1000.0;//B[1][j+1]=-1*((Xo[2]-Xo[1])/(so[1]*so[1])-(Xo[0]-Xo[1])/(so[0]*so[0]))*p/1000.0;B[1, j + 1] = -1 * ((Yo[2] - Yo[1]) / (so[1] * so[1])) * p / 1000.0;B[1, j + 2] = ((Xo[2] - Xo[1]) / (so[1] * so[1])) * p / 1000.0;//i = i + 1;//将第三个角度的系数单独算出i = 2;B[i, j] = ((Yo[i - 1] - Yo[i]) / (so[i - 1] * so[i - 1])) * p / 1000.0;//B[i, j + 1] = -1 * ((Xo[i - 1] - Xo[i]) / (so[i - 1] * so[i - 1])) * p / 1000.0;B[i, j + 1] = ((Yo[i + 1] - Yo[i]) / (so[i] * so[i]) - (Yo[i - 1] - Yo[i]) / (so[i - 1] * so[i - 1])) * p / 1000.0;B[i, j + 2] = -1 * ((Xo[i + 1] - Xo[i]) / (so[i] * so[i]) - (Xo[i - 1] - Xo[i]) / (so[i - 1] * so[i - 1])) * p / 1000.0;B[i, j + 3] = -1 * ((Yo[i + 1] - Yo[i]) / (so[i] * so[i])) * p / 1000.0; B[i, j + 4] = (Xo[i + 1] - Xo[i]) / (so[i] * so[i]) * p / 1000.0;//求其他角度改正的系数for (i = 3; i < 15; i++){if (i < 14){B[i, j + 1] = ((Yo[i - 1] - Yo[i]) / (so[i - 1] * so[i - 1])) * p / 1000.0;B[i, j + 2] = -1 * ((Xo[i - 1] - Xo[i]) / (so[i - 1] * so[i - 1])) * p / 1000.0;B[i, j + 3] = ((Yo[i + 1] - Yo[i]) / (so[i] * so[i]) - (Yo[i - 1] - Yo[i]) / (so[i - 1] * so[i - 1])) * p / 1000.0;B[i, j + 4] = -1 * ((Xo[i + 1] - Xo[i]) / (so[i] * so[i]) - (Xo[i - 1] - Xo[i]) / (so[i - 1] * so[i - 1])) * p / 1000.0;B[i, j + 5] = -1 * ((Yo[i + 1] - Yo[i]) / (so[i] * so[i])) * p / 1000.0;B[i, j + 6] = (Xo[i + 1] - Xo[i]) / (so[i] * so[i]) * p / 1000.0; }else{B[i, j + 1] = (Yo[i - 1] - Yo[i]) / (so[i - 1] * so[i - 1]) * p / 1000.0;B[i, j + 2] = -1 * ((Xo[i - 1] - Xo[i]) / (so[i - 1] * so[i - 1])) * p / 1000.0;B[i, j + 3] = ((Yo[i + 1] - Yo[i]) / (so[i] * so[i]) - (Yo[i - 1] - Yo[i]) / (so[i - 1] * so[i - 1])) * p / 1000.0;B[i, j + 4] = -1 * ((Xo[i + 1] - Xo[i]) / (so[i] * so[i]) - (Xo[i - 1] - Xo[i]) / (so[i - 1] * so[i - 1])) * p / 1000.0;}j = j + 2;}//求边长改正的系数j = 0;//将第一个边长的系数单独算出B[i, 0] = (Xo[1] - Xo[0]) / so[0]; //(Xo[1]-Xo[0])/(so[0]*so[0]);i = i + 1;//将第二个边长的系数单独算出B[i, j] = -1 * (Xo[i - 14] - Xo[i - 15]) / so[i - 15];//B[i, j + 1] = -1 * (Yo[i - 14] - Yo[i - 15]) / so[i - 15];B[i, j + 1] = -1 * B[i, j];B[i, j + 2] = (Yo[i - 14] - Yo[i - 15]) / so[i - 15];for (i = i + 1; i < 30; i++){if (i < 29){B[i, j + 1] = -1 * (Xo[i - 14] - Xo[i - 15]) / so[i - 15];B[i, j + 2] = -1 * (Yo[i - 14] - Yo[i - 15]) / so[i - 15];B[i, j + 3] = -1 * B[i, j + 1];B[i, j + 4] = -1 * B[i, j + 2];}else{B[i, j + 1] = -1 * (Xo[i - 14] - Xo[i - 15]) / so[i - 15];B[i, j + 2] = -1 * (Yo[i - 14] - Yo[i - 15]) / so[i - 15];}j = j + 2;}jisuanjieguo += "B矩阵：" + "\r\n";for (i = 0; i < 30; i++){for (j = 0; j < 27; j++){jisuanjieguo += Math.Round(B[i, j], 6) + "\t";}jisuanjieguo += "\r\n";}//求L矩阵，角度和边长分别求解L[0, 0] = (hudu[0, 0] - (alfo[0] - alfo[14] + Math.PI)) * p;for (i = 1; i < 30; i++){if (i < 15){L[i, 0] = (hudu[i, 0] - (alfo[i] - alfo[i - 1] + Math.PI)) * p;if (L[i, 0] > Math.PI * p){L[i, 0] = (L[i, 0] - 2 * Math.PI * p);}else{L[i, 0] = L[i, 0];}}else{L[i, 0] = (bian[i - 15, 0] - so[i - 15]) * 1000.0;}}jisuanjieguo += "L矩阵(单位：秒和mm)：" + "\r\n";for (i = 0; i < 30; i++){for (j = 0; j < 1; j++){jisuanjieguo += Math.Round(L[i, j], 6) + "\t\t";}jisuanjieguo += "\r\n";}//求权阵Pfor (i = 0; i < 30; i++){if (i < 15){P[i, i] = 1;}else{sigma_s = 5 + 10 * 0.000001 * bian[i - 15, 0] * 1000; //单位为（''/mm）的平方P[i, i] = sigma_hudu * sigma_hudu / (sigma_s * sigma_s);}}jisuanjieguo += "P矩阵：" + "\r\n";for (i = 0; i < 30; i++){for (j = 0; j < 30; j++){jisuanjieguo += Math.Round(P[i, j], 6) + "\t";}jisuanjieguo += "\r\n";}Nbb =matrix_yusuan.matrix_cheng(matrix_yusuan.matrix_cheng(matrix_yusuan.matrix_t(B), P), B); //计算Nbb矩阵double[,] temp2 = new double[27, 30];。

常用测量程序设计代码（1）用全站仪在A点观测了B点斜边和垂直角,求A到B(de)高差.（提示：22sin(1)cos2ABDh D a k a i vR=+-+-,D--斜边,a--垂直角,i--仪器高,v--反光镜高,k--大气折光系数）using System;using ;namespace ConsoleApplication1{class Application{static void Main(string[] args){("请输入斜边=");double D = ());("请输入垂直角[]=");double a = DEG()));("请输入仪器高=");double i = ());("请输入反光镜高=");double v = ());double h = D (a) + (1 - D / D / (a) (a) / +i - v;("高差为{0}",h);}MainMainaMaindd(DEG(str)));}else{break;}} while (true);i = 1;foreach (double a in导线转角集合) {a0 += a + ;if (a0 > 2{a0 -= 2 ;}else if (a0 <{a0 += 2 ;}("第{0}条边(de)方位角为{1}",i++,DMS(a0)); }}Maindd(DEG(str)));}else{break;}("请输入第{0}条边长值=", i++);导线边长集合.Add( ()));} while (true);ount; j++){a0 += 导线转角集合[j] + ;if (a0 > 2 a0 -= 2 ;x0 = x0 + 导线边长集合[j] (a0);y0 = y0 + 导线边长集合[j] (a0);("P{0}点(de)坐标是:{1},{2}",j+2,x0,y0); }}Maindd(str));}else{break;}("请输入第{0}段高差=", i++);路线段高差集合.Add()));} while (true);ount; j++){w += 路线段高差集合[j];总长度 += 路线段长度集合[j];}w = w - (Hb - Ha); ount; j++){H += 路线段高差集合[j] + v 路线段长度集合[j]; ("P{0}点(de)高程是:{1}",j+1,H);}("最后一点就是B点");}}}（7）在如图所示(de)前方交会中,ABJ三点按逆时针方向排列,已知AB两点(de)坐标和两个交会角度,求J点(de)坐标.using System;using ;namespace ConsoleApplication1{class Application{static void Main(string[] args) {("请输入A点(de)x坐标=");double Ax = ());("请输入A点(de)y坐标=");double Ay = ());("请输入B点(de)x坐标=");double Bx = ());("请输入B点(de)y坐标=");double By = ());("请输入α(de)角度值=");doubleα = DEG()));("请输入β(de)角度值=");doubleβ = DEG()));1111()2ni i i i i P x y x y ++==-∑Maindd(str));}else{break ;}("请输入第{0}点(de)y 坐标=",i++);多边形Y 坐标集.Add( ())); } while (true );ount;for (int j = 0; j < Count-1; j++){area += 多边形X 坐标集[j] 多边形Y 坐标集[j + 1] - 多边形X 坐标集[j + 1] 多边形Y 坐标集[j];}22sin (1)cos 2AB D h D a k a i v R=+-+-a i v kMainEG()));差(D, α, 仪器高, 反光镜高, out 水平距离);水平距离平方集合.Add(水平距离 水平距离); 高差集合.Add(h);水平距离平方和 += 水平距离 水平距离;高差和 += h;} while (true );ount; j++){H += 高差集合[j] + 单位长度改正数 水平距离平方集合[j];("第{0}点(de)高程H={1}",j+2,H);}("最后一个点(de)高程应该与B 点(de)高程相同"); }}class 导线{iMain位角(MX, MY, AX, AY);List<double> 坐标增量X集合 = new List<double>();List<double> 坐标增量Y集合 = new List<double>();int i = 1;double坐标增量X和 = 0;double坐标增量Y和 = 0;double坐标增量X绝对值和 = 0;double坐标增量Y绝对值和 = 0;double水平距离和 = 0;do{("请输入{0}点到{1}点(de)水平距离S<直接回车结束输入>=",i,i+1);string str = ();if (str == "") break;double S = (str);("请输入{0}点(de)水平角<左角为正,右角为负>=",i);doubleβ = 导线.DEG()));i++;标(0, 0, β, S, α0, out dx, out dy);坐标增量X集合.Add(dx);坐标增量Y集合.Add(dy);坐标增量X和 += dx; ;坐标增量Y和 += dy; ;坐标增量X绝对值和 += (dx);坐标增量Y绝对值和 += (dy);水平距离和 += S;} while (true);ount; j++){X += 坐标增量X集合[j] + X坐标单位长度改正数 (坐标增量X集合[j]);Y += 坐标增量Y集合[j] + Y坐标单位长度改正数 (坐标增量Y集合[j]);("第{0}点(de)X坐标={1},Y坐标={2}",j+2,X,Y);}("最后一个点(de)坐标应该与B点(de)坐标相同");}}class 导线{i β22sin (1)cos 2AB D h D a k a i v R=+-+-a i v kMain 位角(MX, MY, AX, AY);List <double > 坐标增量X 集合 = new List <double >(); List <double > 坐标增量Y 集合 = new List <double >(); List <double > 高差集合 = new List <double >();List <double > 水平距离平方集合 = new List <double >(); int i = 1;double 坐标增量X 和 = 0;double 坐标增量Y 和 = 0;double 坐标增量X 绝对值和 = 0;double 坐标增量Y 绝对值和 = 0;double 水平距离和 = 0;double 高差和 = 0;double 水平距离平方和 = 0;do{("请输入{0}点到{1}点(de)斜距D<直接回车结束输入>=", i, i + 1);string str = ();if (str == "") break;double D = (str);("请输入{0}点(de)水平角<左角为正,右角为负>=", i);doubleβ = 导线.DEG()));("请输入{0}点到{1}点(de)垂直角=", i, i + 1);doubleα = 导线.DEG()));差(D, α, 仪器高, 反光镜高, out水平距离);水平距离平方集合.Add(水平距离水平距离);高差集合.Add(h);水平距离平方和 += 水平距离水平距离;高差和 += h;标(0, 0, β, 水平距离, α0, out dx, out dy);坐标增量X集合.Add(dx);坐标增量Y集合.Add(dy);坐标增量X和 += dx; ;坐标增量Y和 += dy; ;坐标增量X绝对值和 += (dx);坐标增量Y绝对值和 += (dy);水平距离和 += 水平距离;} while (true);ount; j++){H += 高差集合[j] + 单位长度改正数水平距离平方集合[j];X += 坐标增量X集合[j] + X坐标单位长度改正数 (坐标增量X集合[j]);Y += 坐标增量Y集合[j] + Y坐标单位长度改正数 (坐标增量Y集合[j]);("第{0}点(de)X坐标={1},Y坐标={2},高程H={3}",j+2,X,Y,H);}("最后一个点(de)坐标和高程应该与B点(de)坐标和高程相同");}}class导线{i β22sin (1)cos 2AB D h D a k a i v R=+-+-a i v kMain 位角(MX, MY, AX, AY);位角(BX, BY, NX, NY);List <double > 坐标增量X 集合 = new List <double >(); List <double > 坐标增量Y 集合 = new List <double >(); List <double > 高差集合 = new List <double >();List <double > 水平距离平方集合 = new List <double >(); List <double > 方位角集合 = new List <double >(); int i = 1;double 坐标增量X 和 = 0;double 坐标增量Y 和 = 0;double 坐标增量X 绝对值和 = 0;double 坐标增量Y 绝对值和 = 0;double 水平距离和 = 0;double 高差和 = 0;double 水平距离平方和 = 0;double α = α0;do{("请输入{0}点到{1}点(de)斜距D<直接回车结束输入>=", i, i + 1);string str = ();if (str == "") break;double D = (str);("请输入{0}点(de)水平角<左角为正,右角为负>=", i);doubleβ = 导线.DEG()));("请输入{0}点到{1}点(de)垂直角=", i, i + 1);double垂直角 = 导线.DEG()));差(D, 垂直角, 仪器高, 反光镜高, out水平距离);水平距离平方集合.Add(水平距离水平距离);高差集合.Add(h);水平距离平方和 += 水平距离水平距离;高差和 += h;标(0, 0, β, 水平距离, α, out dx, out dy);方位角集合.Add(α);坐标增量X集合.Add(dx);坐标增量Y集合.Add(dy);坐标增量X和 += dx; ;坐标增量Y和 += dy; ;坐标增量Y绝对值和 += (dy);水平距离和 += 水平距离;} while (true);("请输入最后一个连接角<左角为正,右角为负>=");doubleβn = 导线.DEG()));位角(方位角集合[方位角集合.Count-1],βn);ount + 1);lear();坐标增量Y集合.Clear();坐标增量X和 = 0;坐标增量Y和 = 0;坐标增量X绝对值和 = 0;坐标增量Y绝对值和 = 0;for (int j = 0; j < 方位角集合.Count; j++){方位角集合[j] += 方位角改正数 (j + 1); dd(dx);坐标增量Y集合.Add(dy);坐标增量X和 += dx; ;坐标增量Y和 += dy; ;坐标增量X绝对值和 += (dx);}ount; j++){H += 高差集合[j] + 单位长度改正数水平距离平方集合[j];X += 坐标增量X集合[j] + X坐标单位长度改正数 (坐标增量X集合[j]);Y += 坐标增量Y集合[j] + Y坐标单位长度改正数 (坐标增量Y集合[j]);("第{0}点(de)X坐标={1},Y坐标={2},高程H={3}",j+2,X,Y,H);}("最后一个点(de)坐标和高程应该与B点(de)坐标和高程相同");}}class导线{iMainEG()));("请输入B点到N点(de)坐标方位角=");doubleα1 = 导线.DEG()));List<double> 方位角集合 = new List<double>();int i = 1;doubleα = α0;do{("请输入{0}点(de)水平角[直接回车结束输入]<左角为正,右角为负>=", i);string str = ();if (str == "") break;doubleβ = 导线.DEG(str));i++;α = 导线.方位角(α,β);方位角集合.Add(α);} while (true);ount-1] - α1;double方位角改正数 = -方位角闭合差 / 方位角集合.Count;ount; j++){方位角集合[j] += 方位角改正数 (j + 1);("{0}点到{1}点(de)坐标方位角={2}", j, j + 2, 导线.DMS(方位角集合[j]));}("最后一个坐标方位角应该与B点到N(de)坐标方位角相同");}}class导线{iMain位角(MX, MY, AX, AY);位角(BX, BY, NX, NY);List<double> 坐标增量X集合 = new List<double>();List<double> 坐标增量Y集合 = new List<double>();List<double> 水平距离平方集合 = new List<double>();List<double> 方位角集合 = new List<double>();int i = 1;double坐标增量X和 = 0;double坐标增量Y和 = 0;double坐标增量X绝对值和 = 0;double坐标增量Y绝对值和 = 0;double水平距离和 = 0;double水平距离平方和 = 0;doubleα = α0;do{("请输入{0}点到{1}点(de)平距S<直接回车结束输入>=", i, i + 1);string str = ();if (str == "") break;double S = (str);("请输入{0}点(de)水平角<左角为正,右角为负>=", i);doubleβ = 导线.DEG()));i++;水平距离平方集合.Add(S S);水平距离平方和 += S S;标(0, 0, β, S, α, out dx, out dy);方位角集合.Add(α);坐标增量X集合.Add(dx);坐标增量Y集合.Add(dy);坐标增量X和 += dx; ;坐标增量Y和 += dy; ;坐标增量X绝对值和 += (dx);坐标增量Y绝对值和 += (dy);} while (true);("请输入最后一个连接角<左角为正,右角为负>=");doubleβn = 导线.DEG()));位角(方位角集合[方位角集合.Count-1],βn);ount + 1);lear();坐标增量Y集合.Clear();坐标增量X和 = 0;坐标增量Y和 = 0;坐标增量X绝对值和 = 0;坐标增量Y绝对值和 = 0;for (int j = 0; j < 方位角集合.Count; j++){方位角集合[j] += 方位角改正数 (j + 1);dd(dx);坐标增量Y集合.Add(dy);坐标增量X和 += dx; ;坐标增量Y和 += dy; ;坐标增量X绝对值和 += (dx);坐标增量Y绝对值和 += (dy);}ount; j++){X += 坐标增量X集合[j] + X坐标单位长度改正数 (坐标增量X集合[j]);Y += 坐标增量Y集合[j] + Y坐标单位长度改正数 (坐标增量Y集合[j]);("第{0}点(de)X坐标={1},Y坐标={2}",j+2,X,Y);}("最后一个点(de)坐标应该与B点(de)坐标相同");}}class导线{//将转为弧度static public double DEG(double ang){int fuhao = (int)(ang / (ang));ang = (ang);int d = (int)ang;int m = ((int)(ang 100)) - d 100;double s = ang 10000 - m 100 - d 10000;return ((d + m / + s / fuhao) / ;}//计算方位角,返回弧度值public static double方位角(double x1, double y1, double x2, double y2){double deltaX = x2 - x1;double deltaY = y2 - y1;double angle = ;if (deltaX) >{angle = (deltaY, deltaX);}if (angle < 0){angle += ;}if (deltaY <{angle += ;}return angle;}//计算坐标,返回已知点到计算点(de)方位角public static double坐标(double x0, double y0, double左角, double水平距离,double已知方位角,out double x,out double y) {double方位角 = 已知方位角 + 左角 + ;//将方位角调整到0到2π之间if (方位角 >= 2) 方位角 -= 2;if (方位角 < 方位角 += 2;x = x0 + 水平距离 (方位角);y = y0 + 水平距离 (方位角);return方位角;}//根据后视边(de)方位角与左角,计算前进边(de)方位角public static double方位角(double后视边方位角, double左角){double方位角 = 后视边方位角 + 左角 + ; //将方位角调整到0到2π之间if (方位角 >= 2) 方位角 -= 2;if (方位角 < 方位角 += 2;return方位角;}}}。

误差理论与测量平差上机指导书钱建国张恒憬编写辽宁工程技术大学测绘与地理科学学院测绘工程系目录Visual C++平差编程实现 (2)1矩阵加法 (2)2矩阵乘法 (2)3 矩阵转置 (4)4 矩阵求逆 (4)5 水准网间接平差实例（分组选做） (11)Matlab平差编程实现（分组选做） (19)1 间接平差 (19)Visual C++平差编程实现一、实验名称：解算法方程。

二、实验目的和任务：掌握矩阵加法、乘法与求逆的通用程序的编写。

三、实验要求：1每人独立编写出矩阵加法与乘法的程序，并上机调试通过；2采用VC++6.0开发平台，C或者C++语言编写程序；3写出矩阵运算的结果。

四、实验内容：1矩阵加法矩阵加法的示例函数(C语言)void JZjiafa(double a[15][15],double b[15][15],double c[15][15],intm,int n){for (int i=0;i<=m-1;i++)for(int j=0;j<=n-1;j++){c[i][j]=a[i][j]+b[i][j];}return;}2矩阵乘法矩阵乘法的示例程序(C语言)#include "stdafx.h"void matrixMultiply(double a[14][15],double b[15][13], doublec[14][13],long m,long n,long k){for (long i = 0; i<= m-1; i++){for (long j=0; j<=k-1; j++){c[i][j] =0.0;for (long q=0; q<=n-1;q++){c[i][j] = c[i][j] + a[i][q] * b[q][j];}}}return;}int main(int argc, char* argv[]){long n,m,k,i,j;double a[14][15],c[14][13],b[15][13];FILE *stream;stream = fopen("矩阵输入.txt","r");fscanf(stream,"%ld %ld",&n,&m);for (i=0;i<n;i++){for(j=0;j<m;j++){fscanf(stream,"%lf",&a[i][j]);}}fscanf(stream,"%ld %ld",&m,&k);for(i=0;i<m;i++){for(j=0;j<k;j++){fscanf(stream,"%lf",&b[i][j]);}}fclose(stream);matrixMultiply(a,b,c,4,5,3);stream = fopen("矩阵计算结果.txt","w");for (i=0;i<=3;i++){for(j=0;j<=2;j++)fprintf(stream,"%16.7e ",c[i][j]);fprintf(stream,"\n");}fprintf(stream,"\n");fclose(stream);return 0;}3 矩阵转置矩阵的转置示例函数(C语言)double JZzhuanzhi(double a[15][15], double b[15][15], int m,int n) {{for(int i=0;i<m;i++)for(int j=0;j<n;j++)b[j][i]=a[i][j];}return 0.0;}4 矩阵求逆矩阵求逆的示例函数(C语言)int invGJ(double **a,int n){int *is,*js,i,j,k,l,u,v;double d,p;is=(int *)malloc(n*sizeof(int));js=(int *)malloc(n*sizeof(int));for(k=0;k<=n-1;k++){d=0.0;for(i=k;i<=n-1;i++)for(j=k;j<=n-1;j++){l=i*n+j;p=fabs(a[i][j]);if(p>d){d=p;is[k]=i;js[k]=j;}}if(d+1.0==1.0){free(is);free(js);printf("error not inv\n");return (0);}if(is[k]!=k)for(j=0;j<=n-1;j++){u=k*n+j;v=is[k]*n+j;p=a[k][j];a[k][j]=a[is[k]][j];a[is[k]][j]=p;}if(js[k]!=k)for(i=0;i<=n-1;i++){u=i*n+k;v=i*n+js[k];p=a[i][k];a[i][k]=a[i][js[k]];a[i][js[k]]=p;}l=k*n+k;a[k][k]=1.0/a[k][k];for(j=0;j<=n-1;j++)if(j!=k){u=k*n+j;a[k][j]=a[k][j]*a[k][k];}for(i=0;i<=n-1;i++)if(i!=k)for(j=0;j<=n-1;j++)if(j!=k){u=i*n+j;a[i][j]=a[i][j]-a[i][k]*a[k][j];}for(i=0;i<=n-1;i++)if(i!=k){u=i*n+k;a[i][k]=-a[i][k]*a[k][k];}}for(k=n-1;k>=0;k--){if(js[k]!=k)for(j=0;j<=n-1;j++){u=k*n+j;v=js[k]*n+j;p=a[k][j];a[k][j]=a[js[k]][j];a[js[k]][j]=p;}if(is[k]!=k)for(i=0;i<=n-1;i++){u=i*n+k;v=i*n+is[k];p=a[i][k];a[i][k]=a[i][is[k]];a[i][is[k]]=p;}}free(is);free(js);return (1);} int invGJ(double **a,int n){int *is,*js,i,j,k,l,u,v;double d,p;is=(int *)malloc(n*sizeof(int));js=(int *)malloc(n*sizeof(int));for(k=0;k<=n-1;k++){d=0.0;for(i=k;i<=n-1;i++)for(j=k;j<=n-1;j++){l=i*n+j;p=fabs(a[i][j]);if(p>d){d=p;is[k]=i;js[k]=j;}}if(d+1.0==1.0){free(is);free(js);printf("error not inv\n");return (0);}if(is[k]!=k)for(j=0;j<=n-1;j++){u=k*n+j;v=is[k]*n+j;p=a[k][j];a[k][j]=a[is[k]][j];a[is[k]][j]=p;}if(js[k]!=k)for(i=0;i<=n-1;i++){u=i*n+k;v=i*n+js[k];p=a[i][k];a[i][k]=a[i][js[k]];a[i][js[k]]=p;}l=k*n+k;a[k][k]=1.0/a[k][k];for(j=0;j<=n-1;j++)if(j!=k){u=k*n+j;a[k][j]=a[k][j]*a[k][k];}for(i=0;i<=n-1;i++)if(i!=k)for(j=0;j<=n-1;j++)if(j!=k){u=i*n+j;a[i][j]=a[i][j]-a[i][k]*a[k][j];}for(i=0;i<=n-1;i++)if(i!=k){u=i*n+k;a[i][k]=-a[i][k]*a[k][k];}}for(k=n-1;k>=0;k--){if(js[k]!=k)for(j=0;j<=n-1;j++){u=k*n+j;v=js[k]*n+j;p=a[k][j];a[k][j]=a[js[k]][j];a[js[k]][j]=p;}if(is[k]!=k)for(i=0;i<=n-1;i++){u=i*n+k;v=i*n+is[k];p=a[i][k];a[i][k]=a[i][is[k]];a[i][is[k]]=p;}}free(is);free(js);return (1);}矩阵求逆函数的调用(C语言)#include <stdio.h>#include <stdlib.h>#include <math.h>int invGJ(double **a,int n);void main(){int i,j;double **AA;//首先对二维指针Naa分配内存，采用C语言的方法/* AA=(double **)malloc(sizeof(double)*2);for(i=0;i<2;i++){AA[i]=(double *)mallo(sizeof(double)*2);}*///首先对二维指针Naa分配内存，采用C++语言的方法AA=new double * [2];for(i=0;i<2;i++){AA[i]=new double[2];}double BB[2][2]={1,2,3,4};for(i=0;i<2;i++){for(j=0;j<2;j++){AA[i][j]=BB[i][j];}}//调用矩阵求逆函数invGJ(AA,2);printf("矩阵AA的逆阵如下\n");for(i=0;i<2;i++){for(j=0;j<2;j++){printf("%10.4lf",AA[i][j]);}printf("\n");}double CC[2][2];printf("AA与其逆阵的乘积如下(理论上是单位阵)\n"); for(i=0;i<2;i++){for(j=0;j<2;j++){CC[i][j]=0.0;for(int k=0;k<2;k++){CC[i][j]+=AA[i][k]*BB[k][j];}printf("%10.4lf",CC[i][j]);}printf("\n");}//C 语言释放AA 二维指针的方法 /* for(i=0;i<2;i++){free(AA[i]);}free(AA);*/ //C++语言释放AA 二维指针的方法 for(i=0;i<2;i++) { delete AA[i]; } delete AA;}5 水准网间接平差实例（分组选做）例1：在图1所示的水准网中，已知水准点A 的高程为H A =237.483，为求B 、C 、D 三点的高程，进行了水准测量，测得高差5×1L和水准路线的长度5×1S ，其结果见表1，试按间接平差求定B 、C 、D 三点的高程平差值。

Windows绘图#include "stdafx.h"#include "cbs.h"#include "cbsDoc.h"#include "cbsView.h"#ifdef _DEBUG#define new DEBUG_NEW#undef THIS_FILEstatic char THIS_FILE[] = __FILE__;#endif/////////////////////////////////////////// //////////////////////////////////// CCbsViewIMPLEMENT_DYNCREATE(CCbsView, CView)BEGIN_MESSAGE_MAP(CCbsView, CView)//{{AFX_MSG_MAP(CCbsView)ON_WM_MOUSEMOVE()ON_WM_LBUTTONDOWN()ON_WM_LBUTTONUP()//}}AFX_MSG_MAP// Standard printing commandsON_COMMAND(ID_FILE_PRINT,CView::OnFilePrint)ON_COMMAND(ID_FILE_PRINT_DIRECT, CView::OnFilePrint)ON_COMMAND(ID_FILE_PRINT_PREVIEW, CView::OnFilePrintPreview)END_MESSAGE_MAP()/////////////////////////////////////////// //////////////////////////////////// CCbsView construction/destruction CCbsView::CCbsView(){// TODO: add construction code herethis->m_bDragging=false;}CCbsView::~CCbsView(){}BOOL CCbsView::PreCreateWindow(CREATESTRUCT& cs){// TODO: Modify the Window class or styles here by modifying// the CREATESTRUCT csreturn CView::PreCreateWindow(cs);}/////////////////////////////////////////// //////////////////////////////////// CCbsView drawingvoid CCbsView::OnDraw(CDC* pDC){CCbsDoc* pDoc = GetDocument();ASSERT_VALID(pDoc);// TODO: add draw code for native data here /* pDC->SetMapMode(MM_TEXT);pDC->Rectangle(CRect(50, 50, 100, 100));// 直接基于屏幕绘制pDC->SetMapMode(MM_TEXT);//pDC->SetWindowOrg(50, 50); //屏幕左上角的坐标设置为(50, 50)pDC->Rectangle(CRect(50, 50, 100, 100));pDC->SetMapMode(MM_TEXT);pDC->SetViewportOrg(50,50);//当前原点位置移动到(50, 50)的位置pDC->Rectangle(CRect(50, 50, 100, 100));pDC->SetMapMode(MM_TEXT);pDC->SetViewportOrg(100,100); // 当前原点位置移动到(50, 50)的位置pDC->SetWindowOrg(50, 50); // 当前的原点坐标设置为(50, 50)pDC->Rectangle(CRect(100, 100, 100, 100));pDC->SetWindowOrg(100,100);*/pDC->SetViewportOrg(150,150); // 当前原点位置移动到(50, 50)的位置COLORREF rgbBkClr = RGB(192,192,192);// 定义灰色pDC->SetBkColor(rgbBkClr);// 背景色为灰色pDC->SetTextColor(RGB(0,0,128));// 文本颜色为兰色pDC->TextOut(250,10,"我*****的*****绘*****图*****板！");CPen *pPenOld, PenNew;int nPenStyle[]= { PS_SOLID, // 实线PS_DOT,// 点线PS_DASH,// 虚线PS_DASHDOT,// 点划线PS_DASHDOTDOT, // 双点划线PS_NULL,// 空的边框PS_INSIDEFRAME,// 边框实线};char *strStyle[]={"Solid", "Dot", "Dash", "DashDot","DashDotDot", "Null", "InsideFrame"};pDC->TextOut(50, 65, "用不同样式的画笔绘图"); pDC->TextOut(400, 40, "***制作者***陈斌生***");for(int i=0; i<7; i++){ // 用不同样式的画笔绘图if(PenNew.CreatePen(nPenStyle[i],1,RGB(0,0, 0))){ //创建画笔pPenOld=pDC->SelectObject(&PenNew);// 选择画笔pDC->TextOut(45,90+20*i,strStyle[i]);pDC->MoveTo(150,100+20*i);pDC->LineTo(200,100+20*i);pDC->SelectObject(pPenOld);// 恢复原来的画笔PenNew.DeleteObject();// 删除底层的GDI对象}else {MessageBox("不能创建画笔！"); } }char *strWidth[]={"1", "2", "3", "4", "5", "6", "7"};pDC->TextOut(260, 65, "用不同宽度的画笔绘图");for(i=0; i<7; i++){ // 用不同宽度的画笔绘图if(PenNew.CreatePen(PS_SOLID,i+1,RGB(0, 0,0))){ // 创建画笔pPenOld=pDC->SelectObject(&PenNew); // 选择画笔pDC->TextOut(260,90+20*i,strWidth[i]);pDC->MoveTo(300, 100+20*i);pDC->LineTo(400, 100+20*i);pDC->SelectObject(pPenOld);// 恢复原来的画笔PenNew.DeleteObject();// 删除底层的GDI对象}else { MessageBox("不能创建画笔！"); }}char *strColor[]={"红","绿","蓝","黄","紫","青","灰"};COLORREF rgbPenClr[]={RGB(255,0,0), RGB(0,255,0),RGB(0,0,255), RGB(255,255,0), RGB(255,0,255),RGB(0,255,255),RGB(192,192,192)};pDC->TextOut(460,65,"用不同颜色的画笔绘图");for(i=0; i<7; i++){ // 用不同颜色的画笔绘图CPen *pPenNew=new CPen(PS_INSIDEFRAME,10,rgbPenClr[i]);// 创建画笔的另一种方法pPenOld=pDC->SelectObject(pPenNew);// 选择创建的画笔pDC->TextOut(460,90+20*i, strColor[i]);pDC->MoveTo(500,100+20*i);pDC->LineTo(600,100+20*i);pDC->SelectObject(pPenOld); // 恢复原来的画笔delete pPenNew; // 自动删除底层的GDI对象}}/////////////////////////////////////////// //////////////////////////////////// CCbsView printingBOOL CCbsView::OnPreparePrinting(CPrintInfo* pInfo){// default preparationreturn DoPreparePrinting(pInfo);} void CCbsView::OnBeginPrinting(CDC* /*pDC*/, CPrintInfo* /*pInfo*/){// TODO: add extra initialization before printing}void CCbsView::OnEndPrinting(CDC* /*pDC*/, CPrintInfo* /*pInfo*/){// TODO: add cleanup after printing}/////////////////////////////////////////// //////////////////////////////////// CCbsView diagnostics#ifdef _DEBUGvoid CCbsView::AssertValid() const{CView::AssertValid();}void CCbsView::Dump(CDumpContext& dc) const {CView::Dump(dc);}CCbsDoc* CCbsView::GetDocument() // non-debug version is inline{ASSERT(m_pDocument->IsKindOf(RUNTIME_CL ASS(CCbsDoc)));return (CCbsDoc*)m_pDocument;}#endif //_DEBUG/////////////////////////////////////////// //////////////////////////////////// CCbsView message handlersvoid CCbsView::OnMouseMove(UINT nFlags, CPoint point)// TODO: Add your message handler code here and/or call defaultif(m_bDragging){CClientDC dc(this);OnPrepareDC(&dc);dc.DPtoLP(&point);dc.MoveTo(this->m_ptOrigin);dc.LineTo(point);m_ptOrigin = point;}CView::OnMouseMove(nFlags, point);}void CCbsView::OnLButtonDown(UINT nFlags, CPoint point){// TODO: Add your message handler code here and/or call defaultCClientDC dc(this);OnPrepareDC(&dc); // 调整设备环境的属性dc.DPtoLP(&point); // 将设备坐标转换为逻辑坐标SetCapture(); // 捕捉鼠标//::SetCursor(m_hCross); // 设置十字光标m_ptOrigin=point;m_bDragging=TRUE; // 设置拖拽标记CView::OnLButtonDown(nFlags, point);}void CCbsView::OnLButtonUp(UINT nFlags, CPoint point){// TODO: Add your message handler code here and/or call defaultm_bDragging = false;ReleaseCapture();CView::OnLButtonUp(nFlags, point);}结果水准网平差结果#include<iostream.h>#include<stdlib.h>#include<iomanip.h>#include<math.h>#define max 50class CMatrix{public:CMatrix(){row=0; column=0;}; // 默认构造函数CMatrix(int i, int j){row=i;column=j;} // 构造函数一CMatrix(const CMatrix& m); // 复制构造函数~CMatrix(void){/*cout<<"谢谢使用，矩阵所占空间以释放！"<<endl;*/} // 默认析构函数CMatrix& operator=(const CMatrix& m); // 赋值运算符bool operator==(const CMatrix& m); // 比括较运算符bool operator!=(const CMatrix& m); // 比括较运算符CMatrix operator+(const CMatrix& m); // 加运算符CMatrix operator-(const CMatrix& m); // 减运算符CMatrix& operator+=(const CMatrix& m); // 自加运算符CMatrix& operator-=(const CMatrix& m); // 自减运算符CMatrix operator-();// 取负数CMatrix& operator*(const CMatrix& m); // 乘法运算符void input(); //输入矩阵void outputMatrix(); // 输出该矩阵void setValue(int row, int column, double value) { A[row-1][column-1] = value; }// 设置(i,j)的值double getValue(int row, int column) const { return A[row-1][column-1]; }// 设置行、列的值void setRow(const int i) { row = i; }int getRow() const { return row; }void setColunm(const int i) { column = i; }int getColumn() const { return column; }CMatrix& change(int i, int j);//交换矩阵的行CMatrix& transpose(); // 矩阵转置CMatrix& inverse(); // 矩阵求逆void find(int& f)const;// 判断该矩阵是否可用于迭代求解friend void jocabi(const CMatrix& a) ; //迭代求解void lzys(); //列主元素法求解void solve(); //可逆线性矩阵求解void qxnh(); //曲线拟合private:// 成员变量double A[max][max];int row;// 行int column;// 列};void CMatrix::input() //输入{ cout<<"开始输入矩阵值："<<endl;int i, j;double z;for(i=0;i<row;i++){ cout<<"请输入第"<<i+1<<"行的值："<<endl;for(j=0;j<column;j++){ cin>>z;A[i][j]=z;}}cout<<endl;};CMatrix::CMatrix(const CMatrix& m) // 复制构造函数{ int i, j;for(i=0;i<m.row;i++)for(j=0;j<m.column;j++)this->A[i][j]=m.A[i][j];};CMatrix& CMatrix::operator=(const CMatrix& m) // 赋值运算符{int i,j;for(i=0;i<row;i++){for(j=0;j<column;j++)A[i][j]=m.A[i][j];}return *this;};bool CMatrix::operator ==(const CMatrix& m) // 比括较运算符{ int i,j,k;for(i=0;i<m.row;i++){ for(j=0;j<m.column;j++)if(this->A[i][j]=m.A[i][j])k=1;else k=0;}if(k=1) return true;else return false;};bool CMatrix::operator !=(const CMatrix& m) // 比括较运算符{ int i,j,k;for(i=0;i<m.row;i++){ for(j=0;j<m.column;j++)if(this->A[i][j]=m.A[i][j])k=1;else k=0;}if(k=0) return true;else return false;};CMatrix CMatrix::operator+(const CMatrix& m)// 加运算符{ int i,j;if((this->row==m.row)&&(this->column==m .column)){ for(i=0;i<m.row;i++)for(j=0;j<m.column;j++)this->A[i][j]+=m.A[i][j];}else { cout<<"此两矩阵不能相加，请检查！"<<endl; }return *this;};CMatrix CMatrix::operator-(const CMatrix& m)// 减运算符{ int i,j;if((this->row==m.row)&&(this->column==m.column)){ for(i=0;i<m.row;i++)for(j=0;j<m.column;j++)this->A[i][j]-=m.A[i][j];}else { cout<<"此两矩阵不能相加，请检查！"<<endl; }return *this;};CMatrix& CMatrix::operator+=(const CMatrix& m) //自加运算符{ int i,j;for(i=0;i<m.row;i++)for(j=0;j<m.column;j++)this->A[i][j]=2*m.A[i][j];return *this;};CMatrix& CMatrix::operator-=(const CMatrix& m) //自减运算符{ int i,j;for(i=0;i<m.row;i++)for(j=0;j<m.column;j++)this->A[i][j]=m.A[i][j]-m.A[i][j];return *this;};void CMatrix::find(int& f)const{ int i;for(i=0;i<this->row;i++) if(this->A[i][i]!=0) f=1;else f=0;};CMatrix CMatrix::operator-() // 取负数{ int i,j;for(i=0;i<this->row;i++)for(j=0;j<this->column;j++)this->A[i][j]=-this->A[i][j];return *this;};CMatrix& CMatrix::operator*(const CMatrix& m) // 乘法运算符{ int i,j,t;CMatrix n;if(this->column==m.row){for(i=0;i<this->row;i++)for(j=0;j<m.column;j++){ double sum=0.0;for(t=0;t<m.row;t++)sum+=this->A[i][t]*m.A[t][j];n.A[i][j]=sum;}}else {cerr<<"此两矩阵不能相乘，请检查！"<<endl; exit(1);}return n;};void CMatrix::outputMatrix()// 输出该矩阵{ int i,j;for(i=1;i<=row;i++){ for(j=1;j<=column;j++)cout<<A[i-1][j-1]<<" ";cout<<endl;}};CMatrix& CMatrix::transpose()// 矩阵转置{ int i,j;CMatrix m(this->column,this->row);for(i=0;i<this->row;i++)for(j=0;j<this->column;j++){ m.A[j][i]=this->A[i][j];}return m;};void jocabi(const CMatrix& a) //高斯迭代求解{ int f=1;a.find(f);if(f==0) {cerr<<"该矩阵不满足迭代求解条件！"<<endl; exit(1); }else{CMatrix x,w;x.setColunm(1);x.setRow(a.getColumn());w.setColunm(1);w.setRow(a.getColumn());int i;double z;for(i=1;i<=x.row;i++){ cout<<"请输入等式右侧的第"<<i<<"个值：";cin>>z;w.A[i-1][0]=z;}for(i=1;i<=x.row;i++){ cout<<"请输入X"<<i<<"的近似值：";cin>>z;x.setValue(i, 1, z);}i=1;while(i<=20){ int j, k;for(j=1;j<=x.row;j++){ double sum=0.0;for(k=1;k<j;k++){sum=sum-(a.A[j-1][k-1])*(x.A[k-1][0]);}for(k=j+1;k<=x.row;k++){sum=sum-(a.A[j-1][k-1])*(x.A[k-1][0]); }sum+=w.A[j-1][0];x.A[j-1][0]=sum/a.A[j-1][j-1];}i++;}for(i=1;i<=x.row;i++)cout<<"X"<<i<<" = "<<x.A[i-1][0]<<endl; }};double jdz(double a) //绝对值{ if(a>0.0) return a;else return -a;};CMatrix& CMatrix::change(int i, int j)//交换矩阵的行{ int k;double z;for(k=1;k<=this->column;k++){ z=A[i-1][k-1];A[i-1][k-1]=A[j-1][k-1];A[j-1][k-1]=z;}return *this;};void CMatrix::lzys() //列主元素法求解{ CMatrix x,w;x.setColunm(1);x.setRow(getColumn());w.setColunm(1);w.setRow(getColumn());int i;double z;for(i=0;i<row;i++){for(int j=0;j<row; j++)w.A[i][j]=0.0;}for(i=1;i<=x.row;i++)x.setValue(i, 1, 0.0);for(i=1;i<=x.row;i++){ cout<<"请输入等式右侧的第"<<i<<"个值：";cin>>z;w.A[i-1][0]=z;}i=0;while(i<x.row-1){ int j,t,h=i;for(j=i;j<x.row-1;j++){ if(jdz(A[j][i])<jdz(A[j+1][i]))h=j+1; }this->change(i+1,h+1);w.change(i+1,1);for(j=i+1;j<row;j++){ double k;k=A[j][i]/A[i][i];for(t=i;t<column;t++){A[j][t]=A[j][t]-k*A[i][t];}w.A[j][0]=w.A[j][0]-k*w.A[i][0];}i++;}if(this->A[row-1][column-1]==0){cerr<<"此矩阵对应的方程组有无穷解！"<<endl; exit(1);}i=column-1;while(i>0){ int j, t;for(j=i-1;j>=0;j--){ double k;k=A[j][i]/A[i][i];for(t=i;t>j;t--){ A[j][t]=A[j][t]-k*A[i][t]; }w.A[j][0]=w.A[j][0]-k*w.A[i][0];}i--;}int j;for(j=0;j<x.row;j++){x.A[j][0]=w.A[j][0]/this->A[j][j];}for(i=1;i<=x.row;i++)cout<<"X"<<i<<" = "<<x.A[i-1][0]<<endl; };CMatrix& CMatrix::inverse() //矩阵求逆{ if(this->row!=this->column) {cerr<<"该矩阵不符合求逆条件！"<<endl; exit(1);}else{CMatrix w;w.setColunm(this->getRow());w.setRow(getColumn());int i;for(i=0;i<row;i++){for(int j=0;j<row; j++){w.A[i][j]=0.0;}w.A[i][i]=1;}i=0;while(i<column-1){ int j,t,h=i;for(j=i;j<row-1;j++){if(jdz(A[j][i])<jdz(A[j+1][i]))h=j+1;}this->change(i+1,h+1);w.change(i+1,h+1);for(j=i+1;j<row;j++){ double k;k=A[j][i]/A[i][i];for(t=i;t<column;t++){A[j][t]=A[j][t]-k*A[i][t];w.A[j][t]=w.A[j][t]-k*w.A[i][t];}}i++;}if(this->A[row-1][column-1]==0){cerr<<"此矩阵求逆不成功，其所对应的方程组有无穷解！"<<endl; exit(1);} i=column-1;while(i>0){ int j, t;for(j=i-1;j>=0;j--){ double k;k=A[j][i]/A[i][i];for(t=i;t>j;t--){ A[j][t]=A[j][t]-k*A[i][t]; }for(t=column-1;t>=0;t--){w.A[j][t]=w.A[j][t]-k*w.A[i][t];}}i--;}int j,k;for(j=0;j<row;j++){ for(k=0;k<w.column;k++) {w.A[j][k]=w.A[j][k]/this->A[j][j];}this->A[j][j]=this->A[j][j]/this->A[j][j];}return w;}};void CMatrix::solve() //可逆线性矩阵求解{CMatrix x,w,c;x.setColunm(1);x.setRow(getColumn());w.setColunm(1);w.setRow(getColumn());c.setColunm(getColumn());c.setRow(getRow());int i;double z;for(i=1;i<=x.row;i++){ cout<<"请输入等式右侧的第"<<i<<"个值：";cin>>z;w.A[i-1][0]=z;}c.operator =(this->inverse());x.operator =(c.operator *(w));cout<<"求解结果："<<endl;for(i=1;i<=x.row;i++)cout<<"X"<<i<<" = "<<x.A[i-1][0]<<endl; };void CMatrix::qxnh() //曲线拟合{cout<<"用矩阵曲线拟合："<<endl;int i, j,k,t;double x, y;cout<<"请输入已知坐标点的个数：";cin>>i;CMatrix a(i,1), g(i,i) , w(i,1),c(i,i),q(i,i),p(i,i);for(j=1;j<=i;j++){ cout<<"请输入第"<<j<<"个已知点坐标(先横坐标后纵坐标)：";cin>>x>>y;w.A[j-1][0]=y;g.A[j-1][0]=1.0;for(k=1;k<i;k++){ g.setValue(j,k+1,1.0);for(t=1;t<=k;t++)g.A[j-1][k]*=x;}}c=g.transpose();q=c*g;p=q.inverse();w=c*w;a=p*w;cout<<"Y="<<a.A[0][0];for(j=1;j<i;j++){ cout<<'+'<<'('<<a.A[j][0]<<')';for(k=0;k<j;k++)cout<<"*x";}cout<<endl;cout<<" 1 .求取其它点坐标 0. 退出 "<<endl; while(1){ cout<<"请选择：";cin>>k;switch(k){ case 1:{ double x, y,sum;cout<<"请输入待求点的横坐标（x）：";cin>>x;y=a.A[0][0];for(j=1;j<i;j++){ sum=a.A[j][0];for(k=0;k<j;k++){sum*=x;}y+=sum;}cout<<"纵坐标：y="<<y<<endl; break;}case 0:{cerr<<"谢谢使用！"<<endl; exit(1);break;}default: cout<<"选择有误，请检查！"<<endl; break;}}};struct CElvDif{double value; // 观测值double weight; // 权重long startPoint; // 起始点编号long endPoint; // 终点编号};//水准点类的设计struct CLevelPoint{ // 高程平差值=高程值+高程值改正数long index; // 水准点编号double eleValue; // 高程值double dv; // 高程值改正数（初始化为0）bool isKnown; // 是否为已知点};class CElevationNet{public: // 成员函数CElevationNet(){numElvDif=0;numPoints=0 ;numKnPoint=0;}; // 构造函数~CElevationNet(){}; // 析构函数 void input();void output();void jsgc();int getgz(){return numElvDif;} //返回观测值数目int getzs(){return numPoints;} //返回总点数int getys(){return numKnPoint;} //返回已知点数int getws(){return numPoints-numKnPoint;} //返回未知点数 double geteleValue(int i){return this->lpVec[i-1].eleValue;} //获得高程值 long getindex(int i){return lpVec[i-1].index;} //返回高程点编号void seteleValue(int i, double value){lpVec[i-1].eleValue+=value;} //修改高程值void setdv(int i,double value){this->lpVec[i-1].dv=value;} //修改改正数double getdv(int i){return lpVec[i-1].dv;} //返回改正数friend void xishu(CMatrix& B, CMatrix& X,CElevationNet A); //求取系数矩阵和未知点高程矩阵friend void quanzhen(CMatrix& Q, CElevationNet A); //求取权阵friend void l_zhen( CMatrix& l ,CMatrix& L, CElevationNet A); //求取观测值L矩阵和l阵private: // 成员变量int numElvDif; // 高差总数int numPoints; // 控制网中点的总数目int numKnPoint; //控制网中已知点的数目CElvDif edVec[max]; // 观测值数组CLevelPoint lpVec[max]; // 高程值数组};void CElevationNet::input(){ int a,b,c,i;double z,l,g;cout<<"请对已知点和未知点进行编号，注意已知点的编号应小于未知点的编号！"<<endl<<endl; cout<<"请输入控制网中水准点总数和已知点数：";cin>>a>>b;numPoints=a;numKnPoint=b;for(i=0;i<numPoints;i++){ lpVec[i].index=i+1;lpVec[i].eleValue=0.0;lpVec[i].dv=0.0;lpVec[i].isKnown=0;}for(i=0;i<numKnPoint;i++){ cout<<"请输入第"<<i+1<<"个已知点高程值：";cin>>z;lpVec[i].eleValue=z;lpVec[i].isKnown=1;}cout<<"请输入观测值个数：";cin>>c;numElvDif=c;for(i=0;i<numElvDif;i++){ cout<<"请输入第"<<i+1<<"段观测段的高差值(单位米)、长度(单位千米)、起始点编号和终点编号："<<endl;cin>>g>>l>>a>>b;edVec[i].value=g;edVec[i].weight=l;edVec[i].startPoint=a;edVec[i].endPoint=b;}};void CElevationNet::output(){int i;for(i=0;i<this->numPoints;i++)cout<<"编号为"<<this->lpVec[i].index<<"的点的高程值为："<<this->lpVec[i].eleValue<<endl;for(i=0;i<this->numElvDif;i++)cout<<"观测段"<<i+1<<"的平差后的值为："<<this->edVec[i].value<<endl;};void CElevationNet::jsgc(){ int i,j;for(i=0;i<this->numElvDif;i++){ if((lpVec[edVec[i].startPoint-1].eleValu e!=0.0)&&(lpVec[edVec[i].endPoint-1].eleVal ue==0.0)){lpVec[edVec[i].endPoint-1].eleValue=(l pVec[edVec[i].startPoint-1].eleValue)+(edVe c[i].value);}elseif(lpVec[edVec[i].startPoint-1].eleValue==0 .0&&lpVec[edVec[i].endPoint-1].eleValue!=0.0){lpVec[edVec[i].startPoint-1].eleValue= (lpVec[edVec[i].endPoint-1].eleValue)-(edVe c[i].value);}}for(j=numKnPoint;j<this->numPoints;j++) {if(this->lpVec[j].eleValue==0.0){for(i=0;i<this->numElvDif;i++){ if(this->edVec[i].startPoint==lpVec[ j].index&&this->lpVec[this->edVec[i].endPoi nt-1].eleValue!=0.0)this->lpVec[j].eleValue=this->lpVec[thi s->edVec[i].endPoint-1].eleValue-this->edVe c[i].value;elseif(this->edVec[i].endPoint==lpVec[j].index& &this->lpVec[this->edVec[i].startPoint-1].e leValue!=0.0){this->lpVec[j].eleValue=this->lpVec[this-> edVec[i].startPoint-1].eleValue+this->edVec [i].value;}}}//cout<<"第"<<this->lpVec[j].index<<"个点的高程近似值："<<this->lpVec[j].eleValue<<endl;}};void xishu(CMatrix& B, CMatrix& X,CElevationNet A){int i,j;for(i=1;i<=B.getRow();i++)for(j=1;j<=B.getColumn();j++){ B.setValue(i,j,0);}for(i=0;i<A.numElvDif;i++) { if(A.lpVec[A.edVec[i].startPoint-1].isKn own==1&&A.lpVec[A.edVec[i].endPoint-1].isKn own==0){B.setValue(i+1,A.edVec[i].endPoint-A.numKn Point,1);}elseif(A.lpVec[A.edVec[i].startPoint-1].isKnown ==0&&A.lpVec[A.edVec[i].endPoint-1].isKnown ==1){ B.setValue(i+1,A.edVec[i].startPoint -A.numKnPoint,-1);}elseif(A.lpVec[A.edVec[i].startPoint-1].isKnown ==0&&A.lpVec[A.edVec[i].endPoint-1].isKnown ==0){ B.setValue(i+1,A.edVec[i].endPoint-A. numKnPoint,1);B.setValue(i+1,A.edVec[i].startPoint-A. numKnPoint,-1);}}for(i=0;i<X.getRow();i++)X.setValue(i+1,1,A.lpVec[A.numKnPoint+i ].eleValue);};void quanzhen(CMatrix& Q, CElevationNet A) {int i,j;double sum=0.0;for(i=0;i<A.numElvDif;i++)sum+=A.edVec[i].weight;sum=sum/A.numElvDif;for(i=1;i<=Q.getRow();i++)for(j=1;j<=Q.getColumn();j++){ Q.setValue(i,j,0);}for(i=1;i<=Q.getRow();i++)Q.setValue(i,i,sum/A.edVec[i-1].weight) ;};void l_zhen( CMatrix& l ,CMatrix& L, CElevationNet A){int i;double m;for(i=0;i<A.numElvDif;i++)L.setValue(i+1,1,A.edVec[i].value);for(i=0;i<l.getRow();i++){m=(A.lpVec[A.edVec[i].startPoint-1].ele Value+A.edVec[i].value-A.lpVec[A.edVec[i].e ndPoint-1].eleValue)*1000;l.setValue(i+1,1,m);}};void main(){CElevationNet A;A.input();A.jsgc();CMatrixB(A.getgz(),A.getws()),Q(A.getgz(),A.getgz( )),X(A.getws(),1),l(A.getgz(),1),b(A.getws( ),A.getgz()),V(A.getgz(),1);CMatrixNBB(A.getws(),A.getws()),x(A.getws(),1) ,N(A.getws(),A.getws()),L(A.getgz(),1),M(A. getws(),A.getgz()),R(A.getws(),A.getgz()),W (A.getws(),A.getgz());CMatrixvz(1,A.getgz()),vp(1,A.getgz()),p(1,1);xishu(B,X,A);quanzhen(Q,A);l_zhen(l,L,A);b.operator =(B.transpose());M.operator =(b.operator *(Q));NBB.operator =(M.operator *(B));N.operator =(NBB.inverse());R.operator =(N.operator *(b));W.operator =(R.operator *(Q));x.operator =(W.operator *(l));int i,j;for(i=A.getys()+1;i<=A.getzs();i++){A.setdv(i,(x.getValue(i-A.getys(),1)/1000));A.seteleValue(i,A.getdv(i));cout<<"平差后编号为"<<A.getindex(i)<<"的点高程值："<<A.geteleValue(i)<<"米"<<endl;}V.operator =(B.operator *(x));V.operator -(l);for(j=0;j<V.getRow();j++)cout<<"第"<<j+1<<"段观测高差平差值为："<<(L.getValue(j+1,1)+V.getValue(j+1,1)/100 0)<<"米"<<endl;vz.operator =(V.transpose());vp=(vz.operator *(Q));p=(vp.operator *(V));cout<<"平差后单位权中误差为："<<"+_"<<sqrt(p.getValue(1,1)/A.getws())<<"毫米"<<endl;}。

四等水准测量开源代码namespace测量程序设计{public partial class Form3 : Form{public Form3(){InitializeComponent();}private void button3_Click(object sender, EventArgs e){}private void radioButton1_CheckedChanged(object sender, EventArgs e){//选择路线类型附合有终点高程闭合和支水准都不存在终点高程//选中附合txtz可用而闭合和和支水准不可用choiceStyle();}//定义全局变量方便其他控件调用int istation ;double[] dis;double []detH;double startpoint;double endpoint;//声明泛型集合来存储文件全路径List<double > listdis = new List<double >();List<double> listdetH= new List<double>();private void button1_Click(object sender, EventArgs e)AddDate();}double[] BMPoint;//转点高程数组double closeDetH;//闭合差double totalDetH = 0;//累计高差private void button2_Click(object sender, EventArgs e) {double[] diss = listdis.ToArray();double[] detHh = listdetH.ToArray();CalcDate(diss, detHh);}private void txtc_TextChanged(object sender, EventArgs e) {}//选择类型public void choiceStyle(){if (rb1.Checked == true){txtz.Enabled = true;}else if (rb2.Checked == true){txtz.Enabled = false;}else if (rb3.Checked == true){txtz.Enabled = false;}//读入数据public void AddDate(){if (txtq.Text == ""){MessageBox.Show("还没有输入起始高程");txtq.Focus();}else if (rb1.Checked == true && txtz.Text == ""){MessageBox.Show("还没有输入终点高程");txtz.Focus();}else if (txtc.Text == ""){MessageBox.Show("还没有输入测站数");txtc.Focus();}else{//点击弹出对话框OpenFileDialog ofd = new OpenFileDialog();//设置对话框可以多选ofd.Multiselect = true;//设置对话框的初始目录ofd.InitialDirectory = @"C:\Documents and Settings\Administrator\桌面";//设置对话框的文件类型ofd.Filter = "文本文件|*.txt|所有文件|*.*";//展示对话框ofd.ShowDialog();//获得在打开对话框中选中文件的路径string path = ofd.FileName;if (path == ""){return;}istation = int.Parse(txtc.Text);startpoint = double.Parse(txtq.Text);endpoint = double.Parse(txtz.Text);string[] line = File.ReadAllLines(path);for (int i = 0; i < istation; i++){dis = new double[istation];detH = new double[istation];char[] chs = { ',' };string[] str = line[i].Split(chs, StringSplitOptions.RemoveEmptyEntries);dis[i] = double.Parse(str[0]);//高差中数detH[i] = double.Parse(str[1]);//距离int station = i + 1;//数组转化为集合方便后面数据的使用listdis.Add(dis[i]);listdetH.Add(detH[i]);txtre.Text = txtre.Text + "\r\n" + "第" + station + "站：" + "\r\n" + "距离：" + dis[i].ToString() + " " + "高差中数：" + detH[i].ToString();}}}//平差计算public void CalcDate(double []diss,double []detHh){txtre.Text = txtre.Text + "\r\n" + "计算结果为：" + "\r\n";double tDist = 0;for (int i = 0; i < istation; i++){tDist = tDist + diss[i];//计算距离之和}for (int i = 0; i < istation; i++){totalDetH = totalDetH + detHh[i];//计算累计高差}if (rb1.Checked == true)//附合水准{closeDetH = totalDetH - (endpoint - startpoint);}else{closeDetH = -(totalDetH);//闭合水准和支水准闭合差取负}if (Math.Abs(closeDetH) > 40 * Math.Sqrt(tDist / 1000)){MessageBox.Show("闭合差超限，测量成果不合格！");txtre.Text = txtre.Text + "闭合差超限，成果不合格";return;}else{MessageBox.Show("闭合差合格，继续计算转点高程！");}BMPoint = new double[istation];BMPoint[0] = startpoint;txtre.Text = txtre.Text + "平差后的高程为：" + "\r\n" + "（" + "0" + "）" + startpoint.ToString() + "\r\n"; for (int i = 1; i < istation; i++){BMPoint[i] = BMPoint[i - 1] + detHh[i - 1] + closeDetH * diss[i - 1] / tDist;txtre.Text = txtre.Text + "\r\n" + "（" + i.ToString() + "）" + BMPoint[i].ToString("#0.0000") + "\r\n";}txtre.Text = txtre.Text + "\r\n" + "（" + istation.ToString() + "）" + endpoint.ToString(); ;}}。

实验一一．设计原始资料水准网周密平差及精度评定示例。

如图所示水准网，有2个已知点，3个未知点，7个测段。

各已知数据及观测值见下表（1）已知点高程H1= H2=（2）高差观测值(m)高差观测值(m)端点高差观测值测段距离序号号1-3 11-4 22-3 32-4 43-4 53-5 65-2 7（3）求各待定点的高程；3-4点的高差中误差；3号点、4号点的高程中误差。

（提示，本网可采用以测段的高差为平差元素，采用间接平差法编写程序计算。

）二、水准网间接平差思路⑴．按照网型肯定已知水准点数H1 H2，未知水准点数u ，总点数n ，必要观测数t ，多余观测数r 。

⑵.肯定参数。

为平差后能直接求得待定点高程平差值，以3个待定点高程平差值为参数。

设3,4,5点的高程平差值别离为，， 。

⑶．列立条件方程.左侧为观测值（系数为1），右边为参数和常数项，并进一步改化成误差方程，最终写成矩阵形式。

取得系数矩阵A 和常数项L⑷.列立法方程，并解求法方程。

⑸．精度评定。

计算单位权中误差的估值：评定各待定点的高程中误差： 各待定点的精度即为： 评定高程平差值的精度： 四、平差程序设计思路1、 已知数据的输入需要输入的数据包括水准网中已知点数、未知点数和这些点的点号、已知高程和高差观测值、距离观测值等。

本程序采用文件方式进行输入，文件输入的格式如下： 第一行：已知点个数、未知点个数、观测值个数 第二行：点号 （已知点在前，为支点在后） 第三行：已知高程 （顺序与上一行的点号对应）第四行起：高差观测值，依照“起点点号，终点点号，高差观测值，距离观测值”的顺序输入。

2、 平差计算进程 （1）近似高程计算。

uc PV V r PV V T T -==20ˆσ120ˆˆ20ˆˆˆˆ-==bbx x x x N Q D σσjj X X j X Q ˆˆ0ˆˆσσ±=FN F F Q F Q BB T X X T h h 1ˆˆˆˆ-==X F hT ˆˆ=（2）列立观测值的误差方程。