import type { Meters } from '../units';

/**
 * Umrechnung geografischer Koordinaten in ein metrisches Bezugssystem.
 *
 * WARUM NICHT WEB-MERCATOR
 *
 * Die naheliegende Wahl waere EPSG:3857 (Web-Mercator), das jeder Kartendienst
 * spricht. Dort sind Karteneinheiten aber KEINE Meter: Bei der geografischen
 * Breite phi entspricht eine Karteneinheit nur cos(phi) Metern am Boden. In
 * Deutschland - etwa bei 52 Grad Nord - sind das 0,62 m. Ein daraus
 * abgegriffener Raeumweg waere um mehr als 60 Prozent zu gross.
 *
 * Fuer die Vermassung wird deshalb ETRS89 / UTM Zone 32N (EPSG:25832) genutzt.
 * Das ist zugleich das Bezugssystem der deutschen Landesvermessung; die
 * Karteneinheiten sind Meter, der Massstabsfehler liegt bei hoechstens etwa
 * 0,4 Promille - also unter einem Millimeter je Meter und damit weit unterhalb
 * der Ablesegenauigkeit eines Lageplans.
 *
 * Zone 32N deckt 6 bis 12 Grad oestlicher Laenge ab und damit den groessten Teil
 * Deutschlands. Fuer Gebiete weiter westlich oder oestlich stehen die Zonen 31
 * und 33 zur Verfuegung.
 */

/** Geografische Koordinate in Grad. */
export interface GeoPunkt {
  readonly breite: number;
  readonly laenge: number;
}

/** Koordinate in einem metrischen Bezugssystem. */
export interface UtmPunkt {
  readonly ost: Meters;
  readonly nord: Meters;
  readonly zone: number;
}

// WGS84 bzw. GRS80 - die Unterschiede liegen im Zentimeterbereich und sind
// fuer die Vermassung eines Knotenpunkts ohne Belang.
const A = 6378137.0;
const F = 1 / 298.257223563;
const K0 = 0.9996;
const FALSE_EAST = 500000;

const E2 = F * (2 - F);
const EP2 = E2 / (1 - E2);

/** UTM-Zone zu einer geografischen Laenge. */
export function utmZone(laenge: number): number {
  return Math.floor((laenge + 180) / 6) + 1;
}

/** EPSG-Code der ETRS89/UTM-Zone. */
export function epsgFuerZone(zone: number): number {
  return 25800 + zone;
}

/**
 * Rechnet eine geografische Koordinate in UTM um.
 * Reihenentwicklung nach Krueger; fuer Abstaende bis einige Kilometer vom
 * Mittelmeridian auf Millimeter genau.
 */
export function nachUtm(punkt: GeoPunkt, zone = utmZone(punkt.laenge)): UtmPunkt {
  const phi = (punkt.breite * Math.PI) / 180;
  const lambda = (punkt.laenge * Math.PI) / 180;
  const lambda0 = (((zone - 1) * 6 - 180 + 3) * Math.PI) / 180;

  const sinPhi = Math.sin(phi);
  const cosPhi = Math.cos(phi);
  const tanPhi = Math.tan(phi);

  const N = A / Math.sqrt(1 - E2 * sinPhi * sinPhi);
  const T = tanPhi * tanPhi;
  const C = EP2 * cosPhi * cosPhi;
  const Aa = (lambda - lambda0) * cosPhi;

  const M =
    A *
    ((1 - E2 / 4 - (3 * E2 * E2) / 64 - (5 * E2 * E2 * E2) / 256) * phi -
      ((3 * E2) / 8 + (3 * E2 * E2) / 32 + (45 * E2 * E2 * E2) / 1024) * Math.sin(2 * phi) +
      ((15 * E2 * E2) / 256 + (45 * E2 * E2 * E2) / 1024) * Math.sin(4 * phi) -
      ((35 * E2 * E2 * E2) / 3072) * Math.sin(6 * phi));

  const ost =
    FALSE_EAST +
    K0 *
      N *
      (Aa +
        ((1 - T + C) * Aa ** 3) / 6 +
        ((5 - 18 * T + T * T + 72 * C - 58 * EP2) * Aa ** 5) / 120);

  const nord =
    K0 *
    (M +
      N *
        tanPhi *
        ((Aa * Aa) / 2 +
          ((5 - T + 9 * C + 4 * C * C) * Aa ** 4) / 24 +
          ((61 - 58 * T + T * T + 600 * C - 330 * EP2) * Aa ** 6) / 720));

  return { ost, nord, zone };
}

