From b465ad8b470a5528b4824d95839142ba969581c7 Mon Sep 17 00:00:00 2001
From: csimal <csimal@student.unamur.be>
Date: Wed, 28 Apr 2021 13:42:28 +0200
Subject: [PATCH 1/8] test commit

---
 src/recipes.jl | 1 +
 1 file changed, 1 insertion(+)
 create mode 100644 src/recipes.jl

diff --git a/src/recipes.jl b/src/recipes.jl
new file mode 100644
index 0000000..ad96c72
--- /dev/null
+++ b/src/recipes.jl
@@ -0,0 +1 @@
+# foo
\ No newline at end of file

From 3464d3b1b3f246a627924651dfd786ade5f42246 Mon Sep 17 00:00:00 2001
From: csimal <csimal@student.unamur.be>
Date: Sat, 14 Aug 2021 13:11:26 +0200
Subject: [PATCH 2/8] Proof of Concept Spectral portrait recipe

---
 Project.toml                               |  2 +
 src/plots/plotsrecipes/PSAPlots.jl         | 11 ++++++
 src/plots/plotsrecipes/spectralportrait.jl | 43 ++++++++++++++++++++++
 3 files changed, 56 insertions(+)
 create mode 100644 src/plots/plotsrecipes/PSAPlots.jl
 create mode 100644 src/plots/plotsrecipes/spectralportrait.jl

diff --git a/Project.toml b/Project.toml
index 4cf9c8e..c947f75 100644
--- a/Project.toml
+++ b/Project.toml
@@ -11,6 +11,7 @@ Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 ProgressMeter = "92933f4c-e287-5a05-a399-4b506db050ca"
 Requires = "ae029012-a4dd-5104-9daa-d747884805df"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
+RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
 
 [compat]
 Arpack = "0.3, 0.4, 0.5"
@@ -18,6 +19,7 @@ PlotUtils = "1.0"
 Colors = "0.12"
 ProgressMeter = "1"
 Requires = "1.0"
+RecipesBase = "1.1.1"
 julia = "1"
 
 [extras]
diff --git a/src/plots/plotsrecipes/PSAPlots.jl b/src/plots/plotsrecipes/PSAPlots.jl
new file mode 100644
index 0000000..3ff2ab2
--- /dev/null
+++ b/src/plots/plotsrecipes/PSAPlots.jl
@@ -0,0 +1,11 @@
+module PSAPlots
+    using RecipesBase
+    using Pseudospectra
+    using LinearAlgebra
+
+    import Pseudospectra: vec2ax, _basic_psa_opts, psa_compute
+
+    include("spectralportrait.jl")
+
+    export spectralportrait
+end
\ No newline at end of file
diff --git a/src/plots/plotsrecipes/spectralportrait.jl b/src/plots/plotsrecipes/spectralportrait.jl
new file mode 100644
index 0000000..cf978de
--- /dev/null
+++ b/src/plots/plotsrecipes/spectralportrait.jl
@@ -0,0 +1,43 @@
+@userplot SpectralPortrait
+
+RecipesBase.@recipe function f(p::SpectralPortrait; npts::Integer = 100)
+
+    if length(p.args) != 1 || !(typeof(p.args[1]) <: AbstractMatrix)
+        error("spectralportrait must be given an AbstractMatrix. Got $(typeof(p.args))")
+    end
+    A₀ = p.args[1]
+
+    local ps_data
+    try
+        ps_data = new_matrix(A₀)
+    catch JE
+        @warn "spectralportrait only works for simple cases."
+        rethrow(JE)
+    end
+
+    A = ps_data.matrix
+    ps_dict = ps_data.ps_dict
+    B = get(ps_dict, :matrix2, I)
+    eigA = ps_dict[:ews]
+    zoom = ps_data.zoom_list[ps_data.zoom_pos]
+    if isempty(zoom.ax)
+        zoom.ax = vec2ax(eigA)
+    end
+    psa_opts = _basic_psa_opts(zoom, ps_dict)
+
+    Z, xs, ys, t_levels, err, Tproj, eigAproj, algo = psa_compute(A, npts, zoom.ax, eigA, psa_opts, B)
+
+    @series begin
+        seriestype := :contour
+        xs, ys, log10.(Z)
+    end
+
+    @series begin
+        seriestype := :scatter
+        markercolor := :black
+        markersize := 2
+        label := "eigvals"
+        real.(eigA), imag.(eigA)
+    end
+
+end
\ No newline at end of file

From 5a143fe70f048ab5f3fa33f0d44db7064bb6d479 Mon Sep 17 00:00:00 2001
From: csimal <csimal@student.unamur.be>
Date: Sat, 14 Aug 2021 13:19:22 +0200
Subject: [PATCH 3/8] Added comments

---
 src/plots/plotsrecipes/spectralportrait.jl | 9 +++++++--
 1 file changed, 7 insertions(+), 2 deletions(-)

diff --git a/src/plots/plotsrecipes/spectralportrait.jl b/src/plots/plotsrecipes/spectralportrait.jl
index cf978de..83e5d70 100644
--- a/src/plots/plotsrecipes/spectralportrait.jl
+++ b/src/plots/plotsrecipes/spectralportrait.jl
@@ -1,10 +1,13 @@
 @userplot SpectralPortrait
