Sunpath.cpp


1
/*
2
 * Sunpath.cpp
3
 *
4
 *  Created on: 25.11.2013
5
 *      Author: luoc
6
 */
7
8
#include "Sunpath.h"
9
10
Sunpath::Sunpath(float longitude, float latitude) 
11
{
12
  this->longitude = longitude;
13
  this->latitude = latitude;
14
}
15
16
// private functions
17
18
boolean Sunpath::leapYear(time_t t)
19
{
20
  int y = year(t);
21
  if ( y % 4 == 0 && y % 100 != 0 || y % 400 == 0 ) 
22
  {
23
    return true;
24
  }
25
  else
26
  {
27
    return false;
28
  }
29
}
30
31
// public functions
32
33
double Sunpath::getAzimuth(time_t currentTimeUTC)
34
{
35
  int minutesOfDay = hour(getTimeSolar(currentTimeUTC)) * 60 + minute(getTimeSolar(currentTimeUTC));
36
  
37
  if(minutesOfDay < 0) // make shure that the angle is between 0° and 360°
38
  {
39
    minutesOfDay += 1440;
40
  }
41
  
42
  if(minutesOfDay >= 1440) // make shure that the angle is between 0° and 360°
43
  {
44
    minutesOfDay -= 1440;
45
  }
46
  
47
  return minutesOfDay * 0.25;
48
}
49
50
double Sunpath::getElevationAngle(time_t currentTimeUTC)
51
{
52
  time_t currentTimeSolar = getTimeSolar(currentTimeUTC);
53
  double h; // hours since 12:00
54
  int d; // days since 21. June
55
  
56
  switch(month(currentTimeSolar))
57
  {
58
    case 1: 
59
    d = 193 + day(currentTimeSolar); 
60
    break;
61
    
62
    case 2: 
63
    d = 224 + day(currentTimeSolar); 
64
    break;
65
    
66
    case 3: 
67
    d = 252 + day(currentTimeSolar);
68
    if(leapYear(currentTimeSolar)){ d++; }
69
    break;
70
    
71
    case 4: 
72
    d = 283 + day(currentTimeSolar); 
73
    if(leapYear(currentTimeSolar)){ d++; }
74
    break;
75
    
76
    case 5: 
77
    d = 313 + day(currentTimeSolar);
78
    if(leapYear(currentTimeSolar)){ d++; } 
79
    break;
80
    
81
    case 6: 
82
    if(day(currentTimeSolar) < 21)
83
    {
84
      d = 344 + day(currentTimeSolar);
85
      if(leapYear(currentTimeSolar)){ d++; }
86
    }else
87
    {
88
      d = day(currentTimeSolar) - 21;
89
    } 
90
    break;
91
    
92
    case 7: 
93
    d = 9 + day(currentTimeSolar); 
94
    break;
95
    
96
    case 8: 
97
    d = 40 + day(currentTimeSolar); 
98
    break;
99
    
100
    case 9: 
101
    d = 71 + day(currentTimeSolar); 
102
    break;
103
    
104
    case 10: 
105
    d = 101 + day(currentTimeSolar); 
106
    break;
107
    
108
    case 11: 
109
    d = 132 + day(currentTimeSolar); 
110
    break;
111
    
112
    case 12: 
113
    d = 162 + day(currentTimeSolar); 
114
    break;
115
    
116
    default: d = 0; break;
117
  }
118
  
119
  if(isAM(currentTimeSolar))
120
  {
121
    h = hour(currentTimeSolar) + 12;
122
    if(month(currentTimeSolar) == 6 && day(currentTimeSolar) == 21)
123
    {
124
      if(leapYear(currentTimeSolar))
125
      {
126
        d = 365;
127
      }else
128
      {
129
        d = 364;
130
      }
131
    }else
132
    {
133
      d--;
134
    }
135
  }else
136
  {
137
    h = hour(currentTimeSolar) - 12;
138
  }
139
  
140
   h += double(minute(currentTimeSolar)) / 60.0;
141
142
  // I think it's not necessary to explain :P
143
  return asin(sin(AXIAL_TILT*(PI/180.0))*sin(latitude*(PI/180.0))*cos((2*PI*double(d))/(366))+cos(latitude*(PI/180.0))*sqrt(1.0-square(sin(AXIAL_TILT*(PI/180.0))*cos(sqrt((2*PI*double(d))/(366)))))*cos((2*PI*h)/(24))) * (180/PI); 
144
}
145
146
time_t Sunpath::getTimeSolar(time_t currentTime)
147
{
148
  return long(currentTime + longitude * 4 * 60);
149
}
150
151
Sunpath::~Sunpath() {
152
  // TODO Auto-generated destructor stub
153
}