/**
	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)));
}