-
+# Proof of concept Plot Recipe for spectral portraits.
+# See https://docs.juliaplots.org/latest/recipes/#case-studies for documentation on user recipes for Plots.jl
 RecipesBase.@recipe function f(p::SpectralPortrait; npts::Integer = 100)
-
+    # RecipesBase automatically generates the spectralportrait and spectralportrait! functions, so we have to do the typechecking here
     if length(p.args) != 1 || !(typeof(p.args[1]) <: AbstractMatrix)
         error("spectralportrait must be given an AbstractMatrix. Got $(typeof(p.args))")
     end
+
+    # Below is the computation of of the spectral portrait taken from the base implementation
     A₀ = p.args[1]
 
     local ps_data
@@ -27,6 +30,8 @@ RecipesBase.@recipe function f(p::SpectralPortrait; npts::Integer = 100)
 
     Z, xs, ys, t_levels, err, Tproj, eigAproj, algo = psa_compute(A, npts, zoom.ax, eigA, psa_opts, B)
 
+    # Recipe code below
+
     @series begin
         seriestype := :contour
         xs, ys, log10.(Z)

From 2958545df7c6c61c000ed902fb8674701bfd0c95 Mon Sep 17 00:00:00 2001
From: csimal <csimal@student.unamur.be>
Date: Sat, 14 Aug 2021 13:21:47 +0200
Subject: [PATCH 4/8] Remove unnecessary file

---
 src/recipes.jl | 1 -
 1 file changed, 1 deletion(-)
 delete mode 100644 src/recipes.jl

diff --git a/src/recipes.jl b/src/recipes.jl
deleted file mode 100644
index ad96c72..0000000
--- a/src/recipes.jl
+++ /dev/null
@@ -1 +0,0 @@
-# foo
\ No newline at end of file

From cf52f56fa34d4629c93e2b6624b1d86a92e93b2b Mon Sep 17 00:00:00 2001
From: csimal <csimal@student.unamur.be>
Date: Tue, 17 Aug 2021 22:52:28 +0200
Subject: [PATCH 5/8] [Feature] Added basic recipes for all plots

---
 src/plots/PSAPlots.jl                         | 22 +++++
 src/plots/colormap.jl                         | 88 +++++++++++++++++++
 src/plots/plotsrecipes/PSAPlots.jl            | 11 ---
 src/plots/recipes/mtxexpsplot.jl              | 36 ++++++++
 src/plots/recipes/mtxpowersplot.jl            | 49 +++++++++++
 src/plots/recipes/plotmode.jl                 | 41 +++++++++
 .../spectralportrait.jl                       |  1 +
 src/plots/recipes/surfplot.jl                 | 17 ++++
 src/plots/utils.jl                            | 61 +++++++++++++
 9 files changed, 315 insertions(+), 11 deletions(-)
 create mode 100644 src/plots/PSAPlots.jl
 create mode 100644 src/plots/colormap.jl
 delete mode 100644 src/plots/plotsrecipes/PSAPlots.jl
 create mode 100644 src/plots/recipes/mtxexpsplot.jl
 create mode 100644 src/plots/recipes/mtxpowersplot.jl
 create mode 100644 src/plots/recipes/plotmode.jl
 rename src/plots/{plotsrecipes => recipes}/spectralportrait.jl (97%)
 create mode 100644 src/plots/recipes/surfplot.jl
 create mode 100644 src/plots/utils.jl