/** Umkehrung - fuer die Anzeige und zur Pruefung der Hinrechnung. */
export function nachGeo(punkt: UtmPunkt): GeoPunkt {
  const x = punkt.ost - FALSE_EAST;
  const y = punkt.nord;
  const lambda0 = (((punkt.zone - 1) * 6 - 180 + 3) * Math.PI) / 180;

  const e1 = (1 - Math.sqrt(1 - E2)) / (1 + Math.sqrt(1 - E2));
  const M = y / K0;
  const mu = M / (A * (1 - E2 / 4 - (3 * E2 * E2) / 64 - (5 * E2 * E2 * E2) / 256));

  const phi1 =
    mu +
    ((3 * e1) / 2 - (27 * e1 ** 3) / 32) * Math.sin(2 * mu) +
    ((21 * e1 * e1) / 16 - (55 * e1 ** 4) / 32) * Math.sin(4 * mu) +
    ((151 * e1 ** 3) / 96) * Math.sin(6 * mu) +
    ((1097 * e1 ** 4) / 512) * Math.sin(8 * mu);

  const sinPhi1 = Math.sin(phi1);
  const cosPhi1 = Math.cos(phi1);
  const tanPhi1 = Math.tan(phi1);

  const N1 = A / Math.sqrt(1 - E2 * sinPhi1 * sinPhi1);
  const T1 = tanPhi1 * tanPhi1;
  const C1 = EP2 * cosPhi1 * cosPhi1;
  const R1 = (A * (1 - E2)) / (1 - E2 * sinPhi1 * sinPhi1) ** 1.5;
  const D = x / (N1 * K0);

  const phi =
    phi1 -
    ((N1 * tanPhi1) / R1) *
      ((D * D) / 2 -
        ((5 + 3 * T1 + 10 * C1 - 4 * C1 * C1 - 9 * EP2) * D ** 4) / 24 +
        ((61 + 90 * T1 + 298 * C1 + 45 * T1 * T1 - 252 * EP2 - 3 * C1 * C1) * D ** 6) / 720);

  const lambda =
    lambda0 +
    (D -
      ((1 + 2 * T1 + C1) * D ** 3) / 6 +
      ((5 - 2 * C1 + 28 * T1 - 3 * C1 * C1 + 8 * EP2 + 24 * T1 * T1) * D ** 5) / 120) /
      cosPhi1;

  return { breite: (phi * 180) / Math.PI, laenge: (lambda * 180) / Math.PI };
}

/** Ausschnitt in Kartenkoordinaten - Reihenfolge wie im WMS-Parameter BBOX. */
export interface Ausschnitt {
  readonly minOst: Meters;
  readonly minNord: Meters;
  readonly maxOst: Meters;
  readonly maxNord: Meters;
  readonly zone: number;
}

/**
 * Quadratischer Ausschnitt um einen Punkt mit vorgegebener Kantenlaenge.
 * Weil in UTM gerechnet wird, ist die Kantenlaenge tatsaechlich in Metern am
 * Boden zu verstehen.
 */
export function ausschnittUm(mitte: GeoPunkt, kantenlaenge: Meters): Ausschnitt {
  const utm = nachUtm(mitte);
  const halb = kantenlaenge / 2;
  return {
    minOst: utm.ost - halb,
    minNord: utm.nord - halb,
    maxOst: utm.ost + halb,
    maxNord: utm.nord + halb,
    zone: utm.zone,
  };
}

