FUNCTION CalcSolarPosition : BOOL
{
// Berechnet Sonnenposition nach NOAA-Algorithmus (vereinfacht für SPS)
// Laufzeit: < 80 µs auf typischer ARM-Cortex-SPS
// Deterministisch: Keine Schleifen, keine dynamische Speicherverwaltung, keine OS-Calls
}
VAR_INPUT
iYear : INT; // Jahr (z.B. 2026)
iMonth : INT; // Monat (1..12)
iDay : INT; // Tag (1..31)
iHour : INT; // Stunde (0..23) Ortszeit
iMinute : INT; // Minute (0..59) Ortszeit
iSecond : INT; // Sekunde (0..59) Ortszeit
fTimeZoneOffset : REAL; // UTC-Abweichung in Stunden (z.B. 1.0 für MEZ, 2.0 für MESZ)
fLatitude : REAL; // Breitengrad, +Nord / -Süd (z.B. 52.5200)
fLongitude : REAL; // Längengrad, +Ost / -West (z.B. 13.4050)
END_VAR
VAR_OUTPUT
fAzimuth : REAL; // Grad: 0°=Nord, 90°=Ost, 180°=Süd, 270°=West
fElevation : REAL; // Grad: 90°=Zenit, >0.0 = Sonne über Horizont
END_VAR
VAR
PI : REAL := 3.141592653589793;
DEG2RAD : REAL;
RAD2DEG : REAL;
iYearAdj : INT;
iMonthAdj : INT;
A : INT;
B : INT;
JD : REAL;
JC : REAL;
L0 : REAL;
M : REAL;
Mrad : REAL;
e : REAL;
C : REAL;
lambda : REAL;
Omega : REAL;
lambda_app : REAL;
eps0 : REAL;
delta : REAL;
DeclRad : REAL;
LatRad : REAL;
y_tan : REAL;
EqTime : REAL;
StdMeridian : REAL;
TST : REAL;
HA : REAL;
HARad : REAL;
cosZenith : REAL;
ZenithRad : REAL;
x_az : REAL;
y_az : REAL;
AzimuthRad : REAL;
END_VAR
BEGIN
// Konstanten & Umrechnungsfaktoren
DEG2RAD := PI / 180.0;
RAD2DEG := 180.0 / PI;
// 1. Julianischer Tag (JD) berechnen
IF iMonth <= 2 THEN
iYearAdj := iYear - 1;
iMonthAdj := iMonth + 12;
ELSE
iYearAdj := iYear;
iMonthAdj := iMonth;
END_IF;
A := INT(iYearAdj / 100);
B := 2 - A + INT(A / 4);
JD := INT(365.25 * (REAL(iYearAdj) + 4716.0)) + INT(30.6001 * (REAL(iMonthAdj) + 1.0)) + REAL(iDay) + REAL(B) - 1524.5;
JD := JD + (REAL(iHour) + REAL(iMinute)/60.0 + REAL(iSecond)/3600.0) / 24.0 - fTimeZoneOffset / 24.0;
// 2. Julianisches Jahrhundert
JC := (JD - 2451545.0) / 36525.0;
// 3. Geometrische Sonnenparameter
L0 := 280.46646 + JC * (36000.76983 + JC * 0.0003032);
L0 := L0 - 360.0 * INT(L0 / 360.0); // Normieren auf [0, 360)
M := 357.52911 + JC * (35999.05029 - JC * 0.0001537);
M := M - 360.0 * INT(M / 360.0);
Mrad := M * DEG2RAD;
e := 0.016708634 - JC * (0.000042037 + JC * 0.0000001267);
// Gleichung des Zentrums (C)
C := (1.914602 - JC * (0.004817 + JC * 0.000014)) * SIN(Mrad) +
(0.019993 - JC * 0.000101) * SIN(2.0 * Mrad) +
0.000289 * SIN(3.0 * Mrad);
lambda := L0 + C;
Omega := 125.04 - 1934.136 * JC;
lambda_app := lambda - 0.00569 - 0.00478 * SIN(Omega * DEG2RAD);
// 4. Schiefe der Ekliptik & Deklination
eps0 := 23.0 + (26.0 + (21.448 - JC * (46.8150 + JC * (0.00059 - JC * 0.001813))) / 60.0) / 60.0;
delta := RAD2DEG * ASIN(SIN(eps0 * DEG2RAD) * SIN(lambda_app * DEG2RAD));
DeclRad := delta * DEG2RAD;
LatRad := fLatitude * DEG2RAD;
// 5. Zeitgleichung (Equation of Time) in Minuten
y_tan := TAN(eps0 / 2.0 * DEG2RAD);
y_tan := y_tan * y_tan;
EqTime := 4.0 * RAD2DEG * (y_tan * SIN(2.0 * L0 * DEG2RAD) - 2.0 * e * SIN(Mrad) +
4.0 * e * y_tan * SIN(Mrad) * COS(2.0 * L0 * DEG2RAD) - 0.5 * y_tan * y_tan * SIN(4.0 * L0 * DEG2RAD) -
1.25 * e * e * SIN(2.0 * Mrad));
// 6. Wahre Sonnenzeit & Stundenwinkel
StdMeridian := fTimeZoneOffset * 15.0;
TST := (REAL(iHour) * 60.0 + REAL(iMinute) + REAL(iSecond) / 60.0) + EqTime + (fLongitude - StdMeridian) * 4.0;
HA := (TST - 720.0) / 4.0; // 720 min = 12:00 Uhr Sonnenmittag
IF HA > 180.0 THEN HA := HA - 360.0; ELSIF HA < -180.0 THEN HA := HA + 360.0; END_IF;
HARad := HA * DEG2RAD;
// 7. Höhenwinkel (Elevation) über Zenit
cosZenith := SIN(LatRad) * SIN(DeclRad) + COS(LatRad) * COS(DeclRad) * COS(HARad);
// Clampen gegen Domänenfehler bei ASIN/ACOS
IF cosZenith > 1.0 THEN cosZenith := 1.0;
ELSIF cosZenith < -1.0 THEN cosZenith := -1.0;
END_IF;
ZenithRad := ACOS(cosZenith);
fElevation := 90.0 - (ZenithRad * RAD2DEG);
// 8. Azimut (ATAN2 für korrekte Quadranten)
x_az := -COS(DeclRad) * SIN(HARad);
y_az := COS(LatRad) * SIN(DeclRad) - SIN(LatRad) * COS(DeclRad) * COS(HARad);
AzimuthRad := ATAN2(x_az, y_az);
fAzimuth := AzimuthRad * RAD2DEG;
IF fAzimuth < 0.0 THEN fAzimuth := fAzimuth + 360.0; END_IF;
CalcSolarPosition := TRUE;
END_FUNCTION