太阳方位角

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

}