/**
 * Teilausschnitt eines in Kacheln zerlegten Ausschnitts.
 *
 * WOZU
 *
 * Ein Kartendienst rechnet ein grosses Bild in einem Stueck aus, und der
 * Aufwand waechst mit der Flaeche; jenseits einiger Millionen Bildpunkte
 * bricht der Abruf ab oder laeuft in die Zeitgrenze. Ein grosser Ausschnitt
 * wird deshalb kachelweise geholt und danach zusammengesetzt. Diese Funktion
 * liefert den Kartenausschnitt EINER Kachel.
 *
 * Zerlegt wird in `spalten` mal `zeilen` gleich grosse Kacheln. Gezaehlt wird
 * wie im Bild: `spalte` 0 liegt links, `zeile` 0 liegt OBEN.
 *
 * DIE STELLE, AN DER ES SCHIEFGEHT
 *
 * Bildzeile 0 liegt oben, also am NOERDLICHEN Rand des Ausschnitts. Die
 * Nordwerte laufen der Bildzeile damit ENTGEGEN: Je groesser die Zeile, desto
 * kleiner der Nordwert. Wird das vertauscht, ist das zusammengesetzte Mosaik
 * senkrecht gespiegelt - Zeile fuer Zeile in sich richtig, in der Reihenfolge
 * aber verkehrt.
 *
 * Auffallen wuerde das niemandem: Ein Luftbild hat keine Leserichtung. Der
 * Anwender zeichnet seine Fahrlinien in ein Bild, in dem die noerdliche Zufahrt
 * unten liegt, misst darin voellig plausible Laengen - und vermasst dabei die
 * falschen Wege. Deshalb wird hier von `maxNord` nach unten gerechnet und nicht
 * von `minNord` nach oben.
 *
 * Unbrauchbare Angaben werden begrenzt statt zurueckgewiesen: Der Rueckgabewert
 * ist immer ein echter Teil des Gesamtausschnitts, nie einer daneben. Eine
 * Kachel ausserhalb des Bildes waere ein Bildteil aus einer fremden Gegend -
 * und die faellt im Luftbild ebenso wenig auf wie die Spiegelung.
 */
export function unterAusschnitt(
  gesamt: Ausschnitt,
  spalte: number,
  zeile: number,
  spalten: number,
  zeilen: number,
): Ausschnitt {
  const anzahlSpalten = anzahl(spalten);
  const anzahlZeilen = anzahl(zeilen);
  const s = index(spalte, anzahlSpalten);
  const z = index(zeile, anzahlZeilen);

  return {
    // Ost laeuft mit der Spalte: links ist Westen.
    minOst: teile(gesamt.minOst, gesamt.maxOst, s, anzahlSpalten),
    maxOst: teile(gesamt.minOst, gesamt.maxOst, s + 1, anzahlSpalten),
    // Nord laeuft der Zeile entgegen: Zeile 0 ist der NOERDLICHE Rand.
    maxNord: teile(gesamt.maxNord, gesamt.minNord, z, anzahlZeilen),
    minNord: teile(gesamt.maxNord, gesamt.minNord, z + 1, anzahlZeilen),
    zone: gesamt.zone,
  };
}

/**
 * Kachelgrenze als Anteil zwischen zwei Raendern.
 *
 * Bewusst als Anteil des Ganzen gerechnet und nicht als "Anfang plus i mal
 * Kachelbreite": So trifft die letzte Kante den Rand des Gesamtausschnitts
 * genau. Aufsummierte Kachelbreiten liessen dort einen Spalt von der Groesse
 * des aufgelaufenen Rundungsfehlers - im Mosaik ein Versatz zwischen den
 * Kacheln, der sich als scheinbarer Zeichenfehler auf jeden Weg legt.
 */
function teile(von: Meters, bis: Meters, i: number, kacheln: number): Meters {
  return von + (bis - von) * (i / kacheln);
}

/** Kachelzahl je Richtung; weniger als eine Kachel gibt es nicht. */
function anzahl(wert: number): number {
  if (!Number.isFinite(wert)) return 1;
  return Math.max(1, Math.floor(wert));
}

/** Kachelnummer, auf das vorhandene Raster begrenzt. */
function index(wert: number, kacheln: number): number {
  if (!Number.isFinite(wert)) return 0;
  return Math.min(kacheln - 1, Math.max(0, Math.floor(wert)));
}

/** BBOX-Parameter fuer einen WMS-Aufruf. */
export function bboxParameter(ausschnitt: Ausschnitt): string {
  return [ausschnitt.minOst, ausschnitt.minNord, ausschnitt.maxOst, ausschnitt.maxNord]
    .map((wert) => wert.toFixed(3))
    .join(',');
}

