Full-sized photovoltaic system

1#include <SPI.h>
2#include <PortentaEthernet.h>
3#include <Ethernet.h>
4#include <EthernetUdp.h>
5#include <RTClib.h>
6#include <Arduino.h>
7#include "RealTimeClock.h"
8#include "mbed.h"
9
10#define PI 3.14159265358979323846
11
12// Network related definitions
13IPAddress arduinoIp(192, 168, 3, 130);
14unsigned int localUdpPort = 8888;
15const char ntpServer[] = "europe.pool.ntp.org";
16const int NTP_PACKET_SIZE = 48; // NTP time stamp is in the first 48 bytes of the message
17byte packetBuffer[NTP_PACKET_SIZE]; //buffer to hold incoming and outgoing packets
18EthernetUDP Udp; // A UDP instance to let us send and receive packets over UDP
19const unsigned long ethernetTimeout = 2000;
20const unsigned long ethernetResponseTimeout = 2000;
21
22// Time constants for calculations
23const time_t secondsInSeventyYears = 2208988800;
24const time_t secondsInADay = 86400;
25const double julianDate1970 = 2440587.5;
26
27// Location of PV plant
28const double latitude = 10.0; // + to N
29const double longitude = 10.0; // + to E
30
31// All measurements in mm
32const int elevationMotorLengthTop = 1067;
33const int elevationMotorLengthBeam = 1868;
34const int azimuthMotorLengthTop = 407;
35const int azimuthMotorLengthBeam = 1164;
36const int motorLength = 929;
37const int motorMaxLength = 457;
38const int motorTimeout = 60000;
39
40// Construction limits
41const double angleOffsetSouth = 21.2;
42const double minAngleSouth = 40.4;
43const double maxAngleSouth = 68.4;
44
45const double angleOffsetEast = 9.9;
46const double minAngleEast = 55.5;
47const double maxAngleEast = 125.0;
48
49// Measurement parameters
50const double arduinoReferenceVoltage = 10.0;
51const double potiReferenceVoltage = 10.2;
52const uint8_t analogResolution = 10;
53const int potiOffset = 310;
54const int potiResistanceSum = 9980;
55const double ohmPerMillimeter = 9.84; // Resolution of the motor potentiometer in Ohm/mm
56const double correctionFactor = 1.34;
57
58// Pin mappings
59int elevationMotorExtensionPin = D0;
60int elevationMotorPullPin = D1;
61int azimuthMotorExtensionPin = D2;
62int azimuthMotorPullPin = D3;
63
64int elevationMotorExtensionLedPin = LED_D0;
65int elevationMotorPullLedPin = LED_D1;
66int azimuthMotorExtensionLedPin = LED_D2;
67int azimuthMotorPullLedPin = LED_D3;
68
69int elevationMotorPotiPin = A0;
70int azimuthMotorPotiPin = A1;
71int stormPin = A7;
72
73// Angle array
74#define INDEX_AZIMUTH      0
75#define INDEX_EAST         0
76#define INDEX_ELEVATION    1
77#define INDEX_SOUTH        1
78#define ANGLE_ARR_SIZE     2
79
80double lastSetAngles[ANGLE_ARR_SIZE] = { 0.0 };
81double defaultAngles[ANGLE_ARR_SIZE] = { 90.0, maxAngleSouth };
82#define ANGLES_DELTA               3.0
83#define AZIMUTH_OFFSET_SOUTH     180.0
84#define SWING_ANGLE_OFFSET_SOUTH  90.0
85
86// Time array
87#define INDEX_SUNRISE      0
88#define INDEX_SUNSET       1
89#define INDEX_TIMESTAMP    2
90#define TIME_ARR_SIZE      3
91
92void setup()
93{
94  Serial.begin(9600); 
95
96  // Initialize relay outputs
97  pinMode(elevationMotorExtensionPin, OUTPUT);
98  pinMode(elevationMotorPullPin, OUTPUT);
99  pinMode(azimuthMotorExtensionPin, OUTPUT);
100  pinMode(azimuthMotorPullPin, OUTPUT);
101
102  // Initialize LED outputs
103  pinMode(elevationMotorExtensionLedPin, OUTPUT);
104  pinMode(elevationMotorPullLedPin, OUTPUT);
105  pinMode(azimuthMotorExtensionLedPin, OUTPUT);
106  pinMode(azimuthMotorPullLedPin, OUTPUT);
107
108  // Initialize potentiometer inputs
109  pinMode(elevationMotorPotiPin, INPUT);
110  pinMode(azimuthMotorPotiPin, INPUT);
111  pinMode(stormPin, INPUT);
112
113  analogReadResolution(analogResolution);
114
115  delay(10000);
116}
117
118void loop()
119{
120  double sunAngleResults[ANGLE_ARR_SIZE] = { 0.0, 0.0 };
121  double motorAngleResults[ANGLE_ARR_SIZE] = { 0.0, 0.0 };
122  time_t timeResults[TIME_ARR_SIZE] = { 0 };
123
124  astronomicalCalculations(sunAngleResults, timeResults);
125  calculateMotorAngles(sunAngleResults, motorAngleResults);
126  printResults(sunAngleResults, motorAngleResults, timeResults);
127  followTheSun(motorAngleResults, timeResults);
128  ntpSynchronizeTime();
129  wait();
130}
131
132void calculateMotorAngles(double sunAngleResults[], double motorAngleResults[])
133{
134  double azimuthCorrectedDeg = sunAngleResults[INDEX_EAST] - AZIMUTH_OFFSET_SOUTH;
135  double azimuthCorrectedRad = toRad(azimuthCorrectedDeg);
136  double elevationRad = toRad(sunAngleResults[INDEX_SOUTH]);
137
138  double a = 1000.0;
139  double b = (a / cos(elevationRad));
140  double c = (a * tan(elevationRad));
141  double d = (a * sin(azimuthCorrectedRad));
142  double e = (a * cos(azimuthCorrectedRad));
143  double inclinationRad = atan(c / e);
144  double inclinationDeg = toDeg(inclinationRad);
145
146  if((azimuthCorrectedDeg < -89) || (azimuthCorrectedDeg > 89))
147  {
148    inclinationDeg = maxAngleSouth;
149  }
150
151  if(inclinationDeg < minAngleSouth)
152  {
153    inclinationDeg = minAngleSouth;
154  }
155  else if(inclinationDeg > maxAngleSouth)
156  {
157    inclinationDeg = maxAngleSouth;
158  }
159
160  double f = (e / cos(inclinationRad));
161  double swingAngleDeg = toDeg(atan(d / f));
162  swingAngleDeg += SWING_ANGLE_OFFSET_SOUTH;
163
164  if(azimuthCorrectedDeg < -89)
165  {
166    swingAngleDeg = minAngleEast;
167  }
168  else if(azimuthCorrectedDeg > 89)
169  {
170    swingAngleDeg = maxAngleEast;
171  }
172
173  if(swingAngleDeg < minAngleEast)
174  {
175    swingAngleDeg = minAngleEast;
176  }
177  else if(swingAngleDeg > maxAngleEast)
178  {
179    swingAngleDeg = maxAngleEast;
180  }
181
182  motorAngleResults[INDEX_SOUTH] = inclinationDeg;
183  motorAngleResults[INDEX_EAST] = swingAngleDeg;
184}
185
186void wait()
187{
188  Serial.println("Wait");
189
190  int lightningLed = 0;
191
192  // Wait 5 minutes
193  for (int i = 0; i < 3000; i++)
194  {
195    if (i % 10 == 0)
196    {
197      switch(lightningLed)
198      {
199        case 0:
200          digitalWrite(LED_D3, LOW);
201          digitalWrite(LED_D0, HIGH);
202          lightningLed++;
203          break;
204        case 1:
205          digitalWrite(LED_D0, LOW);
206          digitalWrite(LED_D1, HIGH);
207          lightningLed++;
208          break;
209        case 2:
210          digitalWrite(LED_D1, LOW);
211          digitalWrite(LED_D2, HIGH);
212          lightningLed++;
213          break;
214        case 3:
215          digitalWrite(LED_D2, LOW);
216          digitalWrite(LED_D3, HIGH);
217          lightningLed = 0;
218          break;
219      }
220    }
221    
222    detectStorm();
223    delay(100);
224  }
225
226  digitalWrite(LED_D0, LOW);
227  digitalWrite(LED_D1, LOW);
228  digitalWrite(LED_D2, LOW);
229  digitalWrite(LED_D3, LOW);
230}
231
232void detectStorm()
233{
234  if (digitalRead(stormPin) == HIGH)
235  {
236    Serial.println("Storm detected");
237
238    digitalWrite(LED_D0, LOW);
239    digitalWrite(LED_D1, LOW);
240    digitalWrite(LED_D2, LOW);
241    digitalWrite(LED_D3, LOW);
242
243    setMotorsToDefault();
244    
245    // Wait 30 minutes
246    for (int i = 0; i < 18000; i++)
247    {
248      delay(100);
249
250      if(i % 10 == 0)
251      {
252        digitalWrite(LED_D0, HIGH);
253        digitalWrite(LED_D1, HIGH);
254        digitalWrite(LED_D2, HIGH);
255        digitalWrite(LED_D3, HIGH);
256      }
257      else if(i % 5 == 0)
258      {
259        digitalWrite(LED_D0, LOW);
260        digitalWrite(LED_D1, LOW);
261        digitalWrite(LED_D2, LOW);
262        digitalWrite(LED_D3, LOW);
263      }
264    }
265
266    digitalWrite(LED_D0, LOW);
267    digitalWrite(LED_D1, LOW);
268    digitalWrite(LED_D2, LOW);
269    digitalWrite(LED_D3, LOW);
270  }
271}
272
273void setMotorsToDefault()
274{
275  int azimuthMotorPosition = calcAzimuthMotorPosition(defaultAngles[INDEX_EAST]);
276  setMotor(azimuthMotorPotiPin, azimuthMotorExtensionPin, azimuthMotorPullPin, azimuthMotorExtensionLedPin, azimuthMotorPullLedPin, azimuthMotorPosition);
277  lastSetAngles[INDEX_EAST] = defaultAngles[INDEX_EAST];
278
279  int elevationMotorPosition = calcElevationMotorPosition(defaultAngles[INDEX_SOUTH]);
280  setMotor(elevationMotorPotiPin, elevationMotorExtensionPin, elevationMotorPullPin, elevationMotorExtensionLedPin, elevationMotorPullLedPin, elevationMotorPosition);
281  lastSetAngles[INDEX_SOUTH] = defaultAngles[INDEX_SOUTH];
282}
283
284void followTheSun(double motorAngleResults[], time_t timeResults[])
285{
286  double deltaEast = 0.0;
287  double deltaSouth = 0.0;
288
289  if (isInDayHours(timeResults))
290  {
291    deltaEast = abs(motorAngleResults[INDEX_EAST] - lastSetAngles[INDEX_EAST]);
292    deltaSouth = abs(motorAngleResults[INDEX_SOUTH] - lastSetAngles[INDEX_SOUTH]);
293  }
294  else
295  {
296    Serial.println("Sun not shining - default position");
297    deltaEast = abs(defaultAngles[INDEX_EAST] - lastSetAngles[INDEX_EAST]);
298    motorAngleResults[INDEX_EAST] = defaultAngles[INDEX_EAST];
299    deltaSouth = abs(defaultAngles[INDEX_SOUTH] - lastSetAngles[INDEX_SOUTH]);
300    motorAngleResults[INDEX_SOUTH] = defaultAngles[INDEX_SOUTH];
301  }
302
303  if (deltaEast > ANGLES_DELTA)
304  {
305    Serial.print("Setting east motor - delta is ");
306    Serial.print(deltaEast);
307    Serial.println("°");
308
309    int azimuthMotorPosition = calcAzimuthMotorPosition(motorAngleResults[INDEX_EAST]);
310    setMotor(azimuthMotorPotiPin, azimuthMotorExtensionPin, azimuthMotorPullPin, azimuthMotorExtensionLedPin, azimuthMotorPullLedPin, azimuthMotorPosition);
311    Serial.print("Set east motor to ");
312    Serial.print(azimuthMotorPosition);
313    Serial.println(" mm");
314    lastSetAngles[INDEX_EAST] = motorAngleResults[INDEX_EAST];
315  }
316
317  if (deltaSouth > ANGLES_DELTA)
318  {
319    Serial.print("Setting south motor - delta is ");
320    Serial.print(deltaSouth);
321    Serial.println("°");
322
323    int elevationMotorPosition = calcElevationMotorPosition(motorAngleResults[INDEX_SOUTH]);
324    setMotor(elevationMotorPotiPin, elevationMotorExtensionPin, elevationMotorPullPin, elevationMotorExtensionLedPin, elevationMotorPullLedPin, elevationMotorPosition);
325    Serial.print("Set south motor to ");
326    Serial.print(elevationMotorPosition);
327    Serial.println(" mm");
328    lastSetAngles[INDEX_SOUTH] = motorAngleResults[INDEX_SOUTH];
329  }
330}
331
332void setMotor(int motorPotiPin, int motorExtensionPin, int motorPullPin, int extensionPinLed, int pullPinLed, int targetMotorPosition)
333{
334  int currentMotorPosition = getMotorPosition(motorPotiPin);
335
336  if (targetMotorPosition > currentMotorPosition)
337  {
338    digitalWrite(extensionPinLed, HIGH);
339    digitalWrite(motorExtensionPin, HIGH);
340
341    for (int i = 0; (targetMotorPosition > currentMotorPosition) && (i < (motorTimeout / 10)); i++)
342    {
343      currentMotorPosition = getMotorPosition(motorPotiPin);
344      delay(10);
345    }
346
347    digitalWrite(motorExtensionPin, LOW);
348    digitalWrite(extensionPinLed, LOW);
349
350    return;
351  }
352
353  if (targetMotorPosition < currentMotorPosition)
354  {
355    digitalWrite(pullPinLed, HIGH);
356    digitalWrite(motorPullPin, HIGH);
357
358    for (int i = 0; (targetMotorPosition < currentMotorPosition) && (i < (motorTimeout / 10)); i++)
359    {
360      currentMotorPosition = getMotorPosition(motorPotiPin);
361      delay(10);
362    }
363
364    digitalWrite(motorPullPin, LOW);
365    digitalWrite(pullPinLed, LOW);
366  }
367}
368
369int getMotorPosition(int potiPin)
370{
371  int rawValue = analogRead(potiPin);
372  int resolution = pow(2, analogResolution);
373  double inputVoltage = rawValue * (arduinoReferenceVoltage / resolution);
374  double potiResistance = (inputVoltage / potiReferenceVoltage) * potiResistanceSum;
375  potiResistance *= correctionFactor;
376  int motorPosition = (potiResistance - potiOffset) / ohmPerMillimeter;
377
378  if (motorPosition < 0)
379  {
380    motorPosition = 0;
381  }
382  
383  if (motorPosition > motorMaxLength)
384  {
385    motorPosition = motorMaxLength;
386  }
387
388  return motorPosition;
389}
390
391int calcAzimuthMotorPosition(double azimuthAngle)
392{
393  double a = azimuthMotorLengthTop;
394  double b = azimuthMotorLengthBeam;
395  double realAngle = toRad(azimuthAngle - angleOffsetEast);
396  double c = sqrt(pow(a, 2) + pow(b, 2) - 2*a*b*cos(realAngle));
397  int motorExtension = c - motorLength;
398
399  if (motorExtension < 0)
400  {
401    motorExtension = 0;
402  }
403
404  return motorExtension;
405}
406
407int calcElevationMotorPosition(double elevationAngle)
408{
409  double a = elevationMotorLengthTop;
410  double b = elevationMotorLengthBeam;
411  double realAngle = toRad(elevationAngle - angleOffsetSouth);
412  double c = sqrt(pow(a, 2) + pow(b, 2) - 2*a*b*cos(realAngle));
413  int motorExtension = c - motorLength;
414
415  if (motorExtension < 0)
416  {
417    motorExtension = 0;
418  }
419
420  return motorExtension;
421}
422
423void printResults(double sunAngleResults[], double motorAngleResults[], time_t timeResults[])
424{
425  Serial.println();
426  Serial.println("=== ASTRONOMICAL DATA ===");
427  Serial.println();
428
429  Serial.print("Current time and date: ");
430  DateTime currentTime = DateTime(timeResults[INDEX_TIMESTAMP]);
431  Serial.println(currentTime.timestamp());
432
433  Serial.print("Azimuth: ");
434  Serial.print(sunAngleResults[INDEX_EAST]);
435  Serial.println("°");
436
437  Serial.print("Elevation: ");
438  Serial.print(sunAngleResults[INDEX_SOUTH]);
439  Serial.println("°");
440
441  Serial.print("South motor: ");
442  Serial.print(motorAngleResults[INDEX_SOUTH]);
443  Serial.println("°");
444
445  Serial.print("East motor: ");
446  Serial.print(motorAngleResults[INDEX_EAST]);
447  Serial.println("°");
448
449  Serial.print("Sunrise: ");
450  displayTimeOfDay(timeResults[INDEX_SUNRISE], timeResults[INDEX_TIMESTAMP]);
451  Serial.println(" UTC");
452  
453  Serial.print("Sunset: ");
454  displayTimeOfDay(timeResults[INDEX_SUNSET], timeResults[INDEX_TIMESTAMP]);
455  Serial.println(" UTC");
456}
457
458void displayTimeOfDay(time_t secondsSinceMidnight, time_t anyTimestampOfTheDay)
459{
460  time_t midnight = getMidnightInUnixTime(anyTimestampOfTheDay);
461  time_t timeToDisplay = midnight + secondsSinceMidnight;
462  DateTime dateTime = DateTime(timeToDisplay);
463  Serial.print(dateTime.timestamp(DateTime::TIMESTAMP_TIME));
464}
465
466time_t getMidnightInUnixTime(time_t unixTimestamp)
467{
468  time_t secondsSinceMidnight = unixTimestamp % secondsInADay;
469  time_t midnight = unixTimestamp - secondsSinceMidnight;
470  return midnight;
471}
472
473bool isInDayHours(time_t timeResults[])
474{
475  time_t currentSecondsOfDay = timeResults[INDEX_TIMESTAMP] % secondsInADay;
476  return (currentSecondsOfDay >= timeResults[INDEX_SUNRISE] && currentSecondsOfDay <= timeResults[INDEX_SUNSET]);
477}
478
479void ntpSynchronizeTime()
480{
481  if (Ethernet.begin(nullptr, ethernetTimeout, ethernetResponseTimeout) == 0)
482  {
483    Serial.println("[WARNING]: No network connection available");
484    return;
485  }
486
487  Udp.begin(localUdpPort);
488  sendNTPpacket(ntpServer);
489  delay(1000);
490
491  if (Udp.parsePacket())
492  {
493    // Packet received
494    Udp.read(packetBuffer, NTP_PACKET_SIZE);
495
496    // Extract the 4 byte long timestamp from the buffer, starting at byte 40
497    unsigned long highWord = word(packetBuffer[40], packetBuffer[41]);
498    unsigned long lowWord = word(packetBuffer[42], packetBuffer[43]);
499    
500    // Combine the four bytes (two words)
501    unsigned long secsSince1900 = highWord << 16 | lowWord;
502    time_t secsSince1970 = secsSince1900 - secondsInSeventyYears;
503    
504    // Set RTC
505    set_time(secsSince1970);
506
507    Serial.print("Synchronized RTC to time: ");
508    Serial.print(DateTime(secsSince1970).timestamp());
509    Serial.println(" UTC");
510  }
511  else
512  {
513    Serial.println("[WARNING]: No UDP packet received. Check the network connection!");
514  }
515}
516
517double calcJulianDay(time_t secsSince1970)
518{
519  double daysSince1970 = secsSince1970 / secondsInADay;
520  return (daysSince1970 + julianDate1970);
521}
522
523double calcJulianCentury(time_t secsSince1970)
524{
525  double julianDay = calcJulianDay(secsSince1970);
526  return ((julianDay - 2451545) / 36525);
527}
528
529double toDeg(double val)
530{
531    return (val * 180) / PI;
532}
533
534double toRad(double val)
535{
536    return (val * PI) / 180;
537}
538
539void sendNTPpacket(const char * address) {
540  // Set all bytes in the buffer to 0
541  memset(packetBuffer, 0, NTP_PACKET_SIZE);
542
543  // Initialize values needed to form NTP request
544  packetBuffer[0] = 0b11100011;   // LI, Version, Mode
545  packetBuffer[1] = 0;            // Stratum, or type of clock
546  packetBuffer[2] = 6;            // Polling Interval
547  packetBuffer[3] = 0xEC;         // Peer Clock Precision
548  // 8 bytes of zero for Root Delay & Root Dispersion
549  packetBuffer[12]  = 49;
550  packetBuffer[13]  = 0x4E;
551  packetBuffer[14]  = 49;
552  packetBuffer[15]  = 52;
553
554  Udp.beginPacket(address, 123); // NTP requests use port 123
555  Udp.write(packetBuffer, NTP_PACKET_SIZE);
556  Udp.endPacket();
557}
558
559void astronomicalCalculations(double angleResults[], time_t timeResults[])
560{
561    time_t secsSince1970 = time(NULL);
562    
563    double julianDay = calcJulianDay(secsSince1970);
564    double julianCentury = calcJulianCentury(secsSince1970);
565    DateTime dateTimeInstance = DateTime(secsSince1970);
566    uint8_t hour = dateTimeInstance.hour();
567    uint8_t minute = dateTimeInstance.minute();
568
569    // Geom Mean (deg)
570    double geomMeanLongSun_Deg = fmod(280.46646 + julianCentury * (36000.76983 + julianCentury * 0.0003032), 360);
571    double geomMeanAnomSun_Deg = 357.52911 + julianCentury * (35999.05029 - 0.0001537 * julianCentury);
572
573    // Eccent Earth Orbit
574    double eccentEarthOrbit = 0.016708634 - julianCentury * (0.000042037 + 0.0000001267 * julianCentury);
575
576    // Sun Eq of Ctr
577    double sunEqOfCtr = sin(toRad(geomMeanAnomSun_Deg)) * (1.914602 - julianCentury * (0.004817 + 0.000014 * julianCentury)) + sin(toRad(2 * geomMeanAnomSun_Deg)) * (0.019993 - 0.000101 * julianCentury) + sin(toRad(3 * geomMeanAnomSun_Deg)) * 0.000289;
578
579    // Sun true long (deg)
580    double sunTrueLong_Deg = geomMeanLongSun_Deg + sunEqOfCtr;
581
582    // Sun true mean (deg)
583    double sunTrueAnom_Deg = geomMeanAnomSun_Deg + sunEqOfCtr;
584
585    // Sun Rad Vector (AUs)
586    double sunRadVector = (1.000001018 * (1 - eccentEarthOrbit * eccentEarthOrbit)) / (1 + eccentEarthOrbit * cos(toRad(sunTrueAnom_Deg)));
587
588    // Sun App Long (deg)
589    double sunAppLong_Deg = sunTrueLong_Deg - 0.00569 - 0.00478 * sin(toRad(125.04 - 1934.136 * julianCentury));
590
591    // Mean Obliq Ecliptic (deg)
592    double meanObliqEcliptic_Deg = 23 + (26 + ((21.448 - julianCentury * (46.815 + julianCentury * (0.00059 - julianCentury * 0.001813)))) / 60) / 60;
593
594    // Obliq Corr (deg)
595    double obliqCorr_Deg = meanObliqEcliptic_Deg + 0.00256 * cos(toRad(125.04 - 1934.136 * julianCentury));
596
597    // Sun Rt Ascen (deg)
598    double sunRtAscen_Deg = toDeg(atan2(cos(toRad(obliqCorr_Deg)) * sin(toRad(sunAppLong_Deg)), cos(toRad(sunAppLong_Deg))));
599
600    // Sun Declin (deg)
601    double sunDeclin_Deg = toDeg(asin(sin(toRad(obliqCorr_Deg)) * sin(toRad(sunAppLong_Deg))));
602
603    // var y
604    double y = tan(toRad(obliqCorr_Deg / 2)) * tan(toRad(obliqCorr_Deg / 2));
605
606    // Eq of time (min)
607    double eqOfTime_Min = 4 * toDeg(y * sin(2 * toRad(geomMeanLongSun_Deg)) - 2 * eccentEarthOrbit * sin(toRad(geomMeanAnomSun_Deg)) + 4 * eccentEarthOrbit * y * sin(toRad(geomMeanAnomSun_Deg)) * cos(2 * toRad(geomMeanLongSun_Deg)) - 0.5 * y * y * sin(4 * toRad(geomMeanLongSun_Deg)) - 1.25 * eccentEarthOrbit * eccentEarthOrbit * sin(2 * toRad(geomMeanAnomSun_Deg)));
608
609    // HA Sunrise (Deg)
610    double haSunrise_Deg = toDeg(acos(cos(toRad(90.833)) / (cos(toRad(latitude)) * cos(toRad(sunDeclin_Deg))) - tan(toRad(latitude)) * tan(toRad(sunDeclin_Deg))));
611
612    // Solar Noon (UTC)
613    double solarNoon = (720 - 4 * longitude - eqOfTime_Min) / 1440; // needs a conversion to hh:mm:ss
614
615    // Sunrise Time (UTC)
616    double sunriseTime = solarNoon - haSunrise_Deg * 4 / 1440; // needs a conversion to hh:mm:ss
617
618    // Sunset Time (UTC)
619    double sunsetTime = solarNoon + haSunrise_Deg * 4 / 1440; // needs a conversion to hh:mm:ss
620
621    // Sunlight Duration (min)
622    double sunlightDuration_Min = 8 * haSunrise_Deg;
623
624    // True Solar Time (min)
625    double trueSolarTime_Min = fmod(hour * 60 + minute + eqOfTime_Min + 4 * longitude, 1440);
626
627    // Hour Angle (Deg)
628    double hourAngle_Deg = ((trueSolarTime_Min / 4) < 0) ? (trueSolarTime_Min / 4 + 180) : (trueSolarTime_Min / 4 - 180);
629
630    // Solar Zenith Angle (Deg)
631    double solarZenithAngle_Deg = toDeg(acos(sin(toRad(latitude)) * sin(toRad(sunDeclin_Deg)) + cos(toRad(latitude)) * cos(toRad(sunDeclin_Deg)) * cos(toRad(hourAngle_Deg))));
632
633    // Solar Elevation Angle (Deg)
634    double solarElevationAngle_Deg = 90 - solarZenithAngle_Deg;
635
636    // Approx Atmospheric Refraction (Deg)
637    double approxAtmosphericRefraction_Deg;
638
639    if (solarElevationAngle_Deg > 85)
640    {
641        approxAtmosphericRefraction_Deg = 0;
642    }
643    else
644    {
645        if (solarElevationAngle_Deg > 5)
646        {
647            approxAtmosphericRefraction_Deg = 58.1 / tan(toRad(solarElevationAngle_Deg)) - 0.07 / pow(tan(toRad(solarElevationAngle_Deg)), 3) + 0.000086 / pow(tan(toRad(solarElevationAngle_Deg)), 5);
648        }
649        else
650        {
651            if (solarElevationAngle_Deg > -0.575)
652            {
653                approxAtmosphericRefraction_Deg = 1735 + solarElevationAngle_Deg * (-518.2 + solarElevationAngle_Deg * (103.4 + solarElevationAngle_Deg * (-12.79 + solarElevationAngle_Deg * 0.711)));
654            }
655            else
656            {
657                approxAtmosphericRefraction_Deg = -20.772 / tan(toRad(solarElevationAngle_Deg));
658            }
659        }
660    }
661
662    approxAtmosphericRefraction_Deg = approxAtmosphericRefraction_Deg / 3600;
663
664    // Solar Elevation corrected for atm refraction (Deg)
665    double solarElevationAngleCorrected_Deg = solarElevationAngle_Deg + approxAtmosphericRefraction_Deg;
666
667    // Solar Azimuth Angle (Deg cw from N)
668    double solarAzimuthAngle;
669
670    if (hourAngle_Deg > 0)
671    {
672        solarAzimuthAngle = fmod(toDeg(acos(((sin(toRad(latitude)) * cos(toRad(solarZenithAngle_Deg))) - sin(toRad(sunDeclin_Deg))) / (cos(toRad(latitude)) * sin(toRad(solarZenithAngle_Deg))))) + 180, 360);
673    }
674    else
675    {
676        solarAzimuthAngle = fmod(540 - toDeg(acos(((sin(toRad(latitude)) * cos(toRad(solarZenithAngle_Deg))) - sin(toRad(sunDeclin_Deg))) / (cos(toRad(latitude)) * sin(toRad(solarZenithAngle_Deg))))), 360);
677    }
678
679    angleResults[INDEX_EAST] = solarAzimuthAngle;
680    angleResults[INDEX_SOUTH] = solarElevationAngleCorrected_Deg;
681    timeResults[INDEX_SUNRISE] = (time_t) (sunriseTime * (double) secondsInADay);
682    timeResults[INDEX_SUNSET] = (time_t) (sunsetTime * (double) secondsInADay);
683    timeResults[INDEX_TIMESTAMP] = secsSince1970;
684}