diff --git a/src/plots/PSAPlots.jl b/src/plots/PSAPlots.jl
new file mode 100644
index 0000000..85f375a
--- /dev/null
+++ b/src/plots/PSAPlots.jl
@@ -0,0 +1,22 @@
+module PSAPlots
+    using RecipesBase
+    using Pseudospectra
+    using LinearAlgebra
+
+    import Pseudospectra: vec2ax, _basic_psa_opts, psa_compute
+
+    include("utils.jl")
+    include("colormap.jl")
+    
+    include("recipes/spectralportrait.jl")
+    include("recipes/plotmode.jl")
+    include("recipes/mtxpowersplot.jl")
+    include("recipes/mtxexpsplot.jl")
+    include("recipes/surfplot.jl")
+
+    export spectralportrait, spectralportrait!
+    export plotmode, plotmode!
+    export mtxpowersplot, mtxpowersplot!
+    export mtxexpsplot, mtxexpsplot!
+    export surfplot, surfplot!
+end
\ No newline at end of file
diff --git a/src/plots/colormap.jl b/src/plots/colormap.jl
new file mode 100644
index 0000000..647e9be
--- /dev/null
+++ b/src/plots/colormap.jl
@@ -0,0 +1,88 @@
+
+"""
+    The colormap from EigTool. This provides clearly
+    distinguishable colors for a line-contour plot.
+"""
+function etcgrad()
+    cm = [
+        0   1.00000000000000   0.10000000000000
+        0   0.92500000000000   0.17500000000000
+        0   0.85000000000000   0.25000000000000
+        0   0.77500000000000   0.32500000000000
+        0   0.70000000000000   0.40000000000000
+        0   0.62500000000000   0.47500000000000
+        0   0.55000000000000   0.55000000000000
+        0   0.47500000000000   0.62500000000000
+        0   0.40000000000000   0.70000000000000
+        0   0.32500000000000   0.77500000000000
+        0   0.25000000000000   0.85000000000000
+        0   0.17500000000000   0.92500000000000
+        0   0.10000000000000   1.00000000000000
+        0.06153846153846   0.09230769230769   1.00000000000000
+        0.12307692307692   0.08461538461538   1.00000000000000
+        0.18461538461538   0.07692307692308   1.00000000000000
+        0.24615384615385   0.06923076923077   1.00000000000000
+        0.30769230769231   0.06153846153846   1.00000000000000
+        0.36923076923077   0.05384615384615   1.00000000000000
+        0.43076923076923   0.04615384615385   1.00000000000000
+        0.49230769230769   0.03846153846154   1.00000000000000
+        0.55384615384615   0.03076923076923   1.00000000000000
+        0.61538461538462   0.02307692307692   1.00000000000000
+        0.67692307692308   0.01538461538462   1.00000000000000
+        0.73846153846154   0.00769230769231   1.00000000000000
+        0.80000000000000                  0   1.00000000000000
+        0.80769230769231                  0   0.92307692307692
+        0.81538461538462                  0   0.84615384615385
+        0.82307692307692                  0   0.76923076923077
+        0.83076923076923                  0   0.69230769230769
+        0.83846153846154                  0   0.61538461538462
+        0.84615384615385                  0   0.53846153846154
+        0.85384615384615                  0   0.46153846153846
+        0.86153846153846                  0   0.38461538461538
+        0.86923076923077                  0   0.30769230769231
+        0.87692307692308                  0   0.23076923076923
+        0.88461538461538                  0   0.15384615384615
+        0.89230769230769                  0   0.07692307692308
+        0.90000000000000                  0                  0
+        0.86923076923077                  0                  0
+        0.83846153846154                  0                  0
+        0.80769230769231                  0                  0
+        0.77692307692308                  0                  0
+        0.74615384615385                  0                  0
+        0.71538461538462                  0                  0
+        0.68461538461538                  0                  0
+        0.65384615384615                  0                  0
+        0.62307692307692                  0                  0
+        0.59230769230769                  0                  0
+        0.56153846153846                  0                  0
+        0.53076923076923                  0                  0
+        0.50000000000000                  0                  0
+        0.54166666666667   0.05000000000000                  0
+        0.58333333333333   0.10000000000000                  0
+        0.62500000000000   0.15000000000000                  0
+        0.66666666666667   0.20000000000000                  0
+        0.70833333333333   0.25000000000000                  0
+        0.75000000000000   0.30000000000000                  0
+        0.79166666666667   0.35000000000000                  0
+        0.83333333333333   0.40000000000000                  0
+        0.87500000000000   0.45000000000000                  0
+        0.91666666666667   0.50000000000000                  0
+        0.95833333333333   0.55000000000000                  0
+        1.00000000000000   0.60000000000000                  0
+    ]
+    nc = size(cm,1)
+    cc = Vector{RGBA{Float64}}()
+    vv = zeros(nc)
+    for i in 1:nc
+        col = RGBA(RGB(cm[i,:]...))
+        push!(cc,col)
+        vv[i]=(i-1)/(nc-1)
+    end
+    if :ContinuousColorGradient in names(Plots.PlotUtils, all=true)
+        return PlotUtils.ContinuousColorGradient(cc,vv)
+    else
+        return PlotUtils.ColorGradient(cc,vv)
+    end
+end
+
+const eigtool_cgrad = etcgrad()
\ No newline at end of file
diff --git a/src/plots/plotsrecipes/PSAPlots.jl b/src/plots/plotsrecipes/PSAPlots.jl
deleted file mode 100644
index 3ff2ab2..0000000
--- a/src/plots/plotsrecipes/PSAPlots.jl
+++ /dev/null
@@ -1,11 +0,0 @@
-module PSAPlots
-    using RecipesBase
-    using Pseudospectra
-    using LinearAlgebra
-
-    import Pseudospectra: vec2ax, _basic_psa_opts, psa_compute
-
-    include("spectralportrait.jl")
-
-    export spectralportrait
-end
\ No newline at end of file
diff --git a/src/plots/recipes/mtxexpsplot.jl b/src/plots/recipes/mtxexpsplot.jl
new file mode 100644
index 0000000..6d9cd14
--- /dev/null
+++ b/src/plots/recipes/mtxexpsplot.jl
@@ -0,0 +1,36 @@
+@userplot MtxExpsPlot
+
+
+RecipesBase.@recipe f(p::MtxExpsPlot, dt=0.1, nmax=50, lbmethod=:none)
+    if length(p.args) != 1 || !(typeof(p.args[1]) <: PSAStruct)
+        error("mtxexpsplot must be given a PSAStruct. Got $(typeof(p.args))")
+    end
+    ps_data::PSAStruct = p.args[1] # NB hint the type to the compiler
+    A = ps_data.input_matrix
+
+    ts, transient = get_mtxexpnorm(A, dt, nmax)
+    alp = maximum(real.(ps_data.ps_dict[:ews]))
+
+    layout := @layout [a; b]
+    link := :x
+    seriestype := :path
+    label := ""
+
+    @series begin
+        subplot := 1
+        ts, transient
+    end
+    @series begin
+        subplot := 1
+        ts, exp.(alp*ts)
+    end
+    @series begin
+        subplot := 2
+        yscale := :log10
+        ts, transient
+    end
+    @series begin
+        subplot := 2
+        ts, exp.(alp*ts)
+    end
+end
\ No newline at end of file
diff --git a/src/plots/recipes/mtxpowersplot.jl b/src/plots/recipes/mtxpowersplot.jl
new file mode 100644
index 0000000..e313905
--- /dev/null
+++ b/src/plots/recipes/mtxpowersplot.jl
@@ -0,0 +1,49 @@
+@userplot MtxPowersPlot
+
+
+RecipesBase.@recipe function f(p::MtxPowersPlot, nmax = 50,
+    lbmethod = :none,
+    lbdk = 0.25
+    )
+    if length(p.args) != 1 || !(typeof(p.args[1]) <: PSAStruct)
+        error("mtxpowersplot must be given a PSAStruct. Got $(typeof(p.args))")
+    end
+    ps_data::PSAStruct = p.args[1]
+    A = ps_data.input_matrix
+
+    powers, transient = get_mtxpowersnorm(A, nmax)
+
+    alp = maximum(abs.(ps_data.ps_dict[:ews]))
+
+    layout := @layout [a; b]
+    link := :x # link the x-axis of both subplots
+
+    # normal plot
+    @series begin
+        seriestype := :path
+        subplot := 1
+        label := ""
+        powers, transient
+    end
+    @series begin
+        seriestype := :path
+        subplot := 1
+        label := ""
+        powers, alp.^powers
+    end
+    # log plot
+    @series begin
+        seriestype := :path
+        subplot := 2
+        label := ""
+        yscale := :log10
+        powers, transient
+    end
+    @series begin
+        seriestype := :path
+        subplot := 2
+        label := ""
+        yscale := :log10
+        powers, alp.^powers
+    end
+end
\ No newline at end of file
diff --git a/src/plots/recipes/plotmode.jl b/src/plots/recipes/plotmode.jl
new file mode 100644
index 0000000..27419d0
--- /dev/null
+++ b/src/plots/recipes/plotmode.jl
@@ -0,0 +1,41 @@
+@userplot PlotMode
+
+struct PlotModeTheme
+    linecolor::Symbol
+    prefix::String
+    
+end
+
+
+
+RecipesBase.@recipe f(p::PlotMode; pseudo::Bool, approx::Bool, verbosity)
+
+    z = p.args[1]
+    A = p.args[2]
+    U = p.args[3]
+
+    SS, q = get_mode(A, z, pseudo, verbosity)
+
+    q = U * q
+    
+    #showlog = sum(abs.(diff(q))) * (1/length(q)) >= 10*eps()
+
+    @series begin
+        seriestype := :path
+        label := "realpt"
+        title := ""
+        real.(q)
+    end
+    @series begin
+        seriestype := :path
+        linecolor := :black
+        label := "abs"
+        abs.(q)
+    end
+    @series begin
+        seriestype := :path
+        linecolor := :black
+        label := ""
+        -abs.(q)
+    end
+end
\ No newline at end of file
diff --git a/src/plots/plotsrecipes/spectralportrait.jl b/src/plots/recipes/spectralportrait.jl
similarity index 97%
rename from src/plots/plotsrecipes/spectralportrait.jl
rename to src/plots/recipes/spectralportrait.jl
index 83e5d70..31f9fb6 100644
--- a/src/plots/plotsrecipes/spectralportrait.jl
+++ b/src/plots/recipes/spectralportrait.jl
@@ -34,6 +34,7 @@ RecipesBase.@recipe function f(p::SpectralPortrait; npts::Integer = 100)
 
     @series begin
         seriestype := :contour
