1、需求分析
总体描述:
功能需求:
2.原理与系统设计
原理:
概算原理:
平差原理:
系统设计:
类整体结构设计:
概算:
平差:
具体步骤:
1.列出误差方程和条件方程
2、定权
三、构成方程
4、精度评定
3.测试结果及分析
显示运行结果:
TXT:编辑
鲁棒检验:
4 、主要代码
Dataimport:
.h:
.cpp:
RS326Dlg类:
.h:
.cpp:
CDlgpaint类:
.h:
.cpp:
PaintUI类:
.h:
.cpp:
工具类:
Datapoint.h:
Dataline.h:
Dataangle.h:
矩阵类:-toc" style="margin-left:40px;">Matrix矩阵类:
.h:
.cpp:
能力有限,代码冗杂!
1、需求分析
总体描述:
控制网平差程序控制野外控制网观测数据平差数据处理的目的是基于最小二乘原理,消除网络中的各种几何矛盾,找出全网各待定元素(未知点平面坐标或三维坐标)。
功能需求:
2.原理与系统设计
概算原理:
平差原理:
平差模型:最小二乘平差模型。
间接平差:一个观测值就是一个误差方程,占用内存相对较大,但易于编程解算,容易编制通用的程序。
系统设计:
类整体结构设计:
概算:
逐点解算法:选择已知点或已计算出坐标的点作为起算点,根据具体的图形,逐个推算出各点的坐标。
平差:
具体步骤:
1、误差方程式及条件方程式的列立
2、定权
在等精度方向观测的控制网中,可把方向观测值的权定为1,不同精度的观测网中,则选择其中一种作为单位权中误差u(单位:秒),其它方向观测值 中误差为m(单位:秒),则。
3、组成法方程
4、精度评定
3、测试结果与分析
运行结果展示:
TXT:
鲁棒性检验:
4 、主要代码
Dataimport:
.h:
#pragma once
#include "Datapoint.h"
#include "Dataangle.h"
#include"Dataline.h"
#include "Matrix.h"
#include "CAngle.h"
class Dataimport
{
public:
Dataimport();
~Dataimport();
public:
Point* Psum;//点集
int Pcount;//点数
Line* Lsum;//边集
int Lcount;//边数
Angle* Asum;//测站集
int Acount;//测站数
Ang* Ansum;//观测角集
int ancount;//观测角数
int unk;//未知点数
CMatrix dLe;//B矩阵
CMatrix f;//F矩阵
CMatrix p1;//权矩阵
int t = 0;//递归次数
double m0;//单位权中误差
public:
int SplitStringArray(CString str, char c, CStringArray& array);//字符串分割
CString remove(CString str);//移除字符串中\r\n
double dms2d(double k);//度分秒转度
void read(CString&input,CString &input1);//数据读取和存入
void dealangle();//出来数据获取观测角
double p2azi(Point p1,Point p2);//计算方位角
int azimuth(CString p1);//查找点号
int azimuth1(CString p1, CString p2);//查找边号
void dofdeal();//数据概算
void output1(CString& output);//输出数据概算
void Do_adjustment();//平差计算
void output(CString &output);//输出平差结果
void Save(CString output);//输出到TXT
void Save1(CString output);//输出到EXcel
};
.cpp:
#include "pch.h"
#include "Dataimport.h"
#include <clocale>
#include "math.h"
//构造和析构函数
Dataimport::Dataimport()
{
}
Dataimport::~Dataimport()
{
if (Psum != NULL) {
delete[] Psum;
Psum = NULL;
}
if (Asum != NULL) {
delete[] Asum;
Asum = NULL;
}
if (Lsum != NULL) {
delete[] Lsum;
Lsum = NULL;
}
if (Ansum != NULL) {
delete[] Ansum;
Ansum = NULL;
}
}
//字符串分割
int Dataimport::SplitStringArray(CString str, char c, CStringArray& array)
{
int start = 0;
int idx = str.Find(c, start);
array.RemoveAll();
while (idx != -1) {
CString stmp = str.Mid(start, idx - start);
array.Add(stmp);
start = idx + 1;
idx = str.Find(c, start);
}
CString stmp = str.Right(str.GetLength() - start);
if (!stmp.IsEmpty())
{
array.Add(stmp);
return array.GetSize();
}
}
//字符串处理去除\r\n
CString Dataimport::remove(CString str)
{
str.Remove('\r');
str.Remove('\n');
return str;
}
//度分秒转度
double Dataimport::dms2d(double k)
{
int a4;
double a2;
double a3;
a4 = int(k);
a2 = int((k - a4) * 100);
a3 = ((k - a4) * 100 - a2) * 100;
return double((a4 + a2 / 60 + a3 / 3600));
}
//数据读取分割存储
void Dataimport::read(CString& input, CString& input1)
{
CFileDialog dlgFile(TRUE, _T("txt"), NULL, OFN_EXPLORER, _T("(文本文件)|*.dat"));
if (dlgFile.DoModal() == IDCANCEL) return;
CString trFileName = dlgFile.GetPathName();
CStdioFile rfile;
if (!rfile.Open(trFileName, CFile::modeRead))AfxMessageBox(_T("文件未找到"));
CString buf = _T("");
CString bufer = _T("");
CString str1 = _T("**********");
int ind = 0;
while (rfile.ReadString(buf))
{
if (ind == 1)
{
bufer += buf + _T("\r\n");
}
if (remove(buf) == str1)
{
ind = 1;
}
}
if (ind==0)
{
AfxMessageBox(_T("输入数据格式不符合"));
return;
}
input1 = trFileName;
input = bufer;
CStringArray array;
SplitStringArray(bufer, '\r', array);
int Pc1 = _tstof(array[0]);
int Pc2 = _tstof(array[Pc1 + 1]);
Pcount = Pc1 + Pc2;
unk = Pc2;
Psum = new Point[Pcount];
for (int i = 0; i < Pc1; i++)
{
CStringArray buf1;
SplitStringArray(array[i + 1], ',', buf1);
CString a=buf1[0];
Psum[i].index = buf1[0];
Psum[i].x = _tstof(buf1[1]);
Psum[i].y = _tstof(buf1[2]);
Psum[i].flag = 0;
Psum[i].flag1 = 0;
}
CStringArray buf1;
SplitStringArray(array[Pc1 + 2], ',', buf1);
int m = 0;
for (int i = Pc1; i < Pcount; i++)
{
Psum[i].index = buf1[m];
Psum[i].x = 0;
Psum[i].y = 0;
Psum[i].flag = 1;
m++;
}
Lcount = _tstof(array[Pc1 + 3]);
Lsum = new Line[Lcount];
for (int i=0; i < Lcount; i++) {
CStringArray buf2;
SplitStringArray(array[Pc1 + 4+i], ',', buf2);
Lsum[i].start = buf2[0];
Lsum[i].end = buf2[1];
Lsum[i].length = _tstof(buf2[2]);
Lsum[i].amangle = 0;
}
Acount = _tstof(array[Pc1 + 4 + Lcount]);
Asum = new Angle[Acount];
int a = -1;
for (int i = 0; i < Acount; i++)
{
CStringArray buf2;
SplitStringArray(array[Pc1 + 5 +Lcount+ i], ',', buf2);
Asum[i].startP = buf2[0];
Asum[i].endP = buf2[1];
Asum[i].angle = _tstof(buf2[2]);
if (Asum[i].angle==0)
{
a++;
}
Asum[i].num = a;
}
}
//数据进一步处理
void Dataimport::dealangle()
{
//计算观测角个数实例化
ancount = 0;
for (int i = 0; i < Acount; i++)
{
if (Asum[i].angle != 0)ancount++;
}
Ansum = new Ang[ancount];
//赋值
int m = 0;
for (int i = 0; i < Acount; i++)
{
if (Asum[i].angle == 0)
{
int s = i;
while (Asum[i+1].angle != 0&&i+1<Acount)
{
Ansum[m].startP = Asum[s].endP;
Ansum[m].mid = Asum[s].startP;
Ansum[m].endP = Asum[i + 1].endP;
Ansum[m].angle = Asum[i + 1].angle;
m++;
i++;
}
}
}
}
void Dataimport::dealangle()
{
//计算观测角个数实例化
ancount = 0;
for (int i = 0; i < Acount; i++)
{
if (Asum[i].angle != 0)ancount++;
}
Ansum = new Ang[ancount];
//赋值
int m = 0;
for (int i = 0; i < Acount; i++)
{
if (Asum[i].angle == 0)
{
int s = i;
while (Asum[i+1].angle != 0&&i+1<Acount)
{
Ansum[m].startP = Asum[s].endP;
Ansum[m].mid = Asum[s].startP;
Ansum[m].endP = Asum[i + 1].endP;
Ansum[m].angle = Asum[i + 1].angle;
m++;
i++;
}
}
}
}
//两点求方位角
double Dataimport::p2azi(Point p1,Point p2)
{
int m = 0;
int n = 0;
//记录两点的下标
int a = 0;
int b = 0;
double d[10];
//查找边起始点是否为已知点
for (int j = 0; j < Pcount; j++)
{
if (remove(p1.index) == remove(Psum[j].index))
{
if (Psum[j].flag == 0)
{
m = 1;
a = j;
break;
}
}
}
//查找边终点是否为已知点
for (int j = 0; j < Pcount; j++)
{
if (remove(p2.index) == remove(Psum[j].index))
{
if (Psum[j].flag == 0)
{
n = 1;
b = j;
break;
}
}
}
//两点已知进行方位角计算
double ang;
if (m == n == 1) {
double dx = Psum[b].x - Psum[a].x;
double dy = Psum[b].y - Psum[a].y;
ang = atan(abs(dy / dx));
if (dx ==0&& dy == 0)return 0;
if (dy > 0 && dx < 0) { ang = acos(-1) - ang; }
if (dy < 0 && dx < 0) { ang = acos(-1) + ang; }
if (dy < 0 && dx > 0) { ang = 2 * acos(-1) - ang; }
}
else
{
ang = 0;
}
return ang;
}
//点号查找点下标
int Dataimport::azimuth(CString p1)
{
int a;
for (int j = 0; j < Pcount; j++)
{
if (remove(p1) == remove(Psum[j].index))
{
a = j;
return a;
}
}
return -1;
}
//两点号查找边下标
int Dataimport::azimuth1(CString p1, CString p2)
{
int d;
for (int j = 0; j < Lcount; j++)
{
if ((remove(p1) == remove(Lsum[j].start)) && (remove(p2) == remove(Lsum[j].end))) {
d = j;
return d;
}
if ((remove(p1) == remove(Lsum[j].end)) && (remove(p2) == remove(Lsum[j].start))) {
d = j;
return d;
}
}
return -1;
}
//数据概算
void Dataimport::dofdeal()
{
int sum = 0;
double d1[28];
double d2[28];
while (sum < Pcount)
{
sum = 0;
for (int i = 0; i < ancount; i++)
{
int m = 0; int n = 0;
int a, b,c,d=-1,e=-1;
a = azimuth(Ansum[i].startP);
b = azimuth(Ansum[i].mid);
c = azimuth(Ansum[i].endP);
d = azimuth1(Ansum[i].startP, Ansum[i].mid);
e = azimuth1(Ansum[i].mid, Ansum[i].endP);
int c1 = Psum[a].flag;
int c2 = Psum[b].flag;
int c3 = Psum[c].flag;
if (c2 == 0) {
if (c1==0&&c3==1)
{
double ang1 = p2azi(Psum[b], Psum[a]);
double length = 0;
length = Lsum[e].length;
if (e == -1) { length = 0; }
if (ang1 != 0 && length != 0)
{
double er = (dms2d(Ansum[i].angle) / 180) * acos(-1);
double ang = ang1 + (dms2d(Ansum[i].angle) / 180) * acos(-1) - 2 * acos(-1);
Psum[c].x = Psum[b].x + length * cos(ang);
Psum[c].y = Psum[b].y + length * sin(ang);
Psum[c].flag = 0;
Lsum[e].amangle = ang;
double ad = cos(ang);
d2[c] = Psum[c].x;
}
}
if (c1 == 1&& c3 ==0)
{
double ang1 = p2azi(Psum[b], Psum[c]);
double length = 0;
length = Lsum[d].length;
if (d == -1) { length = 0; }
if (ang1 != 0 && length != 0)
{
double er = (dms2d(Ansum[i].angle) / 180) * acos(-1);
double ang = ang1 - (dms2d(Ansum[i].angle) / 180) * acos(-1) + 2 * acos(-1);
Psum[a].x = Psum[b].x + length * cos(ang);
Psum[a].y = Psum[b].y + length * sin(ang);
Psum[a].flag = 0;
Lsum[d].amangle = ang;
double ad = cos(ang);
d2[a] = Psum[a].x;
}
}
if (c1 == 1 && c3 == 1)break;
}
}
//统计已知坐标点数
for (int k = 0; k < Pcount; k++)
{
if (Psum[k].flag == 0)sum++;
d1[k] = Psum[k].x;
}
}
}
//数据概算输出
void Dataimport::output1(CString& output)
{
output = _T("");
CString strOut;
strOut.Format(_T("\r\n-------------------------------数据概算结果---------------------------------\r"));
output += strOut;
strOut.Format(_T("\r\n概算后待定点坐标:\r\n"));
output += strOut;
strOut.Format(_T("点号\tX坐标\t\tY坐标\r\n"));
output += strOut;
for (int i = 2; i < Pcount; i++)
{
strOut.Format(_T("%s\t%.4f\t%.4f\r\n"), Psum[i].index, Psum[i].x, Psum[i].y);
output += strOut;
}
}
//平差计算
void Dataimport::Do_adjustment()
{
//平差
t = 0;
double sum = 0;//改正数平均差值
int Ncount = -ancount + Acount;
CMatrix x;//改正数矩阵
x.SetSize(unk * 2 + Ncount, 1);
CMatrix Nbb;
CMatrix f1;
/
do
{
//方向观测误差方程
dLe.SetSize(Acount+Lcount, unk*2+Ncount);
f.SetSize(Acount + Lcount, 1);
int p = 206265;
for (int i = 0; i < Acount; i++)
{
int m = 0;
int n = 0;
int a, b,c,d;
a = azimuth(Asum[i].startP);
b = azimuth(Asum[i].endP);
c = azimuth1(Psum[a].index, Psum[b].index);
double length;
if (c==-1) {
length = sqrt((Psum[a].x- Psum[b].x) * (Psum[a].x - Psum[b].x) + (Psum[a].y- Psum[b].y) * (Psum[a].y - Psum[b].y));
}
else
{
length = Lsum[c].length;
}
//查看是否为控制点
if (Psum[a].flag1 == 0){m = 1;}
if (Psum[b].flag1 == 0){n = 1;}
if ((m == 1 && n == 1))
{
dLe(i, unk * 2 + Asum[i].num) = 0;
f(i,0) = 0;
}
if (m == 1 && n == 0)
{
double a1, b1;
double ang = p2azi(Psum[a], Psum[b]);
double ang1=0;
for (int k = 0; k < Acount; k++)
{
if (Asum[i].num == Asum[k].num) {
int po1 = azimuth(Asum[k].startP);
int po2 = azimuth(Asum[k].endP);
ang1= p2azi(Psum[po1],Psum[po2]);
break;
}
}
a1 = p * sin(ang) /length;
b1 =- p * cos(ang) /length;
dLe(i, (b - 2) * 2) = -a1;
dLe(i, (b - 2) * 2+1) = -b1;
dLe(i, unk * 2 + Asum[i].num) = -1;
double ang2 = 0;
ang2 = ang - ang1;
if ((ang - ang1)<0)
{
ang2 = ang2 + 2 * acos(-1);
}
double d_angle = ang2 - (dms2d(Asum[i].angle) / 180) * acos(-1);
double dd = (dms2d(Asum[i].angle) / 180) * acos(-1);
f(i, 0) =d_angle*p;
}
if (m == 0 && n == 1)
{
double a1, b1;
double ang = p2azi(Psum[a], Psum[b]);
double ang1 = 0;
for (int k = 0; k < Acount; k++)
{
if (Asum[i].num == Asum[k].num) {
int po1 = azimuth(Asum[k].startP);
int po2 = azimuth(Asum[k].endP);
ang1 = p2azi(Psum[po1], Psum[po2]);
break;
}
}
a1 = p * sin(ang) / length;
b1 = -p * cos(ang) / length;
dLe(i, (a - 2) * 2) = a1;
dLe(i, (a - 2) * 2 + 1) = b1;
dLe(i, unk * 2 + Asum[i].num) = -1;
double ang2 = 0;
ang2 = ang - ang1;
if ((ang - ang1) < 0)
{
ang2 = ang2 + 2 * acos(-1);
}
double d_angle = ang2 - (dms2d(Asum[i].angle) / 180) * acos(-1);
f(i, 0) = d_angle * p;
}
if (m == 0 && n == 0)
{
double a1, b1;
double ang = p2azi(Psum[a], Psum[b]);
double ang1 = 0;
for (int k = 0; k < Acount; k++)
{
if (Asum[i].num == Asum[k].num) {
int po1 = azimuth(Asum[k].startP);
int po2 = azimuth(Asum[k].endP);
ang1 = p2azi(Psum[po1], Psum[po2]);
break;
}
}
a1 = p * sin(ang) / length;
b1 = -p * cos(ang) / length;
dLe(i, (b - 2) * 2) = -a1;
dLe(i, (b - 2) * 2 + 1) = -b1;
dLe(i, (a - 2) * 2) = a1;
dLe(i, (a - 2) * 2 + 1) = b1;
dLe(i, unk * 2 + Asum[i].num) = -1;
double ang2 = 0;
ang2 = ang - ang1;
if ((ang - ang1) < 0)
{
ang2 = ang2 + 2 * acos(-1);
}
double d_angle = ang2 - (dms2d(Asum[i].angle) / 180) * acos(-1);
f(i, 0) = d_angle * p;
double dc = f(i, 0);
}
}
//边长观测误差
for (int i = Acount; i < Lcount+Acount; i++)
{
int m =0, n = 0;
int a, b;
a = azimuth(Lsum[i-Acount].start);
b = azimuth(Lsum[i-Acount].end);
if (Psum[a].flag1 == 0)
{
m = 1;
}
if (Psum[b].flag1 == 0)
{
n = 1;
}
double a1, b1,s;
double ang = p2azi(Psum[a], Psum[b]);
a1 = sin(ang);
b1 = cos(ang);
s = sqrt((Psum[a].x- Psum[b].x) * (Psum[a].x - Psum[b].x) + (Psum[a].y - Psum[b].y) * (Psum[a].y - Psum[b].y));
//
if (m == 1 && n == 1) {
f( i, 0) = s - Lsum[i-Acount].length;
}
if (m == 0 && n == 0) {
dLe(i, (a - 2) * 2) = -b1;
dLe(i, (a - 2) * 2 + 1) = -a1;
dLe(i, (b - 2) * 2) = b1;
dLe(i, (b - 2) * 2 + 1) = a1;
f(i, 0) = s - Lsum[i - Acount].length;
}
if (m == 1 && n == 0) {
dLe(i, (b - 2) * 2) = b1;
dLe(i, (b - 2) * 2 + 1) = a1;
f( i, 0) = s - Lsum[i - Acount].length;
}
if (m == 0 && n == 1) {
dLe(i, (a - 2) * 2) = -b1;
dLe(i, (a - 2) * 2 + 1) = -a1;
f( i, 0) = s - Lsum[i - Acount].length;
}
}
CString pp1 = dLe.look();
CString pp2 = f.look();
//权矩阵
p1.SetSize(Acount + Lcount, Acount + Lcount);
for (int i = 0; i < Acount; i++)
{
p1(i, i) = 1;
}
for (int i = Acount; i < Acount+Lcount; i++)
{
p1(i, i) = 100000000/Lsum[i-Acount].length;
}
CString pp3 = p1.look();
//
f1 = -1*f;
Nbb= (~dLe * p1 * dLe).Inv();
x=Nbb*(~dLe*p1*f1);
sum = 0;
for (int i = 0; i <(Pcount-2)*2-1; i++)
{
if (abs(x(i,0))<abs(x(i+1,0)))
{
sum =abs( x(i+1, 0));
}
}
for (int k = 2; k < Pcount; k++)
{
Psum[k].x = Psum[k].x+x(2*(k-2), 0);
Psum[k].y = Psum[k].y+x(2*(k-2)+1, 0);
}
double df2[80];
for (int i = 0; i < x.Row(); i++)
{
df2[i] = x(i, 0);
}
for (int i = 0; i < x.Row(); i++)
{
x(i, 0) = 0;
}
t++;
} while (sum>0.0000001);
double df3[28];
for (int i = 0; i <Pcount; i++)
{
df3[i] = Psum[i].x;
}
//精度评定
CMatrix vpv = ~f1 * p1 * f1 - ~f1 * p1 * dLe * x;
double vpv1 = vpv(0, 0);
m0 =sqrt( vpv1 / (Acount + Lcount - 2 * unk));
CMatrix Qx = Nbb;
for (int i = 0; i < Pcount-2; i++)
{
//待定点点位误差
Psum[i+2].mx = m0 * sqrt(Qx(2 * i, 2 * i));
Psum[i+2].my = m0 * sqrt(Qx(2 * i + 1, 2 * i + 1));
Psum[i+2].mk = m0 * sqrt(Psum[i+2].mx * Psum[i+2].mx + Psum[i+2].my * Psum[i+2].my);
//误差椭圆参数计算
Psum[i+2].Q = 0.5 * atan(2 * Qx(2 * i, 2 * i + 1) / (Qx(2 * i, 2 * i) - Qx(2 * i + 1, 2 * i + 1)));
Psum[i+2].E = m0 * sqrt(Qx(2 * i, 2 * i) + Qx(2 * i, 2 * i + 1) * tan(Psum[i+2].Q));
Psum[i+2].F = m0 * sqrt(Qx(2 * i, 2 * i) + Qx(2 * i, 2 * i + 1) * tan(Psum[i+2].Q + acos(-1) / 2));
}
}
//平差结果输出
void Dataimport::output(CString &output)
{
CString strOut;
strOut.Format(_T("\r\n-------------------------------------------------平差结果-------------------------------------------------\r"));
output += strOut;
strOut.Format(_T("\r\n平差迭代次数:%d\r"), t);
output += strOut;
strOut.Format(_T("\r\n平差单位权中误差:%.8f\r"), m0);
output += strOut;
strOut.Format(_T("\r\n平差后系数阵B的元素:\r\n"));
output += strOut;
for (int i = 0; i <dLe.Row(); i++)
{
for (int j = 0; j < dLe.Col(); j++)
{
strOut.Format(_T("%-10.4f\t"), dLe(i, j));
output += strOut;
}
strOut.Format(_T("\r\n"));
output += strOut;
}
strOut.Format(_T("\r\n\n平差后常数阵F的元素:\r\n"));
output += strOut;
for (int i = 0; i < f.Row(); i++)
{
for (int j = 0; j < f.Col(); j++)
{
if (i < Acount)
{
strOut.Format(_T("%.6f\t"), f(i, j));
output += strOut;
}
else
{
strOut.Format(_T("%.6f\t"), f(i, j));
output += strOut;
}
}
strOut.Format(_T("\r\n"));
output += strOut;
}
strOut.Format(_T("\r\n\n平差后常数阵P的元素:\r\n"));
output += strOut;
for (int i = 0; i < p1.Row(); i++)
{
strOut.Format(_T("%.6f\t\t"), p1(i, i));
output += strOut;
strOut.Format(_T(""));
output += strOut;
}
strOut.Format(_T("\r\n平差后待定点坐标:\r\n"));
output += strOut;
strOut.Format(_T("点号\tX坐标\t\tY坐标\r\n"));
output += strOut;
for (int i = 2; i < Pcount; i++)
{
strOut.Format(_T("%s\t%.4f\t%.4f\r\n"), Psum[i].index, Psum[i].x, Psum[i].y);
output += strOut;
}
strOut.Format(_T("\r\n待定点点位误差及误差椭圆参数:\r\n"));
output += strOut;
strOut.Format(_T("点号\tmx\t\tmy\t\tmk\t\t长半径方位角Q\t\t长半径E\t\t短半径F\r\n"));
output += strOut;
for (int i = 2; i < Pcount; i++)
{
strOut.Format(_T("%s\t%.10f\t%.10f\t%.10f\t%.10f\t\t%.10f\t%.10f\r\n"), Psum[i].index, Psum[i].mx, Psum[i].my, Psum[i].mk, Psum[i].Q, Psum[i].E, Psum[i].F);
output += strOut;
}
}
//文件保存txt
void Dataimport::Save(CString output)
{
if (output == _T("")) { AfxMessageBox(_T("请先输入数据!")); }
else {
CFileDialog dlg(false, NULL, NULL, OFN_HIDEREADONLY | OFN_OVERWRITEPROMPT, _T("Txt Files(*.txt)|*.txt|All Files(*.*)|*.*|"), AfxGetMainWnd());
CString strPath;
if (dlg.DoModal() == IDCANCEL) return;
else
{
strPath = dlg.GetPathName();
if (strPath.Find(_T(".txt")) < 0)
{
if (strPath.Find(_T(".dat")) < 0)
{
strPath += _T(".txt");//默认储存为txt模式
}
}
}
TRY
{
CString cs;
CStdioFile file(strPath, CFile::shareExclusive | CFile::modeWrite | CFile::modeCreate);
setlocale(LC_CTYPE, ("chs")); //设置中文输出
file.WriteString(output);
file.Close();
CString temp;
temp.Format(_T("%s%s"), _T("已导出至"),strPath);
AfxMessageBox(temp);
}
CATCH_ALL(e)
{
e->ReportError();
return;
}
END_CATCH_ALL
}
}
//文件保存Excel
void Dataimport::Save1(CString output)
{
if (output == _T("")) { AfxMessageBox(_T("请先输入数据!")); }
else {
CString strPath;
CFileDialog dlg(false, NULL, NULL, OFN_HIDEREADONLY | OFN_OVERWRITEPROMPT, _T("Excel Files(*.xls)|*.xls|All Files(*.*)|*.*|"), AfxGetMainWnd());
if (dlg.DoModal() == IDCANCEL) return;
else
{
strPath = dlg.GetPathName();
if (strPath.Find(_T(".xls")) < 0)
{
strPath += _T(".xls");
}
}
TRY
{
CString cs;
CStdioFile file(strPath, CFile::shareExclusive | CFile::modeWrite | CFile::modeCreate);
setlocale(LC_CTYPE, ("chs")); //设置中文输出
file.WriteString(output);
file.Close();
CString temp;
temp.Format(_T("%s%s"), _T("已导出至"),strPath);
AfxMessageBox(temp);
}
CATCH_ALL(e)
{
e->ReportError();
return;
}
END_CATCH_ALL
}
}
RS326Dlg类:
.h:
........
public:
CString input;
afx_msg void OnBnClickedButton1();
afx_msg void OnBnClickedButton2();
afx_msg void OnBnClickedButton3();
afx_msg void OnBnClickedButton4();
CString output;
CDlgpaint cdlg;
CString input3;
afx_msg void OnBnClickedButton5();
afx_msg void OnBnClickedButton6();
afx_msg void OnBnClickedButton7();
afx_msg void OnBnClickedButton8();
......
.cpp:
Dataimport im;
int a = 0;
void CRS326Dlg::OnBnClickedButton1()
{
UpdateData(true);
im.dLe.SetSize(1, 1);
im.dLe(0, 0) = 0;
im.f.SetSize(1, 1);
im.f(0, 0) = 0;
im.p1.SetSize(1, 1);
im.p1(0, 0) = 0;
im.read(input,input3);
im.dealangle();
a = 1;
output = _T("");
UpdateData(false);
}
void CRS326Dlg::OnBnClickedButton2()
{
UpdateData(true);
if (a>0)
{
im.dofdeal();
im.output1(input);
a = 2;
output = _T("");
im.t = 0;
}
else
{
AfxMessageBox(_T("请读入数据"));
}
UpdateData(false);
}
void CRS326Dlg::OnBnClickedButton3()
{
UpdateData(true);
if (a==2)
{
im.Do_adjustment();
im.output(output);
a = 3;
im.t = 0;
}
else if(a==0)
{
AfxMessageBox(_T("请读入数据"));
}
else
{
AfxMessageBox(_T("请先进行数据概算"));
}
UpdateData(false);
}
void CRS326Dlg::OnBnClickedButton4()
{
UpdateData(true);
if (a==3)
{
cdlg.set(im.Psum, im.Asum, im.Pcount, im.Acount);
cdlg.DoModal();
}
else if(a==0)
{
AfxMessageBox(_T("请读入数据"));
}
else if(a==1)
{
AfxMessageBox(_T("请先进行进行数据概算"));
}
else
{
AfxMessageBox(_T("请先进行进行平差"));
}
UpdateData(false);
}
void CRS326Dlg::OnBnClickedButton5()
{
UpdateData(true);
if (a>2)
{
im.Save(output);
}
else
{
AfxMessageBox(_T("请先进行平差"));
}
UpdateData(false);
}
void CRS326Dlg::OnBnClickedButton6()
{
UpdateData(true);
input= _T("");
output = _T("");
input3= _T("");
a = 0;
UpdateData(false);
}
void CRS326Dlg::OnBnClickedButton7()
{
UpdateData(true);
if (a > 2)
{
im.Save1(output);
}
else
{
AfxMessageBox(_T("请先进行平差"));
}
UpdateData(false);
}
void CRS326Dlg::OnBnClickedButton8()
{
UpdateData(true);
ShellExecute(NULL, _T("open"), _T("C:\\Users\\lenovo\\Desktop\\帮助文档.docx"), NULL, NULL, SW_SHOW);//弹出使用说明
UpdateData(false);
}
CDlgpaint 类:
.h:
public:
Point* PO;//点数组
Angle* Ap;//测站数组
int count;
int count1;
protected:
virtual void DoDataExchange(CDataExchange* pDX); // DDX/DDV 支持
DECLARE_MESSAGE_MAP()
public:
void set(Point* Psum, Angle* Asum,int Pcount,int Acount);//值传递
CStatic Picture;
afx_msg void OnPaint();//窗口刷新函数
.....
.cpp:
void CDlgpaint::set(Point* Psum, Angle* Asum, int Pcount, int Acount)
{
count = Pcount;
count1 = Acount;
PO = new Point[count];
Ap = new Angle[count1];
for (int i = 0; i < Pcount; i++)
{
PO[i].index = Psum[i].index;
PO[i].x = Psum[i].x;
PO[i].y = Psum[i].y;
PO[i].Q = Psum[i].Q;
PO[i].E = Psum[i].E;
PO[i].F = Psum[i].F;
PO[i].flag1 = Psum[i].flag1;
}
for (int i = 0; i < Acount; i++)
{
Ap[i].startP = Asum[i].startP;
Ap[i].endP = Asum[i].endP;
}
}
BEGIN_MESSAGE_MAP(CDlgpaint, CDialogEx)
ON_WM_PAINT()
END_MESSAGE_MAP()
// CDlgpaint 消息处理程序
PaintUI UI;
int t = 0;
void CDlgpaint::OnPaint()
{
CPaintDC dc(this); // device context for painting
// TODO: 在此处添加消息处理程序代码
// 不为绘图消息调用 CDialogEx::OnPaint()
CWnd* Pwin = GetDlgItem(IDC_STATIC);
CRect rect;
Pwin->GetClientRect(rect);
CDC* pDC = Pwin->GetDC();
UI.dofdeal(PO,Ap,count,count1);
UI.XY2po(rect);
UI.DrawUI(pDC, rect);
t++;
}
PaintUI 类:
.h:
#pragma once
#include "CDlgpaint.h"
#include "Datapoint.h"
#include "Dataangle.h"
struct UIpoint //点结构体
{
CString ID;
double x;
double y;
double Q;
double E;
double F;
int flag;
};
struct UIline//边结构体
{
CString start;
CString end;
};
class PaintUI
{
public:
PaintUI();
~PaintUI();
private:
UIpoint* P;//点数组
int count;//点数
UIline* L;//边数组
int count1;//边数
double width;//画布长度
double heigth;//画布宽度
double scale;//比例尺
//实地原点坐标
int Ox;
int Oy;
//栅格网单位长度和实地对于长度
double len;
double len1;
//栅格网长宽对应边数
int wdcount;
int heicount;
//误差椭圆比例尺
int dscale;
int dxw = 0;
int dyw = 0;
public:
CString remove(CString str);//字符串处理
int foundP(CString p1);//点号查找点对应下标
double change(double z);//坐标改正
void dofdeal(Point*Psum,Angle*Asum,int Pcount,int Acount);//值传递
void XY2po(CRect rect);//预处理,包括比例尺计算,坐标赋值
void DrawUI(CDC* pDC, CRect rect);//绘制图像
};
.cpp:
#include "pch.h"
#include "PaintUI.h"
PaintUI::PaintUI()
{
}
PaintUI::~PaintUI()
{
if (P != NULL) {
delete[]P;
}
if (L!=NULL)
{
delete[]L;
}
}
CString PaintUI::remove(CString str)
{
str.Remove('\r');
str.Remove('\n');
return str;
}
int PaintUI::foundP(CString p1)
{
int a;
for (int j = 0; j < count; j++)
{
if (p1 == P[j].ID)
{
a = j;
return a;
}
}
return -1;
}
void PaintUI::dofdeal(Point* Psum, Angle* Asum, int Pcount,int Acount)
{
P = new UIpoint[Pcount];
count = Pcount;
for (int i = 0; i < Pcount; i++)
{
P[i].ID = remove(Psum[i].index);
P[i].x = Psum[i].y;
P[i].y = Psum[i].x;
P[i].Q = Psum[i].Q;
P[i].F = Psum[i].E;
P[i].E = Psum[i].F;
P[i].flag = Psum[i].flag1;
}
L = new UIline[Acount];
count1 = Acount;
for (int i = 0; i < Acount; i++)
{
L[i].start = remove(Asum[i].startP);
L[i].end = remove(Asum[i].endP);
}
}
double PaintUI::change(double z)
{
double a = 0;
a = -z + heigth;
return a;
}
//
void PaintUI::XY2po(CRect rect)
{
//计算比例尺
width = rect.Width();
heigth = rect.Height();
double minx=P[0].x, maxx= P[0].x, miny= P[0].y, maxy=P[0].y;
double dx, dy;
for (int i = 0; i <count; i++)
{
if (P[i].x<minx)
{
minx = P[i].x;
}
if (P[i].x > maxx)
{
maxx = P[i].x;
}
if (P[i].y < miny)
{
miny = P[i].y;
}
if (P[i].y > maxy)
{
maxy = P[i].y;
}
}
dx = maxx - minx;
dy = maxy - miny;
double rx = dx / (heigth * 0.9);
double ry = dy / (width * 0.9);
scale = rx < ry ? ry : rx;
//计算坐标原点实地坐标及栅格大小
Ox = ((int)(minx / 50)) * 50;
Oy = ((int)(miny / 50)) * 50;
double max = dx < dy ? dy : dx;
int nlen = max/ 4;
int dlen = nlen % 100;
if(dlen<50)
{
len = ((int)(nlen / 100)) * 100;
}
else
{
len = ((int)(nlen / 100)) * 100+100;
}
len1 = len;
len = len / scale;
wdcount = (width / len);
heicount = (heigth / len);
for (int i = 0; i < count; i++)
{
P[i].x = (P[i].x - Ox) / scale;
P[i].y = (P[i].y - Oy) / scale;
P[i].y = change(P[i].y);
P[i].x = P[i].x + len;
}
//居中计算
double minx1 = P[0].x, maxx1 = P[0].x, miny1 = P[0].y, maxy1 = P[0].y;
double dx1, dy1;
for (int i = 0; i < count; i++)
{
if (P[i].x < minx1)
{
minx1 = P[i].x;
}
if (P[i].x > maxx1)
{
maxx1 = P[i].x;
}
if (P[i].y < miny1)
{
miny1 = P[i].y;
}
if (P[i].y > maxy1)
{
maxy1 = P[i].y;
}
}
dx1 = (maxx1 + minx1)/2;
dy1 = (maxy1 + miny1)/2;
int xw = width / 2;
int yw = heigth / 2;
dxw = xw - dx1;
dyw = yw - dy1;
for (int i = 0; i < count; i++)
{
P[i].x += dxw;
P[i].y += dyw;
}
dxw = dxw * scale;
dyw = dyw * scale;
}
void PaintUI::DrawUI(CDC* pDC, CRect rect)
{
CPen pen5(PS_SOLID, 2, RGB(255, 0, 0));
CPen* oldpen5 = pDC->SelectObject(&pen5);
dscale = 3000;
for (int i = 2; i < count; i++)
{
double dsx, dsy, dex, dey;
dsx = (P[i].F * sin(P[i].Q)) * dscale + P[i].x;
dsy = (-P[i].F * cos(P[i].Q)) * dscale + P[i].y;
dex = (-P[i].F * sin(P[i].Q)) * dscale + P[i].x;
dey = (P[i].F * cos(P[i].Q)) * dscale + P[i].y;
pDC->MoveTo(dsx, dsy);
pDC->LineTo(dex, dey);
dsx = (-P[i].E * cos(P[i].Q)) * dscale + P[i].x;
dsy = (-P[i].E * sin(P[i].Q)) * dscale + P[i].y;
dex = (P[i].E * cos(P[i].Q)) * dscale + P[i].x;
dey = (P[i].E * sin(P[i].Q)) * dscale + P[i].y;
pDC->MoveTo(dsx, dsy);
pDC->LineTo(dex, dey);
double ex, fy;
ex = P[i].E;
fy = 0;
pDC->MoveTo(dsx, dsy);
for (int j = 6; j <= 360; j += 6)
{
ex = (dscale * P[i].E) * cos((j * acos(-1)) / 180);
fy = (dscale * P[i].F) * sin((j * acos(-1)) / 180);
dex = ex * cos(P[i].Q) - fy * sin(P[i].Q) + P[i].x;
dey = fy * cos(P[i].Q) + ex * sin(P[i].Q) + P[i].y;
pDC->LineTo(dex, dey);
}
}
pDC->SelectObject(oldpen5);
CPen pen(PS_SOLID, 2, RGB(30, 160, 220));
CPen* oldpen = pDC->SelectObject(&pen);
for (int i = 0; i < count1; i++)
{
int a = foundP(L[i].start);
int b = foundP(L[i].end);
pDC->MoveTo(P[a].x, P[a].y);
pDC->LineTo(P[b].x, P[b].y);
}
CPen pen1(PS_SOLID, 2, RGB(223, 113, 0));
CPen* oldpen1 = pDC->SelectObject(&pen1);
for (int i = 0; i < wdcount; i++)
{
int x = (heigth) / 40+1;
for (int j = 0; j < x; j++)
{
pDC->MoveTo(len * (i + 1), 40 * j);
pDC->LineTo(len * (i + 1), 40 * j+25);
}
}
for (int i = 0; i < heicount; i++)
{
int x = (width) / 40 + 1;
for (int j = 0; j < x; j++)
{
pDC->MoveTo(40 * j, change(len * (i + 1)));
pDC->LineTo(40 * j + 25, change(len * (i + 1)));
}
}
CPen pen2(PS_SOLID, 2, RGB(0, 0, 0));
CPen* oldpen2 = pDC->SelectObject(&pen2);
pDC->MoveTo(0,rect.bottom);
pDC->LineTo(0, rect.top);
pDC->MoveTo(0, rect.bottom);
pDC->LineTo(rect.right, rect.bottom);
pDC->MoveTo(0, rect.top);
pDC->LineTo(8, rect.top+10);
pDC->MoveTo(rect.right, rect.bottom);
pDC->LineTo(rect.right-10, rect.bottom-8);
pDC->TextOutW(8, rect.top + 13, _T("X"));
pDC->TextOutW(rect.right - 20, rect.bottom - 20, _T("Y"));
pDC->TextOutW(3, rect.bottom-20, _T("O"));
CPen pen3(PS_SOLID, 2, RGB(255, 0, 0));
CPen* oldpen3 = pDC->SelectObject(&pen3);
CFont font;//创建字体
font.CreateFontW(
20, // nHeight
0, // nWidth
0, // nEscapement
0, // nOrientation*
FW_NORMAL, // nWeight
FALSE, // bItalic
FALSE, // bUnderline
0, // cStrikeOut
ANSI_CHARSET, // nCharSet
OUT_DEFAULT_PRECIS, // nOutPrecision
CLIP_DEFAULT_PRECIS, // nClipPrecision
DEFAULT_QUALITY, // nQuality
DEFAULT_PITCH | FF_SWISS, // nPitchAndFamily
_T("宋体")); // lpszFacename*/)
CFont* pOldFont = pDC->SelectObject(&font);
for (int i = 0; i < count; i++)
{
pDC->TextOutW(P[i].x+3, P[i].y-2,P[i].ID);
}
for (int i = 0; i < wdcount+1; i++)
{
CString str;
str.Format(_T("%d"), (int)(Ox+len1*i-len1-dxw));
pDC->TextOutW(len*i, rect.bottom + 3, str);
}
for (int i = 0; i < heicount+1; i++)
{
CString str;
str.Format(_T("%d"), (int)(Oy + len1 * i-dyw));
pDC->TextOutW(-50, rect.bottom - len*i, str);
}
CPen pen4(PS_SOLID, 2, RGB(0, 255, 64));
CPen* oldpen4 = pDC->SelectObject(&pen4);
for (int i = 0; i < count; i++)
{
if (P[i].flag == 0) {
pDC->MoveTo(P[i].x, P[i].y+10);
pDC->LineTo(P[i].x+7, P[i].y-4);
pDC->LineTo(P[i].x - 7, P[i].y - 4);
pDC->LineTo(P[i].x, P[i].y + 10);
}
else
{
CBrush fillb;
fillb.CreateSolidBrush(RGB(255, 0, 0));
pDC->SelectObject(&fillb);
pDC->Ellipse(P[i].x-4,P[i].y-4, P[i].x+4, P[i].y+4);
fillb.DeleteObject();
}
}
}
工具类:
Datapoint .h:
#pragma once
enum class Type
{
Know,
UnKnow
};
struct Point
{
CString index;//点号
double x;
double y;
int flag=1;//是否为已知点
int flag1=1;//是否为控制点
double mx=0, my=0, mk=0;//点位误差
double Q=0 ,E=0 , F=0;//误差椭圆参数
};
Dataline .h:
#pragma once
#include <afxstr.h>
struct Line
{
CString start;//边的起始点
CString end;//边的终点
double length;//边长
double amangle;//方位角度值
};
Dataangle.h:
#pragma once
struct Angle
{
CString startP;//测站起始点号
CString endP;//测站重点号
double angle;//角度值
int num;//找准测站号
};
struct Ang
{
CString startP;//角度起始点号
CString mid;//角度中点号
CString endP;//角度终点号
double angle;//角度值
};
Matrix 矩阵类:
.h:
#pragma once
class CMatrix
{
public:
CMatrix(int row = 3, int col = 3);
// copy constructor
CMatrix(const CMatrix& m);
~CMatrix(void);
private:
double** dMatData;//保存矩阵元素数据的二维数组
int iRow;//矩阵的行
int iCol;//矩阵的列
public:
int Row() const { return iRow; }//返回行
int Col() const { return iCol; }//返回列
void SetSize(int row, int col);//调整数组的大小,原有数据不变(未测试)
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();//生成单位矩阵
CString look();
};
.cpp:
#include "pch.h"
#include "Matrix.h"
#include "math.h"
CMatrix::CMatrix(int row, int col)
{
iRow = row;
iCol = col;
dMatData = new double* [row];
for (int i = 0; i < row; i++)
{
dMatData[i] = new double[col];
for (int j = 0; j < col; j++)
{
dMatData[i][j] = 0;
}
}
}
// copy constructor,
//拷贝构造函数的作用:
//(1)以类对象作为函数参数传值调用时;
//(2)函数返回值为类对象;
//(3)用一个已定义的对象去初始化一个新对象时;
CMatrix::CMatrix(const CMatrix& m)
{
iRow = m.Row();
iCol = m.Col();
dMatData = new double* [iRow];
for (int i = 0; i < iRow; i++)
{
dMatData[i] = new double[iCol];
// for(int j=0;j<iCol;j++)
{
memcpy(dMatData[i], m.dMatData[i], sizeof(double) * iCol);
}
}
}
CMatrix::~CMatrix(void)
{
for (int i = 0; i < iRow; i++)
{
delete[] dMatData[i];
}
delete[] dMatData;
}
//返回数组元素(引用返回)
double& CMatrix::operator () (int row, int col)
{
if (row >= iRow || col >= iCol)
{
throw("CMatrix::operator(): Index out of range!");
}
return dMatData[row][col];
}
返回数组元素(重载)
double CMatrix::operator () (int row, int col) const
{
if (row >= iRow || col >= iCol)
{
throw("CMatrix::operator(): Index out of range!");
}
return dMatData[row][col];
}
//重载预算符+
CMatrix operator + (const CMatrix& m1, const CMatrix& m2)
{
if ((m1.Col() != m2.Col()) || (m1.Row() != m2.Row()))
{
throw("CMatrix::operator+: The two matrix have different size!");
}
CMatrix matTmp(m1.Row(), m1.Col());
for (int i = 0; i < m1.Row(); i++)
{
for (int j = 0; j < m1.Col(); j++)
{
matTmp(i, j) = m1(i, j) + m2(i, j);
}
}
return matTmp;
}
//重载赋值运算符=,当左右两边矩阵的大小不相等时,
//以右边的大小为基准,调整左边矩阵的大小
CMatrix& CMatrix::operator = (const CMatrix& m)
{
//revised in 2011-4-1, by Daiwujiao
// if(iRow!=m.Row()||iCol!=m.Col())
//{
// throw( "CMatrix::operator=: The two matrix have different size!");
//}
if (iRow != m.Row() || iCol != m.Col())
{
SetSize(m.Row(), m.Col());
}
for (int i = 0; i < iRow; i++)
{
for (int j = 0; j < iCol; j++)
{
dMatData[i][j] = m(i, j);
}
}
return *this;
}
//调整矩阵大小,原有值不变
void CMatrix::SetSize(int row, int col)
{
if (row == iRow && col == iCol)
{
return;
}
double** rsData = new double* [row];
for (int i = 0; i < row; i++)
{
rsData[i] = new double[col];
for (int j = 0; j < col; j++)
{
rsData[i][j] = 0;
}
}
int minRow = (iRow > row) ? row : iRow;
int minCol = (iCol > col) ? col : iCol;
int colSize = minCol * sizeof(double);
for (int i = 0; i < minRow; i++)
{
memcpy(rsData[i], dMatData[i], colSize);
}
for (int i = 0; i < minRow; i++)
{
delete[] dMatData[i];
}
delete[] dMatData;
dMatData = rsData;
iRow = row;
iCol = col;
return;
}
//重载预算符-
CMatrix operator - (const CMatrix& m1, const CMatrix& m2)
{
if ((m1.Col() != m2.Col()) || (m1.Row() != m2.Row()))
{
throw("CMatrix::operator-: The two matrix have different size!");
}
CMatrix matTmp(m1.Row(), m1.Col());
for (int i = 0; i < m1.Row(); i++)
{
for (int j = 0; j < m1.Col(); j++)
{
matTmp(i, j) = m1(i, j) - m2(i, j);
}
}
return matTmp;
}
//重载预算符*,两个矩阵相乘,m1的列要等于m2的行
CMatrix operator * (const CMatrix& m1, const CMatrix& m2)
{
if ((m1.Col() != m2.Row()))
{
throw("CMatrix::operator*: The col of matrix m1 doesn't equ to row of m2 !");
}
CMatrix matTmp(m1.Row(), m2.Col());
for (int i = 0; i < m1.Row(); i++)
{
for (int j = 0; j < m2.Col(); j++)
{
for (int k = 0; k < m2.Row(); k++)
{
matTmp(i, j) += m1(i, k) * m2(k, j);
}
}
}
return matTmp;
}
//重载预算符*,矩阵右乘一个数
CMatrix operator * (const CMatrix& m1, const double& num)
{
CMatrix matTmp(m1.Row(), m1.Col());
for (int i = 0; i < m1.Row(); i++)
{
for (int j = 0; j < m1.Col(); j++)
{
matTmp(i, j) = m1(i, j) * num;
}
}
return matTmp;
}
//重载预算符*,矩阵左乘一个数
CMatrix operator * (const double& num, const CMatrix& m1)
{
CMatrix matTmp(m1.Row(), m1.Col());
for (int i = 0; i < m1.Row(); i++)
{
for (int j = 0; j < m1.Col(); j++)
{
matTmp(i, j) = m1(i, j) * num;
}
}
return matTmp;
}
//矩阵转置
CMatrix operator ~ (const CMatrix& m)
{
CMatrix matTmp(m.Col(), m.Row());
for (int i = 0; i < m.Row(); i++)
for (int j = 0; j < m.Col(); j++)
{
matTmp(j, i) = m(i, j);
}
return matTmp;
}
//矩阵求逆
//采用选全主元法
CMatrix CMatrix::Inv()
{
if (iRow != iCol)
{
throw("待求逆的矩阵行列不相等!");
}
int i, j, k, vv;
CMatrix InvMat(iRow, iRow);
//复制矩阵
InvMat = *this;
int* MainRow = new int[iRow];
int* MainCol = new int[iRow];//用于记录主元素的行和列
double dMainCell;//主元元素的值
double dTemp;//临时变量
for (k = 0; k < iRow; k++)
{
dMainCell = 0;
//选全主元
for (i = k; i < iRow; i++)
{
for (j = k; j < iRow; j++)
{
dTemp = fabs(InvMat(i, j));
if (dTemp > dMainCell)
{
dMainCell = dTemp;
MainRow[k] = i;
MainCol[k] = j;
}
}
}
if (fabs(dMainCell) < 0.0000000000001)//矩阵秩亏,不能求逆
{
throw("矩阵秩亏");
}
if (MainRow[k] != k)//交换行
{
for (j = 0; j < iRow; j++)
{
vv = MainRow[k];
dTemp = InvMat(k, j);
InvMat(k, j) = InvMat(vv, j);
InvMat(vv, j) = dTemp;
}
}
if (MainCol[k] != k)//交换列
{
for (i = 0; i < iRow; i++)
{
vv = MainCol[k];
dTemp = InvMat(i, k);
InvMat(i, k) = InvMat(i, vv);
InvMat(i, vv) = dTemp;
}
}
InvMat(k, k) = 1.0 / InvMat(k, k);//计算乘数
for (j = 0; j < iRow; j++) //计算主行
{
if (j != k)
{
InvMat(k, j) = InvMat(k, j) * InvMat(k, k);
}
}
for (i = 0; i < iRow; i++)//消元
{
if (i != k)
{
for (j = 0; j < iRow; j++)
{
if (j != k)
{
InvMat(i, j) -= InvMat(i, k) * InvMat(k, j);
}
}
}
}
for (i = 0; i < iRow; i++)//计算主列
{
if (i != k)
{
InvMat(i, k) = -InvMat(i, k) * InvMat(k, k);
}
}
}
for (k = iRow - 1; k >= 0; k--)
{
if (MainCol[k] != k)// 交换行
{
for (j = 0; j < iRow; j++)
{
vv = MainCol[k];
dTemp = InvMat(k, j);
InvMat(k, j) = InvMat(vv, j);
InvMat(vv, j) = dTemp;
}
}
if (MainRow[k] != k)//交换列
{
for (i = 0; i < iRow; i++)
{
vv = MainRow[k];
dTemp = InvMat(i, k);
InvMat(i, k) = InvMat(i, vv);
InvMat(i, vv) = dTemp;
}
}
}
delete[] MainRow;
delete[] MainCol;
return InvMat;
}
//单位化矩阵
void CMatrix::Unit()
{
for (int i = 0; i < iRow; i++)
{
for (int j = 0; j < iCol; j++)
{
dMatData[i][j] = (i == j) ? 1 : 0;
}
}
}
CString CMatrix::look()
{
CString result1, result2;
for (int i = 0; i < iRow; i++)
{
for (int j = 0; j < iCol; j++)
{
result1.Format(_T("%lf "), dMatData[i][j]);
result2 = result2 + result1;
if (j == iCol - 1)
{
result1.Format(_T("\r\n"));
result2 = result2 + result1;
}
}
}
return result2;
}