/*
  Copyright (c) 2021 Corinne Diakhoumpa et Erwan Iev-Le Tac

  Permission is hereby granted, free of charge, to any person obtaining a copy
  of this software and associated documentation files (the "Software"), to deal
  in the Software without restriction, including without limitation the rights
  to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
  copies of the Software, and to permit persons to whom the Software is
  furnished to do so, subject to the following conditions:

  The above copyright notice and this permission notice shall be included in all
  copies or substantial portions of the Software.

  THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
  IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
  FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
  AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
  LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
  OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
  SOFTWARE.
*/

const fs = require("fs");
const assert = require('assert').strict;
const utils = require('./utils');

// This program draws the callidiagramme from "Géodésique".
// Definitely not perfect, but hopefully good enough...

const r = 200;
const epsilon = 0.001;
const arcSegmentPixelsLength = 5;

// Determine the fundamental domain of Bolza's surface
// https://en.wikipedia.org/wiki/Bolza_surface#Triangle_surface
let vertices = [];
for (let k = 0; k < 9; k++) {
    let pho = r / Math.pow(2, .25);
    let arg = Math.PI * (1/8 + (k-1)/4);
    vertices.push(scale(pho, {x: Math.cos(arg), y: Math.sin(arg)}));
}
let edges = [];
for (let k = 0; k < 8; k++)
    edges.push(geodesicSegment(vertices[k % 8], vertices[(k+1)% 8]));

function scale(lambda, v) { return { x: lambda * v.x, y: lambda * v.y}; }
function sum(v, w) { return { x: v.x + w.x, y: v.y + w.y }; }
function diff(v, w) { return sum(v, scale(-1, w)); }
function arg(v) { return Math.atan2(v.y, v.x); }
function norm(v) { return Math.hypot(v.x, v.y); }
function dot(v, w) { return v.x * w.x + v.y * w.y; }
function det(v, w) { return v.x * w.y - w.x * v.y; }
function vector(P, Q) { return diff(Q, P); }
function inversion(P) { return scale(r ** 2 / norm(P) ** 2, P); }
function midPoint(P, Q) { return scale(.5, sum(P, Q)); }
function angle(v, w) {
    let c = dot(v, w)/(norm(v) * norm(w));
    c = Math.max(-1, Math.min(1, c));
    return Math.sign(det(v, w)) * Math.acos(c);
}

function geodesicSegment(P, Q) {
    // https://en.wikipedia.org/wiki/Poincar%C3%A9_disk_model#Compass_and_straightedge_construction
    if (Math.abs(det(P, Q)) < epsilon) {
        // Geodesic is a diameter.
        let v = vector(P, Q);
        return {
            diameter: true,
            start: P,
            end: Q,
            vector: scale(r / norm(v), v),
        };
    }

    // Determine the circular arc.
    let P2 = inversion(P);
    let Q2 = inversion(Q);
    let M = midPoint(P, P2);
    let N = midPoint(Q, Q2);

    let Arow1 = vector(P, P2);
    let Arow2 = vector(Q, Q2);
    let B1 = dot(Arow1, M);
    let B2 = dot(Arow2, N);
    let detA = det(Arow1, Arow2);
    let C = scale(1/detA, {
        x: Arow2.y * B1 - Arow1.y * B2,
        y: -Arow2.x * B1 + Arow1.x * B2
    });

    // Geodesic is the shortest path.
    let CP = vector(C, P);
    return {
        diameter: false,
        start: P,
        end: Q,
        center: C,
        radius: norm(CP)
    };
}