+        linecolor --> eigtool_cgrad
         xs, ys, log10.(Z)
     end
 
diff --git a/src/plots/recipes/surfplot.jl b/src/plots/recipes/surfplot.jl
new file mode 100644
index 0000000..aacec80
--- /dev/null
+++ b/src/plots/recipes/surfplot.jl
@@ -0,0 +1,17 @@
+@userplot SurfPlot
+
+RecipesBase.@recipe f(p::SurfPlot)
+    if length(p.args) != 1 || (typeof(p.args[1]) <: PSAStruct)
+        error("surfplot must be given a PSAStruct. Got $(typeof(p.args))")
+    end
+
+    zoom = ps_data.zoom_list[ps_data.zoom_pos]
+    if !zoom.computed
+        return
+    end
+    
+    @series begin
+        seriestype := :surface
+        zoom.x, zoom.y, -log10.(zoom.Z)
+    end
+end
\ No newline at end of file
diff --git a/src/plots/utils.jl b/src/plots/utils.jl
new file mode 100644
index 0000000..d38dd86
--- /dev/null
+++ b/src/plots/utils.jl
@@ -0,0 +1,61 @@
+
+"""
+    get_mode(A, z, pseudo, verbosity)
+
+Compute the (pseudo)mode of `A` associated with `z`.
+"""
+function get_mode(A, z, pseudo, verbosity)
+    if pseudo
+        n = size(A, 2)
+        q = randn(n) + im*randn(n)
+        normalize!(q)
+        niter = 20
+
+        SS, q = psmode_inv_lanczos(A, q, z, 1e-15, niter)   
+    else
+        SS, q = oneeigcond(A,z, verbosity)
+    end
+    return SS, q
+end
+
+"""
+    get_mtxpowersnorm(A, nmax)
+
+Compute ``\|A^k\|`` for k up to `nmax`.
+"""
+function get_mtxpowersnorm(A, nmax)
+    transient = Vector{Float64}(undef, nmax+1)
+    transient[1] = 1.0
+    (min_tp, max_tp) = (-Inf, Inf)
+    mat = copy(A)
+    i = 2
+    transient[i] = norm(mat)
+    stop = false
+    while !stop && i < nmax + 1
+        mat = A * mat
+        i += 1
+        transient[i] = norm(mat)
+        min_tp = min(min_tp, transient[i])
+        max_tp = max(max_tp, transient[i])
+        if max_tp > 1e130 || min_tp < 1z-130 && min_tp != 0.0
+            @warn("stopping: bounds going out of range")
+            stop = true
+        end
+        if min_tp == 0.0
+            @warn("stopping: exactly zero matrix")
+            stop = true
+        end
+    end
+    return 0:(i-1), transient[1:i]
+end
+
+"""
+    get_mtxexpnorm(A, dt, nmax)
+
+Compute ``\|e^{A*k*dt}\|`` for k up to `nmax`.
+"""
+function get_mtxexpnorm(A, dt, nmax)
+    eAdt = exp(dt*A)
+    ns, transient = get_mtxpowersnorm(eAdt, nmax)
+    return dt*ns, transient
+end
\ No newline at end of file

