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.math.BigDecimal;
  16 import java.math.RoundingMode;
  17 import java.util.Calendar;
  18
  19 import org.openhab.binding.astro.internal.model.Eclipse;
  20 import org.openhab.binding.astro.internal.model.EclipseKind;
  21 import org.openhab.binding.astro.internal.model.EclipseType;
  22 import org.openhab.binding.astro.internal.model.Moon;
  23 import org.openhab.binding.astro.internal.model.MoonDistance;
  24 import org.openhab.binding.astro.internal.model.MoonPhase;
  25 import org.openhab.binding.astro.internal.model.MoonPhaseName;
  26 import org.openhab.binding.astro.internal.model.Position;
  27 import org.openhab.binding.astro.internal.model.Range;
  28 import org.openhab.binding.astro.internal.model.Zodiac;
  29 import org.openhab.binding.astro.internal.model.ZodiacSign;
  30 import org.openhab.binding.astro.internal.util.DateTimeUtils;
  31
  32 /**
  33  * Calculates the phase, eclipse, rise, set, distance, illumination and age of
  34  * the moon.
  35  *
  36  * @author Gerhard Riegler - Initial contribution
  37  * @author Christoph Weitkamp - Introduced UoM
  38  * @implNote based on the calculations of
  39  *           http://www.computus.de/mondphase/mondphase.htm azimuth/elevation and
  40  *           zodiac based on http://lexikon.astronomie.info/java/sunmoon/
  41  */
  42 public class MoonCalc {
  43     private static final double NEW_MOON = 0;
  44     private static final double FULL_MOON = 0.5;
  45     private static final double FIRST_QUARTER = 0.25;
  46     private static final double LAST_QUARTER = 0.75;
  47
  48     /**
  49      * Calculates all moon data at the specified coordinates
  50      */
  51     public Moon getMoonInfo(Calendar calendar, double latitude, double longitude) {
  52         Moon moon = new Moon();
  53
  54         double julianDate = DateTimeUtils.dateToJulianDate(calendar);
  55         double julianDateMidnight = DateTimeUtils.midnightDateToJulianDate(calendar);
  56
  57         double[] riseSet = getRiseSet(calendar, latitude, longitude);
  58         Calendar rise = DateTimeUtils.timeToCalendar(calendar, riseSet[0]);
  59         Calendar set = DateTimeUtils.timeToCalendar(calendar, riseSet[1]);
  60
  61         if (rise == null || set == null) {
  62             Calendar tomorrow = (Calendar) calendar.clone();
  63             tomorrow.add(Calendar.DAY_OF_MONTH, 1);
  64
  65             double[] riseSeTomorrow = getRiseSet(tomorrow, latitude, longitude);
  66             if (rise == null) {
  67                 rise = DateTimeUtils.timeToCalendar(tomorrow, riseSeTomorrow[0]);
  68             }
  69             if (set == null) {
  70                 set = DateTimeUtils.timeToCalendar(tomorrow, riseSeTomorrow[1]);
  71             }
  72         }
  73
  74         moon.setRise(new Range(rise, rise));
  75         moon.setSet(new Range(set, set));
  76
  77         MoonPhase phase = moon.getPhase();
  78         phase.setNew(DateTimeUtils.toCalendar(getNextPhase(calendar, julianDateMidnight, NEW_MOON)));
  79         phase.setFirstQuarter(DateTimeUtils.toCalendar(getNextPhase(calendar, julianDateMidnight, FIRST_QUARTER)));
  80         phase.setFull(DateTimeUtils.toCalendar(getNextPhase(calendar, julianDateMidnight, FULL_MOON)));
  81         phase.setThirdQuarter(DateTimeUtils.toCalendar(getNextPhase(calendar, julianDateMidnight, LAST_QUARTER)));
  82
  83         Eclipse eclipse = moon.getEclipse();
  84         eclipse.getKinds().forEach(eclipseKind -> {
  85             double jdate = getEclipse(calendar, EclipseType.MOON, julianDateMidnight, eclipseKind);
  86             eclipse.set(eclipseKind, DateTimeUtils.toCalendar(jdate), new Position());
  87         });
  88
  89         double decimalYear = DateTimeUtils.getDecimalYear(calendar);
  90         MoonDistance apogee = moon.getApogee();
  91         double apogeeJd = getApogee(julianDate, decimalYear);
  92         apogee.setDate(DateTimeUtils.toCalendar(apogeeJd));
  93         apogee.setDistance(getDistance(apogeeJd));
  94
  95         MoonDistance perigee = moon.getPerigee();
  96         double perigeeJd = getPerigee(julianDate, decimalYear);
  97         perigee.setDate(DateTimeUtils.toCalendar(perigeeJd));
  98         perigee.setDistance(getDistance(perigeeJd));
  99
 100         return moon;
 101     }
 102
 103     /**
 104      * Calculates the moon illumination and distance.
 105      */
 106     public void setPositionalInfo(Calendar calendar, double latitude, double longitude, Moon moon) {
 107         double julianDate = DateTimeUtils.dateToJulianDate(calendar);
 108         setMoonPhase(calendar, moon);
 109         setAzimuthElevationZodiac(julianDate, latitude, longitude, moon);
 110
 111         MoonDistance distance = moon.getDistance();
 112         distance.setDate(Calendar.getInstance());
 113         distance.setDistance(getDistance(julianDate));
 114     }
 115
 116     /**
 117      * Calculates the age and the current phase.
 118      */
 119     private void setMoonPhase(Calendar calendar, Moon moon) {
 120         MoonPhase phase = moon.getPhase();
 121         double julianDateEndOfDay = DateTimeUtils.endOfDayDateToJulianDate(calendar);
 122         double parentNewMoon = getPreviousPhase(calendar, julianDateEndOfDay, NEW_MOON);
 123         double age = Math.abs(parentNewMoon - julianDateEndOfDay);
 124         phase.setAge(age);
 125
 126         long parentNewMoonMillis = DateTimeUtils.toCalendar(parentNewMoon).getTimeInMillis();
 127         long ageRangeTimeMillis = phase.getNew().getTimeInMillis() - parentNewMoonMillis;
 128         long ageCurrentMillis = System.currentTimeMillis() - parentNewMoonMillis;
 129         double agePercent = ageRangeTimeMillis != 0 ? ageCurrentMillis * 100.0 / ageRangeTimeMillis : 0;
 130         phase.setAgePercent(agePercent);
 131         phase.setAgeDegree(3.6 * agePercent);
 132         double illumination = getIllumination(DateTimeUtils.dateToJulianDate(calendar));
 133         phase.setIllumination(illumination);
 134         boolean isWaxing = age < (29.530588853 / 2);
 135         if (DateTimeUtils.isSameDay(calendar, phase.getNew())) {
 136             phase.setName(MoonPhaseName.NEW);
 137         } else if (DateTimeUtils.isSameDay(calendar, phase.getFirstQuarter())) {
 138             phase.setName(MoonPhaseName.FIRST_QUARTER);
 139         } else if (DateTimeUtils.isSameDay(calendar, phase.getThirdQuarter())) {
 140             phase.setName(MoonPhaseName.THIRD_QUARTER);
 141         } else if (DateTimeUtils.isSameDay(calendar, phase.getFull())) {
 142             phase.setName(MoonPhaseName.FULL);
 143         } else if (illumination >= 0 && illumination < 50) {
 144             phase.setName(isWaxing ? MoonPhaseName.WAXING_CRESCENT : MoonPhaseName.WANING_CRESCENT);
 145         } else if (illumination >= 50 && illumination < 100) {
 146             phase.setName(isWaxing ? MoonPhaseName.WAXING_GIBBOUS : MoonPhaseName.WANING_GIBBOUS);
 147         }
 148     }
 149
 150     /**
 151      * Calculates moonrise and moonset.
 152      */
 153     private double[] getRiseSet(Calendar calendar, double latitude, double longitude) {
 154         double lambda = prepareCoordinate(longitude, 180);
 155         if (longitude > 0) {
 156             lambda *= -1;
 157         }
 158         double phi = prepareCoordinate(latitude, 90);
 159         if (latitude < 0) {
 160             phi *= -1;
 161         }
 162
 163         double moonJd = Math.floor(DateTimeUtils.midnightDateToJulianDate(calendar)) - 2400000.0;
 164         moonJd -= ((calendar.get(Calendar.ZONE_OFFSET) + calendar.get(Calendar.DST_OFFSET)) / 60000.0) / 1440.0;
 165
 166         double sphi = SN(phi);
 167         double cphi = CS(phi);
 168         double sinho = SN(8.0 / 60.0);
 169
 170         int hour = 1;
 171         double utrise = -1;
 172         double utset = -1;
 173         do {
 174             double yminus = SINALT(moonJd, hour - 1, lambda, cphi, sphi) - sinho;
 175             double yo = SINALT(moonJd, hour, lambda, cphi, sphi) - sinho;
 176             double yplus = SINALT(moonJd, hour + 1, lambda, cphi, sphi) - sinho;
 177             double[] quadRet = QUAD(yminus, yo, yplus);
 178             if (quadRet[3] == 1) {
 179                 if (yminus < 0) {
 180                     utrise = hour + quadRet[1];
 181                 } else {
 182                     utset = hour + quadRet[1];
 183                 }
 184             }
 185             if (quadRet[3] == 2) {
 186                 if (quadRet[0] < 0) {
 187                     utrise = hour + quadRet[2];
 188                     utset = hour + quadRet[1];
 189                 } else {
 190                     utrise = hour + quadRet[1];
 191                     utset = hour + quadRet[2];
 192                 }
 193             }
 194             yminus = yplus;
 195             hour += 2;
 196         } while (hour < 25 && (utrise == -1 || utset == -1));
 197
 198         double rise = prepareTime(utrise);
 199         double set = prepareTime(utset);
 200
 201         return new double[] { rise, set };
 202     }
 203
 204     /**
 205      * Prepares the coordinate for moonrise and moonset calculation.
 206      */
 207     private double prepareCoordinate(double coordinate, double system) {
 208         double c = Math.abs(coordinate);
 209
 210         if (c - Math.floor(c) >= .599) {
 211             c = Math.floor(c) + (c - Math.floor(c)) / 1 * .6;
 212         }
 213         if (c > system) {
 214             c = Math.floor(c) % system + (c - Math.floor(c));
 215         }
 216         return Math.round(c * 100.0) / 100.0;
 217     }
 218
 219     /**
 220      * Prepares a time value for converting to a calendar object.
 221      */
 222     private double prepareTime(double riseSet) {
 223         if (riseSet == -1) {
 224             return riseSet;
 225         }
 226         double riseMinute = (riseSet - Math.floor(riseSet)) * 60.0 / 100.0;
 227         double rounded;
 228         if (riseMinute >= .595) {
 229             riseMinute = 0;
 230             rounded = riseSet + 1;
 231         } else {
 232             rounded = riseSet;
 233         }
 234         rounded = Math.floor(rounded) + riseMinute;
 235
 236         BigDecimal bd = new BigDecimal(Double.toString(rounded));
 237         bd = bd.setScale(2, RoundingMode.HALF_UP);
 238         return bd.doubleValue();
 239     }
 240
 241     /**
 242      * Calculates the moon phase.
 243      */
 244     private double calcMoonPhase(double k, double mode) {
 245         double kMod = Math.floor(k) + mode;
 246         double t = kMod / 1236.85;
 247         double e = var_e(t);
 248         double m = var_m(kMod, t);
 249         double m1 = var_m1(kMod, t);
 250         double f = var_f(kMod, t);
 251         double o = var_o(kMod, t);
 252         double jd = var_jde(kMod, t);
 253         if (mode == NEW_MOON) {
 254             jd += -.4072 * SN(m1) + .17241 * e * SN(m) + .01608 * SN(2 * m1) + .01039 * SN(2 * f)
 255                     + .00739 * e * SN(m1 - m) - .00514 * e * SN(m1 + m) + .00208 * e * e * SN(2 * m)
 256                     - .00111 * SN(m1 - 2 * f) - .00057 * SN(m1 + 2 * f);
 257             jd += .00056 * e * SN(2 * m1 + m) - .00042 * SN(3 * m1) + .00042 * e * SN(m + 2 * f)
 258                     + .00038 * e * SN(m - 2 * f) - .00024 * e * SN(2 * m1 - m) - .00017 * SN(o)
 259                     - .00007 * SN(m1 + 2 * m) + .00004 * SN(2 * m1 - 2 * f);
 260             jd += .00004 * SN(3 * m) + .00003 * SN(m1 + m - 2 * f) + .00003 * SN(2 * m1 + 2 * f)
 261                     - .00003 * SN(m1 + m + 2 * f) + .00003 * SN(m1 - m + 2 * f) - .00002 * SN(m1 - m - 2 * f)
 262                     - .00002 * SN(3 * m1 + m);
 263             jd += .00002 * SN(4 * m1);
 264         } else if (mode == FULL_MOON) {
 265             jd += -.40614 * SN(m1) + .17302 * e * SN(m) + .01614 * SN(2 * m1) + .01043 * SN(2 * f)
 266                     + .00734 * e * SN(m1 - m) - .00515 * e * SN(m1 + m) + .00209 * e * e * SN(2 * m)
 267                     - .00111 * SN(m1 - 2 * f) - .00057 * SN(m1 + 2 * f);
 268             jd += .00056 * e * SN(2 * m1 + m) - .00042 * SN(3 * m1) + .00042 * e * SN(m + 2 * f)
 269                     + .00038 * e * SN(m - 2 * f) - .00024 * e * SN(2 * m1 - m) - .00017 * SN(o)
 270                     - .00007 * SN(m1 + 2 * m) + .00004 * SN(2 * m1 - 2 * f);
 271             jd += .00004 * SN(3 * m) + .00003 * SN(m1 + m - 2 * f) + .00003 * SN(2 * m1 + 2 * f)
 272                     - .00003 * SN(m1 + m + 2 * f) + .00003 * SN(m1 - m + 2 * f) - .00002 * SN(m1 - m - 2 * f)
 273                     - .00002 * SN(3 * m1 + m);
 274             jd += .00002 * SN(4 * m1);
 275         } else {
 276             jd += -.62801 * SN(m1) + .17172 * e * SN(m) - .01183 * e * SN(m1 + m) + .00862 * SN(2 * m1)
 277                     + .00804 * SN(2 * f) + .00454 * e * SN(m1 - m) + .00204 * e * e * SN(2 * m) - .0018 * SN(m1 - 2 * f)
 278                     - .0007 * SN(m1 + 2 * f);
 279             jd += -.0004 * SN(3 * m1) - .00034 * e * SN(2 * m1 - m) + .00032 * e * SN(m + 2 * f)
 280                     + .00032 * e * SN(m - 2 * f) - .00028 * e * e * SN(m1 + 2 * m) + .00027 * e * SN(2 * m1 + m)
 281                     - .00017 * SN(o);
 282             jd += -.00005 * SN(m1 - m - 2 * f) + .00004 * SN(2 * m1 + 2 * f) - .00004 * SN(m1 + m + 2 * f)
 283                     + .00004 * SN(m1 - 2 * m) + .00003 * SN(m1 + m - 2 * f) + .00003 * SN(3 * m)
 284                     + .00002 * SN(2 * m1 - 2 * f);
 285             jd += .00002 * SN(m1 - m + 2 * f) - .00002 * SN(3 * m1 + m);
 286             double w = .00306 - .00038 * e * CS(m) + .00026 * CS(m1) - .00002 * CS(m1 - m) + .00002 * CS(m1 + m)
 287                     + .00002 * CS(2 * f);
 288             jd += (mode == FIRST_QUARTER) ? w : -w;
 289         }
 290         return moonCorrection(jd, t, kMod);
 291     }
 292
 293     /**
 294      * Calculates the eclipse.
 295      */
 296     private double getEclipse(double k, EclipseType typ, EclipseKind eclipse) {
 297         double kMod = Math.floor(k) + ((typ == EclipseType.SUN) ? 0 : 0.5);
 298         double t = kMod / 1236.85;
 299         double f = var_f(kMod, t);
 300         double jd = 0;
 301         double ringTest = 0;
 302         if (SN(Math.abs(f)) <= .36) {
 303             double o = var_o(kMod, t);
 304             double f1 = f - .02665 * SN(o);
 305             double a1 = 299.77 + .107408 * kMod - .009173 * t * t;
 306             double e = var_e(t);
 307             double m = var_m(kMod, t);
 308             double m1 = var_m1(kMod, t);
 309             double p = .207 * e * SN(m) + .0024 * e * SN(2 * m) - .0392 * SN(m1) + .0116 * SN(2 * m1)
 310                     - .0073 * e * SN(m1 + m) + .0067 * e * SN(m1 - m) + .0118 * SN(2 * f1);
 311             double q = 5.2207 - .0048 * e * CS(m) + .002 * e * CS(2 * m) - .3299 * CS(m1) - .006 * e * CS(m1 + m)
 312                     + .0041 * e * CS(m1 - m);
 313             double g = (p * CS(f1) + q * SN(f1)) * (1 - .0048 * CS(Math.abs(f1)));
 314             double u = .0059 + .0046 * e * CS(m) - .0182 * CS(m1) + .0004 * CS(2 * m1) - .0005 * CS(m + m1);
 315             jd = var_jde(kMod, t);
 316             jd += (typ == EclipseType.MOON) ? -.4065 * SN(m1) + .1727 * e * SN(m) : -.4075 * SN(m1) + .1721 * e * SN(m);
 317
 318             jd += .0161 * SN(2 * m1) - .0097 * SN(2 * f1) + .0073 * e * SN(m1 - m) - .005 * e * SN(m1 + m)
 319                     - .0023 * SN(m1 - 2 * f1) + .0021 * e * SN(2 * m);
 320             jd += .0012 * SN(m1 + 2 * f1) + .0006 * e * SN(2 * m1 + m) - .0004 * SN(3 * m1) - .0003 * e * SN(m + 2 * f1)
 321                     + .0003 * SN(a1) - .0002 * e * SN(m - 2 * f1) - .0002 * e * SN(2 * m1 - m) - .0002 * SN(o);
 322             switch (typ) {
 323                 case MOON:
 324                     if ((1.0248 - u - Math.abs(g)) / .545 <= 0) {
 325                         jd = 0; // no moon eclipse
 326                     }
 327                     if (eclipse == EclipseKind.PARTIAL && (1.0128 - u - Math.abs(g)) / .545 > 0
 328                             && (.4678 - u) * (.4678 - u) - g * g > 0) {
 329                         jd = 0; // no partial moon eclipse
 330                     }
 331                     if (eclipse == EclipseKind.TOTAL
 332                             && ((1.0128 - u - Math.abs(g)) / .545 <= 0 != (.4678 - u) * (.4678 - u) - g * g <= 0)) {
 333                         jd = 0; // no total moon eclipse
 334                     }
 335                     break;
 336                 case SUN:
 337                     if (Math.abs(g) > 1.5433 + u) {
 338                         jd = 0; // no sun eclipse
 339                     }
 340                     if (eclipse == EclipseKind.PARTIAL && ((g >= -.9972 && g <= .9972)
 341                             || (Math.abs(g) >= .9972 && Math.abs(g) < .9972 + Math.abs(u)))) {
 342                         jd = 0; // no partial sun eclipse
 343                     }
 344                     if (eclipse != EclipseKind.PARTIAL) {
 345                         if ((g < -.9972 || g > .9972) || (Math.abs(g) < .9972 && Math.abs(g) > .9972 + Math.abs(u))) {
 346                             jd = 0; // no ring or total sun eclipse
 347                         }
 348                         if (u > .0047 || u >= .00464 * Math.sqrt(1 - g * g)) {
 349                             ringTest = 1; // no total sun eclipse
 350                         }
 351                         if (ringTest == 1 && eclipse == EclipseKind.TOTAL) {
 352                             jd = 0;
 353                         }
 354                         if (ringTest == 0 && eclipse == EclipseKind.RING) {
 355                             jd = 0;
 356                         }
 357                     }
 358                     break;
 359             }
 360         }
 361         return jd;
 362     }
 363
 364     /**
 365      * Calculates the illumination.
 366      */
 367     private double getIllumination(double jd) {
 368         double t = (jd - 2451545) / 36525;
 369         double d = 297.8502042 + 445267.11151686 * t - .00163 * t * t + t * t * t / 545868 - t * t * t * t / 113065000;
 370         double m = 357.5291092 + 35999.0502909 * t - .0001536 * t * t + t * t * t / 24490000;
 371         double m1 = 134.9634114 + 477198.8676313 * t + .008997 * t * t + t * t * t / 69699 - t * t * t * t / 14712000;
 372         double i = 180 - d - 6.289 * SN(m1) + 2.1 * SN(m) - 1.274 * SN(2 * d - m1) - .658 * SN(2 * d)
 373                 - .241 * SN(2 * m1) - .110 * SN(d);
 374         return (1 + CS(i)) / 2 * 100.0;
 375     }
 376
 377     /**
 378      * Calculates the next moon phase.
 379      */
 380     private double getNextPhase(Calendar cal, double midnightJd, double mode) {
 381         double tz = 0;
 382         double phaseJd = 0;
 383         do {
 384             double k = var_k(cal, tz);
 385             tz += 1;
 386             phaseJd = calcMoonPhase(k, mode);
 387         } while (phaseJd <= midnightJd);
 388         return phaseJd;
 389     }
 390
 391     /**
 392      * Calculates the previous moon phase.
 393      */
 394     private double getPreviousPhase(Calendar cal, double jd, double mode) {
 395         double tz = 0;
 396         double phaseJd = 0;
 397         do {
 398             double k = var_k(cal, tz);
 399             tz -= 1;
 400             phaseJd = calcMoonPhase(k, mode);
 401         } while (phaseJd > jd);
 402         return phaseJd;
 403     }
 404
 405     /**
 406      * Calculates the next eclipse.
 407      */
 408     protected double getEclipse(Calendar cal, EclipseType type, double midnightJd, EclipseKind eclipse) {
 409         double tz = 0;
 410         double eclipseJd = 0;
 411         do {
 412             double k = var_k(cal, tz);
 413             tz += 1;
 414             eclipseJd = getEclipse(k, type, eclipse);
 415         } while (eclipseJd <= midnightJd);
 416         return eclipseJd;
 417     }
 418
 419     /**
 420      * Calculates the date, where the moon is furthest away from the earth.
 421      */
 422     private double getApogee(double julianDate, double decimalYear) {
 423         double k = Math.floor((decimalYear - 1999.97) * 13.2555) + .5;
 424         double jd = 0;
 425         do {
 426             double t = k / 1325.55;
 427             double d = 171.9179 + 335.9106046 * k - .010025 * t * t - .00001156 * t * t * t
 428                     + .000000055 * t * t * t * t;
 429             double m = 347.3477 + 27.1577721 * k - .0008323 * t * t - .000001 * t * t * t;
 430             double f = 316.6109 + 364.5287911 * k - .0125131 * t * t - .0000148 * t * t * t;
 431             jd = 2451534.6698 + 27.55454988 * k - .0006886 * t * t - .000001098 * t * t * t + .0000000052 * t * t
 432                     + .4392 * SN(2 * d) + .0684 * SN(4 * d) + (.0456 - .00011 * t) * SN(m)
 433                     + (.0426 - .00011 * t) * SN(2 * d - m) + .0212 * SN(2 * f);
 434             jd += -.0189 * SN(d) + .0144 * SN(6 * d) + .0113 * SN(4 * d - m) + .0047 * SN(2 * d + 2 * f)
 435                     + .0036 * SN(d + m) + .0035 * SN(8 * d) + .0034 * SN(6 * d - m) - .0034 * SN(2 * d - 2 * f)
 436                     + .0022 * SN(2 * d - 2 * m) - .0017 * SN(3 * d);
 437             jd += .0013 * SN(4 * d + 2 * f) + .0011 * SN(8 * d - m) + .001 * SN(4 * d - 2 * m) + .0009 * SN(10 * d)
 438                     + .0007 * SN(3 * d + m) + .0006 * SN(2 * m) + .0005 * SN(2 * d + m) + .0005 * SN(2 * d + 2 * m)
 439                     + .0004 * SN(6 * d + 2 * f);
 440             jd += .0004 * SN(6 * d - 2 * m) + .0004 * SN(10 * d - m) - .0004 * SN(5 * d) - .0004 * SN(4 * d - 2 * f)
 441                     + .0003 * SN(2 * f + m) + .0003 * SN(12 * d) + .0003 * SN(2 * d + 2 * f - m) - .0003 * SN(d - m);
 442             k += 1;
 443         } while (jd < julianDate);
 444         return jd;
 445     }
 446
 447     /**
 448      * Calculates the date, where the moon is closest to the earth.
 449      */
 450     private double getPerigee(double julianDate, double decimalYear) {
 451         double k = Math.floor((decimalYear - 1999.97) * 13.2555);
 452         double jd = 0;
 453         do {
 454             double t = k / 1325.55;
 455             double d = 171.9179 + 335.9106046 * k - .010025 * t * t - .00001156 * t * t * t
 456                     + .000000055 * t * t * t * t;
 457             double m = 347.3477 + 27.1577721 * k - .0008323 * t * t - .000001 * t * t * t;
 458             double f = 316.6109 + 364.5287911 * k - .0125131 * t * t - .0000148 * t * t * t;
 459             jd = 2451534.6698 + 27.55454988 * k - .0006886 * t * t - .000001098 * t * t * t + .0000000052 * t * t
 460                     - 1.6769 * SN(2 * d) + .4589 * SN(4 * d) - .1856 * SN(6 * d) + .0883 * SN(8 * d);
 461             jd += -(.0773 + .00019 * t) * SN(2 * d - m) + (.0502 - .00013 * t) * SN(m) - .046 * SN(10 * d)
 462                     + (.0422 - .00011 * t) * SN(4 * d - m) - .0256 * SN(6 * d - m) + .0253 * SN(12 * d) + .0237 * SN(d);
 463             jd += .0162 * SN(8 * d - m) - .0145 * SN(14 * d) + .0129 * SN(2 * f) - .0112 * SN(3 * d)
 464                     - .0104 * SN(10 * d - m) + .0086 * SN(16 * d) + .0069 * SN(12 * d - m) + .0066 * SN(5 * d)
 465                     - .0053 * SN(2 * d + 2 * f);
 466             jd += -.0052 * SN(18 * d) - .0046 * SN(14 * d - m) - .0041 * SN(7 * d) + .004 * SN(2 * d + m)
 467                     + .0032 * SN(20 * d) - .0032 * SN(d + m) + .0031 * SN(16 * d - m);
 468             jd += -.0029 * SN(4 * d + m) - .0027 * SN(2 * d - 2 * m) + .0024 * SN(4 * d - 2 * m)
 469                     - .0021 * SN(6 * d - 2 * m) - .0021 * SN(22 * d) - .0021 * SN(18 * d - m);
 470             jd += .0019 * SN(6 * d + m) - .0018 * SN(11 * d) - .0014 * SN(8 * d + m) - .0014 * SN(4 * d - 2 * f)
 471                     - .0014 * SN(6 * d - 2 * f) + .0014 * SN(3 * d + m) - .0014 * SN(5 * d + m) + .0013 * SN(13 * d);
 472             jd += .0013 * SN(20 * d - m) + .0011 * SN(3 * d + 2 * m) - .0011 * SN(4 * d + 2 * f - 2 * m)
 473                     - .001 * SN(d + 2 * m) - .0009 * SN(22 * d - m) - .0008 * SN(4 * f) + .0008 * SN(6 * d - 2 * f)
 474                     + .0008 * SN(2 * d - 2 * f + m);
 475             jd += .0007 * SN(2 * m) + .0007 * SN(2 * f - m) + .0007 * SN(2 * d + 4 * f) - .0006 * SN(2 * f - 2 * m)
 476                     - .0006 * SN(2 * d - 2 * f + 2 * m) + .0006 * SN(24 * d) + .0005 * SN(4 * d - 4 * f)
 477                     + .0005 * SN(2 * d + 2 * m) - .0004 * SN(d - m) + .0027 * SN(9 * d) + .0027 * SN(4 * d + 2 * f);
 478             k += 1;
 479         } while (jd < julianDate);
 480         return jd;
 481     }
 482
 483     /**
 484      * Calculates the distance from the moon to earth.
 485      */
 486     private double getDistance(double jd) {
 487         double t = (jd - 2451545) / 36525;
 488         double d = 297.8502042 + 445267.11151686 * t - .00163 * t * t + t * t * t / 545868 - t * t * t * t / 113065000;
 489         double m = 357.5291092 + 35999.0502909 * t - .0001536 * t * t + t * t * t / 24490000;
 490         double m1 = 134.9634114 + 477198.8676313 * t + .008997 * t * t + t * t * t / 69699 - t * t * t * t / 14712000;
 491         double f = 93.27209929999999 + 483202.0175273 * t - .0034029 * t * t - t * t * t / 3526000
 492                 + t * t * t * t / 863310000;
 493         return 385000.56 + getCoefficient(d, m, m1, f) / 1000;
 494     }
 495
 496     private double[] calcMoon(double t) {
 497         double p2 = 6.283185307;
 498         double arc = 206264.8062;
 499         double coseps = .91748;
 500         double sineps = .39778;
 501         double lo = FRAK(.606433 + 1336.855225 * t);
 502         double l = p2 * FRAK(.374897 + 1325.55241 * t);
 503         double ls = p2 * FRAK(.993133 + 99.997361 * t);
 504         double d = p2 * FRAK(.827361 + 1236.853086 * t);
 505         double f = p2 * FRAK(.259086 + 1342.227825 * t);
 506         double dl = 22640 * Math.sin(l) - 4586 * Math.sin(l - 2 * d) + 2370 * Math.sin(2 * d) + 769 * Math.sin(2 * l)
 507                 - 668 * Math.sin(ls) - 412 * Math.sin(2 * f) - 212 * Math.sin(2 * l - 2 * d)
 508                 - 206 * Math.sin(l + ls - 2 * d) + 192 * Math.sin(l + 2 * d) - 165 * Math.sin(ls - 2 * d)
 509                 - 125 * Math.sin(d) - 110 * Math.sin(l + ls) + 148 * Math.sin(l - ls) - 55 * Math.sin(2 * f - 2 * d);
 510         double s = f + (dl + 412 * Math.sin(2 * f) + 541 * Math.sin(ls)) / arc;
 511         double h = f - 2 * d;
 512         double n = -526 * Math.sin(h) + 44 * Math.sin(l + h) - 31 * Math.sin(-l + h) - 23 * Math.sin(ls + h)
 513                 + 11 * Math.sin(-ls + h) - 25 * Math.sin(-2 * l + f) + 21 * Math.sin(-l + f);
 514         double lmoon = p2 * FRAK(lo + dl / 1296000);
 515         double bmoon = (18520 * Math.sin(s) + n) / arc;
 516         double cb = Math.cos(bmoon);
 517         double x = cb * Math.cos(lmoon);
 518         double v = cb * Math.sin(lmoon);
 519         double w = Math.sin(bmoon);
 520         double y = coseps * v - sineps * w;
 521         double z = sineps * v + coseps * w;
 522         double rho = Math.sqrt(1 - z * z);
 523         double dec = (360 / p2) * Math.atan(z / rho);
 524         double ra = (48 / p2) * Math.atan(y / (x + rho));
 525         if (ra < 0) {
 526             ra += 24;
 527         }
 528         return new double[] { dec, ra };
 529     }
 530
 531     private double CS(double x) {
 532         return Math.cos(x * SunCalc.DEG2RAD);
 533     }
 534
 535     private double SN(double x) {
 536         return Math.sin(x * SunCalc.DEG2RAD);
 537     }
 538
 539     private double SINALT(double moonJd, int hour, double lambda, double cphi, double sphi) {
 540         double jdo = moonJd + hour / 24.0;
 541         double t = (jdo - 51544.5) / 36525.0;
 542         double[] decra = calcMoon(t);
 543         double tau = 15.0 * (LMST(jdo, lambda) - decra[1]);
 544         return sphi * SN(decra[0]) + cphi * CS(decra[0]) * CS(tau);
 545     }
 546
 547     private double LMST(double moonJd, double lambda) {
 548         double moonJdo = Math.floor(moonJd);
 549         double ut = (moonJd - moonJdo) * 24.0;
 550         double t = (moonJdo - 51544.5) / 36525.0;
 551         double gmst = 6.697374558 + 1.0027379093 * ut + (8640184.812866 + (.093104 - .0000062 * t) * t) * t / 3600.0;
 552         return 24.0 * FRAK((gmst - lambda / 15.0) / 24.0);
 553     }
 554
 555     private double FRAK(double x) {
 556         double ret = x - (int) (x);
 557         if (ret < 0) {
 558             ret += 1;
 559         }
 560         return ret;
 561     }
 562
 563     private double[] QUAD(double yminus, double yo, double yplus) {
 564         double nz = 0;
 565         double a = .5 * (yminus + yplus) - yo;
 566         double b = .5 * (yplus - yminus);
 567         double c = yo;
 568         double xe = -b / (2 * a);
 569         double ye = (a * xe + b) * xe + c;
 570         double dis = b * b - 4 * a * c;
 571         double zero1 = 0;
 572         double zero2 = 0;
 573         if (dis >= 0) {
 574             double dx = .5 * Math.sqrt(dis) / Math.abs(a);
 575             zero1 = xe - dx;
 576             zero2 = xe + dx;
 577             if (Math.abs(zero1) <= 1) {
 578                 nz += 1;
 579             }
 580             if (Math.abs(zero2) <= 1) {
 581                 nz += 1;
 582             }
 583             if (zero1 < -1) {
 584                 zero1 = zero2;
 585             }
 586         }
 587         return new double[] { ye, zero1, zero2, nz };
 588     }
 589
 590     private double var_o(double k, double t) {
 591         return 124.7746 - 1.5637558 * k + .0020691 * t * t + .00000215 * t * t * t;
 592     }
 593
 594     private double var_f(double k, double t) {
 595         return 160.7108 + 390.67050274 * k - .0016341 * t * t - .00000227 * t * t * t + .000000011 * t * t * t * t;
 596     }
 597
 598     private double var_m1(double k, double t) {
 599         return 201.5643 + 385.81693528 * k + .1017438 * t * t + .00001239 * t * t * t - .000000058 * t * t * t * t;
 600     }
 601
 602     private double var_m(double k, double t) {
 603         return 2.5534 + 29.10535669 * k - .0000218 * t * t - .00000011 * t * t * t;
 604     }
 605
 606     private double var_e(double t) {
 607         return 1 - .002516 * t - .0000074 * t * t;
 608     }
 609
 610     private double var_jde(double k, double t) {
 611         return 2451550.09765 + 29.530588853 * k + .0001337 * t * t - .00000015 * t * t * t
 612                 + .00000000073 * t * t * t * t;
 613     }
 614
 615     private double var_k(Calendar cal, double tz) {
 616         return (cal.get(Calendar.YEAR) + (cal.get(Calendar.DAY_OF_YEAR) + tz) / 365 - 2000) * 12.3685;
 617     }
 618
 619     private double moonCorrection(double jd, double t, double k) {
 620         double ret = jd;
 621         ret += .000325 * SN(299.77 + .107408 * k - .009173 * t * t) + .000165 * SN(251.88 + .016321 * k)
 622                 + .000164 * SN(251.83 + 26.651886 * k) + .000126 * SN(349.42 + 36.412478 * k)
 623                 + .00011 * SN(84.66 + 18.206239 * k);
 624         ret += .000062 * SN(141.74 + 53.303771 * k) + .00006 * SN(207.14 + 2.453732 * k)
 625                 + .000056 * SN(154.84 + 7.30686 * k) + .000047 * SN(34.52 + 27.261239 * k)
 626                 + .000042 * SN(207.19 + .121824 * k) + .00004 * SN(291.34 + 1.844379 * k);
 627         ret += .000037 * SN(161.72 + 24.198154 * k) + .000035 * SN(239.56 + 25.513099 * k)
 628                 + .000023 * SN(331.55 + 3.592518 * k);
 629         return ret;
 630     }
 631
 632     private double getCoefficient(double d, double m, double m1, double f) {
 633         int[] kd = new int[] { 0, 2, 2, 0, 0, 0, 2, 2, 2, 2, 0, 1, 0, 2, 0, 0, 4, 0, 4, 2, 2, 1, 1, 2, 2, 4, 2, 0, 2, 2,
 634                 1, 2, 0, 0, 2, 2, 2, 4, 0, 3, 2, 4, 0, 2, 2, 2, 4, 0, 4, 1, 2, 0, 1, 3, 4, 2, 0, 1, 2, 2 };
 635         int[] km = new int[] { 0, 0, 0, 0, 1, 0, 0, -1, 0, -1, 1, 0, 1, 0, 0, 0, 0, 0, 0, 1, 1, 0, 1, -1, 0, 0, 0, 1, 0,
 636                 -1, 0, -2, 1, 2, -2, 0, 0, -1, 0, 0, 1, -1, 2, 2, 1, -1, 0, 0, -1, 0, 1, 0, 1, 0, 0, -1, 2, 1, 0, 0 };
 637         int[] km1 = new int[] { 1, -1, 0, 2, 0, 0, -2, -1, 1, 0, -1, 0, 1, 0, 1, 1, -1, 3, -2, -1, 0, -1, 0, 1, 2, 0,
 638                 -3, -2, -1, -2, 1, 0, 2, 0, -1, 1, 0, -1, 2, -1, 1, -2, -1, -1, -2, 0, 1, 4, 0, -2, 0, 2, 1, -2, -3, 2,
 639                 1, -1, 3, -1 };
 640         int[] kf = new int[] { 0, 0, 0, 0, 0, 2, 0, 0, 0, 0, 0, 0, 0, -2, 2, -2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2,
 641                 0, 0, 0, 0, 0, 0, -2, 2, 0, 2, 0, 0, 0, 0, 0, 0, -2, 0, 0, 0, 0, -2, -2, 0, 0, 0, 0, 0, 0, 0, -2 };
 642         int[] kr = new int[] { -20905355, -3699111, -2955968, -569925, 48888, -3149, 246158, -152138, -170733, -204586,
 643                 -129620, 108743, 104755, 10321, 0, 79661, -34782, -23210, -21636, 24208, 30824, -8379, -16675, -12831,
 644                 -10445, -11650, 14403, -7003, 0, 10056, 6322, -9884, 5751, 0, -4950, 4130, 0, -3958, 0, 3258, 2616,
 645                 -1897, -2117, 2354, 0, 0, -1423, -1117, -1571, -1739, 0, -4421, 0, 0, 0, 0, 1165, 0, 0, 8752 };
 646         double sr = 0;
 647         for (int t = 0; t < 60; t++) {
 648             sr += kr[t] * CS(kd[t] * d + km[t] * m + km1[t] * m1 + kf[t] * f);
 649         }
 650         return sr;
 651     }
 652
 653     /**
 654      * Sets the azimuth, elevation and zodiac in the moon object.
 655      */
 656     private void setAzimuthElevationZodiac(double julianDate, double latitude, double longitude, Moon moon) {
 657         double lat = latitude * SunCalc.DEG2RAD;
 658         double lon = longitude * SunCalc.DEG2RAD;
 659
 660         double gmst = toGMST(julianDate);
 661         double lmst = toLMST(gmst, lon) * 15. * SunCalc.DEG2RAD;
 662
 663         double d = julianDate - 2447891.5;
 664         double anomalyMean = 360 * SunCalc.DEG2RAD / 365.242191 * d + 4.87650757829735 - 4.935239984568769;
 665         double nu = anomalyMean + 360.0 * SunCalc.DEG2RAD / Math.PI * 0.016713 * Math.sin(anomalyMean);
 666         double sunLon = mod2Pi(nu + 4.935239984568769);
 667
 668         double l0 = 318.351648 * SunCalc.DEG2RAD;
 669         double p0 = 36.340410 * SunCalc.DEG2RAD;
 670         double n0 = 318.510107 * SunCalc.DEG2RAD;
 671         double i = 5.145396 * SunCalc.DEG2RAD;
 672         double l = 13.1763966 * SunCalc.DEG2RAD * d + l0;
 673         double mMoon = l - 0.1114041 * SunCalc.DEG2RAD * d - p0;
 674         double n = n0 - 0.0529539 * SunCalc.DEG2RAD * d;
 675         double c = l - sunLon;
 676         double ev = 1.2739 * SunCalc.DEG2RAD * Math.sin(2 * c - mMoon);
 677         double ae = 0.1858 * SunCalc.DEG2RAD * Math.sin(anomalyMean);
 678         double a3 = 0.37 * SunCalc.DEG2RAD * Math.sin(anomalyMean);
 679         double mMoon2 = mMoon + ev - ae - a3;
 680         double ec = 6.2886 * SunCalc.DEG2RAD * Math.sin(mMoon2);
 681         double a4 = 0.214 * SunCalc.DEG2RAD * Math.sin(2 * mMoon2);
 682         double l2 = l + ev + ec - ae + a4;
 683         double v = 0.6583 * SunCalc.DEG2RAD * Math.sin(2 * (l2 - sunLon));
 684         double l3 = l2 + v;
 685         double n2 = n - 0.16 * SunCalc.DEG2RAD * Math.sin(anomalyMean);
 686
 687         double moonLon = mod2Pi(n2 + Math.atan2(Math.sin(l3 - n2) * Math.cos(i), Math.cos(l3 - n2)));
 688         double moonLat = Math.asin(Math.sin(l3 - n2) * Math.sin(i));
 689
 690         double[] raDec = ecl2Equ(moonLat, moonLon, julianDate);
 691
 692         double distance = (1 - 0.00301401) / (1 + 0.054900 * Math.cos(mMoon2 + ec)) * 384401;
 693
 694         double[] raDecTopo = geoEqu2TopoEqu(raDec, distance, lat, lmst);
 695         double[] azAlt = equ2AzAlt(raDecTopo[0], raDecTopo[1], lat, lmst);
 696
 697         Position position = moon.getPosition();
 698         position.setAzimuth(azAlt[0] * SunCalc.RAD2DEG);
 699         position.setElevation(azAlt[1] * SunCalc.RAD2DEG + refraction(azAlt[1]));
 700
 701         // zodiac
 702         double idxd = Math.floor(moonLon * SunCalc.RAD2DEG / 30);
 703         int idx = 0;
 704         if (idxd < 0) {
 705             idx = (int) (Math.ceil(idxd));
 706         } else {
 707             idx = (int) (Math.floor(idxd));
 708         }
 709
 710         if (idx >= 0 || idx <= ZodiacSign.values().length) {
 711             moon.setZodiac(new Zodiac(ZodiacSign.values()[idx]));
 712         }
 713     }
 714
 715     private double mod2Pi(double x) {
 716         return (mod(x, 2. * Math.PI));
 717     }
 718
 719     private double mod(double a, double b) {
 720         return (a - Math.floor(a / b) * b);
 721     }
 722
 723     /**
 724      * Transform equatorial coordinates (ra/dec) to horizonal coordinates
 725      * (azimuth/altitude).
 726      */
 727     private double[] equ2AzAlt(double ra, double dec, double geolat, double lmst) {
 728         double cosdec = Math.cos(dec);
 729         double sindec = Math.sin(dec);
 730         double lha = lmst - ra;
 731         double coslha = Math.cos(lha);
 732         double sinlha = Math.sin(lha);
 733         double coslat = Math.cos(geolat);
 734         double sinlat = Math.sin(geolat);
 735
 736         double n = -cosdec * sinlha;
 737         double d = sindec * coslat - cosdec * coslha * sinlat;
 738         double az = mod2Pi(Math.atan2(n, d));
 739         double alt = Math.asin(sindec * sinlat + cosdec * coslha * coslat);
 740
 741         return new double[] { az, alt };
 742     }
 743
 744     /**
 745      * Transform ecliptical coordinates (lon/lat) to equatorial coordinates
 746      * (ra/dec)
 747      */
 748     private double[] ecl2Equ(double lat, double lon, double jd) {
 749         double t = (jd - 2451545.0) / 36525.0;
 750         double eps = (23. + (26 + 21.45 / 60.) / 60. + t * (-46.815 + t * (-0.0006 + t * 0.00181)) / 3600.)
 751                 * SunCalc.DEG2RAD;
 752         double coseps = Math.cos(eps);
 753         double sineps = Math.sin(eps);
 754
 755         double sinlon = Math.sin(lon);
 756         double ra = mod2Pi(Math.atan2((sinlon * coseps - Math.tan(lat) * sineps), Math.cos(lon)));
 757         double dec = Math.asin(Math.sin(lat) * coseps + Math.cos(lat) * sineps * sinlon);
 758
 759         return new double[] { ra, dec };
 760     }
 761
 762     /**
 763      * Transform geocentric equatorial coordinates (rA/dec) to topocentric
 764      * equatorial coordinates.
 765      */
 766     private double[] geoEqu2TopoEqu(double[] raDec, double distance, double observerLat, double lmst) {
 767         double cosdec = Math.cos(raDec[1]);
 768         double sindec = Math.sin(raDec[1]);
 769         double coslst = Math.cos(lmst);
 770         double sinlst = Math.sin(lmst);
 771         double coslat = Math.cos(observerLat);
 772         double sinlat = Math.sin(observerLat);
 773         double rho = getCenterDistance(observerLat);
 774
 775         double x = distance * cosdec * Math.cos(raDec[0]) - rho * coslat * coslst;
 776         double y = distance * cosdec * Math.sin(raDec[0]) - rho * coslat * sinlst;
 777         double z = distance * sindec - rho * sinlat;
 778
 779         double distanceTopocentric = Math.sqrt(x * x + y * y + z * z);
 780         double raTopo = mod2Pi(Math.atan2(y, x));
 781         double decTopo = Math.asin(z / distanceTopocentric);
 782
 783         return new double[] { raTopo, decTopo };
 784     }
 785
 786     /**
 787      * Convert julian date to greenwich mean sidereal time.
 788      */
 789     private double toGMST(double jd) {
 790         double ut = (jd - 0.5 - Math.floor(jd - 0.5)) * 24.;
 791         double jdMod = Math.floor(jd - 0.5) + 0.5;
 792         double t = (jdMod - 2451545.0) / 36525.0;
 793         double t0 = 6.697374558 + t * (2400.051336 + t * 0.000025862);
 794         return (mod(t0 + ut * 1.002737909, 24.));
 795     }
 796
 797     /**
 798      * Convert greenwich mean sidereal time to local mean sidereal time.
 799      */
 800     private double toLMST(double gmst, double lon) {
 801         return mod(gmst + SunCalc.RAD2DEG * lon / 15., 24.);
 802     }
 803
 804     /**
 805      * Returns geocentric distance from earth center.
 806      */
 807     private double getCenterDistance(double lat) {
 808         double co = Math.cos(lat);
 809         double si = Math.sin(lat);
 810         double fl = 1.0 - 1.0 / 298.257223563;
 811         fl = fl * fl;
 812         si = si * si;
 813         double u = 1.0 / Math.sqrt(co * co + fl * si);
 814         double a = 6378.137 * u;
 815         double b = 6378.137 * fl * u;
 816         return Math.sqrt(a * a * co * co + b * b * si);
 817     }
 818
 819     /**
 820      * Returns altitude increase in altitude in degrees. Rough refraction
 821      * formula using standard atmosphere: 1015 mbar and 10°C.
 822      */
 823     private double refraction(double alt) {
 824         int pressure = 1015;
 825         int temperature = 10;
 826         double altdeg = alt * SunCalc.RAD2DEG;
 827
 828         if (altdeg < -2 || altdeg >= 90) {
 829             return 0;
 830         }
 831
 832         if (altdeg > 15) {
 833             return 0.00452 * pressure / ((273 + temperature) * Math.tan(alt));
 834         }
 835
 836         double y = alt;
 837         double d = 0.0;
 838         double p = (pressure - 80.0) / 930.0;
 839         double q = 0.0048 * (temperature - 10.0);
 840         double y0 = y;
 841         double d0 = d;
 842         double n = 0.0;
 843
 844         for (int i = 0; i < 3; i++) {
 845             n = y + (7.31 / (y + 4.4));
 846             n = 1.0 / Math.tan(n * SunCalc.DEG2RAD);
 847             d = n * p / (60.0 + q * (n + 39.0));
 848             n = y - y0;
 849             y0 = d - d0 - n;
 850             n = ((n != 0.0) && (y0 != 0.0)) ? y - n * (alt + d - y) / y0 : alt + d;
 851             y0 = y;
 852             d0 = d;
 853             y = n;
 854         }
 855         return d;
 856     }
 857 }