meteoriccorrect.cpp

卫星单点定位程序

#include "MeteoricCorrect.h"

//计算对流层延迟量
double GetTropDelay (double *TropDelay,			//对流层延迟量
					 PMETEODATA pmd,			//气象数据
					 PCRDCARTESIAN pcrdSite,	//测站坐标
					 PCRDCARTESIAN pcrdSat)		//卫星坐标
{
/*	double	RH0	= 50;					//参考面的相对湿度
	double	p0	= 1013.25;				//参考面的气压
	double	T0	= 0;					//参考面的干温
	double	H0	= 0;					//参考面的高度
	double	T, p, h, RH, t;				//测站点的气象数据
	double	e;							//水气压
	double	Kd, Kw, hd, hw;
	CRDGEODETIC			cg;				//测站点的大地坐标
*/
	CRDTOPOCENTRIC		ct;				//卫星在测站的站心地平坐标系中的坐标
	CRDTOPOCENTRICPOLAR	ctp;			//卫星在测站的站心地平极坐标系中的坐标
	
	//计算高度角
	CartesianToTopocentric(&ct, pcrdSat, pcrdSite,
		6378137.0000000000, 1/298.257223563);		//以 WGS-84 参考椭球参数计算
	TopocentricToTopocentricPolar(&ctp, &ct);

	//判断 pmd 是否为参考面上的数据，是则计算测站点气象参数
/*	if (pmd->height == 0 && pmd->pressure == 1013.25 && pmd->RH == 50
		&& pmd->temperature == 0)	{
		CartesianToGeodetic(&cg, pcrdSite, 6378137.0000000000, 1/298.257223563);
		h = cg.height;		//取大地高作为正高计算
		RH	= 50 * exp(-0.0006396 * h);
		p	= 1013.25 * pow(1 - 0.0000226 * h, 5.225);
		t	= -0.0065 * h;
	}
	else	{
		h	= pmd->height;
		p	= pmd->pressure;
		RH	= pmd->RH;
		t	= pmd->temperature;
	}	

	//计算延迟量
	T	= t + 273.16;
	e	= RH * exp(-37.2465 + 0.213166 * T - 0.000256908 * T * T);
	hw	= 11000;
	hd	= 40136 + 148.72 * t;
	Kd	= 155.2 * 0.0000001 * p * (hd - h) / T;
	Kw	= 155.2 * 0.0000001 * 4810 * (hd - h) * e / (T * T);
	TropDelay	= Kd / sin(sqrt(ctp.elevation * ctp.elevation + 6.25))
		+ Kw / sin(sqrt(ctp.elevation * ctp.elevation + 2.25));
*/
	*TropDelay	= 2.47 / (sin(ctp.elevation) + 0.0121);

	return ctp.elevation;
}