#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;
}