/*
Sonnenstandberechnung und Solartracking,
jetzt mit hochpräziser DS3231RTC Echtzeituhr
Sonnenstandsberechnung geht konform mit:
https://lexikon.astronomie.info/zeitgleichung/
- Mit LCD -Ausgabe -
- mit Spannungswert PV-Modul -
- Neu: mit Hall-Winkelsensor statt Drehpotentiometer -
Positionsangaben für Ort: Rönkhausen
Version 10.2, M. Schulte, 12. Juli 2019
Änderung: Variable Grenzwert für Rücklauf eingefügt
Änderung: Neuberechnung der Servoposition
Änderung: Steppermotor vor und zurück von 25 auf 32 steps pro loop
Änderung: Azimutpositionen neu kalibriert
Änderung: Adafruit-Motorshield Version 2.3 ausgetauscht gegen Version 1
Änderung: SaintSmart-20x4LCD ausgetauscht gegen Sunfounder LCD Display
Änderung: Datalogger-Shield mit DS1307RTC fällt weg, DS3231RTC-Uhr hinzugefügt
Änderung: Zeitbasis jetzt hochpräzise mit wenigen Sekunden Abweichung pro Jahr
Änderung: Neue Variablen: AziH, AziB
Änderung: I2C Adresse des Displays per Lötbrücke (A0) von 0x27 auf 0x26 geändert:
nun besteht die Möglichkeit, ein 2. Display (0x27) am I2C-Bus anzuschliessen
*/

#include <Servo.h>
#include <Wire.h>
#include <DS3231.h>
#include <LiquidCrystal_I2C.h>
#include <AFMotor.h>
#define RTC_I2C_ADDRESS 0x68                                    // I2C Adresse des RTC DS3231

AF_Stepper stepper(200, 2);
LiquidCrystal_I2C lcd(0x26, 20, 4);                             // I2C Adresse des Displays von 0x27 auf 0x26 geändert
RTCDateTime dt;
DS3231 clock;
Servo servo2;                                                   // Servomotor für Elevation

double pi = 3.1415926;
double kwert = pi / 180;
double latitude = 51.222191;
double longitude = 7.954169;
double Elevation;
double elevation;
double Aufgang;
double Untergang;

void setup() {

pinMode(A2, INPUT);                                             // Hall-Sensor
pinMode(A3, INPUT);                                             // PV-Panel
pinMode(14, OUTPUT);                                            // grüne LED, Anschluss an A0
pinMode(15, OUTPUT);                                            // rote LED, Anschluss an A1
servo2.attach(9);
servo2.write(5);                                                // Servo Grundstellung
stepper.setSpeed(15);                                           // Stepper 15 U/min
lcd.init();                                                     // initialisiere 20x4 LCD
lcd.backlight();
lcd.clear();
clock.begin();
//clock.setDateTime(2019, 06, 29, 8, 43, 00);                   // Hier Datum und aktuelle Uhrzeit (MEZ) setzen
}


// Beginn der Schleife

