position angle of axis and sub-sun point (feat. in Chapter53 Meeus

[YingChangJ]

Hello,

When I was trying to make this animation display how Moon looks like from earth, it seems that I need more info than just Libration. At least, where the pole points to is a degrees of freedom. (yes, it keeps a small degree away from the normal of ecliptic plane, but I am not very sure how to rotate this pole)

After reading the Chapter53 Meeus, I found the calculation of sub-sun and position angle is not expensive, especially when the libration is alrighdy available.

I think it might help someone have the same problem.

btw, sub-sun is not really needed if the position of Moon and Sun is available, I added just in convenience to draw. (librations and position angle of axis could give us a Rotation matrix of Moon's Selenographic coordinate system, that could be totally enough for all other things.)

Here is the implementation, beyond the old Libration function, (it have the same result as the book, but not test with a longer time range):

export function Libration(date) {
  const time = MakeTime(date);
  const t = time.tt / 36525.0;
  const t2 = t * t;
  const t3 = t2 * t;
  const t4 = t2 * t2;
  const moon = CalcMoon(time); //should be include procession be not include nutation
  const mlon = moon.geo_eclip_lon;
  const mlat = moon.geo_eclip_lat;
  const dist_km = moon.distance_au * KM_PER_AU;
  const tilt = e_tilt(time);
  const dpsi = tilt.dpsi / 3600;
  // Inclination angle
  const I = DEG2RAD * 1.54242;
  // Moon's argument of latitude in radians.
  const f =
    DEG2RAD *
    NormalizeLongitude(
      93.272095 +
        483202.0175233 * t -
        0.0036539 * t2 -
        t3 / 3526000 +
        t4 / 863310000
    );
  // Moon's ascending node's mean longitude in radians.
  const Omega =
    DEG2RAD *
    NormalizeLongitude(
      125.0445479 -
        1934.1362891 * t +
        0.0020754 * t2 +
        t3 / 467441 -
        t4 / 60616000
    );
  // Sun's mean anomaly.
  const m =
    DEG2RAD *
    NormalizeLongitude(
      357.5291092 + 35999.0502909 * t - 0.0001536 * t2 + t3 / 24490000
    );
  // Moon's mean anomaly.
  const mdash =
    DEG2RAD *
    NormalizeLongitude(
      134.9633964 +
        477198.8675055 * t +
        0.0087414 * t2 +
        t3 / 69699 -
        t4 / 14712000
    );
  // Moon's mean elongation.
  const d =
    DEG2RAD *
    NormalizeLongitude(
      297.8501921 +
        445267.1114034 * t -
        0.0018819 * t2 +
        t3 / 545868 -
        t4 / 113065000
    );
  // Eccentricity of the Earth's orbit.
  const e = 1.0 - 0.002516 * t - 0.0000074 * t2;
  const k1 = DEG2RAD * (119.75 + 131.849 * t);
  const k2 = DEG2RAD * (72.56 + 20.186 * t);
  let sigma = 0;
  let rho = 0;
  let tau = 0;
  function sub_planet(lambda, beta) {
    // Optical librations
    const w = lambda - DEG2RAD * dpsi - Omega;
    const a = Math.atan2(
      Math.sin(w) * Math.cos(beta) * Math.cos(I) - Math.sin(beta) * Math.sin(I),
      Math.cos(w) * Math.cos(beta)
    );
    const ldash = LongitudeOffset(RAD2DEG * (a - f));
    const bdash = Math.asin(
      -Math.sin(w) * Math.cos(beta) * Math.sin(I) - Math.sin(beta) * Math.cos(I)
    );
    // Physical librations
    rho =
      -0.02752 * Math.cos(mdash) +
      -0.02245 * Math.sin(f) +
      +0.00684 * Math.cos(mdash - 2 * f) +
      -0.00293 * Math.cos(2 * f) +
      -0.00085 * Math.cos(2 * f - 2 * d) +
      -0.00054 * Math.cos(mdash - 2 * d) +
      -0.0002 * Math.sin(mdash + f) +
      -0.0002 * Math.cos(mdash + 2 * f) +
      -0.0002 * Math.cos(mdash - f) +
      +0.00014 * Math.cos(mdash + 2 * f - 2 * d);
    sigma =
      -0.02816 * Math.sin(mdash) +
      +0.02244 * Math.cos(f) +
      -0.00682 * Math.sin(mdash - 2 * f) +
      -0.00279 * Math.sin(2 * f) +
      -0.00083 * Math.sin(2 * f - 2 * d) +
      +0.00069 * Math.sin(mdash - 2 * d) +
      +0.0004 * Math.cos(mdash + f) +
      -0.00025 * Math.sin(2 * mdash) +
      -0.00023 * Math.sin(mdash + 2 * f) +
      +0.0002 * Math.cos(mdash - f) +
      +0.00019 * Math.sin(mdash - f) +
      +0.00013 * Math.sin(mdash + 2 * f - 2 * d) +
      -0.0001 * Math.cos(mdash - 3 * f);
    tau =
      +0.0252 * e * Math.sin(m) +
      +0.00473 * Math.sin(2 * mdash - 2 * f) +
      -0.00467 * Math.sin(mdash) +
      +0.00396 * Math.sin(k1) +
      +0.00276 * Math.sin(2 * mdash - 2 * d) +
      +0.00196 * Math.sin(Omega) +
      -0.00183 * Math.cos(mdash - f) +
      +0.00115 * Math.sin(mdash - 2 * d) +
      -0.00096 * Math.sin(mdash - d) +
      +0.00046 * Math.sin(2 * f - 2 * d) +
      -0.00039 * Math.sin(mdash - f) +
      -0.00032 * Math.sin(mdash - m - d) +
      +0.00027 * Math.sin(2 * mdash - m - 2 * d) +
      +0.00023 * Math.sin(k2) +
      -0.00014 * Math.sin(2 * d) +
      +0.00014 * Math.cos(2 * mdash - 2 * f) +
      -0.00012 * Math.sin(mdash - 2 * f) +
      -0.00012 * Math.sin(2 * mdash) +
      +0.00011 * Math.sin(2 * mdash - 2 * m - 2 * d);
    const ldash2 =
      -tau + (rho * Math.cos(a) + sigma * Math.sin(a)) * Math.tan(bdash);
    const bdash2 = sigma * Math.cos(a) - rho * Math.sin(a);
    const b = RAD2DEG * bdash + bdash2;
    const l = ldash + ldash2;
    return { b, l };
  }

  const helio_moon_vect = HelioVector(Body.Moon, time);
  const helio_moon_ecl = Ecliptic(helio_moon_vect);

  const sub_sun = sub_planet(
    helio_moon_ecl.elon * DEG2RAD,
    helio_moon_ecl.elat * DEG2RAD
  );
  const { b, l } = sub_planet(mlon + dpsi, mlat);
  const diam_deg =
    2.0 *
    RAD2DEG *
    Math.atan(
      MOON_MEAN_RADIUS_KM /
        Math.sqrt(dist_km * dist_km - MOON_MEAN_RADIUS_KM * MOON_MEAN_RADIUS_KM)
    );

  const v = Omega + dpsi * DEG2RAD + (sigma / Math.sin(I)) * DEG2RAD;
  const x = Math.sin(I + rho * DEG2RAD) * Math.sin(v);
  const y =
    Math.sin(I + rho * DEG2RAD) * Math.cos(v) * Math.cos(tilt.tobl * DEG2RAD) -
    Math.cos(I + rho * DEG2RAD) * Math.sin(tilt.tobl * DEG2RAD);
  const omega = Math.atan2(x, y);

  // Convert spherical coordinates to a vector.
  // The MoonResult angles are already expressed in radians.
  const dist_cos_lat = moon.distance_au * Math.cos(moon.geo_eclip_lat);
  const ecm = [
    dist_cos_lat * Math.cos(moon.geo_eclip_lon),
    dist_cos_lat * Math.sin(moon.geo_eclip_lon),
    moon.distance_au * Math.sin(moon.geo_eclip_lat),
  ];
  // Obtain true and mean obliquity angles for the given time.
  // This serves to pre-calculate the nutation also, and cache it in `time`.
  // Convert ecliptic coordinates to equatorial coordinates, both in mean equinox of date.
  const eqm = obl_ecl2equ_vec(tilt.mobl, ecm);
  // Add nutation to convert ECM to true equatorial coordinates of date (EQD).
  const eqd = nutation(eqm, time, PrecessDirection.From2000);
  const eqd_vec = new VectorFromArray(eqd, time);
  const eqd_sphe = SphereFromVector(eqd_vec);
  const P = Math.asin(
    (Math.sqrt(x * x + y * y) * Math.cos(eqd_sphe.lon * DEG2RAD - omega)) /
      Math.cos(b * DEG2RAD)
  );

  return new LibrationInfo(
    b,
    l,
    RAD2DEG * mlat,
    RAD2DEG * mlon,
    dist_km,
    diam_deg,
    P * RAD2DEG,
    sub_sun.b,
    sub_sun.l
  );
}

[cosinekitty]

This is very cool! Thank you for posting it. I'm not in a position right now to study this carefully, but I feel @ebraminio might find this of interest, especially for rendering images of the Moon in the Persian calendar.

[ebraminio]

It's a fantastic work @YingChangJ and thanks for the ping @cosinekitty it's indeed very educational for me.