coordinate.cpp

读取GPS的Rinex观测文件

#include "coordinate.h"

void CartesianToGeodetic (PCRDGEODETIC pcg, PCRDCARTESIAN pcc,double dSemiMajorAxis, double dFlattening)//由笛卡尔坐标转换为大地坐标
//pcg：指向所转换出的大地坐标的指针
//pcc：指向待转换的笛卡尔坐标的指针
//dSemiMajorAxis：参考椭球的长半轴
//dFlattening：参考椭球的扁率
{
	double e2=1-pow(1-dFlattening,2);
    if(pcc->x>=0&&pcc->y>=0)
    pcg->longitude=atan(pcc->y/pcc->x);
    if(pcc->x>=0&&pcc->y<0)
    pcg->longitude=atan(pcc->y/pcc->x);
    if(pcc->x<0&&pcc->y>=0)
    pcg->longitude=PI+atan(pcc->y/pcc->x);
    if(pcc->x<0&&pcc->y<0)
    pcg->longitude=-PI+atan(pcc->y/pcc->x);

    double m_d_tan_B=pcc->z/sqrt(pow(pcc->x,2)+pow(pcc->y,2));
    double B_temp;//在迭代式时， 存储纬度的中间变量
    double N;
    pcg->latitude=atan(m_d_tan_B);

do
{
	B_temp=pcg->latitude;
	N=dSemiMajorAxis/sqrt(1-e2*pow(sin(pcg->latitude),2));
	m_d_tan_B=pcc->z/sqrt(pow(pcc->x,2)+pow(pcc->y,2))+
		N*e2*sin(pcg->latitude)/sqrt(pow(pcc->x,2)+pow(pcc->y,2));
	pcg->latitude=atan(m_d_tan_B);
}
while(fabs(pcg->latitude-B_temp)>0.0001);

N=dSemiMajorAxis/sqrt(1-e2*pow(sin(pcg->latitude),2));
pcg->height=sqrt(pow(pcc->x,2)+pow(pcc->y,2))/cos(pcg->latitude)-N;
}
void GeodeticToCartesian (PCRDCARTESIAN pcc, PCRDGEODETIC pcg,double dSemiMajorAxis, double dFlattening)//由大地坐标转换为笛卡尔坐标
//pcc：指向所转换出的笛卡尔坐标的指针
//pcg：指向待转换的大地坐标的指针
//dSemiMajorAxis：参考椭球的长半轴
//dFlattening：参考椭球的扁率
{
	double e2=1-pow(1-dFlattening,2);//第一偏心率的平方
	double N=dSemiMajorAxis/sqrt(1-e2*pow(sin(pcg->latitude),2));
	pcc->x=(N+pcg->height)*cos(pcg->latitude)*cos(pcg->longitude);
	pcc->y=(N+pcg->height)*cos(pcg->latitude)*sin(pcg->longitude);
	pcc->z=(N*(1-e2)+pcg->height)*sin(pcg->latitude);
}

void CartesianToTopocentric (PCRDTOPOCENTRIC pct,PCRDCARTESIAN pcc,PCRDCARTESIAN pccCenter,double dSemiMajorAxis,double dFlattening)
//由笛卡尔坐标转换为站心地平坐标

//pct：指向所转换出的站心地平坐标的指针
//pcc：指向待转换的笛卡尔坐标的指针
//pccCenter：指向站心的笛卡尔坐标的指针
//dSemiMajorAxis：参考椭球的长半轴
//dFlattening：参考椭球的扁率
{
	double delta_X,delta_Y,delta_Z;
	delta_X=pcc->x-pccCenter->x;
	delta_Y=pcc->y-pccCenter->y;
	delta_Z=pcc->z-pccCenter->z;
	PCRDGEODETIC pcg=new CRDGEODETIC;
	CartesianToGeodetic ( pcg,  pccCenter, dSemiMajorAxis,  dFlattening);//(X,Y,Z)->(B,L,H)站心转换
	pct->northing= -sin(pcg->latitude)*cos(pcg->longitude)*delta_X - sin(pcg->latitude)*sin(pcg->longitude)*delta_Y + cos(pcg->latitude)*delta_Z;
	pct->easting = -sin(pcg->longitude)*delta_X + cos(pcg->longitude)*delta_Y;
	pct->upping  =  cos(pcg->latitude)*cos(pcg->longitude)*delta_X + cos(pcg->latitude)*sin(pcg->longitude)*delta_Y + sin(pcg->latitude)*delta_Z;
}

void TopocentricToTopocentricPolar (PCRDTOPOCENTRICPOLAR pctp,PCRDTOPOCENTRIC pct)//由站心地平直角坐标转换为站心地平极坐标
//pct ：站心地平  直角坐标    的指针
//pctp：站心地平  极坐标      的指针
{
	pctp->range=sqrt(pow(pct->easting,2)+pow(pct->northing,2)+pow(pct->upping,2));
	pctp->azimuth=atan(pct->easting/pct->northing);
	pctp->elevation=atan(pct->upping/pctp->range);
}