void loop() {

dt = clock.getDateTime();                                       // Zeitwert holen von RTC Time Echtzeituhr

String hour = clock.dateFormat("H", dt);
String minute = clock.dateFormat("i", dt);
String year = clock.dateFormat("z", dt);
String second = clock.dateFormat("s", dt);

int yearday = year.toInt() + 1;
int Jahr = (dt.year);
int Monat = (dt.month);
int Tag = (dt.day);
double Stunde = hour.toInt();
double Minute = minute.toInt();
double Sekunde = second.toInt();

double Deklination;
double Azimut;
double Zeitgleichung;
double Zeitdifferenz;
double AufgangOrtszeit;
double UntergangOrtszeit;
double Tageslaenge;
double Refraktion;
double MEZ;                                                     // Mitteleuropäische Zeit
double MOZ;                                                     // Mittlere Ortszeit
double WOZ;                                                     // Wahre Ortszeit
double ZeitSeitMittag;
double B = latitude * kwert;                                    // geogr. Breite in Radians
double R = 0.00;                                                // Refraktion Anfangswert
double h = -(50.0 / 60.0) * kwert;                              // Höhe des Sonnenmittelpunkts bei Aufgang: Radius + Refraktion
int Grenzwert;
int Zeitzone = 1;                                               // Zeitzone Greenwich + 1 (z.B. Berlin)
float P = 1013.25;                                              // Luftdruck der Standard-Atmosphäre in hPa
float T = 15.0;                                                 // Temperatur der Standard-Atmosphäre in °C
int zeitschritt = 1000;                                         // Zeitschritt loop = 1 sec
int sensorValue;                                                // Analoger Ist-Wert (Hall-Sensor)
int azimutHall;                                                 // Gradwert Ist-Azimutwert (Hall-Sensor)

MEZ = Stunde + Minute / 60 + Sekunde / 3600;

Zeitgleichung = BerechneZeitgleichung (yearday);                             // Aufruf der Funktion: BerechneZeitgleichung
Deklination = BerechneDeklination (yearday);                                 // Aufruf der Funktion: BerechneDeklination
Zeitdifferenz = BerechneZeitdifferenz (h, latitude, kwert, Deklination, pi); // Aufruf der Funktion: BerechneZeitdifferenz
ZeitSeitMittag = MEZ + longitude / 15.0 - Zeitzone - 12 + Zeitgleichung;     // Berechne ZeitSeitMittag
Azimut = SunAngles(Deklination, B, ZeitSeitMittag, Elevation);               // Aufruf der Funktion: SunAngles
Refraktion = BerechneRwert(Elevation, kwert, P, T);                          // Aufruf der Funktion: BerechneRwert
Elevation = (Elevation + Refraktion) / kwert;                                // Höhe mit Refraktionskorrektur in Grad
Tageslaenge = BerechneTageslaenge (Zeitdifferenz, Zeitgleichung, longitude, Zeitzone, Aufgang, Untergang);   // Aufruf der Funktion: BerechneTageslaenge

MOZ = MEZ - (-longitude / 15) - 1;                                           // Berechne Mittlere Ortszeit aus MEZ
WOZ = MOZ + Zeitgleichung;                                                   // Berechne Wahre Ortszeit

int WOZdez = (int)WOZ;
int WOZmin = (WOZ - WOZdez) * 60;
int WOZsec = (WOZ - WOZdez) * 60;

int MOZdez = (int)MOZ;
int MOZmin = (MOZ - MOZdez) * 60;
int MOZsec = (MOZ - MOZdez) * 60;

int StdL = (int)Tageslaenge;
int MinL = (int)((Tageslaenge - StdL) * 60);

int StdA = (int)Aufgang;
int MinA = (int)((Aufgang - StdA) * 60);

int StdU = (int)Untergang;
int MinU = (int)((Untergang - StdU) * 60);

if (StdL >= 7 && StdL < 9) {
Grenzwert = 135;
};                                                                           // Grenzwerte für Rücklaufposition
if (StdL >= 9 && StdL < 10) {
Grenzwert = 120;
};                                                                           // abhängig von der Tageslänge (nur Stunde)
if (StdL >= 10 && StdL < 11) {
Grenzwert = 105;
};
if (StdL >= 11 && StdL < 13) {
Grenzwert = 90;
};
if (StdL >= 13 && StdL < 14) {
Grenzwert = 75;
};
if (StdL >= 14 && StdL < 15) {
Grenzwert = 60;
};
if (StdL >= 15 && StdL < 17) {
Grenzwert = 45;
};


// Servomotor für Elevation


elevation = map(Elevation, 90, 0, 105, 3);               // Kalibrierung der Servoposition (elevation)
servo2.write(elevation);                                 // Schreiben der Servoposition (Nachführung)


// Steppermotor für Azimutposition

sensorValue = analogRead(A2);                            // Auslesen Hall-Sensorwert an Pin A2

if (sensorValue >= 125 && sensorValue < 506) {           // Skalierung der Hall-Sensorwerte in Gradwerte 45-180 Grad
azimutHall = map(sensorValue, 125, 506, 45, 180);
}
if (sensorValue >= 506 && sensorValue <= 886) {          // Skalierung der Hall-Sensorwerte in Gradwerte 180-315 Grad
azimutHall = map(sensorValue, 506, 886, 180, 315);
}

int AziH = azimutHall;
int AziB = (int)Azimut;

if (Elevation > 0) {                                     // Solange Elevation > 0 --> Sonne steht über dem Horizont

if (azimutHall < (int)Azimut) {                          // Azimut-Nachführung: Wenn Hall-Sensorwert kleiner als Azimut, dann:
digitalWrite(14, HIGH);                                  // LED grün an
stepper.step(32, FORWARD, MICROSTEP);                    // Steppermotor dreht 32 steps vorwärts (entspricht +1 Grad des Drehkranzes)
digitalWrite(15, LOW);                                   // LED rot aus
}
if (azimutHall > (int)Azimut) {                          // Azimut-Nachführung: Wenn Hall-Sensorwert grösser als Azimut, dann:
digitalWrite(15, HIGH);                                  // LED rot an
stepper.step(32, BACKWARD, MICROSTEP);                   // Steppermotor dreht 32 steps rückwärts (entspricht -1 Grad des Drehkranzes)
digitalWrite(14, LOW);                                   // LED grün aus
}
}
else {                                                   // Sonne ist untergegangen --> Rückführung auf Anfangsposition
if (Elevation < 0 && (azimutHall - 1) >= Grenzwert) {    // Azimut-Rückführung: Wenn Elevation < -1 (Sonne ist untergegangen), dann:
digitalWrite(15, HIGH);                                  // LED rot an: Stepper dreht rückwärts
stepper.setSpeed(30);                                    // Stepper 30 U/min
stepper.step(200, BACKWARD, INTERLEAVE);                 // Schneller Rücklauf mit 200 steps/loop solange Hall-Sensorwert >= Grenzwert
digitalWrite(14, LOW);                                   // LED grün aus
}
}

if (AziH = AziB) {                                       // Damit alle LEDs aus sind beim Erreichen der Zielposition (azimutHall = Azimut)
digitalWrite(14, LOW);
digitalWrite(15, LOW);
}


// Deklarationen und Typumwandlungen für LCD-Anzeige

const char* windrichtung[16] = {"N ", "NNO", "NO ", "ONO", "O ", "OSO", "SO ", "SSO", "S ", "SSW", "SW ", "WSW", "W ", "WNW", "NW ", "NNW"};

int i = Windrichtung(azimutHall);                        // Aufruf Funktion: Windrichtung
String wr = windrichtung[i];

int analogPin = A3;                                      // Solarpanel Pluspol an A3
float pvmod;
float volt;
byte stellen = 2;
pvmod = analogRead(analogPin);                           // Abfrage Spannung Solar Panel
volt = pvmod / 1023 * 4.75;                              // Umrechnung Analogwert in Spannungswert (Volt)

int Std = (int)Tageslaenge;
int Min = (int)((Tageslaenge - Std) * 60);


// Ausgabe der Werte auf 20x4 LCD-Display

lcd.setCursor(0, 0);
lcd.print(clock.dateFormat("d.M.Y H:i:s", dt));

lcd.setCursor(0, 1);
lcd.print("AzS");
lcd.setCursor(4, 1);
lcd.print(Azimut);
lcd.write(0xDF);
lcd.print(" ");

lcd.setCursor(12, 1);
lcd.print("WZ ");
if (WOZdez < 10) lcd.print("0");
lcd.print(WOZdez);
lcd.print(":");
if (WOZmin < 10) lcd.print("0");
lcd.print(WOZmin);

lcd.setCursor(0, 2);
lcd.print("AzI ");
lcd.print(azimutHall);
lcd.write(0xDF);
lcd.print(wr);
lcd.print(" ");

lcd.setCursor(12, 2);
lcd.print("TL ");
if (Std < 10) lcd.print("0");
lcd.print(Std);
lcd.print(":");
if (Min < 10) lcd.print("0");
lcd.print(Min);

lcd.setCursor(0, 3);
lcd.print("Elv");
lcd.setCursor(12, 3);
lcd.print("PV");
lcd.setCursor(3, 3);
lcd.print(" ");
lcd.setCursor(4, 3);
lcd.print(Elevation, 2);
lcd.write(0xDF);
lcd.setCursor(15, 3);
lcd.print(volt, 1);
lcd.print(" V");

delay(zeitschritt); // Programmpause
}


