clcx.cpp

该程序是我自己为测量外业中常用的导线简易平差程序

// aaa.cpp : Defines the entry point for the console application.
//
//
//#include "stdafx.h"
#include <stdio.h>
#include <math.h>
#include <string.h>
#define PI  3.1415926535897932384626433832795
#define exit
struct Pnt
{   
    char name[255]; 
    double x;
    double y;
} pnts[100];

struct Observe
{   
    Pnt p;
    char hname[255]; //后视
    char cname[255]; //测站
    char qname[255]; //前视
	double front_ang;//前视方向
    double back_ang;
    double front_dist;//前视距离
	double back_dist;
	double ang;    
    double dist;  
}obs[100];

double Azimuth(double x0, double y0, double x1, double y1)  //根据坐标计算坐标方位角
{
    double dx = x1 - x0;
    double dy = y1 - y0;
    return atan2(dy, dx) + (dy < 0 ? 1 : 0) *2*PI;
}
double DmsToRad(double dms)    //角度化弧度
{
    dms+=(dms>=0?1:-1)*1e-7;
    int deg=(int)dms;
	dms=(dms-deg)*100;
    int min=(int)dms;
    double sec=(dms-min)*100;
    return(deg+min/60.0+sec/3600.0)*PI/180;
}
double RadToDms(double rad)   //弧度化角度
{
    int f=(rad >=0 ?1:-1);
    rad = (rad+f*0.001/206265) * 180/PI;
    int deg = (int)rad;
    int min = (int)((rad - deg) * 60.0);
    double sec = (rad - deg - min/60.0) * 3600.0;
    return (deg + min/100.0 + sec/10000.0-(1e-7*f));
}
void dxdy(double azimuth, double dist,double &dx, double &dy)   //计算坐标增量
{
    dx = cos(azimuth) * dist;
    dy = sin(azimuth) * dist;
}
void Coordinate(double x0, double y0,double dx, double dy,   
                double vdx, double vdy,double &xi, double &yi)   //计算坐标
{
    xi = x0 + dx +vdx*0.001;
    yi = y0 + dy +vdy*0.001;
}
int shj_sheru(double shuju)//数据舍入函数
{ int A;
A = ( shuju >= 0 ? int (shuju + 0.5 )  : int (shuju - 0.5 ) );
return A;
}
int main(int argc, char* argv[])
{
    FILE * in, * out;
    char inFile[256], outFile[256];
	//char str1 = {"******"};
	double dms_front[100],dms_back[100];
    double dx[100], dy[100];
    double azimuth[100];
    double fB,fB_limit,VBi;
    double EBi=0,Es=0,Edx=0,Edy=0;
    double vdx[100],vdy[100];
    double fs,fx,fy;//vv = 0;
    int i,t,n,m=8,K;     
    if(argc == 1)
    {
        printf (" 请输入文件名: ");
        scanf  ("%s", inFile );
        printf (" 请输入cg名: " );
        scanf  ("%s", outFile );
    }
    else
        strcpy ( inFile , argv[1] );
    in = fopen ( inFile, "r" );
    if( in == NULL)
    {
        printf ( "不能创建你指定的文件！" );
        exit(-1);
    }
    out = fopen( outFile, "w" );
    if( out == NULL )
    {
        printf( "不能创建你指定的文件！" );
        exit(-1);
    }
    i = 0;
	while (true) 
	{
	    fscanf( in, "%s %lf %lf", pnts[i].name, &pnts[i].x, &pnts[i].y );
		if ( pnts[i].x  == 0 || pnts[i].y == 0 ) break; 
		i ++ ; 
	}
      t =  i +1;
	   i = 0;
    while (true) 
	{
        fscanf(in, "%s %s %s %lf %lf %lf %lf %lf %lf",
			obs[i].hname, obs[i].cname, obs[i].qname, &dms_front[i], &dms_back[i], &obs[i].front_dist, &obs[i].back_dist);
        obs[i].ang = DmsToRad ( dms_back[i] ) - DmsToRad ( dms_front[i] );
		if ( obs[i].ang > ( 2 * PI) )   obs[i].ang -= (2 * PI);
		if ( obs[i].ang < 0 )           obs[i].ang += (2 * PI); 
		EBi += obs[i].ang;
        if (obs[i].back_dist == 0) break;
        i++;
    }
	  n = i + 1;//n=5
	 for ( i = 0 ; i< (n-1); i++)
	   {    obs[i].dist = ( obs[i].back_dist + obs[i+1].front_dist ) / 2; }
	   azimuth[0] = Azimuth(pnts[0].x, pnts[0].y, pnts[1].x, pnts[1].y) ;  //已知边的坐标方位角
	   azimuth[n] = Azimuth(pnts[2].x, pnts[2].y, pnts[3].x, pnts[3].y) ;
	   fB = EBi + azimuth[0] - azimuth[n];
	   if( fB + n*PI >=0) fB = fB-n*PI;
	   else fB = fB + n*PI;
	   fB = RadToDms(fB) * 1e4;
	   fB_limit = 2 * m * sqrt(n);
	   VBi = -fB/n; 
       fprintf ( out, "                 ");
	   fprintf ( out, "                附和导线近似平差结果如下：\n\n");
	   fprintf ( out, " ***************************************************************************\n");
	   fprintf ( out, "                角度闭合差      fB = %-8d \n",shj_sheru (fB) );
	   fprintf ( out, "                角度允许闭合差  fB_limit = %-8d\n",shj_sheru (fB_limit) );
	   if ( fB < fB_limit )  fprintf( out, "               方位角闭合差未超限！\n");
	   else fprintf( out, "        方位角闭合差超限！\n");
	   fprintf ( out, "                角度改正数      VBi = %-8d\n\n",shj_sheru(VBi));
	   fprintf ( out, " ***************************************************************************\n");
	   fprintf ( out, "                 未知边的坐标方位角为(度分秒)：\n");
	   for ( i = 1; i < n; i++) 
	   {
		   azimuth[i] = azimuth[i-1] + obs[i-1].ang + DmsToRad(VBi*1e-4)-PI;
		   fprintf ( out, "         %s-%s: %10.4lf\n",obs[i-1].cname, obs[i-1].qname, RadToDms(azimuth[i]));
	   }
	   fprintf ( out, " ***************************************************************************\n");
	   fprintf ( out, "\n                        坐标增量：\n");
	   for ( i = 1; i < n; i++)                              //x,y的坐标增量
	   {
		   dxdy ( azimuth[i], obs[i-1].dist, dx[i], dy[i]);
		   fprintf ( out, " 测站-前视点 ： %s-%s: dx=%.3lf,dy=%.3lf\n", obs[i-1].cname, obs[i-1].qname, dx[i], dy[i]);
		   Edx += dx[i];
		   Edy += dy[i];
		   Es += obs[i-1].dist;
	   }
	   fprintf ( out, " ***************************************************************************\n");
	   fprintf ( out, "                 Edx=%.3lf\n                 Edy=%.3lf\n                 Es=%.3lf\n",Edx,Edy,Es);
	   fx = ( pnts[1].x + Edx-pnts[2].x) * 1000;
	   fy = ( pnts[1].y + Edy-pnts[2].y) * 1000;
	   fs = sqrt( fx*fx + fy*fy );
	   K = int ( Es / (fs*0.001) ); 
	   fprintf ( out, " ***************************************************************************\n");
	   fprintf ( out, "                     坐标增量闭合差：fx=%d,fy=%d\n",shj_sheru (fx),shj_sheru(fy) );
	   fprintf ( out, "                     全长闭合差：fs=%d\n",shj_sheru(fs) );
	   fprintf ( out, "                     全长相对闭合差：K=1/");
	   fprintf ( out, "%d\n",K);
	   fprintf ( out, "***************************************************************************\n");
	   fprintf ( out, "                      各坐标增量改正数：\n");
	   for ( i = 1; i < n; i++)
	   {
		   vdx[i] = ( -fx / Es * obs[i-1].dist);
		   vdy[i] = ( -fy / Es * obs[i-1].dist); 
		  
		   fprintf ( out, "                  %s-%s: vdx = %d, vdy = %d\n", 
			                     obs[i-1].cname, obs[i-1].qname, shj_sheru (vdx[i]),shj_sheru (vdy[i]) );
	   }
	   obs[0].p.x = pnts[1].x;
	   obs[0].p.y = pnts[1].y;
	   
	   fprintf ( out, " ***************************************************************************\n");
	   fprintf ( out, "                       待测点的坐标：\n");
	   for ( i = 1; i < n-1; i++)                   
	   {   
		   Coordinate(obs[i-1].p.x, obs[i-1].p.y, dx[i], dy[i],
			   vdx[i], vdy[i], obs[i].p.x, obs[i].p.y);
		   fprintf ( out, "              %s: x = %.3lf, y = %.3 lf\n",obs[i].cname, obs[i].p.x, obs[i].p.y);
	   }
     fprintf(out, " ***************************************************************************\n");
	   /*for (i =1;i<n-1;i++)
	 { vv += ( vdx[i] * vdx[i] + vdy[i]*vdy[i] );
	 
	 }
	 vv = sqrt( vv / (2*(n-1)) );
		 fprintf(out, "                   测边中误差为：%.3lf\n",vv);*/

	   fclose ( in );
	   fclose ( out );
	   return 0;
}