From 2a80b996d6af9f878bf4280756ac284a2acab287 Mon Sep 17 00:00:00 2001
From: csimal <csimal@student.unamur.be>
Date: Wed, 18 Aug 2021 16:34:06 +0200
Subject: [PATCH 6/8] [Feature] Added title in plotmode

---
 src/plots/recipes/plotmode.jl | 18 ++++++++++--------
 1 file changed, 10 insertions(+), 8 deletions(-)

diff --git a/src/plots/recipes/plotmode.jl b/src/plots/recipes/plotmode.jl
index 27419d0..46a53b2 100644
--- a/src/plots/recipes/plotmode.jl
+++ b/src/plots/recipes/plotmode.jl
@@ -1,13 +1,14 @@
 @userplot PlotMode
 
-struct PlotModeTheme
-    linecolor::Symbol
-    prefix::String
-    
+function get_mode_title(pseudo, approx, z, s)
+    prefix = pseudo ? "Pseudomode: " : "Eigenmode: "
+    rel = approx ? "\\approx" : "="
+    formatted = @sprintf("%15.2e", s)
+    infix = pseudo ? ("\$\\epsilon = " * formatted * "\$") : ("\$\\kappa(\\lambda) " * rel * formatted * "\$")
+    λ_str = @sprintf("%12g%+12gi", real(z), imag(z))
+    return prefix * infix * "\n  \$\\lambda = " * λ_str * "\$"
 end
 
-
-
 RecipesBase.@recipe f(p::PlotMode; pseudo::Bool, approx::Bool, verbosity)
 
     z = p.args[1]
@@ -19,11 +20,12 @@ RecipesBase.@recipe f(p::PlotMode; pseudo::Bool, approx::Bool, verbosity)
     q = U * q
     
     #showlog = sum(abs.(diff(q))) * (1/length(q)) >= 10*eps()
-
+    title_str = get_mode_title(pseudo, approx, z, SS[end])
     @series begin
         seriestype := :path
+        linecolor := (pseudo ? :magenta : :cyan)
         label := "realpt"
-        title := ""
+        title := title_str
         real.(q)
     end
     @series begin

From 35ae8e17ed120a519a68bbe1efec30351cc0cb89 Mon Sep 17 00:00:00 2001
From: csimal <csimal@student.unamur.be>
Date: Sun, 22 Aug 2021 18:18:38 +0200
Subject: [PATCH 7/8] [Fix] Compiles, and mostly works

---
 src/Pseudospectra.jl                  |  2 ++
 src/plots/PSAPlots.jl                 |  8 ++++++--
 src/plots/defaults.jl                 |  8 ++++++++
 src/plots/recipes/mtxexpsplot.jl      | 14 +++++++++-----
 src/plots/recipes/mtxpowersplot.jl    |  9 ++++++---
 src/plots/recipes/plotmode.jl         | 27 +++++++++++++++++----------
 src/plots/recipes/spectralportrait.jl |  4 ++--
 src/plots/recipes/surfplot.jl         |  6 +++---
 src/plots/utils.jl                    |  4 ++--
 9 files changed, 55 insertions(+), 27 deletions(-)
 create mode 100644 src/plots/defaults.jl

diff --git a/src/Pseudospectra.jl b/src/Pseudospectra.jl
index 9c7928b..6dedf18 100644
--- a/src/Pseudospectra.jl
+++ b/src/Pseudospectra.jl
@@ -116,6 +116,8 @@ include("transients.jl")
 include("plotter.jl")
 include("zooming.jl")
 
