/*
输入参数：
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;
}