function geodesicWith(P, v) {
    function intersectionsWith(g) {
        let result = [];
        for (let k = 0; k < 8; k++) {
            let I = intersection(g, edges[k]);
            if (I) result.push(I);
        }
        return result;
    }
    function isFarEnough(g, e) {
        // If g.start is close to the bound, choose another end. This is to
        // workaround rounding errors....
        return norm(vector(g.start, e)) > epsilon;
    }
    if (Math.abs(det(P, v)) < epsilon) {
        let g = {
            diameter: true,
            vector: scale(r / norm(v), v),
            start: P
        };
        let I = intersectionsWith(g);
        let i = I.findIndex(e => {
            return dot(vector(P, e), v) >= 0 && isFarEnough(g, e);
        });
        assert(i >= 0);
        g.end = I[i];
        return g;
    }

    // Determine center of circular arc.
    let Arow1 = scale(2, P);
    let Arow2 = v;
    let B1 = r ** 2 + dot(P, P);
    let B2 = dot(P, v);
    let detA = det(Arow1, Arow2);
    let C = scale(1/detA, {
        x: Arow2.y * B1 - Arow1.y * B2,
        y: -Arow2.x * B1 + Arow1.x * B2
    });

    let CP = vector(C, P);
    let g = {
        diameter: false,
        center: C,
        radius: norm(CP),
        start: P,
    };
    let I = intersectionsWith(g);
    let i = I.findIndex(e => {
        let theta_delta = angle(CP, vector(C, e));
        return det(CP, v) * theta_delta >= 0 && isFarEnough(g, e);
    });
    assert(i >= 0);
    g.end = I[i];
    return g;
}

function tangentVector(geodesic, P) {
    if (geodesic.diameter)
        return geodesic.vector;

    // Geodesic parametrization: center + radius * e^{sign * i t}
    var sign = geodesic.end > geodesic.start ? 1 : -1;
    var CP = vector(geodesic.center, P);
    var t = Math.acos(CP.x / geodesic.radius);
    if (CP.y * Math.sin(sign * t) < 0)
        t = -t;
    return scale(sign * r,
                 {x: -Math.sin(sign * t), y: Math.cos(sign * t)});
}

function intersection(g1, g2) {
    function intersectionsLineAndCircle(a, b, c, center, r) {
        // https://en.wikipedia.org/wiki/Intersection_%28Euclidean_geometry%29#A_line_and_a_circle
        c -= a*center.x + b*center.y;
        let D = r**2 * (a**2 + b**2) - c**2;
        if (D < 0)
            return [];
        let root = Math.sqrt(D);
        return [
            sum(center,
                scale(1/(a**2 + b**2), { x: a*c + b*root, y : b*c - a*root })),
            sum(center,
                scale(1/(a**2 + b**2), { x: a*c - b*root, y : b*c + a*root }))
        ];
    }

    if (g1.diameter && g2.diameter) {
        // Diameters intersect at the center.
        return {x: 0, y: 0};
    }

    if (g2.diameter) {
        // Invert geodesics to ensure the second is not a diameter.
        return getIntersection(g2, g1);
    }

    let I;
    if (g1.diameter) {
        // Diameter and circle.
        I = intersectionsLineAndCircle(g1.vector.y,
                                       -g1.vector.x,
                                       0,
                                       g2.center,
                                       g2.radius);
    } else {
        // Two circles.
        // https://en.wikipedia.org/wiki/Intersection_%28Euclidean_geometry%29#Two_circles
        let C1C2 = vector(g1.center, g2.center);
        let c = g1.radius ** 2 - norm(g1.center)**2 -
            (g2.radius ** 2 - norm(g2.center)**2);
        I = intersectionsLineAndCircle(2 * C1C2.x,
                                       2 * C1C2.y,
                                       c,
                                       g1.center,
                                       g1.radius);
    }
    return I.length ? (norm(I[0]) < r ? I[0] : I[1]) : null;
}

function pathForSegment(g)
{
    if (g.diameter)
        return `M ${g.start.x},${g.start.y} L${g.end.x},${g.end.y}`;

    let CP = vector(g.center, g.start);
    let CQ = vector(g.center, g.end);
    let theta_start = arg(CP);
    let theta_delta = angle(CP, CQ);
    let arcLength = Math.abs(theta_delta) * g.radius;
    let segmentCount = 1 + Math.ceil(arcLength / arcSegmentPixelsLength);
    let path = "";
    for (let i = 0; i < segmentCount; i++) {
        let theta_i = theta_start + theta_delta * i / (segmentCount - 1);
        let x = g.center.x + g.radius * Math.cos(theta_i);
        let y = g.center.y + g.radius * Math.sin(theta_i);
        path += `${i ? 'L' : 'M'}  ${x},${y} `;
    }
    return path;
}