+include("plots/PSAPlots.jl")
+
 """
     new_matrix(A::AbstractMatrix, opts::Dict{Symbol,Any}=()) -> ps_data
 
diff --git a/src/plots/PSAPlots.jl b/src/plots/PSAPlots.jl
index 85f375a..96bb7b8 100644
--- a/src/plots/PSAPlots.jl
+++ b/src/plots/PSAPlots.jl
@@ -2,11 +2,15 @@ module PSAPlots
     using RecipesBase
     using Pseudospectra
     using LinearAlgebra
+    using PlotUtils
+    using Printf
 
-    import Pseudospectra: vec2ax, _basic_psa_opts, psa_compute
+    import Pseudospectra: vec2ax, psa_compute
+    import Pseudospectra: psmode_inv_lanczos, oneeigcond
 
     include("utils.jl")
-    include("colormap.jl")
+    include("defaults.jl")
+    #include("colormap.jl")
     
     include("recipes/spectralportrait.jl")
     include("recipes/plotmode.jl")
diff --git a/src/plots/defaults.jl b/src/plots/defaults.jl
new file mode 100644
index 0000000..9afd5d5
--- /dev/null
+++ b/src/plots/defaults.jl
@@ -0,0 +1,8 @@
+
+_basic_psa_opts(zoom,ps_dict) = Dict{Symbol,Any}(
+    #:levels=>expandlevels(zoom.levels),
+    :recompute_levels=>zoom.autolev,
+    :proj_lev=>zoom.proj_lev,
+    :scale_equal=>zoom.scale_equal,
+    :real_matrix=>ps_dict[:Aisreal],
+    :verbosity=>0)
\ No newline at end of file
diff --git a/src/plots/recipes/mtxexpsplot.jl b/src/plots/recipes/mtxexpsplot.jl
index 6d9cd14..b5094e6 100644
--- a/src/plots/recipes/mtxexpsplot.jl
+++ b/src/plots/recipes/mtxexpsplot.jl
@@ -1,7 +1,7 @@
 @userplot MtxExpsPlot
 
 
-RecipesBase.@recipe f(p::MtxExpsPlot, dt=0.1, nmax=50, lbmethod=:none)
+RecipesBase.@recipe function f(p::MtxExpsPlot; dt=0.1, nmax=50, lbmethod=:none)
     if length(p.args) != 1 || !(typeof(p.args[1]) <: PSAStruct)
         error("mtxexpsplot must be given a PSAStruct. Got $(typeof(p.args))")
     end
@@ -11,19 +11,22 @@ RecipesBase.@recipe f(p::MtxExpsPlot, dt=0.1, nmax=50, lbmethod=:none)
     ts, transient = get_mtxexpnorm(A, dt, nmax)
     alp = maximum(real.(ps_data.ps_dict[:ews]))
 
-    layout := @layout [a; b]
-    link := :x
+    #layout := @layout [a; b]
+    #link := :x
     seriestype := :path
     label := ""
 
     @series begin
-        subplot := 1
+        #subplot := 1
+        label := "Matrix norm"
         ts, transient
     end
     @series begin
-        subplot := 1
+        #subplot := 1
+        label := "Spectral Radius"
         ts, exp.(alp*ts)
     end
+    #=
     @series begin
         subplot := 2
         yscale := :log10
@@ -33,4 +36,5 @@ RecipesBase.@recipe f(p::MtxExpsPlot, dt=0.1, nmax=50, lbmethod=:none)
         subplot := 2
         ts, exp.(alp*ts)
     end
+    =#
 end
\ No newline at end of file
diff --git a/src/plots/recipes/mtxpowersplot.jl b/src/plots/recipes/mtxpowersplot.jl
index e313905..1f90451 100644
--- a/src/plots/recipes/mtxpowersplot.jl
+++ b/src/plots/recipes/mtxpowersplot.jl
@@ -1,7 +1,7 @@
 @userplot MtxPowersPlot
 
 
-RecipesBase.@recipe function f(p::MtxPowersPlot, nmax = 50,
+RecipesBase.@recipe function f(p::MtxPowersPlot, nmax::Integer = 50;
     lbmethod = :none,
     lbdk = 0.25
     )
@@ -15,8 +15,9 @@ RecipesBase.@recipe function f(p::MtxPowersPlot, nmax = 50,
 
     alp = maximum(abs.(ps_data.ps_dict[:ews]))
 
-    layout := @layout [a; b]
-    link := :x # link the x-axis of both subplots
+    # NB. We need Plots to use @layout
+    #layout := @layout [a; b]
+    #link := :x # link the x-axis of both subplots
 
     # normal plot
     @series begin
@@ -32,6 +33,7 @@ RecipesBase.@recipe function f(p::MtxPowersPlot, nmax = 50,
         powers, alp.^powers
     end
     # log plot
+    #=
     @series begin
         seriestype := :path
         subplot := 2
@@ -46,4 +48,5 @@ RecipesBase.@recipe function f(p::MtxPowersPlot, nmax = 50,
         yscale := :log10
         powers, alp.^powers
     end
+    =#
 end
\ No newline at end of file
diff --git a/src/plots/recipes/plotmode.jl b/src/plots/recipes/plotmode.jl
index 46a53b2..24d56c9 100644
--- a/src/plots/recipes/plotmode.jl
+++ b/src/plots/recipes/plotmode.jl
@@ -9,8 +9,10 @@ function get_mode_title(pseudo, approx, z, s)
     return prefix * infix * "\n  \$\\lambda = " * λ_str * "\$"
 end
 
-RecipesBase.@recipe f(p::PlotMode; pseudo::Bool, approx::Bool, verbosity)
-
+RecipesBase.@recipe function f(p::PlotMode; pseudo=false, approx=false, verbosity=0)
+    if length(p.args) != 3
+        error("plotmode must be given 3 arguments")
+    end
     z = p.args[1]
     A = p.args[2]
     U = p.args[3]
@@ -20,23 +22,28 @@ RecipesBase.@recipe f(p::PlotMode; pseudo::Bool, approx::Bool, verbosity)
     q = U * q
     
     #showlog = sum(abs.(diff(q))) * (1/length(q)) >= 10*eps()
-    title_str = get_mode_title(pseudo, approx, z, SS[end])
+
+    # NB. The title is broken
+    #title_str = get_mode_title(pseudo, approx, z, SS[end])
+    #title := title_str
+    
     @series begin
         seriestype := :path
         linecolor := (pseudo ? :magenta : :cyan)
-        label := "realpt"
-        title := title_str
+        label := "Real part"
         real.(q)
     end
+
+    seriestype := :path
+    linecolor := :black
+    linestyle := :dash
+
     @series begin
-        seriestype := :path
-        linecolor := :black
-        label := "abs"
+        label := "Modulus"
         abs.(q)
     end
+
     @series begin
-        seriestype := :path
-        linecolor := :black
         label := ""
         -abs.(q)
     end
diff --git a/src/plots/recipes/spectralportrait.jl b/src/plots/recipes/spectralportrait.jl
index 31f9fb6..c453b0e 100644
--- a/src/plots/recipes/spectralportrait.jl
+++ b/src/plots/recipes/spectralportrait.jl
@@ -34,7 +34,7 @@ RecipesBase.@recipe function f(p::SpectralPortrait; npts::Integer = 100)
 
     @series begin
         seriestype := :contour
-        linecolor --> eigtool_cgrad
+        #linecolor --> eigtool_cgrad
         xs, ys, log10.(Z)
     end
 
@@ -42,7 +42,7 @@ RecipesBase.@recipe function f(p::SpectralPortrait; npts::Integer = 100)
         seriestype := :scatter
         markercolor := :black
         markersize := 2
-        label := "eigvals"
+        label --> "eigvals"
         real.(eigA), imag.(eigA)
     end
 
diff --git a/src/plots/recipes/surfplot.jl b/src/plots/recipes/surfplot.jl
index aacec80..9d70261 100644
--- a/src/plots/recipes/surfplot.jl
+++ b/src/plots/recipes/surfplot.jl
@@ -1,10 +1,10 @@
 @userplot SurfPlot
 
-RecipesBase.@recipe f(p::SurfPlot)
-    if length(p.args) != 1 || (typeof(p.args[1]) <: PSAStruct)
+RecipesBase.@recipe function f(p::SurfPlot)
+    if length(p.args) != 1 || !(typeof(p.args[1]) <: PSAStruct)
         error("surfplot must be given a PSAStruct. Got $(typeof(p.args))")
     end
-
+    ps_data = p.args[1]
     zoom = ps_data.zoom_list[ps_data.zoom_pos]
     if !zoom.computed
         return
diff --git a/src/plots/utils.jl b/src/plots/utils.jl
index d38dd86..9bf3be4 100644
--- a/src/plots/utils.jl
+++ b/src/plots/utils.jl
@@ -21,7 +21,7 @@ end
 """
     get_mtxpowersnorm(A, nmax)
 
