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  * @see 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         double sr = 385000.56 + getCoefficient(d, m, m1, f) / 1000;
 494         return sr;
 495     }
 496
 497     private double[] calcMoon(double t) {
 498         double p2 = 6.283185307;
 499         double arc = 206264.8062;
 500         double coseps = .91748;
 501         double sineps = .39778;
 502         double lo = FRAK(.606433 + 1336.855225 * t);
 503         double l = p2 * FRAK(.374897 + 1325.55241 * t);
 504         double ls = p2 * FRAK(.993133 + 99.997361 * t);
 505         double d = p2 * FRAK(.827361 + 1236.853086 * t);
 506         double f = p2 * FRAK(.259086 + 1342.227825 * t);
 507         double dl = 22640 * Math.sin(l) - 4586 * Math.sin(l - 2 * d) + 2370 * Math.sin(2 * d) + 769 * Math.sin(2 * l)
 508                 - 668 * Math.sin(ls) - 412 * Math.sin(2 * f) - 212 * Math.sin(2 * l - 2 * d)
 509                 - 206 * Math.sin(l + ls - 2 * d) + 192 * Math.sin(l + 2 * d) - 165 * Math.sin(ls - 2 * d)
 510                 - 125 * Math.sin(d) - 110 * Math.sin(l + ls) + 148 * Math.sin(l - ls) - 55 * Math.sin(2 * f - 2 * d);
 511         double s = f + (dl + 412 * Math.sin(2 * f) + 541 * Math.sin(ls)) / arc;
 512         double h = f - 2 * d;
 513         double n = -526 * Math.sin(h) + 44 * Math.sin(l + h) - 31 * Math.sin(-l + h) - 23 * Math.sin(ls + h)
 514                 + 11 * Math.sin(-ls + h) - 25 * Math.sin(-2 * l + f) + 21 * Math.sin(-l + f);
 515         double lmoon = p2 * FRAK(lo + dl / 1296000);
 516         double bmoon = (18520 * Math.sin(s) + n) / arc;
 517         double cb = Math.cos(bmoon);
 518         double x = cb * Math.cos(lmoon);
 519         double v = cb * Math.sin(lmoon);
 520         double w = Math.sin(bmoon);
 521         double y = coseps * v - sineps * w;
 522         double z = sineps * v + coseps * w;
 523         double rho = Math.sqrt(1 - z * z);
 524         double dec = (360 / p2) * Math.atan(z / rho);
 525         double ra = (48 / p2) * Math.atan(y / (x + rho));
 526         if (ra < 0) {
 527             ra += 24;
 528         }
 529         return new double[] { dec, ra };
 530     }
 531
 532     private double CS(double x) {
 533         return Math.cos(x * SunCalc.DEG2RAD);
 534     }
 535
 536     private double SN(double x) {
 537         return Math.sin(x * SunCalc.DEG2RAD);
 538     }
 539
 540     private double SINALT(double moonJd, int hour, double lambda, double cphi, double sphi) {
 541         double jdo = moonJd + hour / 24.0;
 542         double t = (jdo - 51544.5) / 36525.0;
 543         double decra[] = calcMoon(t);
 544         double tau = 15.0 * (LMST(jdo, lambda) - decra[1]);
 545         return sphi * SN(decra[0]) + cphi * CS(decra[0]) * CS(tau);
 546     }
 547
 548     private double LMST(double moonJd, double lambda) {
 549         double moonJdo = Math.floor(moonJd);
 550         double ut = (moonJd - moonJdo) * 24.0;
 551         double t = (moonJdo - 51544.5) / 36525.0;
 552         double gmst = 6.697374558 + 1.0027379093 * ut + (8640184.812866 + (.093104 - .0000062 * t) * t) * t / 3600.0;
 553         return 24.0 * FRAK((gmst - lambda / 15.0) / 24.0);
 554     }
 555
 556     private double FRAK(double x) {
 557         double ret = x - (int) (x);
 558         if (ret < 0) {
 559             ret += 1;
 560         }
 561         return ret;
 562     }
 563
 564     private double[] QUAD(double yminus, double yo, double yplus) {
 565         double nz = 0;
 566         double a = .5 * (yminus + yplus) - yo;
 567         double b = .5 * (yplus - yminus);
 568         double c = yo;
 569         double xe = -b / (2 * a);
 570         double ye = (a * xe + b) * xe + c;
 571         double dis = b * b - 4 * a * c;
 572         double zero1 = 0;
 573         double zero2 = 0;
 574         if (dis >= 0) {
 575             double dx = .5 * Math.sqrt(dis) / Math.abs(a);
 576             zero1 = xe - dx;
 577             zero2 = xe + dx;
 578             if (Math.abs(zero1) <= 1) {
 579                 nz += 1;
 580             }
 581             if (Math.abs(zero2) <= 1) {
 582                 nz += 1;
 583             }
 584             if (zero1 < -1) {
 585                 zero1 = zero2;
 586             }
 587         }
 588         return new double[] { ye, zero1, zero2, nz };
 589     }
 590
 591     private double var_o(double k, double t) {
 592         return 124.7746 - 1.5637558 * k + .0020691 * t * t + .00000215 * t * t * t;
 593     }
 594
 595     private double var_f(double k, double t) {
 596         return 160.7108 + 390.67050274 * k - .0016341 * t * t - .00000227 * t * t * t + .000000011 * t * t * t * t;
 597     }
 598
 599     private double var_m1(double k, double t) {
 600         return 201.5643 + 385.81693528 * k + .1017438 * t * t + .00001239 * t * t * t - .000000058 * t * t * t * t;
 601     }
 602
 603     private double var_m(double k, double t) {
 604         return 2.5534 + 29.10535669 * k - .0000218 * t * t - .00000011 * t * t * t;
 605     }
 606
 607     private double var_e(double t) {
 608         return 1 - .002516 * t - .0000074 * t * t;
 609     }
 610
 611     private double var_jde(double k, double t) {
 612         return 2451550.09765 + 29.530588853 * k + .0001337 * t * t - .00000015 * t * t * t
 613                 + .00000000073 * t * t * t * t;
 614     }
 615
 616     private double var_k(Calendar cal, double tz) {
 617         return (cal.get(Calendar.YEAR) + (cal.get(Calendar.DAY_OF_YEAR) + tz) / 365 - 2000) * 12.3685;
 618     }
 619
 620     private double moonCorrection(double jd, double t, double k) {
 621         double ret = jd;
 622         ret += .000325 * SN(299.77 + .107408 * k - .009173 * t * t) + .000165 * SN(251.88 + .016321 * k)
 623                 + .000164 * SN(251.83 + 26.651886 * k) + .000126 * SN(349.42 + 36.412478 * k)
 624                 + .00011 * SN(84.66 + 18.206239 * k);
 625         ret += .000062 * SN(141.74 + 53.303771 * k) + .00006 * SN(207.14 + 2.453732 * k)
 626                 + .000056 * SN(154.84 + 7.30686 * k) + .000047 * SN(34.52 + 27.261239 * k)
 627                 + .000042 * SN(207.19 + .121824 * k) + .00004 * SN(291.34 + 1.844379 * k);
 628         ret += .000037 * SN(161.72 + 24.198154 * k) + .000035 * SN(239.56 + 25.513099 * k)
 629                 + .000023 * SN(331.55 + 3.592518 * k);
 630         return ret;
 631     }
 632
 633     private double getCoefficient(double d, double m, double m1, double f) {
 634         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,
 635                 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 };
 636         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,
 637                 -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 };
 638         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,
 639                 -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,
 640                 1, -1, 3, -1 };
 641         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,
 642                 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 };
 643         int[] kr = new int[] { -20905355, -3699111, -2955968, -569925, 48888, -3149, 246158, -152138, -170733, -204586,
 644                 -129620, 108743, 104755, 10321, 0, 79661, -34782, -23210, -21636, 24208, 30824, -8379, -16675, -12831,
 645                 -10445, -11650, 14403, -7003, 0, 10056, 6322, -9884, 5751, 0, -4950, 4130, 0, -3958, 0, 3258, 2616,
 646                 -1897, -2117, 2354, 0, 0, -1423, -1117, -1571, -1739, 0, -4421, 0, 0, 0, 0, 1165, 0, 0, 8752 };
 647         double sr = 0;
 648         for (int t = 0; t < 60; t++) {
 649             sr += kr[t] * CS(kd[t] * d + km[t] * m + km1[t] * m1 + kf[t] * f);
 650         }
 651         return sr;
 652     }
 653
 654     /**
 655      * Sets the azimuth, elevation and zodiac in the moon object.
 656      */
 657     private void setAzimuthElevationZodiac(double julianDate, double latitude, double longitude, Moon moon) {
 658         double lat = latitude * SunCalc.DEG2RAD;
 659         double lon = longitude * SunCalc.DEG2RAD;
 660
 661         double gmst = toGMST(julianDate);
 662         double lmst = toLMST(gmst, lon) * 15. * SunCalc.DEG2RAD;
 663
 664         double d = julianDate - 2447891.5;
 665         double anomalyMean = 360 * SunCalc.DEG2RAD / 365.242191 * d + 4.87650757829735 - 4.935239984568769;
 666         double nu = anomalyMean + 360.0 * SunCalc.DEG2RAD / Math.PI * 0.016713 * Math.sin(anomalyMean);
 667         double sunLon = mod2Pi(nu + 4.935239984568769);
 668
 669         double l0 = 318.351648 * SunCalc.DEG2RAD;
 670         double p0 = 36.340410 * SunCalc.DEG2RAD;
 671         double n0 = 318.510107 * SunCalc.DEG2RAD;
 672         double i = 5.145396 * SunCalc.DEG2RAD;
 673         double l = 13.1763966 * SunCalc.DEG2RAD * d + l0;
 674         double mMoon = l - 0.1114041 * SunCalc.DEG2RAD * d - p0;
 675         double n = n0 - 0.0529539 * SunCalc.DEG2RAD * d;
 676         double c = l - sunLon;
 677         double ev = 1.2739 * SunCalc.DEG2RAD * Math.sin(2 * c - mMoon);
 678         double ae = 0.1858 * SunCalc.DEG2RAD * Math.sin(anomalyMean);
 679         double a3 = 0.37 * SunCalc.DEG2RAD * Math.sin(anomalyMean);
 680         double mMoon2 = mMoon + ev - ae - a3;
 681         double ec = 6.2886 * SunCalc.DEG2RAD * Math.sin(mMoon2);
 682         double a4 = 0.214 * SunCalc.DEG2RAD * Math.sin(2 * mMoon2);
 683         double l2 = l + ev + ec - ae + a4;
 684         double v = 0.6583 * SunCalc.DEG2RAD * Math.sin(2 * (l2 - sunLon));
 685         double l3 = l2 + v;
 686         double n2 = n - 0.16 * SunCalc.DEG2RAD * Math.sin(anomalyMean);
 687
 688         double moonLon = mod2Pi(n2 + Math.atan2(Math.sin(l3 - n2) * Math.cos(i), Math.cos(l3 - n2)));
 689         double moonLat = Math.asin(Math.sin(l3 - n2) * Math.sin(i));
 690
 691         double raDec[] = ecl2Equ(moonLat, moonLon, julianDate);
 692
 693         double distance = (1 - 0.00301401) / (1 + 0.054900 * Math.cos(mMoon2 + ec)) * 384401;
 694
 695         double raDecTopo[] = geoEqu2TopoEqu(raDec, distance, lat, lmst);
 696         double azAlt[] = equ2AzAlt(raDecTopo[0], raDecTopo[1], lat, lmst);
 697
 698         Position position = moon.getPosition();
 699         position.setAzimuth(azAlt[0] * SunCalc.RAD2DEG);
 700         position.setElevation(azAlt[1] * SunCalc.RAD2DEG + refraction(azAlt[1]));
 701
 702         // zodiac
 703         double idxd = Math.floor(moonLon * SunCalc.RAD2DEG / 30);
 704         int idx = 0;
 705         if (idxd < 0) {
 706             idx = (int) (Math.ceil(idxd));
 707         } else {
 708             idx = (int) (Math.floor(idxd));
 709         }
 710
 711         if (idx >= 0 || idx <= ZodiacSign.values().length) {
 712             moon.setZodiac(new Zodiac(ZodiacSign.values()[idx]));
 713         }
 714     }
 715
 716     private double mod2Pi(double x) {
 717         return (mod(x, 2. * Math.PI));
 718     }
 719
 720     private double mod(double a, double b) {
 721         return (a - Math.floor(a / b) * b);
 722     }
 723
 724     /**
 725      * Transform equatorial coordinates (ra/dec) to horizonal coordinates
 726      * (azimuth/altitude).
 727      */
 728     private double[] equ2AzAlt(double ra, double dec, double geolat, double lmst) {
 729         double cosdec = Math.cos(dec);
 730         double sindec = Math.sin(dec);
 731         double lha = lmst - ra;
 732         double coslha = Math.cos(lha);
 733         double sinlha = Math.sin(lha);
 734         double coslat = Math.cos(geolat);
 735         double sinlat = Math.sin(geolat);
 736
 737         double n = -cosdec * sinlha;
 738         double d = sindec * coslat - cosdec * coslha * sinlat;
 739         double az = mod2Pi(Math.atan2(n, d));
 740         double alt = Math.asin(sindec * sinlat + cosdec * coslha * coslat);
 741
 742         return new double[] { az, alt };
 743     }
 744
 745     /**
 746      * Transform ecliptical coordinates (lon/lat) to equatorial coordinates
 747      * (ra/dec)
 748      */
 749     private double[] ecl2Equ(double lat, double lon, double jd) {
 750         double t = (jd - 2451545.0) / 36525.0;
 751         double eps = (23. + (26 + 21.45 / 60.) / 60. + t * (-46.815 + t * (-0.0006 + t * 0.00181)) / 3600.)
 752                 * SunCalc.DEG2RAD;
 753         double coseps = Math.cos(eps);
 754         double sineps = Math.sin(eps);
 755
 756         double sinlon = Math.sin(lon);
 757         double ra = mod2Pi(Math.atan2((sinlon * coseps - Math.tan(lat) * sineps), Math.cos(lon)));
 758         double dec = Math.asin(Math.sin(lat) * coseps + Math.cos(lat) * sineps * sinlon);
 759
 760         return new double[] { ra, dec };
 761     }
 762
 763     /**
 764      * Transform geocentric equatorial coordinates (rA/dec) to topocentric
 765      * equatorial coordinates.
 766      */
 767     private double[] geoEqu2TopoEqu(double[] raDec, double distance, double observerLat, double lmst) {
 768         double cosdec = Math.cos(raDec[1]);
 769         double sindec = Math.sin(raDec[1]);
 770         double coslst = Math.cos(lmst);
 771         double sinlst = Math.sin(lmst);
 772         double coslat = Math.cos(observerLat);
 773         double sinlat = Math.sin(observerLat);
 774         double rho = getCenterDistance(observerLat);
 775
 776         double x = distance * cosdec * Math.cos(raDec[0]) - rho * coslat * coslst;
 777         double y = distance * cosdec * Math.sin(raDec[0]) - rho * coslat * sinlst;
 778         double z = distance * sindec - rho * sinlat;
 779
 780         double distanceTopocentric = Math.sqrt(x * x + y * y + z * z);
 781         double raTopo = mod2Pi(Math.atan2(y, x));
 782         double decTopo = Math.asin(z / distanceTopocentric);
 783
 784         return new double[] { raTopo, decTopo };
 785     }
 786
 787     /**
 788      * Convert julian date to greenwich mean sidereal time.
 789      */
 790     private double toGMST(double jd) {
 791         double ut = (jd - 0.5 - Math.floor(jd - 0.5)) * 24.;
 792         double jdMod = Math.floor(jd - 0.5) + 0.5;
 793         double t = (jdMod - 2451545.0) / 36525.0;
 794         double t0 = 6.697374558 + t * (2400.051336 + t * 0.000025862);
 795         return (mod(t0 + ut * 1.002737909, 24.));
 796     }
 797
 798     /**
 799      * Convert greenwich mean sidereal time to local mean sidereal time.
 800      */
 801     private double toLMST(double gmst, double lon) {
 802         return mod(gmst + SunCalc.RAD2DEG * lon / 15., 24.);
 803     }
 804
 805     /**
 806      * Returns geocentric distance from earth center.
 807      */
 808     private double getCenterDistance(double lat) {
 809         double co = Math.cos(lat);
 810         double si = Math.sin(lat);
 811         double fl = 1.0 - 1.0 / 298.257223563;
 812         fl = fl * fl;
 813         si = si * si;
 814         double u = 1.0 / Math.sqrt(co * co + fl * si);
 815         double a = 6378.137 * u;
 816         double b = 6378.137 * fl * u;
 817         return Math.sqrt(a * a * co * co + b * b * si);
 818     }
 819
 820     /**
 821      * Returns altitude increase in altitude in degrees. Rough refraction
 822      * formula using standard atmosphere: 1015 mbar and 10°C.
 823      */
 824     private double refraction(double alt) {
 825         int pressure = 1015;
 826         int temperature = 10;
 827         double altdeg = alt * SunCalc.RAD2DEG;
 828
 829         if (altdeg < -2 || altdeg >= 90) {
 830             return 0;
 831         }
 832
 833         if (altdeg > 15) {
 834             return 0.00452 * pressure / ((273 + temperature) * Math.tan(alt));
 835         }
 836
 837         double y = alt;
 838         double d = 0.0;
 839         double p = (pressure - 80.0) / 930.0;
 840         double q = 0.0048 * (temperature - 10.0);
 841         double y0 = y;
 842         double d0 = d;
 843         double n = 0.0;
 844
 845         for (int i = 0; i < 3; i++) {
 846             n = y + (7.31 / (y + 4.4));
 847             n = 1.0 / Math.tan(n * SunCalc.DEG2RAD);
 848             d = n * p / (60.0 + q * (n + 39.0));
 849             n = y - y0;
 850             y0 = d - d0 - n;
 851             n = ((n != 0.0) && (y0 != 0.0)) ? y - n * (alt + d - y) / y0 : alt + d;
 852             y0 = y;
 853             d0 = d;
 854             y = n;
 855         }
 856         return d;
 857     }
 858 }