git.basschouten.com Git - openhab-addons.git/blob

   1 /**
   2  * Copyright (c) 2010-2022 Contributors to the openHAB project
   3  *
   4  * See the NOTICE file(s) distributed with this work for additional
   5  * information.
   6  *
   7  * This program and the accompanying materials are made available under the
   8  * terms of the Eclipse Public License 2.0 which is available at
   9  * http://www.eclipse.org/legal/epl-2.0
  10  *
  11  * SPDX-License-Identifier: EPL-2.0
  12  */
  13 package org.openhab.binding.astro.internal.calc;
  14
  15 import java.util.Calendar;
  16 import java.util.Map.Entry;
  17
  18 import org.openhab.binding.astro.internal.model.Eclipse;
  19 import org.openhab.binding.astro.internal.model.EclipseType;
  20 import org.openhab.binding.astro.internal.model.Position;
  21 import org.openhab.binding.astro.internal.model.Radiation;
  22 import org.openhab.binding.astro.internal.model.Range;
  23 import org.openhab.binding.astro.internal.model.Sun;
  24 import org.openhab.binding.astro.internal.model.SunPhaseName;
  25 import org.openhab.binding.astro.internal.util.DateTimeUtils;
  26
  27 /**
  28  * Calculates the SunPosition (azimuth, elevation) and Sun data.
  29  *
  30  * @author Gerhard Riegler - Initial contribution
  31  * @author Christoph Weitkamp - Introduced UoM
  32  * @see based on the calculations of http://www.suncalc.net
  33  */
  34 public class SunCalc {
  35     private static final double J2000 = 2451545.0;
  36     private static final double SC = 1367; // Solar constant in W/m²
  37     public static final double DEG2RAD = Math.PI / 180;
  38     public static final double RAD2DEG = 180. / Math.PI;
  39
  40     private static final double M0 = 357.5291 * DEG2RAD;
  41     private static final double M1 = 0.98560028 * DEG2RAD;
  42     private static final double J0 = 0.0009;
  43     private static final double J1 = 0.0053;
  44     private static final double J2 = -0.0069;
  45     private static final double C1 = 1.9148 * DEG2RAD;
  46     private static final double C2 = 0.0200 * DEG2RAD;
  47     private static final double C3 = 0.0003 * DEG2RAD;
  48     private static final double P = 102.9372 * DEG2RAD;
  49     private static final double E = 23.45 * DEG2RAD;
  50     private static final double TH0 = 280.1600 * DEG2RAD;
  51     private static final double TH1 = 360.9856235 * DEG2RAD;
  52     private static final double SUN_ANGLE = -0.83;
  53     private static final double SUN_DIAMETER = 0.53 * DEG2RAD; // sun diameter
  54     private static final double H0 = SUN_ANGLE * DEG2RAD;
  55     private static final double H1 = -6.0 * DEG2RAD; // nautical twilight angle
  56     private static final double H2 = -12.0 * DEG2RAD; // astronomical twilight
  57                                                       // angle
  58     private static final double H3 = -18.0 * DEG2RAD; // darkness angle
  59     private static final double MINUTES_PER_DAY = 60 * 24;
  60     private static final int CURVE_TIME_INTERVAL = 20; // 20 minutes
  61     private static final double JD_ONE_MINUTE_FRACTION = 1.0 / 60 / 24;
  62
  63     /**
  64      * Calculates the sun position (azimuth and elevation).
  65      */
  66     public void setPositionalInfo(Calendar calendar, double latitude, double longitude, Double altitude, Sun sun) {
  67         double lw = -longitude * DEG2RAD;
  68         double phi = latitude * DEG2RAD;
  69
  70         double j = DateTimeUtils.dateToJulianDate(calendar);
  71         double m = getSolarMeanAnomaly(j);
  72         double c = getEquationOfCenter(m);
  73         double lsun = getEclipticLongitude(m, c);
  74         double d = getSunDeclination(lsun);
  75         double a = getRightAscension(lsun);
  76         double th = getSiderealTime(j, lw);
  77
  78         double azimuth = getAzimuth(th, a, phi, d) / DEG2RAD;
  79         double elevation = getElevation(th, a, phi, d) / DEG2RAD;
  80         double shadeLength = getShadeLength(elevation);
  81
  82         Position position = sun.getPosition();
  83         position.setAzimuth(azimuth + 180);
  84         position.setElevation(elevation);
  85         position.setShadeLength(shadeLength);
  86
  87         setRadiationInfo(calendar, elevation, altitude, sun);
  88     }
  89
  90     /**
  91      * Calculates sun radiation data.
  92      */
  93     private void setRadiationInfo(Calendar calendar, double elevation, Double altitude, Sun sun) {
  94         double sinAlpha = Math.sin(DEG2RAD * elevation);
  95
  96         int dayOfYear = calendar.get(Calendar.DAY_OF_YEAR);
  97         int daysInYear = calendar.getActualMaximum(Calendar.DAY_OF_YEAR);
  98
  99         // Direct Solar Radiation (in W/m²) at the atmosphere entry
 100         // At sunrise/sunset - calculations limits are reached
 101         double rOut = (elevation > 3) ? SC * (0.034 * Math.cos(DEG2RAD * (360 * dayOfYear / daysInYear)) + 1) : 0;
 102         double altitudeRatio = (altitude != null) ? 1 / Math.pow((1 - (6.5 / 288) * (altitude / 1000.0)), 5.256) : 1;
 103         double m = (Math.sqrt(1229 + Math.pow(614 * sinAlpha, 2)) - 614 * sinAlpha) * altitudeRatio;
 104
 105         // Direct radiation after atmospheric layer
 106         // 0.6 = Coefficient de transmissivité
 107         double rDir = rOut * Math.pow(0.6, m) * sinAlpha;
 108
 109         // Diffuse Radiation
 110         double rDiff = rOut * (0.271 - 0.294 * Math.pow(0.6, m)) * sinAlpha;
 111         double rTot = rDir + rDiff;
 112
 113         Radiation radiation = sun.getRadiation();
 114         radiation.setDirect(rDir);
 115         radiation.setDiffuse(rDiff);
 116         radiation.setTotal(rTot);
 117     }
 118
 119     /**
 120      * Returns true, if the sun is up all day (no rise and set).
 121      */
 122     private boolean isSunUpAllDay(Calendar calendar, double latitude, double longitude, Double altitude) {
 123         Calendar cal = DateTimeUtils.truncateToMidnight(calendar);
 124         Sun sun = new Sun();
 125         for (int minutes = 0; minutes <= MINUTES_PER_DAY; minutes += CURVE_TIME_INTERVAL) {
 126             setPositionalInfo(cal, latitude, longitude, altitude, sun);
 127             if (sun.getPosition().getElevationAsDouble() < SUN_ANGLE) {
 128                 return false;
 129             }
 130             cal.add(Calendar.MINUTE, CURVE_TIME_INTERVAL);
 131         }
 132         return true;
 133     }
 134
 135     /**
 136      * Calculates all sun rise and sets at the specified coordinates.
 137      */
 138     public Sun getSunInfo(Calendar calendar, double latitude, double longitude, Double altitude,
 139             boolean useMeteorologicalSeason) {
 140         return getSunInfo(calendar, latitude, longitude, altitude, false, useMeteorologicalSeason);
 141     }
 142
 143     private Sun getSunInfo(Calendar calendar, double latitude, double longitude, Double altitude, boolean onlyAstro,
 144             boolean useMeteorologicalSeason) {
 145         double lw = -longitude * DEG2RAD;
 146         double phi = latitude * DEG2RAD;
 147         double j = DateTimeUtils.midnightDateToJulianDate(calendar) + 0.5;
 148         double n = getJulianCycle(j, lw);
 149         double js = getApproxSolarTransit(0, lw, n);
 150         double m = getSolarMeanAnomaly(js);
 151         double c = getEquationOfCenter(m);
 152         double lsun = getEclipticLongitude(m, c);
 153         double d = getSunDeclination(lsun);
 154         double jtransit = getSolarTransit(js, m, lsun);
 155         double w0 = getHourAngle(H0, phi, d);
 156         double w1 = getHourAngle(H0 + SUN_DIAMETER, phi, d);
 157         double jset = getSunsetJulianDate(w0, m, lsun, lw, n);
 158         double jsetstart = getSunsetJulianDate(w1, m, lsun, lw, n);
 159         double jrise = getSunriseJulianDate(jtransit, jset);
 160         double jriseend = getSunriseJulianDate(jtransit, jsetstart);
 161         double w2 = getHourAngle(H1, phi, d);
 162         double jnau = getSunsetJulianDate(w2, m, lsun, lw, n);
 163         double jciv2 = getSunriseJulianDate(jtransit, jnau);
 164
 165         double w3 = getHourAngle(H2, phi, d);
 166         double w4 = getHourAngle(H3, phi, d);
 167         double jastro = getSunsetJulianDate(w3, m, lsun, lw, n);
 168         double jdark = getSunsetJulianDate(w4, m, lsun, lw, n);
 169         double jnau2 = getSunriseJulianDate(jtransit, jastro);
 170         double jastro2 = getSunriseJulianDate(jtransit, jdark);
 171
 172         Sun sun = new Sun();
 173         sun.setAstroDawn(new Range(DateTimeUtils.toCalendar(jastro2), DateTimeUtils.toCalendar(jnau2)));
 174         sun.setAstroDusk(new Range(DateTimeUtils.toCalendar(jastro), DateTimeUtils.toCalendar(jdark)));
 175
 176         if (onlyAstro) {
 177             return sun;
 178         }
 179
 180         sun.setNoon(new Range(DateTimeUtils.toCalendar(jtransit),
 181                 DateTimeUtils.toCalendar(jtransit + JD_ONE_MINUTE_FRACTION)));
 182         sun.setRise(new Range(DateTimeUtils.toCalendar(jrise), DateTimeUtils.toCalendar(jriseend)));
 183         sun.setSet(new Range(DateTimeUtils.toCalendar(jsetstart), DateTimeUtils.toCalendar(jset)));
 184
 185         sun.setCivilDawn(new Range(DateTimeUtils.toCalendar(jciv2), DateTimeUtils.toCalendar(jrise)));
 186         sun.setCivilDusk(new Range(DateTimeUtils.toCalendar(jset), DateTimeUtils.toCalendar(jnau)));
 187
 188         sun.setNauticDawn(new Range(DateTimeUtils.toCalendar(jnau2), DateTimeUtils.toCalendar(jciv2)));
 189         sun.setNauticDusk(new Range(DateTimeUtils.toCalendar(jnau), DateTimeUtils.toCalendar(jastro)));
 190
 191         boolean isSunUpAllDay = isSunUpAllDay(calendar, latitude, longitude, altitude);
 192
 193         // daylight
 194         Range daylightRange = new Range();
 195         if (sun.getRise().getStart() == null && sun.getRise().getEnd() == null) {
 196             if (isSunUpAllDay) {
 197                 daylightRange = new Range(DateTimeUtils.truncateToMidnight(calendar),
 198                         DateTimeUtils.truncateToMidnight(addDays(calendar, 1)));
 199             }
 200         } else {
 201             daylightRange = new Range(sun.getRise().getEnd(), sun.getSet().getStart());
 202         }
 203         sun.setDaylight(daylightRange);
 204
 205         // morning night
 206         Sun sunYesterday = getSunInfo(addDays(calendar, -1), latitude, longitude, altitude, true,
 207                 useMeteorologicalSeason);
 208         Range morningNightRange = null;
 209         if (sunYesterday.getAstroDusk().getEnd() != null
 210                 && DateTimeUtils.isSameDay(sunYesterday.getAstroDusk().getEnd(), calendar)) {
 211             morningNightRange = new Range(sunYesterday.getAstroDusk().getEnd(), sun.getAstroDawn().getStart());
 212         } else if (isSunUpAllDay || sun.getAstroDawn().getStart() == null) {
 213             morningNightRange = new Range();
 214         } else {
 215             morningNightRange = new Range(DateTimeUtils.truncateToMidnight(calendar), sun.getAstroDawn().getStart());
 216         }
 217         sun.setMorningNight(morningNightRange);
 218
 219         // evening night
 220         Range eveningNightRange = null;
 221         if (sun.getAstroDusk().getEnd() != null && DateTimeUtils.isSameDay(sun.getAstroDusk().getEnd(), calendar)) {
 222             eveningNightRange = new Range(sun.getAstroDusk().getEnd(),
 223                     DateTimeUtils.truncateToMidnight(addDays(calendar, 1)));
 224         } else {
 225             eveningNightRange = new Range();
 226         }
 227         sun.setEveningNight(eveningNightRange);
 228
 229         // night
 230         if (isSunUpAllDay) {
 231             sun.setNight(new Range());
 232         } else {
 233             Sun sunTomorrow = getSunInfo(addDays(calendar, 1), latitude, longitude, altitude, true,
 234                     useMeteorologicalSeason);
 235             sun.setNight(new Range(sun.getAstroDusk().getEnd(), sunTomorrow.getAstroDawn().getStart()));
 236         }
 237
 238         // eclipse
 239         Eclipse eclipse = sun.getEclipse();
 240         MoonCalc mc = new MoonCalc();
 241
 242         eclipse.getKinds().forEach(eclipseKind -> {
 243             double jdate = mc.getEclipse(calendar, EclipseType.SUN, j, eclipseKind);
 244             eclipse.set(eclipseKind, DateTimeUtils.toCalendar(jdate), new Position());
 245         });
 246
 247         SunZodiacCalc zodiacCalc = new SunZodiacCalc();
 248         zodiacCalc.getZodiac(calendar).ifPresent(z -> sun.setZodiac(z));
 249
 250         SeasonCalc seasonCalc = new SeasonCalc();
 251         sun.setSeason(seasonCalc.getSeason(calendar, latitude, useMeteorologicalSeason));
 252
 253         // phase
 254         for (Entry<SunPhaseName, Range> rangeEntry : sun.getAllRanges().entrySet()) {
 255             SunPhaseName entryPhase = rangeEntry.getKey();
 256             if (rangeEntry.getValue().matches(Calendar.getInstance())) {
 257                 if (entryPhase == SunPhaseName.MORNING_NIGHT || entryPhase == SunPhaseName.EVENING_NIGHT) {
 258                     sun.getPhase().setName(SunPhaseName.NIGHT);
 259                 } else {
 260                     sun.getPhase().setName(entryPhase);
 261                 }
 262             }
 263         }
 264
 265         return sun;
 266     }
 267
 268     /**
 269      * Adds the specified days to the calendar.
 270      */
 271     private Calendar addDays(Calendar calendar, int days) {
 272         Calendar cal = (Calendar) calendar.clone();
 273         cal.add(Calendar.DAY_OF_MONTH, days);
 274         return cal;
 275     }
 276
 277     // all the following methods are translated to java based on the javascript
 278     // calculations of http://www.suncalc.net
 279     private double getJulianCycle(double j, double lw) {
 280         return Math.round(j - J2000 - J0 - lw / (2 * Math.PI));
 281     }
 282
 283     private double getApproxSolarTransit(double ht, double lw, double n) {
 284         return J2000 + J0 + (ht + lw) / (2 * Math.PI) + n;
 285     }
 286
 287     private double getSolarMeanAnomaly(double js) {
 288         return M0 + M1 * (js - J2000);
 289     }
 290
 291     private double getEquationOfCenter(double m) {
 292         return C1 * Math.sin(m) + C2 * Math.sin(2 * m) + C3 * Math.sin(3 * m);
 293     }
 294
 295     private double getEclipticLongitude(double m, double c) {
 296         return m + P + c + Math.PI;
 297     }
 298
 299     private double getSolarTransit(double js, double m, double lsun) {
 300         return js + (J1 * Math.sin(m)) + (J2 * Math.sin(2 * lsun));
 301     }
 302
 303     private double getSunDeclination(double lsun) {
 304         return Math.asin(Math.sin(lsun) * Math.sin(E));
 305     }
 306
 307     private double getRightAscension(double lsun) {
 308         return Math.atan2(Math.sin(lsun) * Math.cos(E), Math.cos(lsun));
 309     }
 310
 311     private double getSiderealTime(double j, double lw) {
 312         return TH0 + TH1 * (j - J2000) - lw;
 313     }
 314
 315     private double getAzimuth(double th, double a, double phi, double d) {
 316         double h = th - a;
 317         return Math.atan2(Math.sin(h), Math.cos(h) * Math.sin(phi) - Math.tan(d) * Math.cos(phi));
 318     }
 319
 320     private double getElevation(double th, double a, double phi, double d) {
 321         return Math.asin(Math.sin(phi) * Math.sin(d) + Math.cos(phi) * Math.cos(d) * Math.cos(th - a));
 322     }
 323
 324     private double getShadeLength(double elevation) {
 325         return 1 / Math.tan(elevation * DEG2RAD);
 326     }
 327
 328     private double getHourAngle(double h, double phi, double d) {
 329         return Math.acos((Math.sin(h) - Math.sin(phi) * Math.sin(d)) / (Math.cos(phi) * Math.cos(d)));
 330     }
 331
 332     private double getSunsetJulianDate(double w0, double m, double Lsun, double lw, double n) {
 333         return getSolarTransit(getApproxSolarTransit(w0, lw, n), m, Lsun);
 334     }
 335
 336     private double getSunriseJulianDate(double jtransit, double jset) {
 337         return jtransit - (jset - jtransit);
 338     }
 339 }