-Compute ``\|A^k\|`` for k up to `nmax`.
+Compute ``|A^k|`` for k up to `nmax`.
 """
 function get_mtxpowersnorm(A, nmax)
     transient = Vector{Float64}(undef, nmax+1)
@@ -52,7 +52,7 @@ end
 """
     get_mtxexpnorm(A, dt, nmax)
 
-Compute ``\|e^{A*k*dt}\|`` for k up to `nmax`.
+Compute ``|e^{A*k*dt}|`` for k up to `nmax`.
 """
 function get_mtxexpnorm(A, dt, nmax)
     eAdt = exp(dt*A)

From 1cbbdd17fb84c1bdf604d5d122778df2664ca8ca Mon Sep 17 00:00:00 2001
From: csimal <csimal@student.unamur.be>
Date: Mon, 13 Sep 2021 08:54:14 +0200
Subject: [PATCH 8/8] [Feature] Added Documentation

---
 src/plots/recipes/mtxpowersplot.jl    | 29 ++++++++++++++++++++---
 src/plots/recipes/plotmode.jl         | 33 ++++++++++++++++++++++-----
 src/plots/recipes/spectralportrait.jl | 19 ++++++++++++---
 src/plots/utils.jl                    |  4 ++--
 4 files changed, 71 insertions(+), 14 deletions(-)

diff --git a/src/plots/recipes/mtxpowersplot.jl b/src/plots/recipes/mtxpowersplot.jl
index 1f90451..64357ae 100644
--- a/src/plots/recipes/mtxpowersplot.jl
+++ b/src/plots/recipes/mtxpowersplot.jl
@@ -1,7 +1,30 @@
 @userplot MtxPowersPlot
 
+"""
+    mtxpowersplot(A, nmax)
 
