egnos.c

使用欧洲航空航天局推荐的算法估计电波在大气对流层中的传播延迟

/**
	egnos.c
	欧洲航空局推荐的电波在大气传播延迟的估计算法

	更新日志：
	20080903:1.将EGNOS模型算法整理形成egnos.c
 */
#include "comnavpp.h"
#include "mathspt.h"

/*
	Egnos_TropDelay
	用EGNOS模型计算对流层延迟
	《中国地区GPS中性大气天顶延迟研究及应用》，曲伟箐，中国科学院上海天文台
@para	elev,	卫星观测仰角(°)
@para	doy,	年积日(天)
@para	prcl,	测站大地坐标(rad,rad,m)
@return	F4型，对流层延迟估计量(m)
 */
F4 Egnos_TropDelay(	register F4 elev,
					register I4 doy,
					register COORDBLH * prcl)
{
	static const F4 __egnos_avrg[5][5] = {
	//	P0			T0	e0		beta0	lamda0	 latitude
		{1013.25,299.65,26.31,6.30e-3,2.77},	//15
		{1017.25,294.15,21.79,6.05e-3,3.15},	//30
		{1015.75,283.15,11.66,5.58e-3,2.57},	//45
		{1011.75,272.15,06.78,5.39e-3,1.81},	//60
		{1013.00,263.65,04.11,4.53e-3,1.55}		//75
	};
	static const F4 __egnos_dpara[5][5] = {
	//	  P0   	T0	  e0  beta0	lamda0	 latitude
		{-0.00,00.00,0.00,0.00e-3,0.00},	//15
		{-3.75,07.00,8.85,0.25e-3,0.33},	//30
		{-2.25,11.00,7.24,0.32e-3,0.46},	//45
		{-1.75,15.00,5.36,0.81e-3,0.74},	//60
		{-0.50,14.50,3.39,0.62e-3,0.30}		//75
	};
#define __egnos_p_i		0
#define __egnos_t_i		1
#define __egnos_e_i		2
#define __egnos_b_i		3
#define __egnos_l_i		4
#define __egnos_rd		287.054
#define __egnos_gm		9.784
#define __egnos_k1		77.604
#define __egnos_g		9.80665
#define __egnos_k2		3.82e5
	register I4 ilatbot,ilattop;
	register const F4 * egnos_ptop,* egnos_pbot;
	register F8 botscale,topscale,dayscale,dry,wet,beta,lambda,t,p,e,tmp,tmp1;
	register F8 lat = prcl->B*(180/PI);
	if(lat < 0) doy -= 211;
	else doy -= 28;
	dayscale = cos(2*PI*doy*(1.0/365.25));
	lat = fabs(lat);
	ilatbot = (I4)((lat + 00)*(1.0/15));
	ilattop = (I4)((lat + 15)*(1.0/15));
	botscale = (ilattop*15 - lat)*(1.0/15);
	topscale = (lat - ilatbot*15)*(1.0/15);
	if(ilatbot < 0) ilatbot = 0;
	else if(ilatbot > 4) ilatbot = 4;
	if(ilattop < 0) ilattop = 0;
	else if(ilattop > 4) ilattop = 4;
	//气象参数年平均值
	egnos_ptop = __egnos_avrg[ilattop];
	egnos_pbot = __egnos_avrg[ilatbot];
	beta = egnos_ptop[__egnos_b_i]*topscale + egnos_pbot[__egnos_b_i]*botscale;
	lambda = egnos_ptop[__egnos_l_i]*topscale + egnos_pbot[__egnos_l_i]*botscale;
	t = egnos_ptop[__egnos_t_i]*topscale + egnos_pbot[__egnos_t_i]*botscale;
	p = egnos_ptop[__egnos_p_i]*topscale + egnos_pbot[__egnos_p_i]*botscale;
	e = egnos_ptop[__egnos_e_i]*topscale + egnos_pbot[__egnos_e_i]*botscale;
	//加上气象参数的季节变化拟合
	egnos_ptop = __egnos_dpara[ilattop];
	egnos_pbot = __egnos_dpara[ilatbot];
	beta -= (egnos_ptop[__egnos_b_i]*topscale + egnos_pbot[__egnos_b_i]*botscale)*dayscale;
	lambda -= (egnos_ptop[__egnos_l_i]*topscale + egnos_pbot[__egnos_l_i]*botscale)*dayscale;
	t -= (egnos_ptop[__egnos_t_i]*topscale + egnos_pbot[__egnos_t_i]*botscale)*dayscale;
	p -= (egnos_ptop[__egnos_p_i]*topscale + egnos_pbot[__egnos_p_i]*botscale)*dayscale;
	e -= (egnos_ptop[__egnos_e_i]*topscale + egnos_pbot[__egnos_e_i]*botscale)*dayscale;
	t = recip(t);			//公式中仅用到t的倒数，所以在这里先求倒
	lambda += 1;			//公式中仅用到lambda+1，所以在这里直接先加1
	//zdry，平均海平面的干天顶延迟
	dry = p*(1e-6*__egnos_k1*__egnos_rd/__egnos_gm);
	//zwet，平均海平面的湿天顶延迟
	wet = e*t*(1e-6*__egnos_k2*__egnos_rd)
		* recip(__egnos_gm*lambda - beta*__egnos_rd);
	//ddry，接收机高度处的干天顶延迟
	tmp = 1.0 - beta*prcl->H*t;
	tmp1 = beta*(__egnos_g/__egnos_rd);
	dry *= pow(tmp,tmp1);
	//dwet，接收机高度处的湿天顶延迟
	wet *= pow(tmp,tmp1*lambda-1);

	tmp = elev*elev;
	return dry*recip(sin(sqrt(tmp + 6.25)*(PI/180)))
		+  wet*recip(sin(sqrt(tmp + 2.25)*(PI/180)));
}