[c++] 计算太阳高度角

/*
输入参数：                            
Longitude - 经度（单位"度"）
Latitude  - 纬度（单位"度"）
Year      - 年
Month     - 月
Day       - 日
Hour      - 时
Minute    - 分
Second    - 秒

输出参数：
Height    - 太阳高度角   （单位"弧度"）
Direction - 太阳方位角   （单位"弧度"）
*/
void SunPositionNR(double Longitude, double Latitude, int year, int month, int day, int hour, int minute, int second, double& height, double& direction)
{
	double A = year/4;   
	double B = A-floor(A);   
	double C = 32.8;   
	if(month <= 2)   
		C = 30.6;   
	if(B==0 && month>2)   
		C = 31.8;   
	double G = floor(30.6*month-C+0.5)+day;   
	double L = Longitude/15.0;                    //经度修正    
	double H = hour-8 + minute/60.0;              //时刻修正    
	double N = G+(H-L)/24.0;                      //计算积日 
	double theta, NO;     //theta为日角
	NO = 79.6764 + 0.2422*(year-1985) - floor((year-1985)/4.0);
	theta = 2*PI*(N-NO)/365.2422;

	double lat;//赤纬角
	lat = 0.3723+23.2567*sin(theta)+0.1149*sin(2*theta)-0.1712*sin(3*theta)-0.758*cos(theta)+0.3656*cos(2*theta)+0.0201*cos(3*theta);
	lat = lat*3.14159265358979/180.0;

	double Et;//时差（以分钟为单位）
	Et = 0.0028-1.9857*sin(theta)+9.9059*sin(2*theta)-7.0924*cos(theta)-0.6882*cos(2*theta);

	double timeAngle;//时角
	timeAngle = ( (hour-12) + (minute - (120.0-Longitude) * 4.0 + Et)/60.0 ) * 15.0;
	timeAngle = timeAngle * 3.14159265358979 / 180.0;

	// 将纬度由角度转换为弧度
	Latitude = Latitude * 3.14159265358979 / 180.0;

	// 太阳高度角
	height = sin(Latitude)*sin(lat)+cos(Latitude)*cos(lat)*cos(timeAngle);
	height = asin(height);

	// 太阳方位角
	direction = (sin(height)*sin(Latitude)-sin(lat))/(cos(height)*cos(Latitude));
	direction = acos(direction);

	if (timeAngle<0)  //(以正南方向为0度，顺时针为正)
		direction *= -1.0;
}