-RecipesBase.@recipe function f(p::MtxPowersPlot, nmax::Integer = 50;
+Plot the matrix norm of `A^k` for `k` up to `nmax`
+
+# Examples
+```jldoctest
+julia> opts = Dict{Symbol,Any}(:npts => 40,:ax => [-110,110,-110,110])
+julia> A = diagm(1 => 100.0*ones(4)) + diagm(-2 => 1e-4*ones(3))
+5×5 Matrix{Float64}:
+ 0.0     100.0       0.0       0.0    0.0
+ 0.0       0.0     100.0       0.0    0.0
+ 0.0001    0.0       0.0     100.0    0.0
+ 0.0       0.0001    0.0       0.0  100.0
+ 0.0       0.0       0.0001    0.0    0.0
+ julia> ps_data = new_matrix(A,opts)
+ [...]
+ julia> driver!(ps_data, opts)
+ julia> mtxpowersplot(ps_data, nmax=20)
+```
+"""
+function mtxpowersplot end
+
+RecipesBase.@recipe function f(p::MtxPowersPlot;
+    nmax::Integer = 50,
     lbmethod = :none,
     lbdk = 0.25
     )
@@ -23,13 +46,13 @@ RecipesBase.@recipe function f(p::MtxPowersPlot, nmax::Integer = 50;
     @series begin
         seriestype := :path
         subplot := 1
-        label := ""
+        label := "Matrix norm"
         powers, transient
     end
     @series begin
         seriestype := :path
         subplot := 1
-        label := ""
+        label := "Spectral Radius"
         powers, alp.^powers
     end
     # log plot
diff --git a/src/plots/recipes/plotmode.jl b/src/plots/recipes/plotmode.jl
index 24d56c9..d437e80 100644
--- a/src/plots/recipes/plotmode.jl
+++ b/src/plots/recipes/plotmode.jl
@@ -9,13 +9,34 @@ function get_mode_title(pseudo, approx, z, s)
     return prefix * infix * "\n  \$\\lambda = " * λ_str * "\$"
 end
 
-RecipesBase.@recipe function f(p::PlotMode; pseudo=false, approx=false, verbosity=0)
-    if length(p.args) != 3
-        error("plotmode must be given 3 arguments")
+"""
+    plotmode(A, z, U=I; ...)
+
+Plot the (pseudo)eigenmode of `A` associated with (pseudo)eigenvalue `z`.
+
+# Arguments
+* `A` an arbitrary square matrix
+* `z` the real or complex (pseudo)eigenvalue whose eigenmode is to be plotted.
+* `U` optional Schur factor. The eigenmode is premultiplied by `U` before plotting.
+# Keyword arguments
+* `pseudo::Bool = false` Whether to plot the eigenmode or pseudomode.
+* `approx::Bool = false`
+* `verbosity::Int = 0` Verbosity parameter passed when computing the eigenmode.
+
+# Examples
+```jldoctest
+julia> A = Pseudospectra.grcar(40)
+julia> plotmode(A, 0.5+2.0im)
+```
+"""
+function plotmode end
+
+RecipesBase.@recipe function f(p::PlotMode, U = I; pseudo=false, approx=false, verbosity=0)
+    if length(p.args) != 2
+        error("plotmode must be given 2 arguments")
     end
-    z = p.args[1]
-    A = p.args[2]
-    U = p.args[3]
+    A = p.args[1]
+    z = p.args[2]
 
     SS, q = get_mode(A, z, pseudo, verbosity)
 
diff --git a/src/plots/recipes/spectralportrait.jl b/src/plots/recipes/spectralportrait.jl
index c453b0e..a690626 100644
--- a/src/plots/recipes/spectralportrait.jl
+++ b/src/plots/recipes/spectralportrait.jl
@@ -1,8 +1,21 @@
 @userplot SpectralPortrait
-# Proof of concept Plot Recipe for spectral portraits.
-# See https://docs.juliaplots.org/latest/recipes/#case-studies for documentation on user recipes for Plots.jl
+
+"""
+    spectralportrait(A; npts=100)
+
+Plot the spectral portrait of the square matrix `A`.
+
+# Keyword arguments
+* `npts::Integer = 100` The size of the grid used for computing the levels of the pseudospectrum
+# Examples
+```jldoctest
+julia> A = Pseudospectra.grcar(40)
+julia> spectralportrait(A)
+```
+"""
+function spectralportrait end
+
 RecipesBase.@recipe function f(p::SpectralPortrait; npts::Integer = 100)
-    # RecipesBase automatically generates the spectralportrait and spectralportrait! functions, so we have to do the typechecking here
     if length(p.args) != 1 || !(typeof(p.args[1]) <: AbstractMatrix)
         error("spectralportrait must be given an AbstractMatrix. Got $(typeof(p.args))")
     end
diff --git a/src/plots/utils.jl b/src/plots/utils.jl
index 9bf3be4..c64f762 100644
--- a/src/plots/utils.jl
+++ b/src/plots/utils.jl
@@ -26,7 +26,7 @@ Compute ``|A^k|`` for k up to `nmax`.
 function get_mtxpowersnorm(A, nmax)
     transient = Vector{Float64}(undef, nmax+1)
     transient[1] = 1.0
-    (min_tp, max_tp) = (-Inf, Inf)
+    (min_tp, max_tp) = (Inf, -Inf)
     mat = copy(A)
     i = 2
     transient[i] = norm(mat)
@@ -37,7 +37,7 @@ function get_mtxpowersnorm(A, nmax)
         transient[i] = norm(mat)
         min_tp = min(min_tp, transient[i])
         max_tp = max(max_tp, transient[i])
-        if max_tp > 1e130 || min_tp < 1z-130 && min_tp != 0.0
+        if (max_tp > 1e130) || ((min_tp < 1e-130) && (min_tp != 0.0))
             @warn("stopping: bounds going out of range")
             stop = true
         end