Solar.h

Go to the documentation of this file.

#pragma once

#include <math.h>

#define PI 3.1415926f
// official is 90 at 50'
// civil is 96
// nautical is 102
// astronomical is 108

namespace Library
{
        class Solar
        {
                float lat = 0.0f, lng = 0.0f, ALTITUDE = 0.0f, ZENITH = 0.0f;

        public:
                Solar() {};
                Solar(float prmLatitude, float prmLongitude, float prmAltitudeFeet, float prmZenithDegrees)
                {
                        lat = prmLatitude;
                        lng = prmLongitude;
                        ALTITUDE = prmAltitudeFeet;
                        ZENITH = prmZenithDegrees;
                }

                float calculateSunrise(int year, int month, int day, int localOffset, int daylightSavings) const
                {
                        /*
                        localOffset will be <0 for western hemisphere and >0 for eastern hemisphere
                        daylightSavings should be 1 if it is in effect during the summer otherwise it should be 0
                        */
                        //1. first calculate the day of the year
                        float N1 = floor(275.0f * month / 9.0f);
                        float N2 = floor((month + 9.0f) / 12.0f);
                        float N3 = 1.0f + floor((year - 4.0f * floor(year / 4.0f) + 2.0f) / 3.0f);
                        float N = N1 - (N2 * N3) + day - 30.0f;

                        //2. convert the longitude to hour value and calculate an approximate time
                        float lngHour = lng / 15.0f;
                        float t = N + ((6.0f - lngHour) / 24.0f);   //if rising time is desired:

                        //3. calculate the Sun's mean anomaly   
                        float M = (0.9856f * t) - 3.289f;

                        //4. calculate the Sun's true longitude
                        float L = fmod(M + (1.916f * sin((PI / 180.0f)*M)) + (0.020f * sin(2.0f * (PI / 180.0f) * M)) + 282.634f, 360.0f);

                        //5a. calculate the Sun's right ascension      
                        float RA = fmod(180.0f / PI*atan(0.91764f * tan((PI / 180.0f)*L)), 360.0f);

                        //5b. right ascension value needs to be in the same quadrant as L   
                        float Lquadrant = floor(L / 90.0f) * 90.0f;
                        float RAquadrant = floor(RA / 90.0f) * 90.0f;
                        RA = RA + (Lquadrant - RAquadrant);

                        //5c. right ascension value needs to be converted into hours   
                        RA = RA / 15.0f;

                        //6. calculate the Sun's declination
                        float sinDec = 0.39782f * sin((PI / 180.0f)*L);
                        float cosDec = cos(asin(sinDec));

                        //7a. calculate the Sun's local hour angle
                        float altAdj = 1.15f * sqrtf(ALTITUDE) / 60.0f;
                        float cosH = (cos((PI / 180.0f)*(ZENITH + altAdj)) - (sinDec * sin((PI / 180.0f)*lat))) / (cosDec * cos((PI / 180.0f)*lat));

                        assert(cosH <= 1.0f && cosH >= -1.0f);

                        //7b. finish calculating H and convert into hours
                        float H = 360.0f - (180.0f / PI)*acos(cosH);   //   if if rising time is desired:
                        H = H / 15.0f;

                        //8. calculate local mean time of rising/setting      
                        float T = H + RA - (0.06571f * t) - 6.622f;

                        //9. adjust back to UTC
                        float UT = T - lngHour;

                        //10. convert UT value to local time zone of latitude/longitude
                        float LT = UT + localOffset + daylightSavings;
                        if (LT < 0.0f) LT += 24.0f;
                        if (LT >= 24.0f) LT -= 24.0f;
                        return LT;
                }

                float calculateSunset(int year, int month, int day, int localOffset, int daylightSavings) const
                {
                        /*
                        localOffset will be <0 for western hemisphere and >0 for eastern hemisphere
                        daylightSavings should be 1 if it is in effect during the summer otherwise it should be 0
                        */
                        //1. first calculate the day of the year
                        float N1 = floor(275.0f * month / 9.0f);
                        float N2 = floor((month + 9.0f) / 12.0f);
                        float N3 = (1.0f + floor((year - 4.0f * floor(year / 4.0f) + 2.0f) / 3.0f));
                        float N = N1 - (N2 * N3) + day - 30.0f;

                        //2. convert the longitude to hour value and calculate an approximate time
                        float lngHour = lng / 15.0f;
                        float t = N + ((18 - lngHour) / 24.0f);   //if setting time is desired:

                        //3. calculate the Sun's mean anomaly   
                        float M = (0.9856f * t) - 3.289f;

                        //4. calculate the Sun's true longitude
                        float L = fmod(M + (1.916f * sin((PI / 180.0f)*M)) + (0.020f * sin(2.0f * (PI / 180.0f) * M)) + 282.634f, 360.0f);

                        //5a. calculate the Sun's right ascension      
                        float RA = fmod(180.0f / PI*atan(0.91764f * tan((PI / 180.0f)*L)), 360.0f);

                        //5b. right ascension value needs to be in the same quadrant as L   
                        float Lquadrant = floor(L / 90.0f) * 90.0f;
                        float RAquadrant = floor(RA / 90.0f) * 90.0f;
                        RA = RA + (Lquadrant - RAquadrant);

                        //5c. right ascension value needs to be converted into hours   
                        RA = RA / 15.0f;

                        //6. calculate the Sun's declination
                        float sinDec = 0.39782f * sin(PI / 180.0f*L);
                        float cosDec = cos(asin(sinDec));

                        //7a. calculate the Sun's local hour angle
                        float altAdj = 1.15f * sqrtf(ALTITUDE) / 60.0f;
                        float cosH = (cos(PI / 180.0f*(ZENITH + altAdj)) - (sinDec * sin(PI / 180.0f*lat))) / (cosDec * cos(PI / 180.0f*lat));

                        assert(cosH <= 1.0f && cosH >= -1.0f);

                        //7b. finish calculating H and convert into hours
                        float H = (180.0f / PI)*acos(cosH);
                        H = H / 15.0f;

                        //8. calculate local mean time of rising/setting      
                        float T = H + RA - (0.06571f * t) - 6.622f;

                        //9. adjust back to UTC
                        float UT = T - lngHour;

                        //10. convert UT value to local time zone of latitude/longitude
                        float LT = UT + localOffset + daylightSavings;
                        if (LT < 0.0f) LT += 24.0f;
                        if (LT >= 24.0f) LT -= 24.0f;
                        return LT + 12.0f; // PM
                }
        };
}