coordinate_system.c

本程序适用于Gps坐标系和北京54或西安坐标系的转换

/*****************coordinate_system.c ***********************/
/**********坐标系统统一函数(coordinate_system.c)*************/
/*  编写人：吴光跃 时间：20060316****************************/
/*  修改人：吴光跃 时间：20060323****************************/
/*  测试人：       时间：        ****************************/


/**************************************************************************/
#include <stdio.h>

#include "coordinate_system.h"

/**************************************************************************/
/*布尔莎(Bursa)模型
  X_wgs84=X0  +(1+u)*X_cgs2000     +ex*Y_cgs2000      -ey*Z_cgs2000;
  Y_wgs84=Y0     -ex*X_cgs2000  +(1+u)*Y_cgs2000      +ez*Z_cgs2000;
  Z_wgs84=Z0     +ey*X_cgs2000     -ez*Y_cgs2000   +(1+u)*Z_cgs2000;
其中(X_wgs84,Y_wgs84,Z_wgs84)为WGS-84坐标系下坐标，
    (X_cgs2000,Y_cgs2000,Z_cgs2000)为CGS2000坐标系下坐标*/

double X0,Y0,Z0,U,ex,ey,ez;  /*全局变量，布尔莎(Bursa)模型待标定的方程系数*/

/**************************************************************************/
/*标定方程系数函数(seek_coefficient.c)
  通过在两个坐标系中都具有精确已知坐标的公共点加以确定，
  分别取CGS2000坐标系下如下P1—P7点，及其对应在WGS-84坐标系下的点进行标定。
    公共点： CGS2000坐标系下坐标：  WGS-84坐标系下坐标：
      P1：   (0 ,0 ,0)              (X_wgs84[0],Y_wgs84[0],Z_wgs84[0])
      P2：   (X_cgs2000[1],0,0)     (X_wgs84[1],Y_wgs84[1],Z_wgs84[1])
      P3：   (-X_cgs2000[1],0,0)    (X_wgs84[2],Y_wgs84[2],Z_wgs84[2]) 
      P4：   (0,Y_cgs2000[3],0)     (X_wgs84[3],Y_wgs84[3],Z_wgs84[3])
      P5：   (0,-Y_cgs2000[3],0)    (X_wgs84[4],Y_wgs84[4],Z_wgs84[4])
      P6：   (0,0,Z_cgs2000[5])     (X_wgs84[5],Y_wgs84[5],Z_wgs84[5])
      P7：   (0,0,-Z_cgs2000[5])     (X_wgs84[6],Y_wgs84[6],Z_wgs84[6])
  以上两坐标系下7个公共点可以标定出模型方程的全部系数*/

void seek_coefficient()   /*标定模型方程系数函数*/
{
	float X_cgs2000[7]={0,1,-1,0,0,0,0}; /*7个公共点在CGS2000坐标系下X轴坐标值(待定值)*/
	float Y_cgs2000[7]={0,0,0,1,-1,0,0}; 
	float Z_cgs2000[7]={0,0,0,0,0,1,-1};
	float X_wgs84[7]={1,2,3,4,5,6,7};    /*7个公共点在WGS-84坐标系下X轴坐标值(待定值)*/
	float Y_wgs84[7]={1,7,6,5,4,3,2};
	float Z_wgs84[7]={1,3,4,5,6,7,8};

	X0=X_wgs84[0];  /*X0,Y0,Z0为平移参数*/
	Y0=Y_wgs84[0];
	Z0=Z_wgs84[0];

	U=((X_wgs84[1]-X_wgs84[2])/(2*X_cgs2000[1])+   /*U为尺度参数*/
	   (Y_wgs84[3]-Y_wgs84[4])/(2*Y_cgs2000[3])+
	   (Z_wgs84[5]-Z_wgs84[6])/(2*Z_cgs2000[5])-
	   3)/3;

	ex=(-(Y_wgs84[1]-Y_wgs84[2])/(2*X_cgs2000[1])+    /*ex,ey,ez为旋转参数*/
		(X_wgs84[3]-X_wgs84[4])/(2*Y_cgs2000[3]))/2;

	ey=((Z_wgs84[1]-Z_wgs84[2])/(2*X_cgs2000[1])-
		(X_wgs84[5]-X_wgs84[6])/(2*Z_cgs2000[5]))/2;

	ez=(-(Z_wgs84[3]-Z_wgs84[4])/(2*Y_cgs2000[3])+
		(Y_wgs84[5]-Y_wgs84[6])/(2*Z_cgs2000[5]))/2;

	/*printf("标定出模型方程的系数如下：\n");*/
	/*printf("X0=%f,Y0=%f,Z0=%f,\n U=%f,ex=%f,ey=%f,ez=%f\n",X0,Y0,Z0,U,ex,ey,ez);*/
}


/*求矩阵乘A[3][3]*B[3]函数(Matrix3_Multiply)*/
void Matrix3_Multiply(double A[3][3],double B[3],double Y[3])  
{
	int i;

	for(i=0;i<=2;i++)
		Y[i]=A[i][0]*B[0]+A[i][1]*B[1]+A[i][2]*B[2]; /*矩阵乘结果存入Y[3]*/
}

/*double X_WGS84=0,Y_WGS84=0,Z_WGS84=0;*/
/*将CGS2000坐标转换为WGS-84坐标函数*/
void coordinate_CGS2000_to_WGS84(double X_CGS2000,double Y_CGS2000,double Z_CGS2000,
								 double *X_WGS84,double *Y_WGS84,double *Z_WGS84)
{
	double A[3][3]={{1+U,ex,-ey},{-ex,1+U,ez},{ey,-ez,1+U}};
	double B[3]={X_CGS2000,Y_CGS2000,Z_CGS2000};
	double Y[3];

    Matrix3_Multiply(A,B,Y);  /*矩阵乘结果存入矩阵Y[3]*/

    *X_WGS84=X0+Y[0];   /*根据布尔莎(Bursa)模型转换*/
    *Y_WGS84=Y0+Y[1];
    *Z_WGS84=Z0+Y[2];
	/*printf("(%f,  %f,  %f)\n",*X_WGS84,*Y_WGS84,*Z_WGS84);*/
}

/**************************************************************************/