// Funktion: Windrichtung

int Windrichtung (float azimutHall) {

int i;

if (azimutHall >= 0 && azimutHall < 11) {
i = 0;
}
if (azimutHall >= 11 && azimutHall < 34) {
i = 1;
}
if (azimutHall >= 34 && azimutHall < 56) {
i = 2;
}
if (azimutHall >= 56 && azimutHall < 79) {
i = 3;
}
if (azimutHall >= 79 && azimutHall < 101) {
i = 4;
}
if (azimutHall >= 101 && azimutHall < 124) {
i = 5;
}
if (azimutHall >= 124 && azimutHall < 146) {
i = 6;
}
if (azimutHall >= 146 && azimutHall < 169) {
i = 7;
}
if (azimutHall >= 169 && azimutHall < 191) {
i = 8;
}
if (azimutHall >= 191 && azimutHall < 214) {
i = 9;
}
if (azimutHall >= 214 && azimutHall < 236) {
i = 10;
}
if (azimutHall >= 236 && azimutHall < 259) {
i = 11;
}
if (azimutHall >= 259 && azimutHall < 281) {
i = 12;
}
if (azimutHall >= 281 && azimutHall < 304) {
i = 13;
}
if (azimutHall >= 304 && azimutHall < 326) {
i = 14;
}
if (azimutHall >= 326 && azimutHall < 349) {
i = 15;
}
return i;
}


