travernet.cpp

利用已知控制点坐标和方位角

// Travernet.cpp: implementation of the Travernet class.
//
//////////////////////////////////////////////////////////////////////

#include "stdafx.h"
#include "xinyan.h"
#include "Travernet.h"
#include "math.h"
#include "constfile.h"

#ifdef _DEBUG
#undef THIS_FILE
static char THIS_FILE[]=__FILE__;
#define new DEBUG_NEW
#endif

//////////////////////////////////////////////////////////////////////
// Construction/Destruction
//////////////////////////////////////////////////////////////////////

Travernet::Travernet()
{

}

Travernet::~Travernet()
{

}
bool Travernet::ReadData(CString filename)
{
	CStdioFile fp;
	CFileException e;
	angobs *p1; 	
	disobs *p3;

	int i;

	if(!fp.Open( filename,
		CFile::modeRead|CFile::typeText,NULL))
	{
		AfxGetApp()->m_pMainWnd->MessageBox(
			"数据文件不存在或数据文件错!",
			"提示",MB_OK|MB_ICONWARNING);
			return FALSE;
	}
    int MAXLINE=254;
	char strBuffer[255],buff[255];
	fp.ReadString(buff,MAXLINE);		
	sscanf(buff,"%s",strBuffer);
	if(strcmpi(strBuffer,"TRAVERNET"))
	{
		AfxGetApp()->m_pMainWnd->MessageBox(
			"非测量数据处理文件!",
			"提示",MB_OK|MB_ICONWARNING);
		fp.Close();
		return FALSE;
	}
	fp.ReadString(buff,MAXLINE);
	sscanf(buff,"%d %d %d %d",
		&totalN,&knownN,&azi_N,&dis_N);
	 unknownN= totalN-knownN;
	fp.ReadString(buff,MAXLINE);
	sscanf(buff,"%lf %lf %lf",
		&mb,&md,&ms);
	if(!Allocmem())
	{
		AfxMessageBox("内存申请失败",MB_OK|MB_ICONERROR);
		return false;
	}
	for(i=0;i<knownN;i++)//读入已知点坐标
	{
		fp.ReadString(buff,MAXLINE);
		sscanf(buff,"%lf %lf",x+i,y+i);
	}
	p1=anglevalue;//读入角度观测值
	for(i=0;i<azi_N;i++)
	{
		fp.ReadString(buff,MAXLINE);
		sscanf(buff,"%d %d %lf\n",&p1->from,&p1->to,&p1->angle);
		p1->angle=angle_r(p1->angle);
		p1=p1->next;
	}
    p3=distancevalue;//读入距离观测值
	for(i=0;i<dis_N;i++)
	{
		fp.ReadString(buff,MAXLINE);
		sscanf(buff,"%d%d%lf\n",&p3->from,&p3->to,&p3->distance);
		p3=p3->next;
	}
	fp.Close();
	return true;

}
bool Travernet::Allocmem()
{
	int i;
	x=(double*)calloc(totalN,sizeof(double));if(!x) return false;//近似坐标
	y=(double*)calloc(totalN,sizeof(double));if(!y) return false;
	ed=(int*)calloc(totalN,sizeof(int));if(!ed) return false;
	v=(double*)calloc(azi_N+dis_N,sizeof(double));if(!v) return false;
	l=(double*)calloc(azi_N+dis_N,sizeof(double));if(!l) return false;
	X=(double*)calloc(totalN,sizeof(double));if(!X) return false;//改正后坐标
	Y=(double*)calloc(totalN,sizeof(double));if(!Y) return false;

	mx=(double*)calloc(unknownN,sizeof(double));if(!mx) return false;
	my=(double*)calloc(unknownN,sizeof(double));if(!my) return false;
	m=(double*)calloc(unknownN,sizeof(double));if(!m) return false;
    
	
	angobs *p1,*p2;//为方向值链表申请空间
	anglevalue=(angobs*)calloc(1,sizeof(angobs));
	anglevalue->next=NULL;
	p1=anglevalue;
	for(i=0;i<azi_N-1;i++)
	{
		p2=(angobs*)calloc(1,sizeof(angobs));
		p2->next=NULL;
		p1->next=p2;
		p1=p2;
	}
    disobs *p3,*p4;//为距离值链表申请空间
	distancevalue=(disobs*)calloc(1,sizeof(disobs));
	distancevalue->next=NULL;
	p3=distancevalue;
	for(i=0;i<dis_N-1;i++)
	{
		p4=(disobs*)calloc(1,sizeof(disobs));
		p4->next=NULL;
		p3->next=p4;
		p3=p4;
	}
    return true;
}
void Travernet::Precalcoor()//计算近似坐标主函数
{
	int i;
	for(i=0;i<knownN;i++)
		ed[i]=1;
	caled=knownN;
	while(caled<totalN)
	{
		for(i=knownN;i<totalN;i++)
		{
			if(ed[i]!=1)
				if(Canbecal(i)==1)
				{ed[i]=1;
				caled++;}
		}
	}
}
int Travernet::Canbecal(int i)//计算近似坐标
{
	angobs*p1;
	disobs*p3;
    int temp1,temp2;
	int cond;
	double ang1,ang2,dtemp;
	p1=anglevalue;
	cond=4;
		while(p1!=NULL)
		{
			if(p1->to==i)
			{
				if(ed[p1->from]==1)
				{
					cond-=1;
					temp1=p1->from;
					ang1=p1->angle;
					break;
				}

			}
			p1=p1->next;
		}
		if(cond==4)
			return 4;
		p1=anglevalue;
		while(p1!=NULL)
		{
			if(p1->from==temp1)
			{
				if(ed[p1->to]==1)
				{
					cond-=1;
					ang2=p1->angle;
					temp2=p1->to;
					break;
				}
			}
			p1=p1->next;
		}
		if(cond==3)
			return 3;
		p3=distancevalue;
		while(p3!=NULL)
		{
			if(p3->from==temp1||p3->to==temp1)
				if(p3->to==i||p3->from==i)
			{
				cond-=1;dtemp=p3->distance;
				break;
			}
			p3=p3->next;
		}


		if(cond=1)
		{
			double dx,dy,b;
			dx=x[temp2]-x[temp1];
			dy=y[temp2]-y[temp1];
			b=dxy_az(dx,dy)+(ang1-ang2)+PI2;
			x[i]=x[temp1]+dtemp*cos(b);
			y[i]=y[temp1]+dtemp*sin(b);
		} 
		return 1;
}
void Travernet::CalunsN()//计算未设站数
{
  unsN=knownN;
  angobs*p1;
  int i;
  for(i=0;i<knownN;i++)
  {
    p1=anglevalue; 
	  while(p1!=NULL)
	  {
		  if(p1->from==i)
		  {unsN-=1;
		  break;
          }
		  else
		  p1=p1->next;
	  }
  }
   BN=2*unknownN+totalN-unsN;
	 
}



