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

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