/* Compute the longitude of the sun at any time.
Parameter: floating number jdn, the number of days since 1/1/4713 BC noon
Algorithm from: "Astronomical Algorithms" by Jean Meeus, 1998
347712
function SunLongitude (jdn) {
  var T, T2, dr, M, L0, DL, lambda, theta, omega
  T = (jdn - 2451545.0) / 36525 // Time in Julian centuries from 2000-01-01 12:00:00 GMT
  T2 = T * T
  dr = Math.PI / 180 // degree to radian
  M = 357.52910 + 35999.05030 * T - 0.0001559 * T2 - 0.00000048 * T * T2 // mean anomaly, degree
  L0 = 280.46645 + 36000.76983 * T + 0.0003032 * T2 // mean longitude, degree
  DL = (1.914600 - 0.004817 * T - 0.000014 * T2) * Math.sin(dr * M)
  DL = DL + (0.019993 - 0.000101 * T) * Math.sin(dr * 2 * M) + 0.000290 * Math.sin(dr * 3 * M)
  theta = L0 + DL // true longitude, degree
  // obtain apparent longitude by correcting for nutation and aberration
  omega = 125.04 - 1934.136 * T
  lambda = theta - 0.00569 - 0.00478 * Math.sin(omega * dr)
  // Convert to radians
  lambda = lambda * dr
  lambda = lambda - Math.PI * 2 * (Math.floor(lambda / (Math.PI * 2))) // Normalize to (0, 2*PI)
  return lambda
}
1343359
def to_jd(year, week, day):
    ''''''
    return day + n_weeks(SUN, gregorian.to_jd(year - 1, 12, 28), week)
1180275
def pfLon(jd, lat, lon):
    
    sun = swe.sweObjectLon(const.SUN, jd)
    moon = swe.sweObjectLon(const.MOON, jd)
    asc = swe.sweHousesLon(jd, lat, lon,
                           const.HOUSES_DEFAULT)[1][0]
    
    if isDiurnal(jd, lat, lon):
        return angle.norm(asc + moon - sun)
    else:
        return angle.norm(asc + sun - moon)
316722
function jdnToEthiopic(jdn, era = JD_EPOCH_OFFSET_AMETE_MIHRET) {
  const r = mod((jdn - era), 1461);
  const n = mod(r, 365) + 365 * Math.floor(r / 1460);

  const year = 4 * Math.floor((jdn - era) / 1461) + Math.floor(r / 365) - Math.floor(r / 1460);
  const month = Math.floor(n / 30) + 1;
  const day = mod(n, 30) + 1;

  return { year, month, day };
}
16669
def compute_moon
      # === 月が光を発した時刻 t1(jd) の計算
      t_1_jd = calc_t1(:moon, @tdb)
      # === 時刻 t1(= TDB) におけるの位置・速度（ICRS 座標）の計算 (地球, 月, 太陽)
      Const::BODIES.each { |k, v| @icrs_1[k] = get_icrs(v, t_1_jd) }
      # === 時刻 t2 における地球重心から時刻 t1 における月への方向ベクトルの計算
      v_12 = calc_unit_vector(@icrs_2[:earth][0,3], @icrs_1[:moon][0,3])
      # === GCRS 座標系: 光行差の補正（方向ベクトルの Lorentz 変換）
      dd = conv_lorentz(v_12)
      pos_moon = calc_pos(dd, @r_e[:moon])
      # === 瞬時の真座標系: GCRS への bias & precession（歳差） & nutation（章動）の適用
      bpn = EphBpn.new(@tdb.strftime("%Y%m%d%H%M%S"))
      pos_moon_bpn = bpn.apply_bpn(pos_moon)
      # === 座標変換
      eq_pol_m  = MkCoord.rect2pol(pos_moon_bpn)
      ec_rect_m = MkCoord.rect_eq2ec(pos_moon_bpn, bpn.eps)
      ec_pol_m  = MkCoord.rect2pol(ec_rect_m)
      # === 視半径／（地平）視差計算
      radius = Math.asin(@am / (eq_pol_m[2] * Const::AU / 1000))
      radius *= 180 / Math::PI * 3600
      parallax = Math.asin(@re / (eq_pol_m[2] * Const::AU / 1000))
      parallax *= 180 / Math::PI * 3600
      return [eq_pol_m, ec_pol_m, [radius, parallax]]
    rescue => e
      raise
    end
246110
function () {
            if (this.time === undefined) {
                this.computeTime();
            }
            var cal = new GregorianCalendar(this.timezoneOffset, this.locale);
            cal.setTime(this.time);
            return cal;
        }
461017
public Equatorial getMoonPosition()
    {
        //
        // See page 142 of "Practial Astronomy with your Calculator",
        // by Peter Duffet-Smith, for details on the algorithm.
        //
        if (moonPosition == null) {
            // Calculate the solar longitude.  Has the side effect of
            // filling in "meanAnomalySun" as well.
            double sunLong = getSunLongitude();

            //
            // Find the # of days since the epoch of our orbital parameters.
            // TODO: Convert the time of day portion into ephemeris time
            //
            double day = getJulianDay() - JD_EPOCH;       // Days since epoch

            // Calculate the mean longitude and anomaly of the moon, based on
            // a circular orbit.  Similar to the corresponding solar calculation.
            double meanLongitude = norm2PI(13.1763966*PI/180*day + moonL0);
            double meanAnomalyMoon = norm2PI(meanLongitude - 0.1114041*PI/180 * day - moonP0);

            //
            // Calculate the following corrections:
            //  Evection:   the sun's gravity affects the moon's eccentricity
            //  Annual Eqn: variation in the effect due to earth-sun distance
            //  A3:         correction factor (for ???)
            //
            double evection = 1.2739*PI/180 * Math.sin(2 * (meanLongitude - sunLong)
                                                - meanAnomalyMoon);
            double annual   = 0.1858*PI/180 * Math.sin(meanAnomalySun);
            double a3       = 0.3700*PI/180 * Math.sin(meanAnomalySun);

            meanAnomalyMoon += evection - annual - a3;

            //
            // More correction factors:
            //  center  equation of the center correction
            //  a4      yet another error correction (???)
            //
            // TODO: Skip the equation of the center correction and solve Kepler's eqn?
            //
            double center = 6.2886*PI/180 * Math.sin(meanAnomalyMoon);
            double a4 =     0.2140*PI/180 * Math.sin(2 * meanAnomalyMoon);

            // Now find the moon's corrected longitude
            moonLongitude = meanLongitude + evection + center - annual + a4;

            //
            // And finally, find the variation, caused by the fact that the sun's
            // gravitational pull on the moon varies depending on which side of
            // the earth the moon is on
            //
            double variation = 0.6583*PI/180 * Math.sin(2*(moonLongitude - sunLong));

            moonLongitude += variation;

            //
            // What we've calculated so far is the moon's longitude in the plane
            // of its own orbit.  Now map to the ecliptic to get the latitude
            // and longitude.  First we need to find the longitude of the ascending
            // node, the position on the ecliptic where it is crossed by the moon's
            // orbit as it crosses from the southern to the northern hemisphere.
            //
            double nodeLongitude = norm2PI(moonN0 - 0.0529539*PI/180 * day);

            nodeLongitude -= 0.16*PI/180 * Math.sin(meanAnomalySun);

            double y = Math.sin(moonLongitude - nodeLongitude);
            double x = Math.cos(moonLongitude - nodeLongitude);

            moonEclipLong = Math.atan2(y*Math.cos(moonI), x) + nodeLongitude;
            double moonEclipLat = Math.asin(y * Math.sin(moonI));

            moonPosition = eclipticToEquatorial(moonEclipLong, moonEclipLat);
        }
        return moonPosition;
    }
4764
def to_jd(year, month, day):
    ''''''
    return (day + ceil(29.5 * (month - 1)) + (year - 1) * 354 + trunc((3 + (11 * year)) / 30) + EPOCH) - 1
790368
final boolean hasNoMajorSolarTerm(long utcDays) {
        double jd0 = JulianDay.ofEphemerisTime(this.midnight(utcDays)).getValue();
        int index0 = (2 + (int) Math.floor(SolarTerm.solarLongitude(jd0) / 30)) % 12;
        double jd1 = JulianDay.ofEphemerisTime(this.midnight(this.newMoonOnOrAfter(utcDays + 1))).getValue();
        int index1 = (2 + (int) Math.floor(SolarTerm.solarLongitude(jd1) / 30)) % 12;
        return (index0 == index1);
    }
1097608
def to_lat_lon
      k0 = 0.9996
      x  = @e - 500000 # Remove Longitude offset
      y  = @n

      # Set hemisphere (1=Northern, 0=Southern)
      y    -= 10000000.0 unless UTMZones::northern_hemisphere? @zone

      ecc = @ellipsoid.eccentricity
      eccPrimeSquared = (ecc)/(1-ecc)
      m  = y / k0
      mu = m / (@ellipsoid.radius * (1 - ecc / 4 - 3 * ecc ** 2 / 64 - 5 * ecc ** 3 / 256))
      e1 = (1 - Math::sqrt(1 - ecc)) / (1 + Math::sqrt(1 - ecc))
      phi1rad = mu +
        (3 * e1 / 2 - 27 * e1 ** 3 / 32) * Math::sin(2 * mu) +
        (21 * e1 ** 2 / 16 - 55 * e1 ** 4 / 32) * Math::sin(4 * mu) +
        (151 * e1 ** 3 / 96) * Math::sin(6 * mu)
      n1 = @ellipsoid.radius / Math::sqrt(1 - ecc * Math::sin(phi1rad) ** 2)
      t1 = Math::tan(phi1rad) ** 2
      c1 = ecc * Math::cos(phi1rad)**2
      r1 = @ellipsoid.radius * (1 - ecc) / ((1 - ecc * Math::sin(phi1rad) ** 2) ** 1.5)
      d = x / (n1 * k0)
      latitude_rad = phi1rad-
        (n1 * Math::tan(phi1rad) / r1) * (d * d / 2 - 
                                    (5 + 3 * t1 + 10 * c1 - 4 * c1 * c1 - 9 * eccPrimeSquared) * d ** 4 / 24 +
                                    (61 + 90 * t1 + 298 * c1 + 45 * t1 * t1 - 
                                     252 * eccPrimeSquared - 3 * c1 * c1) * d ** 6 / 720)
      latitude_deg = latitude_rad * 180.0 / Math::PI
      lon_tmp = (d - (1 + 2 * t1 + c1) * d ** 3 / 6 +
                   (5 - 2 * c1 + 28 * t1 - 3 * c1 * c1 + 8 * eccPrimeSquared +
                    24 * t1 * t1) * d ** 5 / 120) / Math::cos(phi1rad)
      longitude_deg = UTMZones::lon_origin(@zone) + lon_tmp * 180.0 / Math::PI
      LatLon.new latitude_deg, longitude_deg
    end
576209
function d2j(jdn) {
    var gy = d2g(jdn).gy // Calculate Gregorian year (gy).
        , jy = gy - 621
        , r = jalCal(jy)
        , jdn1f = g2d(gy, 3, r.march)
        , jd
        , jm
        , k;

    // Find number of days that passed since 1 Farvardin.
    k = jdn - jdn1f;
    if (k >= 0) {
        if (k <= 185) {
            // The first 6 months.
            jm = 1 + div(k, 31);
            jd = mod(k, 31) + 1;
            return  { jy: jy
                , jm: jm
                , jd: jd
            };
        } else {
            // The remaining months.
            k -= 186;
        }
    } else {
        // Previous Jalali year.
        jy -= 1;
        k += 179;
        if (r.leap === 1)
            k += 1;
    }
    jm = 7 + div(k, 30);
    jd = mod(k, 30) + 1;
    return  { jy: jy
        , jm: jm
        , jd: jd
    };
}
312719
def getLunarMonth11(yy, timeZone):
    ''''''
    # off = jdFromDate(31, 12, yy) \
    #            - 2415021.076998695
    off = jdFromDate(31, 12, yy) - 2415021.
    k = int(off / 29.530588853)
    nm = getNewMoonDay(k, timeZone)
    sunLong = getSunLongitude(nm, timeZone)
    # sun longitude at local midnight
    if (sunLong >= 9):
        nm = getNewMoonDay(k - 1, timeZone)
    return nm
193542
def sun_ecliptic_longitude(utc_time):
    
    jdate = jdays2000(utc_time) / 36525.0
    # mean anomaly, rad
    m_a = np.deg2rad(357.52910 +
                     35999.05030 * jdate -
                     0.0001559 * jdate * jdate -
                     0.00000048 * jdate * jdate * jdate)
    # mean longitude, deg
    l_0 = 280.46645 + 36000.76983 * jdate + 0.0003032 * jdate * jdate
    d_l = ((1.914600 - 0.004817 * jdate - 0.000014 * jdate * jdate) * np.sin(m_a) +
           (0.019993 - 0.000101 * jdate) * np.sin(2 * m_a) + 0.000290 * np.sin(3 * m_a))
    # true longitude, deg
    l__ = l_0 + d_l
    return np.deg2rad(l__)
557068
def nextSunrise(date, pos):
    
    jd = eph.nextSunrise(date.jd, pos.lat, pos.lon)
    return Datetime.fromJD(jd, date.utcoffset)
996624
function gregorianToJDN(year = 1, month = 1, day = 1, JD_OFFSET = JD_EPOCH_OFFSET_GREGORIAN) {
  const s = Math.floor(year / 4)
          - Math.floor((year - 1) / 4)
          - Math.floor(year / 100)
          + Math.floor((year - 1) / 100)
          + Math.floor(year / 400)
          - Math.floor((year - 1) / 400);
  const t = Math.floor((14 - month) / 12);
  const n = 31 * t * (month - 1)
          + (1 - t) * (59 + s + 30 * (month - 3) + Math.floor((3 * month - 7) / 5))
          + day - 1;
  const j = JD_OFFSET
          + 365 * (year - 1)
          + Math.floor((year - 1) / 4)
          - Math.floor((year - 1) / 100)
          + Math.floor((year - 1) / 400)
          + n;

  return j;
}
1322351
long calculateFirstDayOfYearMillis(int year) {
        if (year > MAX_YEAR) {
            throw new ArithmeticException("Year is too large: " + year + " > " + MAX_YEAR);
        }
        if (year < MIN_YEAR) {
            throw new ArithmeticException("Year is too small: " + year + " < " + MIN_YEAR);
        }

        // Java epoch is 1970-01-01 Gregorian which is 0622-07-16 Islamic.
        // 0001-01-01 Islamic is -42520809600000L
        // would prefer to calculate against year zero, but leap year
        // can be in that year so it doesn't work
        year--;
        long cycle = year / CYCLE;
        long millis = MILLIS_YEAR_1 + cycle * MILLIS_PER_CYCLE;
        int cycleRemainder = (year % CYCLE) + 1;
        
        for (int i = 1; i < cycleRemainder; i++) {
            millis += (isLeapYear(i) ? MILLIS_PER_LONG_YEAR : MILLIS_PER_SHORT_YEAR);
        }
        
        return millis;
    }
1261802
private static int[] jdToGregorian(double jd) {
       double JULIAN_EPOCH = 1721425.5;
       double wjd, depoch, quadricent, dqc, cent, dcent, quad, dquad, yindex, yearday, leapadj;
       int year, month, day;
       
       wjd = Math.floor(jd - 0.5) + 0.5;
       depoch = wjd - JULIAN_EPOCH;
       quadricent = Math.floor(depoch / 146097);
       dqc = depoch % 146097;
       cent = Math.floor(dqc / 36524);
       dcent = dqc % 36524;
       quad = Math.floor(dcent / 1461);
       dquad = dcent % 1461;
       yindex = Math.floor(dquad / 365);
       year = (int)((quadricent * 400) + (cent * 100) + (quad * 4) + yindex);
       
       if (!((cent == 4) || (yindex == 4))) {
          year++;
       }
       
       yearday = wjd - gregorianToJD(year, 1, 1);
       leapadj = ((wjd < gregorianToJD(year, 3, 1)) ? 0
             :
             (isGregorianLeap(year) ? 1 : 2)
             );
       
       month = (int)Math.floor((((yearday + leapadj) * 12) + 373) / 367);
       day = (int)(wjd - gregorianToJD(year, month, 1)) + 1;

       int[] julianDate = new int[3];
       
       julianDate[0] = year;
       julianDate[1] = month;
       julianDate[2] = day;
       
       return julianDate;
    }