bool Travernet::Allocmem1()
{
	int i;
	xx=(double*)calloc(BN,sizeof(double));if(!xx) return false;
	U=(double*)calloc(BN,sizeof(double));if(!U) return false;
	p=(double*)calloc(azi_N+dis_N,sizeof(double));if(!p) return false;
    z=(double*)calloc(totalN-unsN,sizeof(double));if(!z) return false;
	B=(double**)calloc(azi_N+dis_N,sizeof(double));
	for(i=0;i<azi_N+dis_N;i++)
		B[i]=(double*)calloc(BN,sizeof(double));

    
	N=(double**)calloc(BN,sizeof(double));
	for(i=0;i<BN;i++)
		N[i]=(double*)calloc(BN,sizeof(double));

    return true;
}
void Travernet::Calz()//计算定向角平均值
{
	int i;
	for(i=unsN;i<totalN;i++)
	{
		int temp=0;
		z[i]=0;
		angobs*p1=anglevalue;
		while(p1!=NULL)
		{ 
			if(p1->from==i)
			{
				double dtemp,dx,dy;
				dx=x[p1->to]-x[p1->from];
				dy=y[p1->to]-y[p1->from];

				dtemp=dxy_az(dx,dy)-p1->angle;
				if(dtemp<0)
					dtemp+=PI2;
				z[i]+=dtemp;
				temp+=1;
			}
			p1=p1->next;
		}
		if(temp!=0)
			z[i]/=temp;
	}
}


