trails/packages/gpx/src/compute-days.ts

import type { Waypoint } from "@trails-cool/types";
import type { TrackPoint } from "./types.ts";
import { despike, cumulativeFilteredTotals, type ElevPoint } from "./elevation-clean.ts";

export interface DayStage {
  dayNumber: number;
  startWaypointIndex: number;
  endWaypointIndex: number;
  startName?: string;
  endName?: string;
  /** Distance in meters */
  distance: number;
  /** Ascent in meters */
  ascent: number;
  /** Descent in meters */
  descent: number;
}

/**
 * Split a route into day stages based on waypoints with isDayBreak.
 *
 * Day boundaries are defined by waypoints where isDayBreak is true.
 * The first waypoint is the implicit start of Day 1, the last waypoint
 * is the implicit end of the final day, and each isDayBreak waypoint
 * marks the end of one day and the start of the next.
 *
 * If no waypoints have isDayBreak, returns a single day covering the
 * entire route.
 */
export function computeDays(
  waypoints: Waypoint[],
  tracks: TrackPoint[][],
): DayStage[] {
  if (waypoints.length === 0 || tracks.length === 0) return [];

  // Flatten all tracks into a single point array
  const allPoints: TrackPoint[] = tracks.flat();
  if (allPoints.length === 0) return [];

  // Find waypoint indices that are day breaks
  const breakIndices: number[] = [];
  for (let i = 0; i < waypoints.length; i++) {
    if (waypoints[i]!.isDayBreak) {
      breakIndices.push(i);
    }
  }

  // Build day boundaries: [0, break1, break2, ..., lastWaypointIndex]
  const boundaries = [0, ...breakIndices];
  if (boundaries[boundaries.length - 1] !== waypoints.length - 1) {
    boundaries.push(waypoints.length - 1);
  }

  // For each waypoint, find the closest point in the track
  const waypointTrackIndices = waypoints.map((wp) => {
    let bestIdx = 0;
    let bestDist = Infinity;
    for (let i = 0; i < allPoints.length; i++) {
      const dx = allPoints[i]!.lat - wp.lat;
      const dy = allPoints[i]!.lon - wp.lon;
      const d = dx * dx + dy * dy;
      if (d < bestDist) {
        bestDist = d;
        bestIdx = i;
      }
    }
    return bestIdx;
  });

  // Cumulative distance per flat point index.
  const cumDist: number[] = [0];
  for (let i = 1; i < allPoints.length; i++) {
    const prev = allPoints[i - 1]!;
    const curr = allPoints[i]!;
    cumDist.push(cumDist[i - 1]! + haversine(prev.lat, prev.lon, curr.lat, curr.lon));
  }

  // Cumulative ascent/descent via the shared cleaning path (spec:
  // elevation-profile-hardening) instead of summing every raw delta: despike
  // each segment, then run the hysteresis accumulator once over the despiked
  // track. Day totals (cumulative[end] − cumulative[start]) therefore match
  // the filtered route total and won't overstate by 20–50%.
  const eleSeq: ElevPoint[] = [];
  const eleFlatIdx: number[] = [];
  let flat = 0;
  for (const seg of tracks) {
    let segCum = 0;
    const segPts: ElevPoint[] = [];
    const segIdx: number[] = [];
    for (let i = 0; i < seg.length; i++) {
      const p = seg[i]!;
      if (i > 0) segCum += haversine(seg[i - 1]!.lat, seg[i - 1]!.lon, p.lat, p.lon);
      if (p.ele !== undefined) {
        segPts.push({ distance: segCum, elevation: p.ele });
        segIdx.push(flat);
      }
      flat++;
    }
    const cleaned = despike(segPts);
    for (let i = 0; i < cleaned.length; i++) {
      eleSeq.push(cleaned[i]!);
      eleFlatIdx.push(segIdx[i]!);
    }
  }
  const { ascent, descent } = cumulativeFilteredTotals(eleSeq);
  // Map running totals onto every flat point index, carrying the last value
  // forward across points that carry no elevation.
  const cumAscent: number[] = new Array(allPoints.length);
  const cumDescent: number[] = new Array(allPoints.length);
  let sp = 0;
  let curA = 0;
  let curD = 0;
  for (let i = 0; i < allPoints.length; i++) {
    while (sp < eleFlatIdx.length && eleFlatIdx[sp] === i) {
      curA = ascent[sp]!;
      curD = descent[sp]!;
      sp++;
    }
    cumAscent[i] = curA;
    cumDescent[i] = curD;
  }

  // Build day stages
  const days: DayStage[] = [];
  for (let d = 0; d < boundaries.length - 1; d++) {
    const startWpIdx = boundaries[d]!;
    const endWpIdx = boundaries[d + 1]!;
    const startTrackIdx = waypointTrackIndices[startWpIdx]!;
    const endTrackIdx = waypointTrackIndices[endWpIdx]!;

    days.push({
      dayNumber: d + 1,
      startWaypointIndex: startWpIdx,
      endWaypointIndex: endWpIdx,
      startName: waypoints[startWpIdx]!.name,
      endName: waypoints[endWpIdx]!.name,
      distance: Math.round(cumDist[endTrackIdx]! - cumDist[startTrackIdx]!),
      ascent: Math.round(cumAscent[endTrackIdx]! - cumAscent[startTrackIdx]!),
      descent: Math.round(cumDescent[endTrackIdx]! - cumDescent[startTrackIdx]!),
    });
  }

  return days;
}

function haversine(lat1: number, lon1: number, lat2: number, lon2: number): number {
  const R = 6371000;
  const toRad = (d: number) => (d * Math.PI) / 180;
  const dLat = toRad(lat2 - lat1);
  const dLon = toRad(lon2 - lon1);
  const a =
    Math.sin(dLat / 2) ** 2 +
    Math.cos(toRad(lat1)) * Math.cos(toRad(lat2)) * Math.sin(dLon / 2) ** 2;
  return R * 2 * Math.atan2(Math.sqrt(a), Math.sqrt(1 - a));
}