function moveToOtherSide(stream, g) {
    let P = g.end;
    let k = edges.findIndex(e => {
        return Math.abs(norm(vector(e.center, P)) - e.radius) < epsilon;
    });
    let alpha = angle(tangentVector(edges[k], P), tangentVector(g, P));

    k = (k + 4) % 8;
    if (k % 4 == 0) {
        // Symmetry with respect to x = 0
        P.x = -P.x;
    } else if (k % 4 == 1) {
        // Symmetriy with respect to y = -x.
        let tmp = P.x;
        P.x = -P.y;
        P.y = -tmp;
    } else if (k % 4  == 2) {
        // Symmetry with respect to y = 0
        P.y = -P.y;
    } else {
        assert(k % 4 == 3);
        // Symmetry with respect to x = y.
        let tmp = P.x;
        P.x = P.y
        P.y = tmp;
    }
    let v = tangentVector(edges[k], P);
    v = {
         x: Math.cos(alpha) * v.x - Math.sin(alpha) * v.y,
         y: Math.sin(alpha) * v.x + Math.cos(alpha) * v.y
    };
    if (norm(vector(edges[k].center, sum(P, scale(1/norm(v), v)))) <=
        edges[k].radius + epsilon)
        v = scale(-1, v);
    return geodesicWith(P, v);
}

const filename = "geodesique.svg"
console.log(`Generate ${filename}...`)

let text = "Suivre son chemin malgré les obstacles puis laisser son empreinte dans l’Histoire.";

let stream = fs.createWriteStream(filename);
stream.write(utils.svgProlog(`width="8.5cm" height="8.5cm" viewBox="0 0 ${2*r+10} ${2*r+10}"`));
stream.write(utils.svgHeaderWithCopyLeft("Géodésique", `Un cercle en pointillé représentant le bord du disque de Poincaré. Un octogone régulier représentant le domaine fondamental de la surface de Bolza (visuellement il ressemble à une étoile). Le texte suivant est tracé à l'intérieur, le long d'une géodésique: '${text}'. Lorsque le texte sort par un coté de l'octogone, il rentre à nouveau par le coté opposé.`));

// Draw the path.
stream.write("<defs>");
let position = {x: -r / 2, y: 0}
let theta = Math.PI/4;
let direction = {x: r * Math.cos(theta), y: r * Math.sin(theta)};
g = geodesicWith(position, direction);
g = moveToOtherSide(stream, g);
stream.write(`<path id="path" stroke="black" fill="none" d="${pathForSegment(g)} `);
g = moveToOtherSide(stream, g);
stream.write(`${pathForSegment(g)} `);
g = moveToOtherSide(stream, g);
stream.write(`${pathForSegment(g)} `);
g = moveToOtherSide(stream, g);
stream.write(`${pathForSegment(g)} `);
g = moveToOtherSide(stream, g);
stream.write(`${pathForSegment(g)} `);
g = moveToOtherSide(stream, g);
stream.write(`${pathForSegment(g)} `);
g = moveToOtherSide(stream, g);
stream.write(`${pathForSegment(g)} `);
stream.write(`"/>`);
stream.write("</defs>");

stream.write(`<g style="font-family: \'EB Garamond\'; font-size: 24px" transform="translate(5, 5) rotate(180, ${r}, ${r}) translate(${r}, ${r})">`);

// Draw boundary of Poincaré's disk
stream.write(`<circle stroke="black" fill="none" stroke-dasharray="5" r="${r}"/>`);

// Draw octogon.
for (let k = 0; k < 8; k++) {
    stream.write(`<path stroke="black" fill="none" d="${pathForSegment(edges[k])}"/>`);
}

// Draw the text on the path.
stream.write(`<text><textPath xlink:href="#path">${text}</textPath></text>`)

stream.write("</g>")
stream.write("</svg>")
stream.end();
console.log("done")