void Travernet::Constitute()//组成误差方程式阵，常数项阵，权阵
{
	angobs*p1=anglevalue;
	disobs*p3=distancevalue;
	
	int i,j,k;
	double dx,dy,T,D,a,b;


	for(i=0;i<azi_N;i++)
	{
		dx=x[p1->to]-x[p1->from];
		dy=y[p1->to]-y[p1->from];
		T=dxy_az(dx,dy);
		D=dxy_s(dx,dy);
		a=206265*sin(T)/(100*D);
		b=-206265*cos(T)/(100*D);

		if(p1->from>=knownN)
		{
			B[i][2*(p1->from-knownN)]=a;
			B[i][2*(p1->from-knownN)+1]=b;
			
		}
		if(p1->to>=knownN)
		{
			B[i][2*(p1->to-knownN)]=-a;
			B[i][2*(p1->to-knownN)+1]=-b;
		}
	
		B[i][2*unknownN+p1->from-unsN]=-1;
		l[i]=p1->angle+z[p1->from];
		if(l[i]>=PI2)
			l[i]-=PI2;
		l[i]-=T;
		l[i]*=(180/PI*3600);
		p[i]=1;

		
		p1=p1->next;
		
		
	}
	for(i=azi_N;i<azi_N+dis_N;i++)
	{
		double T1;
		dx=x[p3->to]-x[p3->from];
		dy=y[p3->to]-y[p3->from];
		T1=dxy_az(dx,dy);
	
		if(p3->from>=knownN)
		{
			B[i][2*(p3->from-knownN)]=-cos(T1);
			B[i][2*(p3->from-knownN)+1]=-sin(T1);
		}
		if(p3->to>=knownN)
		{
			B[i][2*(p3->to-knownN)]=cos(T1);
			B[i][2*(p3->to-knownN)+1]=sin(T1);
		}
		l[i]=(p3->distance-dxy_s(dx,dy))*100.0;
		p[i]=pow(mb,2)/pow((md/10+ms*p3->distance/10000),2);
	
		p3=p3->next;
	
	
	}
	//组法方程阵及自由项阵
	for(i=0;i<BN;i++)
	{
		for(j=0;j<BN;j++)
		{
			for(k=0;k<azi_N+dis_N;k++)
			{	N[i][j]+=B[k][i]*B[k][j]*p[k];
			}
		}
	}
	for(i=0;i<BN;i++)
	{
		for(j=0;j<azi_N+dis_N;j++)
		{
		  U[i]+=p[j]*B[j][i]*l[j];
		}
	}
}
void Travernet::inverse(double**N,int BN)//矩阵求逆
{
	double**NC, temp;
	int i,j,k;
    NC=(double**)calloc(BN,sizeof(double));
	for(i=0;i<BN;i++)
		NC[i]=(double*)calloc(BN,sizeof(double));
    for(i=0;i<BN;i++)
	NC[i][i]=1;
    for(i=0;i<BN-1;i++)
	{
		for(j=i+1;j<BN;j++)
		{
			temp=N[j][i]/N[i][i];
			for(k=0;k<BN;k++)
			{N[j][k]-=temp*N[i][k];
			NC[j][k]-=temp*NC[i][k];
			}
		}
	}
	for(i=BN-1;i>0;i--)
	{
		for(j=i-1;j>=0;j--)
		{
			temp=N[j][i]/N[i][i];
			for(k=0;k<BN;k++)
			{N[j][k]-=temp*N[i][k];
			NC[j][k]-=temp*NC[i][k];
			}
		}
	}
	for(i=0;i<BN;i++)
	{
		temp=N[i][i];
		for(j=0;j<BN;j++)
			N[i][j]=NC[i][j]/temp;
	}
	for(i=0;i<BN;i++)
	free(NC[i]);
}
void Travernet::Process()//精度评定
{
	int i,j;
	inverse(N,BN);
	for(i=0;i<BN;i++)
	{
		for(j=0;j<BN;j++)
		{xx[i]+=N[i][j]*U[j];
		}
	}
	for(i=knownN;i<totalN;i++)
	{
		X[i]=x[i]+xx[2*(i-knownN)]/100.;
		Y[i]=y[i]+xx[2*(i-knownN)+1]/100.;
	
	}
	pvv=0;
	for(i=0;i<azi_N+dis_N;i++)
	{
		 double temp=0;
		for(j=0;j<BN;j++)
			temp+=B[i][j]*xx[j];
		temp-=l[i];
		pvv+=temp*temp*p[i];
	}
	m0=sqrt(pvv/(azi_N+dis_N-BN));
		for(i=0;i<2*unknownN;i+=2)
		{
			mx[i/2+knownN]=m0*sqrt(N[i][i]);
			my[i/2+knownN]=m0*sqrt(N[i+1][i+1]);
			m[i/2+knownN]=sqrt(mx[i/2+knownN]*mx[i/2+knownN]+my[i/2+knownN]*my[i/2+knownN]);
		}
		prompt();
}

 

