module BoundaryCondition
using CairoMakie
using StaticArrays
using ForwardDiff
using LinearAlgebra

function ϕ(x, R, ε, m)
    r = norm(x)
    θ = atan(x[2], x[1])
    return r - rstar(θ, R, ε, m)
end

function iso_d_point(θ, R, ε, m, d)
    x = rstar(θ, R, ε, m) .* SA[cos(θ), sin(θ)]
    return iso_d_point(x, R, ε, m, d)
end

function iso_d_point(x::AbstractVector, R, ε, m, d)
    ν = compute_ν(x, R, ε, m, d)
    return x + d * ν
end

function compute_ν(θ, R, ε, m, d)
    x = rstar(θ, R, ε, m) .* SA[cos(θ), sin(θ)]
    return compute_ν(x, R, ε, m, d)
end

function compute_ν(x::AbstractVector, R, ε, m, d)
    ν = ForwardDiff.gradient(_x -> ϕ(_x, R, ε, m), x)
    ν = normalize(ν)
    return ν
end

rstar(θ, R, ε, m) = R * (1 + ε * cos(m * θ))

function boundary_condition()

    # Settings
    R = 1
    ε = 0.044
    m = 7
    xmin = 0
    xmax = 2R
    ymin = 0.0
    ymax = xmax

    fontsize = 28
    fig = Figure(; size=(800, 500), fontsize)


    ax = Axis(
        fig[1, 1],
        limits=(xmin, xmax, ymin, ymax),
        leftspinevisible=false,
        rightspinevisible=false,
        bottomspinevisible=false,
        topspinevisible=false,
        aspect=1,
        # backgroundcolor=:transparent
    )
    hidedecorations!(ax)
    deregister_interaction!(ax, :rectanglezoom)
    deregister_interaction!(ax, :dragpan)
    deregister_interaction!(ax, :scrollzoom)
    deregister_interaction!(ax, :limitreset)

    # The mesh
    nx = ny = 10
    for i in range(xmin, xmax, length=nx)
        lines!(ax, [i, i], [ymin, ymax], color=:gray, linestyle=:dot)
    end
    for j in range(ymin, ymax, length=ny)
        lines!(ax, [xmin, xmax], [j, j], color=:gray, linestyle=:dot)
    end

    # The star
    θ = LinRange(0, 1.5 * π / 2, 100)
    color = colorant"#1154a6"
    points_d0 = [Point2f(rstar(_θ, R, ε, m) .* [cos(_θ), sin(_θ)]...) for _θ in θ]
    lines!(points_d0; color, linewidth=3)

    # Draw some iso-d
    color = :red
    linestyle = :dash
    linewidth = 2
    for d in (-0.3, 0.4, 0.8, 1.2)
        points = [Point2f(iso_d_point(_θ, R, ε, m, d)...) for _θ in θ]
        lines!(points; color, linestyle, linewidth)
    end

    # Projections
    linestyle = (:dot, :dense)
    θ_proj = (deg2rad(14), deg2rad(58))
    color = :green
    linewidth = 2
    for (i, _θ) in enumerate(θ_proj)
        ν = compute_ν(_θ, R, ε, m, 0.0)
        A = Point2f(rstar(_θ, R, ε, m) .* [cos(_θ), sin(_θ)])
        points = [A, Point2f(iso_d_point(_θ, R, ε, m, 0.6xmax))]
        lines!(points; linestyle, color, linewidth)
        l = 0.05
        B = Point2f(A.data .+ l * ν)
        u = SA[-ν[2], ν[1]]
        C = Point2f(B.data .+ l * normalize(u))
        D = Point2f(C.data .- l * ν)
        lines!([A, B, C, D], color=:black, linewidth=0.7)
        text!(points[2], text=L"x_{%$i}", align=(:center, :top), offset=(-15, 0))
        text!(points[1], text=L"p(x_{%$i})", align=(:right, :top))
    end

    # The rectangle
    color = :black
    linewidth = 8
    lines!([xmax, xmax], [ymin, ymax]; color, linewidth)
    lines!([xmax, xmin], [ymax, ymax]; color, linewidth)
    # poly!(ax, [xmin, xmax, xmax, xmin], [ymin, ymin, ymax, ymax], strokewidth=2, color=:transparent, strokecolor=:black)


    display(fig)
    foreach(ext -> save(joinpath(@__DIR__, "boundary-condition.$ext"), fig), ("png", "svg", "pdf"))

end

boundary_condition()

end