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