// Funktion: SunAngles für die Berechnung von Azimut und Elevation

double SunAngles (double Deklination, double B, double ZeitSeitMittag, double & Elevation) {
double DK = Deklination;
double cosdec = cos(DK);
double sindec = sin(DK);
double lha = ZeitSeitMittag * (1.0027379 - 1 / 365.25) * 15 * kwert;
double coslha = cos(lha);
double sinlha = sin(lha);
double coslat = cos(B);
double sinlat = sin(B);
double N = -cosdec * sinlha;
double D = sindec * coslat - cosdec * coslha * sinlat;
Elevation = asin(sindec * sinlat + cosdec * coslha * coslat);             // Höhe des Sonnenmittelpunktes über dem Horizont
double Azimut = atan2(N, D);
if (Azimut < 0) {
Azimut += 2 * pi;
}
Azimut = Azimut / kwert;
return Azimut;
}


// Funktion: BerechneZeitgleichung

double BerechneZeitgleichung (int yearday) {
double Zeitgleichung = -0.170869921174742 * sin(0.0336997028793971 * yearday + 0.465419984181394) - 0.129890681040717 * sin(0.0178674832556871 * yearday - 0.167936777524864);
return Zeitgleichung;
}


// Funktion: BerechneDeklination

double BerechneDeklination(int yearday) {
double Deklination = 0.409526325277017 * sin(0.0169060504029192 * (yearday - 80.0856919827619));
return Deklination;
}


// Funktion: BerechneZeitdifferenz

double BerechneZeitdifferenz(double h, double latitude, double kwert, double Deklination, double pi) {
double Zeitdifferenz;
Zeitdifferenz = 12 * acos((sin(h) - sin(latitude * kwert) * sin(Deklination)) / (cos(latitude * kwert) * cos(Deklination))) / pi;
return Zeitdifferenz;
}


// Funktion: BerechneRwert
// Näherungslösung für die Refraktion für ein Objekt bei geringer Höhe über dem mathematischen Horizont
// Refraktion beträgt bei Sonnenaufgang 34 Bogenminuten = 0.56667°
// Falls die Höhe der Sonne nicht genauer als auf 0.5° gewünscht ist, kann diese Funktion ignoriert werden

double BerechneRwert (double Elevation, double kwert, float P, float T) {
double R;
if (Elevation >= 15 * kwert) {
R = 0.00452 * kwert * P / tan(Elevation) / (273 + T); // über 15° - einfachere Formel
}
else if (Elevation > -1 * kwert) {
R = kwert * P * (0.1594 + 0.0196 * Elevation + 0.00002 * (Elevation * Elevation)) / ((273 + T) * (1 + 0.505 * Elevation + 0.0845 * (Elevation * Elevation)));
}
return R;
}


// Funktion: BerechneTageslaenge

double BerechneTageslaenge (double Zeitdifferenz, double Zeitgleichung, double longitude, int Zeitzone, double & Aufgang, double & Untergang) {
double AufgangOrtszeit = 12 - Zeitdifferenz - Zeitgleichung;
double UntergangOrtszeit = 12 + Zeitdifferenz - Zeitgleichung;
Aufgang = AufgangOrtszeit - longitude / 15 + Zeitzone;
Untergang = UntergangOrtszeit - longitude / 15 + Zeitzone;
double Tageslaenge = Untergang - Aufgang;
return Tageslaenge;
}