bool Travernet::WriteData(CString filename)
{
	CStdioFile f32;
	CFileException e32;

	if(!f32.Open(filename,CFile::modeCreate|CFile::modeWrite|CFile::typeText,&e32))
	{
		CString errorMessage="文件"+filename+"不能创建!";
		AfxMessageBox(errorMessage,MB_OK|MB_ICONSTOP,0);
		f32.Close();
		return FALSE;
	}
	CString strBuffer;
	strBuffer.Format("%s","导线网平差结果\n");
	f32.WriteString(strBuffer);

	strBuffer.Format("%s","近似坐标\n");
    f32.WriteString(strBuffer);
	for(int i=knownN;i<totalN;i++)
	{
		strBuffer.Format("x[%2d]=%10.5lf y[%2d]=%10.5lf\n",i,x[i],i,y[i]);
        f32.WriteString(strBuffer);
	}
   
	strBuffer.Format("%s","近似值改正数\n");
	f32.WriteString(strBuffer);
	for( i=0;i<BN;i++)
	{
		strBuffer.Format("xx[%2d]=%2.6lf\n",i,xx[i]);
		f32.WriteString(strBuffer);
	}
	strBuffer.Format("%s","改正后坐标\n");
	f32.WriteString(strBuffer);
	for(i=knownN;i<totalN;i++)
	{
		strBuffer.Format("X[%2d]=%10.5lf Y[%2d]=%10.5lf\n",i,X[i],i,Y[i]);
		f32.WriteString(strBuffer);
	}
	strBuffer.Format("[pvv]=%3.6lf\n",pvv);
	f32.WriteString(strBuffer);
	strBuffer.Format("m0=+ -%3.6lf\n",m0);
	f32.WriteString(strBuffer);
	for(i=knownN;i<totalN;i++)
	{
		strBuffer.Format("mx[%2d]=+ -%3.6lf my[%2d]=+ -%3.6lf m[%2d]=+ -%3.6lf\n",i,mx[i],i,my[i],i,m[i]);
		f32.WriteString(strBuffer);
	}
	f32.Close();
	return TRUE;
}