#include "math.h"
int GetDayOfYear(int Year,int Month,int Day)
{
//功能：给定年月日得到此日在一年中的天数，1月1日为第1天，分闰年和平年
//作者：汪帮主 2009-6-6 wzj08@mails.jlu.edu.cn
int DayOfYear;
int FebDays;
if ((Year%4==0 && Year%100!=0) || (Year%400==0))
{
FebDays = 29;
}
else
{
FebDays = 28;
}
DayOfYear = Day;
if (Month > 1) DayOfYear = DayOfYear + 31;
if (Month > 2) DayOfYear = DayOfYear + FebDays;
if (Month > 3) DayOfYear = DayOfYear + 31 ;
if (Month > 4) DayOfYear = DayOfYear + 30 ;
if (Month > 5) DayOfYear = DayOfYear + 31 ;
if (Month > 6) DayOfYear = DayOfYear + 30 ;
if (Month > 7) DayOfYear = DayOfYear + 31 ;
if (Month > 8) DayOfYear = DayOfYear + 31 ;
if (Month > 9) DayOfYear = DayOfYear + 30 ;
if (Month > 10) DayOfYear = DayOfYear + 31 ;
if (Month > 11) DayOfYear = DayOfYear + 30 ;
return DayOfYear;
}
void GetSolarAngle(int year,int month,int day,int hour,int min,double sec,double lat,double lon,int TimeZone,int N ,double& AzimuthAngle,double& ZenithAngle)
{
//功能：计算太阳天顶角和太阳方位角
//作者：汪帮主 2009-6-6 wzj08@mails.jlu.edu.cn
double N0, sitar, ED, dLon, Et, Ho, viewAng;
const double PI = 3.14159265;
double dTimeAngle, gtdt, latitudeArc, latitude, HeightAngleArc;
double AzimuthAngleArc, HANoon, dTA, TimeNoon, CosAzimuthAngle, HeightAngle;
N0 = 79.6764 + 0.2422*(year-1985) - floor((year-1985)/4.0);
sitar = 2*PI*(N-N0)/365.2422;
ED = 0.3723 + 23.2567*sin(sitar) + 0.1149*sin(2*sitar) - 0.1712*sin(3*sitar)- 0.758*cos(sitar) + 0.3656*cos(2*sitar) + 0.0201*cos(3*sitar);
ED = ED*PI/180; //ED本身有符号，无需判断正负。
//ED = (m_northOrSouth==0)? ED:(-ED);
dLon = 0.0;
//地球上某一点与其所在时区中心的经度差
if (lon >= 0)
{
if (TimeZone == -1)
{
dLon = lon - (floor((lon*10-75)/150)+1)*15.0;
}
else
{
dLon = lon - TimeZone*15.0;
}
}
else
{ if (TimeZone == -1)
{
dLon = (floor((lon*10-75)/150)+1)*15.0- lon;
}
else
{
dLon = TimeZone*15.0- lon;
}
}
Et = 0.0028 - 1.9857*sin(sitar) + 9.9059*sin(2*sitar) - 7.0924*cos(sitar)- 0.6882*cos(2*sitar); //视差
gtdt = hour + min/60.0 + sec/3600.0 + dLon/15; //地方时
gtdt = gtdt + Et/60.0;
dTimeAngle = 15.0*(gtdt-12);
dTimeAngle = dTimeAngle*PI/180;
latitudeArc = lat*PI/180;
HeightAngleArc = asin(sin(latitudeArc)*sin(ED)+cos(latitudeArc)*cos(ED)*cos(dTimeAngle));
CosAzimuthAngle = (sin(HeightAngleArc)*sin(latitudeArc)-sin(ED))/cos(HeightAngleArc)/cos(latitudeArc);
AzimuthAngleArc = acos(CosAzimuthAngle);
HeightAngle = HeightAngleArc*180/PI;
AzimuthAngle = AzimuthAngleArc *180/PI;
if( dTimeAngle < 0 )
{
AzimuthAngle = 180 - AzimuthAngle;
}
else
{
AzimuthAngle = 180 + AzimuthAngle;
}
ZenithAngle = 90-HeightAngle;
}