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

   1 /**
   2  * Copyright (c) 2010-2020 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.util.Calendar;
  17
  18 import org.openhab.binding.astro.internal.model.Eclipse;
  19 import org.openhab.binding.astro.internal.model.EclipseKind;
  20 import org.openhab.binding.astro.internal.model.EclipseType;
  21 import org.openhab.binding.astro.internal.model.Moon;
  22 import org.openhab.binding.astro.internal.model.MoonDistance;
  23 import org.openhab.binding.astro.internal.model.MoonPhase;
  24 import org.openhab.binding.astro.internal.model.MoonPhaseName;
  25 import org.openhab.binding.astro.internal.model.Position;
  26 import org.openhab.binding.astro.internal.model.Range;
  27 import org.openhab.binding.astro.internal.model.Zodiac;
  28 import org.openhab.binding.astro.internal.model.ZodiacSign;
  29 import org.openhab.binding.astro.internal.util.DateTimeUtils;
  30
  31 /**
  32  * Calculates the phase, eclipse, rise, set, distance, illumination and age of
  33  * the moon.
  34  *
  35  * @author Gerhard Riegler - Initial contribution
  36  * @author Christoph Weitkamp - Introduced UoM
  37  * @see based on the calculations of
  38  *      http://www.computus.de/mondphase/mondphase.htm azimuth/elevation and
  39  *      zodiac based on http://lexikon.astronomie.info/java/sunmoon/
  40  */
  41 public class MoonCalc {
  42     private static final double NEW_MOON = 0;
  43     private static final double FULL_MOON = 0.5;
  44     private static final double FIRST_QUARTER = 0.25;
  45     private static final double LAST_QUARTER = 0.75;
  46
  47     /**
  48      * Calculates all moon data at the specified coordinates
  49      */
  50     public Moon getMoonInfo(Calendar calendar, double latitude, double longitude) {
  51         Moon moon = new Moon();
  52
  53         double julianDate = DateTimeUtils.dateToJulianDate(calendar);
  54         double julianDateMidnight = DateTimeUtils.midnightDateToJulianDate(calendar);
  55
  56         double[] riseSet = getRiseSet(calendar, latitude, longitude);
  57         Calendar rise = DateTimeUtils.timeToCalendar(calendar, riseSet[0]);
  58         Calendar set = DateTimeUtils.timeToCalendar(calendar, riseSet[1]);
  59
  60         if (rise == null || set == null) {
  61             Calendar tomorrow = (Calendar) calendar.clone();
  62             tomorrow.add(Calendar.DAY_OF_MONTH, 1);
  63
  64             double[] riseSeTomorrow = getRiseSet(tomorrow, latitude, longitude);
  65             if (rise == null) {
  66                 rise = DateTimeUtils.timeToCalendar(tomorrow, riseSeTomorrow[0]);
  67             }
  68             if (set == null) {
  69                 set = DateTimeUtils.timeToCalendar(tomorrow, riseSeTomorrow[1]);
  70             }
  71         }
  72
  73         moon.setRise(new Range(rise, rise));
  74         moon.setSet(new Range(set, set));
  75
  76         MoonPhase phase = moon.getPhase();
  77         phase.setNew(DateTimeUtils.toCalendar(getNextPhase(calendar, julianDateMidnight, NEW_MOON)));
  78         phase.setFirstQuarter(DateTimeUtils.toCalendar(getNextPhase(calendar, julianDateMidnight, FIRST_QUARTER)));
  79         phase.setFull(DateTimeUtils.toCalendar(getNextPhase(calendar, julianDateMidnight, FULL_MOON)));
  80         phase.setThirdQuarter(DateTimeUtils.toCalendar(getNextPhase(calendar, julianDateMidnight, LAST_QUARTER)));
  81
  82         Eclipse eclipse = moon.getEclipse();
  83         eclipse.getKinds().forEach(eclipseKind -> {
  84             double jdate = getEclipse(calendar, EclipseType.MOON, julianDateMidnight, eclipseKind);
  85             eclipse.set(eclipseKind, DateTimeUtils.toCalendar(jdate), new Position());
  86         });
  87
  88         double decimalYear = DateTimeUtils.getDecimalYear(calendar);
  89         MoonDistance apogee = moon.getApogee();
  90         double apogeeJd = getApogee(julianDate, decimalYear);
  91         apogee.setDate(DateTimeUtils.toCalendar(apogeeJd));
  92         apogee.setDistance(getDistance(apogeeJd));
  93
  94         MoonDistance perigee = moon.getPerigee();
  95         double perigeeJd = getPerigee(julianDate, decimalYear);
  96         perigee.setDate(DateTimeUtils.toCalendar(perigeeJd));
  97         perigee.setDistance(getDistance(perigeeJd));
  98
  99         return moon;
 100     }
 101
 102     /**
 103      * Calculates the moon illumination and distance.
 104      */
 105     public void setPositionalInfo(Calendar calendar, double latitude, double longitude, Moon moon) {
 106         double julianDate = DateTimeUtils.dateToJulianDate(calendar);
 107         setMoonPhase(calendar, moon);
 108         setAzimuthElevationZodiac(julianDate, latitude, longitude, moon);
 109
 110         MoonDistance distance = moon.getDistance();
 111         distance.setDate(Calendar.getInstance());
 112         distance.setDistance(getDistance(julianDate));
 113     }
 114
 115     /**
 116      * Calculates the age and the current phase.
 117      */
 118     private void setMoonPhase(Calendar calendar, Moon moon) {
 119         MoonPhase phase = moon.getPhase();
 120         double julianDateEndOfDay = DateTimeUtils.endOfDayDateToJulianDate(calendar);
 121         double parentNewMoon = getPreviousPhase(calendar, julianDateEndOfDay, NEW_MOON);
 122         double age = Math.abs(parentNewMoon - julianDateEndOfDay);
 123         phase.setAge((int) age);
 124
 125         long parentNewMoonMillis = DateTimeUtils.toCalendar(parentNewMoon).getTimeInMillis();
 126         long ageRangeTimeMillis = phase.getNew().getTimeInMillis() - parentNewMoonMillis;
 127         long ageCurrentMillis = System.currentTimeMillis() - parentNewMoonMillis;
 128         double agePercent = ageRangeTimeMillis != 0 ? ageCurrentMillis * 100.0 / ageRangeTimeMillis : 0;
 129         phase.setAgePercent(agePercent);
 130         phase.setAgeDegree(3.6 * agePercent);
 131         double illumination = getIllumination(DateTimeUtils.dateToJulianDate(calendar));
 132         phase.setIllumination(illumination);
 133         boolean isWaxing = age < (29.530588853 / 2);
 134         if (DateTimeUtils.isSameDay(calendar, phase.getNew())) {
 135             phase.setName(MoonPhaseName.NEW);
 136         } else if (DateTimeUtils.isSameDay(calendar, phase.getFirstQuarter())) {
 137             phase.setName(MoonPhaseName.FIRST_QUARTER);
 138         } else if (DateTimeUtils.isSameDay(calendar, phase.getThirdQuarter())) {
 139             phase.setName(MoonPhaseName.THIRD_QUARTER);
 140         } else if (DateTimeUtils.isSameDay(calendar, phase.getFull())) {
 141             phase.setName(MoonPhaseName.FULL);
 142         } else if (illumination >= 0 && illumination < 50) {
 143             phase.setName(isWaxing ? MoonPhaseName.WAXING_CRESCENT : MoonPhaseName.WANING_CRESCENT);
 144         } else if (illumination >= 50 && illumination < 100) {
 145             phase.setName(isWaxing ? MoonPhaseName.WAXING_GIBBOUS : MoonPhaseName.WANING_GIBBOUS);
 146         }
 147     }
 148
 149     /**
 150      * Calculates moonrise and moonset.
 151      */
 152     private double[] getRiseSet(Calendar calendar, double latitude, double longitude) {
 153         double lambda = prepareCoordinate(longitude, 180);
 154         if (longitude > 0) {
 155             lambda *= -1;
 156         }
 157         double phi = prepareCoordinate(latitude, 90);
 158         if (latitude < 0) {
 159             phi *= -1;
 160         }
 161
 162         double moonJd = Math.floor(DateTimeUtils.midnightDateToJulianDate(calendar)) - 2400000.0;
 163         moonJd -= ((calendar.get(Calendar.ZONE_OFFSET) + calendar.get(Calendar.DST_OFFSET)) / 60000.0) / 1440.0;
 164
 165         double sphi = SN(phi);
 166         double cphi = CS(phi);
 167         double sinho = SN(8.0 / 60.0);
 168
 169         int hour = 1;
 170         double utrise = -1;
 171         double utset = -1;
 172         do {
 173             double yminus = SINALT(moonJd, hour - 1, lambda, cphi, sphi) - sinho;
 174             double yo = SINALT(moonJd, hour, lambda, cphi, sphi) - sinho;
 175             double yplus = SINALT(moonJd, hour + 1, lambda, cphi, sphi) - sinho;
 176             double[] quadRet = QUAD(yminus, yo, yplus);
 177             if (quadRet[3] == 1) {
 178                 if (yminus < 0) {
 179                     utrise = hour + quadRet[1];
 180                 } else {
 181                     utset = hour + quadRet[1];
 182                 }
 183             }
 184             if (quadRet[3] == 2) {
 185                 if (quadRet[0] < 0) {
 186                     utrise = hour + quadRet[2];
 187                     utset = hour + quadRet[1];
 188                 } else {
 189                     utrise = hour + quadRet[1];
 190                     utset = hour + quadRet[2];
 191                 }
 192             }
 193             yminus = yplus;
 194             hour += 2;
 195         } while (hour < 25 && (utrise == -1 || utset == -1));
 196
 197         double rise = prepareTime(utrise);
 198         double set = prepareTime(utset);
 199
 200         return new double[] { rise, set };
 201     }
 202
 203     /**
 204      * Prepares the coordinate for moonrise and moonset calculation.
 205      */
 206     private double prepareCoordinate(double coordinate, double system) {
 207         double c = Math.abs(coordinate);
 208
 209         if (c - Math.floor(c) >= .599) {
 210             c = Math.floor(c) + (c - Math.floor(c)) / 1 * .6;
 211         }
 212         if (c > system) {
 213             c = Math.floor(c) % system + (c - Math.floor(c));
 214         }
 215         return Math.round(c * 100.0) / 100.0;
 216     }
 217
 218     /**
 219      * Prepares a time value for converting to a calendar object.
 220      */
 221     private double prepareTime(double riseSet) {
 222         if (riseSet == -1) {
 223             return riseSet;
 224         }
 225         double riseMinute = (riseSet - Math.floor(riseSet)) * 60.0 / 100.0;
 226         double rounded;
 227         if (riseMinute >= .595) {
 228             riseMinute = 0;
 229             rounded = riseSet + 1;
 230         } else {
 231             rounded = riseSet;
 232         }
 233         rounded = Math.floor(rounded) + riseMinute;
 234
 235         BigDecimal bd = new BigDecimal(Double.toString(rounded));
 236         bd = bd.setScale(2, BigDecimal.ROUND_HALF_UP);
 237         return bd.doubleValue();
 238     }
 239
 240     /**
 241      * Calculates the moon phase.
 242      */
 243     private double calcMoonPhase(double k, double mode) {
 244         double kMod = Math.floor(k) + mode;
 245         double t = kMod / 1236.85;
 246         double e = var_e(t);
 247         double m = var_m(kMod, t);
 248         double m1 = var_m1(kMod, t);
 249         double f = var_f(kMod, t);
 250         double o = var_o(kMod, t);
 251         double jd = var_jde(kMod, t);
 252         if (mode == NEW_MOON) {
 253             jd += -.4072 * SN(m1) + .17241 * e * SN(m) + .01608 * SN(2 * m1) + .01039 * SN(2 * f)
 254                     + .00739 * e * SN(m1 - m) - .00514 * e * SN(m1 + m) + .00208 * e * e * SN(2 * m)
 255                     - .00111 * SN(m1 - 2 * f) - .00057 * SN(m1 + 2 * f);
 256             jd += .00056 * e * SN(2 * m1 + m) - .00042 * SN(3 * m1) + .00042 * e * SN(m + 2 * f)
 257                     + .00038 * e * SN(m - 2 * f) - .00024 * e * SN(2 * m1 - m) - .00017 * SN(o)
 258                     - .00007 * SN(m1 + 2 * m) + .00004 * SN(2 * m1 - 2 * f);
 259             jd += .00004 * SN(3 * m) + .00003 * SN(m1 + m - 2 * f) + .00003 * SN(2 * m1 + 2 * f)
 260                     - .00003 * SN(m1 + m + 2 * f) + .00003 * SN(m1 - m + 2 * f) - .00002 * SN(m1 - m - 2 * f)
 261                     - .00002 * SN(3 * m1 + m);
 262             jd += .00002 * SN(4 * m1);
 263         } else if (mode == FULL_MOON) {
 264             jd += -.40614 * SN(m1) + .17302 * e * SN(m) + .01614 * SN(2 * m1) + .01043 * SN(2 * f)
 265                     + .00734 * e * SN(m1 - m) - .00515 * e * SN(m1 + m) + .00209 * e * e * SN(2 * m)
 266                     - .00111 * SN(m1 - 2 * f) - .00057 * SN(m1 + 2 * f);
 267             jd += .00056 * e * SN(2 * m1 + m) - .00042 * SN(3 * m1) + .00042 * e * SN(m + 2 * f)
 268                     + .00038 * e * SN(m - 2 * f) - .00024 * e * SN(2 * m1 - m) - .00017 * SN(o)
 269                     - .00007 * SN(m1 + 2 * m) + .00004 * SN(2 * m1 - 2 * f);
 270             jd += .00004 * SN(3 * m) + .00003 * SN(m1 + m - 2 * f) + .00003 * SN(2 * m1 + 2 * f)
 271                     - .00003 * SN(m1 + m + 2 * f) + .00003 * SN(m1 - m + 2 * f) - .00002 * SN(m1 - m - 2 * f)
 272                     - .00002 * SN(3 * m1 + m);
 273             jd += .00002 * SN(4 * m1);
 274         } else {
 275             jd += -.62801 * SN(m1) + .17172 * e * SN(m) - .01183 * e * SN(m1 + m) + .00862 * SN(2 * m1)
 276                     + .00804 * SN(2 * f) + .00454 * e * SN(m1 - m) + .00204 * e * e * SN(2 * m) - .0018 * SN(m1 - 2 * f)
 277                     - .0007 * SN(m1 + 2 * f);
 278             jd += -.0004 * SN(3 * m1) - .00034 * e * SN(2 * m1 - m) + .00032 * e * SN(m + 2 * f)
 279                     + .00032 * e * SN(m - 2 * f) - .00028 * e * e * SN(m1 + 2 * m) + .00027 * e * SN(2 * m1 + m)
 280                     - .00017 * SN(o);
 281             jd += -.00005 * SN(m1 - m - 2 * f) + .00004 * SN(2 * m1 + 2 * f) - .00004 * SN(m1 + m + 2 * f)
 282                     + .00004 * SN(m1 - 2 * m) + .00003 * SN(m1 + m - 2 * f) + .00003 * SN(3 * m)
 283                     + .00002 * SN(2 * m1 - 2 * f);
 284             jd += .00002 * SN(m1 - m + 2 * f) - .00002 * SN(3 * m1 + m);
 285             double w = .00306 - .00038 * e * CS(m) + .00026 * CS(m1) - .00002 * CS(m1 - m) + .00002 * CS(m1 + m)
 286                     + .00002 * CS(2 * f);
 287             jd += (mode == FIRST_QUARTER) ? w : -w;
 288         }
 289         return moonCorrection(jd, t, kMod);
 290     }
 291
 292     /**
 293      * Calculates the eclipse.
 294      */
 295     private double getEclipse(double k, EclipseType typ, EclipseKind eclipse) {
 296         double kMod = Math.floor(k) + ((typ == EclipseType.SUN) ? 0 : 0.5);
 297         double t = kMod / 1236.85;
 298         double f = var_f(kMod, t);
 299         double jd = 0;
 300         double ringTest = 0;
 301         if (SN(Math.abs(f)) <= .36) {
 302             double o = var_o(kMod, t);
 303             double f1 = f - .02665 * SN(o);
 304             double a1 = 299.77 + .107408 * kMod - .009173 * t * t;
 305             double e = var_e(t);
 306             double m = var_m(kMod, t);
 307             double m1 = var_m1(kMod, t);
 308             double p = .207 * e * SN(m) + .0024 * e * SN(2 * m) - .0392 * SN(m1) + .0116 * SN(2 * m1)
 309                     - .0073 * e * SN(m1 + m) + .0067 * e * SN(m1 - m) + .0118 * SN(2 * f1);
 310             double q = 5.2207 - .0048 * e * CS(m) + .002 * e * CS(2 * m) - .3299 * CS(m1) - .006 * e * CS(m1 + m)
 311                     + .0041 * e * CS(m1 - m);
 312             double g = (p * CS(f1) + q * SN(f1)) * (1 - .0048 * CS(Math.abs(f1)));
 313             double u = .0059 + .0046 * e * CS(m) - .0182 * CS(m1) + .0004 * CS(2 * m1) - .0005 * CS(m + m1);
 314             jd = var_jde(kMod, t);
 315             jd += (typ == EclipseType.MOON) ? -.4065 * SN(m1) + .1727 * e * SN(m) : -.4075 * SN(m1) + .1721 * e * SN(m);
 316
 317             jd += .0161 * SN(2 * m1) - .0097 * SN(2 * f1) + .0073 * e * SN(m1 - m) - .005 * e * SN(m1 + m)
 318                     - .0023 * SN(m1 - 2 * f1) + .0021 * e * SN(2 * m);
 319             jd += .0012 * SN(m1 + 2 * f1) + .0006 * e * SN(2 * m1 + m) - .0004 * SN(3 * m1) - .0003 * e * SN(m + 2 * f1)
 320                     + .0003 * SN(a1) - .0002 * e * SN(m - 2 * f1) - .0002 * e * SN(2 * m1 - m) - .0002 * SN(o);
 321             switch (typ) {
 322                 case MOON:
 323                     if ((1.0248 - u - Math.abs(g)) / .545 <= 0) {
 324                         jd = 0; // no moon eclipse
 325                     }
 326                     if (eclipse == EclipseKind.PARTIAL && (1.0128 - u - Math.abs(g)) / .545 > 0
 327                             && (.4678 - u) * (.4678 - u) - g * g > 0) {
 328                         jd = 0; // no partial moon eclipse
 329                     }
 330                     if (eclipse == EclipseKind.TOTAL
 331                             && ((1.0128 - u - Math.abs(g)) / .545 <= 0 != (.4678 - u) * (.4678 - u) - g * g <= 0)) {
 332                         jd = 0; // no total moon eclipse
 333                     }
 334                     break;
 335                 case SUN:
 336                     if (Math.abs(g) > 1.5433 + u) {
 337                         jd = 0; // no sun eclipse
 338                     }
 339                     if (eclipse == EclipseKind.PARTIAL && ((g >= -.9972 && g <= .9972)
 340                             || (Math.abs(g) >= .9972 && Math.abs(g) < .9972 + Math.abs(u)))) {
 341                         jd = 0; // no partial sun eclipse
 342                     }
 343                     if (eclipse != EclipseKind.PARTIAL) {
 344                         if ((g < -.9972 || g > .9972) || (Math.abs(g) < .9972 && Math.abs(g) > .9972 + Math.abs(u))) {
 345                             jd = 0; // no ring or total sun eclipse
 346                         }
 347                         if (u > .0047 || u >= .00464 * Math.sqrt(1 - g * g)) {
 348                             ringTest = 1; // no total sun eclipse
 349                         }
 350                         if (ringTest == 1 && eclipse == EclipseKind.TOTAL) {
 351                             jd = 0;
 352                         }
 353                         if (ringTest == 0 && eclipse == EclipseKind.RING) {
 354                             jd = 0;
 355                         }
 356                     }
 357                     break;
 358             }
 359         }
 360         return jd;
 361     }
 362
 363     /**
 364      * Calculates the illumination.
 365      */
 366     private double getIllumination(double jd) {
 367         double t = (jd - 2451545) / 36525;
 368         double d = 297.8502042 + 445267.11151686 * t - .00163 * t * t + t * t * t / 545868 - t * t * t * t / 113065000;
 369         double m = 357.5291092 + 35999.0502909 * t - .0001536 * t * t + t * t * t / 24490000;
 370         double m1 = 134.9634114 + 477198.8676313 * t + .008997 * t * t + t * t * t / 69699 - t * t * t * t / 14712000;
 371         double i = 180 - d - 6.289 * SN(m1) + 2.1 * SN(m) - 1.274 * SN(2 * d - m1) - .658 * SN(2 * d)
 372                 - .241 * SN(2 * m1) - .110 * SN(d);
 373         return (1 + CS(i)) / 2 * 100.0;
 374     }
 375
 376     /**
 377      * Calculates the next moon phase.
 378      */
 379     private double getNextPhase(Calendar cal, double midnightJd, double mode) {
 380         double tz = 0;
 381         double phaseJd = 0;
 382         do {
 383             double k = var_k(cal, tz);
 384             tz += 1;
 385             phaseJd = calcMoonPhase(k, mode);
 386         } while (phaseJd <= midnightJd);
 387         return phaseJd;
 388     }
 389
 390     /**
 391      * Calculates the previous moon phase.
 392      */
 393     public double getPreviousPhase(Calendar cal, double jd, double mode) {
 394         double tz = 0;
 395         double phaseJd = 0;
 396         do {
 397             double k = var_k(cal, tz);
 398             tz -= 1;
 399             phaseJd = calcMoonPhase(k, mode);
 400         } while (phaseJd > jd);
 401         return phaseJd;
 402     }
 403
 404     /**
 405      * Calculates the next eclipse.
 406      */
 407     protected double getEclipse(Calendar cal, EclipseType type, double midnightJd, EclipseKind eclipse) {
 408         double tz = 0;
 409         double eclipseJd = 0;
 410         do {
 411             double k = var_k(cal, tz);
 412             tz += 1;
 413             eclipseJd = getEclipse(k, type, eclipse);
 414         } while (eclipseJd <= midnightJd);
 415         return eclipseJd;
 416     }
 417
 418     /**
 419      * Calculates the date, where the moon is furthest away from the earth.
 420      */
 421     private double getApogee(double julianDate, double decimalYear) {
 422         double k = Math.floor((decimalYear - 1999.97) * 13.2555) + .5;
 423         double jd = 0;
 424         do {
 425             double t = k / 1325.55;
 426             double d = 171.9179 + 335.9106046 * k - .010025 * t * t - .00001156 * t * t * t
 427                     + .000000055 * t * t * t * t;
 428             double m = 347.3477 + 27.1577721 * k - .0008323 * t * t - .000001 * t * t * t;
 429             double f = 316.6109 + 364.5287911 * k - .0125131 * t * t - .0000148 * t * t * t;
 430             jd = 2451534.6698 + 27.55454988 * k - .0006886 * t * t - .000001098 * t * t * t + .0000000052 * t * t
 431                     + .4392 * SN(2 * d) + .0684 * SN(4 * d) + (.0456 - .00011 * t) * SN(m)
 432                     + (.0426 - .00011 * t) * SN(2 * d - m) + .0212 * SN(2 * f);
 433             jd += -.0189 * SN(d) + .0144 * SN(6 * d) + .0113 * SN(4 * d - m) + .0047 * SN(2 * d + 2 * f)
 434                     + .0036 * SN(d + m) + .0035 * SN(8 * d) + .0034 * SN(6 * d - m) - .0034 * SN(2 * d - 2 * f)
 435                     + .0022 * SN(2 * d - 2 * m) - .0017 * SN(3 * d);
 436             jd += .0013 * SN(4 * d + 2 * f) + .0011 * SN(8 * d - m) + .001 * SN(4 * d - 2 * m) + .0009 * SN(10 * d)
 437                     + .0007 * SN(3 * d + m) + .0006 * SN(2 * m) + .0005 * SN(2 * d + m) + .0005 * SN(2 * d + 2 * m)
 438                     + .0004 * SN(6 * d + 2 * f);
 439             jd += .0004 * SN(6 * d - 2 * m) + .0004 * SN(10 * d - m) - .0004 * SN(5 * d) - .0004 * SN(4 * d - 2 * f)
 440                     + .0003 * SN(2 * f + m) + .0003 * SN(12 * d) + .0003 * SN(2 * d + 2 * f - m) - .0003 * SN(d - m);
 441             k += 1;
 442         } while (jd < julianDate);
 443         return jd;
 444     }
 445
 446     /**
 447      * Calculates the date, where the moon is closest to the earth.
 448      */
 449     private double getPerigee(double julianDate, double decimalYear) {
 450         double k = Math.floor((decimalYear - 1999.97) * 13.2555);
 451         double jd = 0;
 452         do {
 453             double t = k / 1325.55;
 454             double d = 171.9179 + 335.9106046 * k - .010025 * t * t - .00001156 * t * t * t
 455                     + .000000055 * t * t * t * t;
 456             double m = 347.3477 + 27.1577721 * k - .0008323 * t * t - .000001 * t * t * t;
 457             double f = 316.6109 + 364.5287911 * k - .0125131 * t * t - .0000148 * t * t * t;
 458             jd = 2451534.6698 + 27.55454988 * k - .0006886 * t * t - .000001098 * t * t * t + .0000000052 * t * t
 459                     - 1.6769 * SN(2 * d) + .4589 * SN(4 * d) - .1856 * SN(6 * d) + .0883 * SN(8 * d);
 460             jd += -(.0773 + .00019 * t) * SN(2 * d - m) + (.0502 - .00013 * t) * SN(m) - .046 * SN(10 * d)
 461                     + (.0422 - .00011 * t) * SN(4 * d - m) - .0256 * SN(6 * d - m) + .0253 * SN(12 * d) + .0237 * SN(d);
 462             jd += .0162 * SN(8 * d - m) - .0145 * SN(14 * d) + .0129 * SN(2 * f) - .0112 * SN(3 * d)
 463                     - .0104 * SN(10 * d - m) + .0086 * SN(16 * d) + .0069 * SN(12 * d - m) + .0066 * SN(5 * d)
 464                     - .0053 * SN(2 * d + 2 * f);
 465             jd += -.0052 * SN(18 * d) - .0046 * SN(14 * d - m) - .0041 * SN(7 * d) + .004 * SN(2 * d + m)
 466                     + .0032 * SN(20 * d) - .0032 * SN(d + m) + .0031 * SN(16 * d - m);
 467             jd += -.0029 * SN(4 * d + m) - .0027 * SN(2 * d - 2 * m) + .0024 * SN(4 * d - 2 * m)
 468                     - .0021 * SN(6 * d - 2 * m) - .0021 * SN(22 * d) - .0021 * SN(18 * d - m);
 469             jd += .0019 * SN(6 * d + m) - .0018 * SN(11 * d) - .0014 * SN(8 * d + m) - .0014 * SN(4 * d - 2 * f)
 470                     - .0014 * SN(6 * d - 2 * f) + .0014 * SN(3 * d + m) - .0014 * SN(5 * d + m) + .0013 * SN(13 * d);
 471             jd += .0013 * SN(20 * d - m) + .0011 * SN(3 * d + 2 * m) - .0011 * SN(4 * d + 2 * f - 2 * m)
 472                     - .001 * SN(d + 2 * m) - .0009 * SN(22 * d - m) - .0008 * SN(4 * f) + .0008 * SN(6 * d - 2 * f)
 473                     + .0008 * SN(2 * d - 2 * f + m);
 474             jd += .0007 * SN(2 * m) + .0007 * SN(2 * f - m) + .0007 * SN(2 * d + 4 * f) - .0006 * SN(2 * f - 2 * m)
 475                     - .0006 * SN(2 * d - 2 * f + 2 * m) + .0006 * SN(24 * d) + .0005 * SN(4 * d - 4 * f)
 476                     + .0005 * SN(2 * d + 2 * m) - .0004 * SN(d - m) + .0027 * SN(9 * d) + .0027 * SN(4 * d + 2 * f);
 477             k += 1;
 478         } while (jd < julianDate);
 479         return jd;
 480     }
 481
 482     /**
 483      * Calculates the distance from the moon to earth.
 484      */
 485     private double getDistance(double jd) {
 486         double t = (jd - 2451545) / 36525;
 487         double d = 297.8502042 + 445267.11151686 * t - .00163 * t * t + t * t * t / 545868 - t * t * t * t / 113065000;
 488         double m = 357.5291092 + 35999.0502909 * t - .0001536 * t * t + t * t * t / 24490000;
 489         double m1 = 134.9634114 + 477198.8676313 * t + .008997 * t * t + t * t * t / 69699 - t * t * t * t / 14712000;
 490         double f = 93.27209929999999 + 483202.0175273 * t - .0034029 * t * t - t * t * t / 3526000
 491                 + t * t * t * t / 863310000;
 492         double sr = 385000.56 + getCoefficient(d, m, m1, f) / 1000;
 493         return sr;
 494     }
 495
 496     public 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 }