geoidIntersect function

GlobalSolarEclipseInfo geoidIntersect(ShadowInfo shadow)

Implementation

GlobalSolarEclipseInfo geoidIntersect(ShadowInfo shadow) {
  var kind = EclipseKind.Partial;
  var peak = shadow.time;
  var distance = shadow.r;
  dynamic latitude; // left undefined for partial eclipses
  dynamic longitude; // left undefined for partial eclipses

  // We want to calculate the intersection of the shadow axis with the Earth's geoid.
  // First we must convert EQJ (equator of J2000) coordinates to EQD (equator of date)
  // coordinates that are perfectly aligned with the Earth's equator at this
  // moment in time.
  final rot = RotationMatrix.rotationEQJtoEQD(shadow.time);
  final v = AstroVector.rotateVector(
      rot, shadow.dir); // shadow-axis vector in equator-of-date coordinates
  final e = AstroVector.rotateVector(
      rot, shadow.target); // lunacentric Earth in equator-of-date coordinates

  // Convert all distances from AU to km.
  // But dilate the z-coordinates so that the Earth becomes a perfect sphere.
  // Then find the intersection of the vector with the sphere.
  // See p 184 in Montenbruck & Pfleger's "Astronomy on the Personal Computer", second edition.
  v.x *= KM_PER_AU;
  v.y *= KM_PER_AU;
  v.z *= KM_PER_AU / EARTH_FLATTENING;
  e.x *= KM_PER_AU;
  e.y *= KM_PER_AU;
  e.z *= KM_PER_AU / EARTH_FLATTENING;

  // Solve the quadratic equation that finds whether and where
  // the shadow axis intersects with the Earth in the dilated coordinate system.
  final R = EARTH_EQUATORIAL_RADIUS_KM;
  final A = v.x * v.x + v.y * v.y + v.z * v.z;
  final B = -2.0 * (v.x * e.x + v.y * e.y + v.z * e.z);
  final C = (e.x * e.x + e.y * e.y + e.z * e.z) - R * R;
  final raDic = B * B - 4 * A * C;

  double? obscuration;

  if (raDic > 0.0) {
    // Calculate the closer of the two intersection points.
    // This will be on the day side of the Earth.
    final u = (-B - sqrt(raDic)) / (2 * A);

    // Convert lunacentric dilated coordinates to geocentric coordinates.
    final px = u * v.x - e.x;
    final py = u * v.y - e.y;
    final pz = (u * v.z - e.z) * EARTH_FLATTENING;

    // Convert cartesian coordinates into geodetic latitude/longitude.
    final proj = hypot(px, py) * EARTH_FLATTENING_SQUARED;
    if (proj == 0.0) {
      latitude = (pz > 0.0) ? 90.0 : -90.0;
    } else {
      latitude = RAD2DEG * atan(pz / proj);
    }

    // Adjust longitude for Earth's rotation at the given UT.
    final gast = _siderealTime(peak);
    longitude = (RAD2DEG * atan2(py, px) - (15 * gast)) % 360.0;
    if (longitude <= -180.0) {
      longitude += 360.0;
    } else if (longitude > 180.0) {
      longitude -= 360.0;
    }

    // We want to determine whether the observer sees a total eclipse or an annular eclipse.
    // We need to perform a series of vector calculations...
    // Calculate the inverse rotation matrix, so we can convert EQD to EQJ.
    final inv = RotationMatrix.inverseRotation(rot);

    // Put the EQD geocentric coordinates of the observer into the vector 'o'.
    // Also convert back from kilometers to astronomical units.
    var o = AstroVector(
        px / KM_PER_AU, py / KM_PER_AU, pz / KM_PER_AU, shadow.time);

    // Rotate the observer's geocentric EQD back to the EQJ system.
    o = AstroVector.rotateVector(inv, o);

    // Convert geocentric vector to lunacentric vector.
    o.x += shadow.target.x;
    o.y += shadow.target.y;
    o.z += shadow.target.z;

    // Recalculate the shadow using a vector from the Moon's center toward the observer.
    final surface =
        ShadowInfo.calcShadow(MOON_POLAR_RADIUS_KM, shadow.time, o, shadow.dir);

    // If we did everything right, the shadow distance should be very close to zero.
    // That's because we already determined the observer 'o' is on the shadow axis!
    if (surface.r > 1.0e-9 || surface.r < 0.0) {
      throw 'Unexpected shadow distance from geoid intersection = ${surface.r}';
    }

    kind = eclipseKindFromUmbra(surface.k);
    obscuration = (kind == EclipseKind.Total)
        ? 1.0
        : solarEclipseObscuration(shadow.dir, o);
  } else {
    // This is a partial solar eclipse. It does not make practical sense to calculate obscuration.
    // Anyone who wants obscuration should use Astronomy.searchLocalSolarEclipse for a specific location on the Earth.
    obscuration = null;
  }

  return GlobalSolarEclipseInfo(
      kind, obscuration, peak, distance, latitude, longitude);
}