/** Kantenlaenge des Ausschnitts in Metern. */
export function kantenlaenge(ausschnitt: Ausschnitt): Meters {
  return ausschnitt.maxOst - ausschnitt.minOst;
}

/** Eine Zahl mit oder ohne Nachkommastellen, Punkt oder Komma. */
const ZAHL = String.raw`-?\d+(?:[.,]\d+)?`;

/**
 * Liest eine Koordinate aus einer eingefuegten Zeichenkette.
 *
 * Erkannt werden:
 *   - "50.9412784, 6.9582814"
 *   - "50,9412784 6,9582814"                     (deutsches Dezimalkomma)
 *   - Google Maps      "@50.9412784,6.9582814"
 *   - Apple Karten     "?ll=50.9412784,6.9582814"
 *   - OpenStreetMap    "#map=19/50.9412784/6.9582814"
 *   - OSM mit Marke    "?mlat=50.9412784&mlon=6.9582814"
 *   - Bing Maps        "?cp=50.9412784~6.9582814"
 *   - "geo:50.9412784,6.9582814"
 *
 * Der Weg ueber den eingefuegten Kartenverweis ist bewusst vorgesehen: Wer eine
 * Kreuzung sucht, findet sie meist zuerst in einem Kartendienst und kann den
 * Verweis dann einfach uebernehmen.
 *
 * OpenStreetMap kam zunaechst nicht vor. Fuer eine Anwendung, deren Ergebnisse
 * in behoerdliche Verfahren gehen, ist gerade das der naheliegende Dienst - und
 * seine Adresse traegt die Koordinate im Fragment ("#map="), das keiner der
 * uebrigen Muster erfasste. Die Zahlen dort haben nicht zwingend
 * Nachkommastellen: "#map=5/51.5/11" ist eine gueltige Adresse.
 */
export function leseKoordinate(text: string): GeoPunkt | null {
  const roh = text.trim();
  if (roh === '') return null;

  /*
   * Reihenfolge mit Absicht: Erst die Muster, die ein Schluesselwort tragen,
   * dann die beiden blanken Zahlen. Ein Verweis enthaelt viele Zahlenpaare -
   * Zoomstufe, Bildgroesse, Zeitstempel -, und ohne Schluesselwort waere nicht
   * zu entscheiden, welches die Koordinate ist.
   */
  const muster: readonly RegExp[] = [
    // Google Maps: der Teil nach dem At-Zeichen.
    new RegExp(String.raw`@(${ZAHL}),\s*(${ZAHL})`),
    // OpenStreetMap und Abkoemmlinge: Zoomstufe, Breite, Laenge im Fragment.
    new RegExp(String.raw`#map=[\d.]+/(${ZAHL})/(${ZAHL})`, 'i'),
    // OpenStreetMap mit gesetzter Marke.
    new RegExp(String.raw`[?&]mlat=(${ZAHL})&mlon=(${ZAHL})`, 'i'),
    // Bing Maps: Mittelpunkt mit Tilde getrennt.
    new RegExp(String.raw`[?&]cp=(${ZAHL})~(${ZAHL})`, 'i'),
    // Ausdrueckliche Parameter und geo:
    new RegExp(String.raw`(?:^geo:|[?&](?:q|ll|center)=)(${ZAHL})[,;]\s*(${ZAHL})`, 'i'),
    // Zwei blanke Zahlen - nur, wenn sonst nichts dasteht.
    new RegExp(String.raw`^(${ZAHL})[\s,;]+(${ZAHL})$`),
  ];

  for (const regel of muster) {
    const treffer = regel.exec(roh);
    if (treffer) return pruefe(zahl(treffer[1]), zahl(treffer[2]));
  }

  return null;
}

function zahl(wert: string | undefined): number {
  return Number((wert ?? '').replace(',', '.'));
}

function pruefe(breite: number, laenge: number): GeoPunkt | null {
  if (!Number.isFinite(breite) || !Number.isFinite(laenge)) return null;
  if (Math.abs(breite) > 90 || Math.abs(laenge) > 180) return null;
  return { breite, laenge };
}