text:
#
# plotting.jl --
#
# Implements high-level interface for plotting with the PGPlot library.
#

module Plotting

export
    Figure,
    colorbar,
    curve!,
    curve,
    figure,
    heatmap!,
    heatmap,
    hist!,
    hist,
    mouse,
    palette,
    scatter!,
    scatter

using Colors, ArrayTools, TwoDimensional
using ..PGPlot.Bindings
using ..PGPlot.Colormaps
using .Bindings: pgarray, DEFAULT_XFORM, get_own_xform, throw_bad_xform_size

# FIXME: type piracy
Base.Tuple(A::AffineTransform) = (A.x, A.xx, A.xy, A.y, A.yx, A.yy)

# Structure used to wrap an integer so that some basic function can be extended
# while avoiding type-piracy.
struct Figure
    num::Int   # Figure number.
end
number(fig::Figure) = fig.num

const ExtentOption = Union{Nothing,NTuple{4,<:Union{Nothing,Real}}}
const AxisStyleOption = Union{Nothing,Symbol,Integer}
const ColorOption = Union{Nothing,AbstractString,Symbol,Integer,Colorant,RGBVec}
const FontOption = Union{Nothing,AbstractString,Symbol,Integer}
const LabelOption = Union{Nothing,AbstractString,Symbol}
const XFormOption = Union{Nothing,AffineTransform,NTuple{6,Real},
                          AbstractVector{<:Real},AbstractMatrix{<:Real}}
const HeightOption = Union{Nothing,Real}
const FigureOption = Union{Nothing,Figure,Integer}

const Marker = Union{Integer,Char} # FIXME: Symbol?
const Markers = Union{Marker,AbstractVector{<:Marker}}

# Color index to use for color set by value/name.
const USER_COLOR_INDEX = PGInt(9)

# Constants for available fonts.
const NORMAL_FONT = PGInt(1)
const ROMAN_FONT  = PGInt(2)
const ITALIC_FONT = PGInt(3)
const SCRIPT_FONT = PGInt(4)

# Transparent text background color.
const TRANSPARENT_BACKGROUND = PGInt(-1)
const BACKGROUND_COLOR_INDEX = PGInt(0)
const FOREGROUND_COLOR_INDEX = PGInt(1)
const DEFAULT_FRAME_COLOR = FOREGROUND_COLOR_INDEX
const DEFAULT_TITLE_COLOR = FOREGROUND_COLOR_INDEX
const DEFAULT_DRAWING_COLOR = FOREGROUND_COLOR_INDEX
const DEFAULT_TEXT_BACKGROUND = TRANSPARENT_BACKGROUND
const DEFAULT_TEXT_COLOR = DEFAULT_DRAWING_COLOR # for annotations
const DEFAULT_TEXT_FONT = NORMAL_FONT
const DEFAULT_TEXT_HEIGHT = PGFloat(1.0)

"""

`Box{T}` represents a rectangular region whose side are aligned with the
coordinate axes.

"""
struct Box{T<:Real}
    x1::T
    x2::T
    y1::T
    y2::T
end

Box(A::NTuple{4,Real}) = Box(A[1], A[2], A[3], A[4])
Box(x1, x2, y1, y2) = Box(promote(x1, x2, y1, y2)...)
Box(B::Box) = B
Box{T}(B::Box{T}) where {T<:Real} = B
Box{T}(B::Box) where {T<:Real} = Box{T}(B.x1, B.x2, B.y1, B.y2)
Base.Tuple(B::Box) = (B.x1, B.x2, B.y1, B.y2)

const BAD_COLOR = colorant"orange"

abstract type AbstractScaling{T<:Integer} end

struct LinearScaling{T} <: AbstractScaling{T}
    lo::T   # first color index
    hi::T   # last color index
    bad::T  # color index for bad values
end

Base.first(scl::AbstractScaling) = scl.lo
Base.last(scl::AbstractScaling) = scl.hi
bad(scl::AbstractScaling) = scl.bad

LinearScaling(lo::Integer, hi::Integer) =
    LinearScaling(lo, hi, get_color(BAD_COLOR))

LinearScaling{T}(lo::Integer, hi::Integer) where {T<:Integer} =
    LinearScaling{T}(lo, hi, get_color(BAD_COLOR))

LinearScaling(; badcolor::ColorOption = BAD_COLOR) =
    LinearScaling{PGInt}(pgqcir()..., get_color(badcolor))

(scl::LinearScaling)(A::AbstractArray, vmin, vmax) =
    rescale(scl, A, vmin, vmax)

have_bad_values(A::AbstractArray{<:Integer}) = false

function have_bad_values(A::AbstractArray{<:Real})
    flag = true
    @inbounds @simd for i in eachindex(A)
        flag &= isfinite(A[i])
    end
    return !flag
end

check(A::AbstractArray) =
    (have_bad_values(A)
     && throw(ArgumentError("non-finite value(s) in argument")))

function check(x::AbstractVector, y::AbstractVector)
    length(x) == length(y) ||
        throw(DimensionMismatch("x and y must have the same length"))
    check(x)
    check(y)
end

"""

```julia
twice(x) -> (x, x)
```

yeilds a 2-tuple with the same value repeated twice.

"""
twice(x) = (x, x)

"""

```julia
confirm(flag)
```

specifies whether or not interactive devices should ask for user
confirmation before advancing to next plot.

"""
confirm(flag::Bool) = pgask(flag)

# PGPlot has a rather small maximum number (GRIMAX = 8) of open devices.
# We keep the list of open devices/figures in a vector.
const DEFAULT_DEVICE = "/XSERVE"
const DEVICES = Vector{PGInt}(undef, 0)
const BAD_DEVICE = eltype(DEVICES)(-1)
const FIGCNT = Ref{Int}(0)
const MAXFIGCNT = 100 # just used to detect wrong figure numbers

"""

```julia
device(fig)
```

yields the PGPlot device number associated with figure `fig`.  The result is a
positive integer if there is an open PGPlot device for that figure, -1
otherwise.

"""
device(fig::Figure) = device(number(fig))
device(fig::Integer) = (1 ≤ fig ≤ length(DEVICES) ? DEVICES[fig] : BAD_DEVICE)

"""

```julia
figure() -> fig
figure(nothing) -> fig
```

create a new figure using the default plotting device and selects it as the
current plotting device.  The figure is returned.

```julia
figure(dev::AbstractString) -> fig
```

opens the plotting device `dev`, selects it as the current plotting device and
associates a new figure with it.  The figure is returned.

```julia
figure(num::Integer) -> fig
```

if figure number `num` exists, makes it associated plotting device the current
one; otherwise opens the default plotting device, selects it as the current
plotting device and associates it with the figure number `num`.  The figure is
returned.

```julia
figure(fig::Figure) -> fig
```

checks that figure `fig` is associated with an open plotting device and makes
it the current one.  The figure is returned.  An exception is thrown if `fig`
is not associated with an open plotting device.

""" figure

# Create a new figure using the default device.
figure() = figure(DEFAULT_DEVICE)
figure(::Nothing) = figure()

# Create a new figure for the given device string.
figure(dev::AbstractString) = register_device(open_device(dev))

function figure(num::Integer)
    1 ≤ num ≤ MAXFIGCNT || throw_invalid_figure_number(num)
    if num ≤ length(DEVICES)
        dev = DEVICES[num]
        if dev ≥ 1
            pgslct(dev)
            return Figure(num)
        end
    end
    while length(DEVICES) < num
        push!(DEVICES, BAD_DEVICE)
    end
    DEVICES[num] = open_device(DEFAULT_DEVICE)
    return Figure(num)
end

# Use an existing figure.
function figure(fig::Figure)
    dev = device(fig)
    if dev > 0
        pgslct(dev)
        return fig
    end
    1 ≤ number(fig) ≤ length(DEVICES) && throw_device_closed(fig)
    throw_nonexisting_figure(fig)
end

Base.isopen(fig::Figure) = device(fig) > 0

function Base.close(fig::Figure)
    dev = device(fig)
    if dev > 0
        DEVICES[dev] = BAD_DEVICE
        pgclos(dev)
        return nothing
    end
    1 ≤ number(fig) ≤ length(DEVICES) && throw_device_closed(fig)
    throw_nonexisting_figure(fig)
end

@noinline throw_invalid_figure_number(num::Integer) =
    throw(ArgumentError(string("invalid figure number: ", num)))

@noinline throw_invalid_device(dev::Integer) =
    throw(ArgumentError(string("invalid PGPlot device number: ", dev)))

@noinline throw_open_device_failure(dev::AbstractString) =
    error(string("cannot open PGPlot device: ", dev))

throw_nonexisting_figure(fig::Figure) = throw_nonexisting_figure(number(fig))
@noinline throw_nonexisting_figure(num::Integer) =
    throw(ArgumentError(string("Figure #", num, " does not exist")))

throw_device_closed(fig::Figure) = throw_device_closed(number(fig))
@noinline throw_device_closed(num::Integer) =
    throw(ArgumentError(string("Plotting device of Figure #", num,
                               " has been closed")))

function current_figure()
    dev = pgqid()
    dev > 0 || error("there is no open plotting devices")
    return register_device(dev)
end

select_figure(fig::Figure) = figure(fig)
select_figure(num::Integer) = figure(num)
select_figure(::Nothing) = select_figure()
select_figure() = (if pgqid() < 1; figure(); end; nothing)

function open_device(name::AbstractString)
    dev = pgopen(name)
    dev > 0 || throw_open_device_failure(name)
    pgask(false) # do not wait for user input
    if pgqinf("HARDCOPY") == "YES"
        pgscr(0, 1,1,1) # set background color to white
        pgscr(1, 0,0,0) # set foreground color to black
    else
        pgscr(0, 0,0,0) # set background color to black
        pgscr(1, 1,1,1) # set foreground color to white
    end
    return dev
end

function register_device(dev::Integer)
    dev > 0 || throw_invalid_device(dev)
    j = 0
    @inbounds for i in 1:length(DEVICES)
        if DEVICES[i] == dev
            return Figure(i)
        elseif j == 0 && DEVICES[i] < 1
            j = i
        end
    end
    if j != 0
        DEVICES[j] = dev
        return Figure(j)
    end
    if length(DEVICES) < MAXFIGCNT
        push!(DEVICES, dev)
        return Figure(length(DEVICES))
    end
    error("too many registered figures")
end

function forget_device(dev::Integer)
    if dev > 0
        @inbounds for i in 1:length(DEVICES)
            if DEVICES[i] == dev
                DEVICES[i] = BAD_DEVICE
                break
            end
        end
    end
    nothing
end

"""

```julia
newpage()
```

starts a new page with the standard viewport.  It is possible to specify
the viewport bounds in normalized device coordinate units (NDC):

```julia
newpage(B::Box)
newpage(x1, x2, y1, y2)
```

"""
newpage() = (pgpage(); setviewport())
newpage(B::Box) = (pgpage(); setviewport(B))
newpage(x1::Real, x2::Real, y1::Real, y2::Real) =
    (pgpage(); setviewport(x1, x2, y1, y2))

"""

```julia
getviewport(units=0) -> x1, x2, y1, y2
```

yields the viewport limits in requested units, normalized device coordinate
units (NDC) by default.  To retrieve the viewport as a box:

```julia
getviewport(Box{T}, units=0) -> box
```

"""
getviewport(units::Integer=0) = pgqvp(units)
getviewport(::Type{T}, units::Integer=0) where {T<:Box} = T(getviewport(units))

"""

```julia
setviewport()
```

sets the standard viewport.  It is possible to specify the viewport
bounds in normalized device coordinate units (NDC):

```julia
setviewport(B::Box)
setviewport(x1, x2, y1, y2)
```

"""
setviewport() = pgvstd()
setviewport(B::Box) = setviewport(B.x1, B.x2, B.y1, B.y2)
setviewport(x1::Real, x2::Real, y1::Real, y2::Real) =
    pgsvp(x1, x2, y1, y2)

"""

```julia
getworld() -> x1, x2, y1, y2
```

yields the current world coordinate limits.  To retrieve these limits as a
box:

```julia
getworld(Box{T}) -> box
```

"""
getworld() = Box(pgqwin()...)
getworld(::Type{T}) where {T<:Box} = T(getworld())

"""

```julia
setworld(B::Box, just=false)
setworld(x1, x2, y1, y2, just=false)
```

both set the world coordinate limits.  If `just` is true, the viewport is
adjusted so that both axes have the same scaling.

"""
setworld(B::Box, just::Bool=false) = setworld(B.x1, B.x2, B.y1, B.y2, just)
setworld(x1::Real, x2::Real, y1::Real, y2::Real, just::Bool=false) =
    (just ? pgwnad(x1, x2, y1, y2) : pgswin(x1, x2, y1, y2))

"""

```julia
frame(..., adv=false)
```

draws the environment of subsequent plots.

If `adv` is true, the graphic page is erased and a new page is created
befroe drawing the frame.

"""
frame(B::Box, adv::Bool = false; kwds...) =
    frame(B.x1, B.x2, B.y1, B.y2, adv; kwds...)

function frame(x1::Real, x2::Real, y1::Real, y2::Real,
               adv::Bool = false;
               fig::FigureOption = nothing,
               wait::Bool = false,
               just::Bool = false,
               style::AxisStyleOption = nothing,
               xlabel::LabelOption = nothing,
               ylabel::LabelOption = nothing,
               title::LabelOption = nothing,
               textcolor::ColorOption = nothing,
               framecolor::ColorOption = nothing)

    xopts, yopts = get_axes_style(style)
    select_figure(fig)
    pgsci(get_color(framecolor, DEFAULT_FRAME_COLOR))
    adv && newpage()
    setworld(x1, x2, y1, y2, just)
    pgbox(xopts, 0.0, 0, yopts, 0.0, 0)
    #pgenv(x1, x2, y1, y2, just, get_axis_style(style))
    if xlabel !== nothing || ylabel !== nothing || title !== nothing
        # deal with colors, size, text style
        pgsci(get_color(textcolor, DEFAULT_TITLE_COLOR))
        pglab(get_label(xlabel), get_label(ylabel), get_label(title))
    end
end

function frame(x::AbstractVector, y::AbstractVector, adv::Bool = false;
               extent::ExtentOption = nothing,
               kwds...)
    check(x, y)
    x1, x2 = get_vrange(PGFloat, x, extent, (1,2))
    y1, y2 = get_vrange(PGFloat, y, extent, (3,4))
    frame(x1, x2, y1, y2, adv; kwds...)
end

get_axes_style(opt) = get_axes_style(get_axis_style(opt))

get_axes_style(opt::Integer) =
    (opt == -2 ? twice(" ") :
     opt == -1 ? twice("BC") :
     opt ==  0 ? twice("BCNST") :
     opt ==  1 ? twice("ABCNST") :
     opt ==  2 ? twice("ABCGNST") :
     opt == 10 ? ("BCNSTL", "BCNST") :
     opt == 20 ? ("BCNST",  "BCNSTL") :
     opt == 30 ? ("BCNSTL", "BCNSTL") :
     throw(ArgumentError("illegal AXIS argument")))

"""

```julia
curve([x,] y; kwds...)
```

draws a frame box (with [`frame`](@ref)) and a curve (as a polyline) joining
the points defined by abscissae `x` and ordinates `y` (both being vectors in
that case).  If `x` is missing, the vector `y` is plotted against its indices.

To plot functions, do one of:

```julia
curve(x, fy; kwds...)      # plot fy(x) against x
curve(fx, y; kwds...)      # plot y against fx(y)
curve(t, fx, fy; kwds...)  # plot fy(t) against fx(t) (parametric plot)
```

To add a curve to an existing plot, call [`curve!`](@ref) instead.

Possible keywords:

- `color` specifies the color of the lines.

Other keywods are transmitted to [`frame`](@ref) used to draw the frame box
(e.g., `fig`, `xlabel`, `ylabel`, `title`, `extent`, ...).

"""
function curve(x::AbstractVector, y::AbstractVector;
              color::ColorOption = nothing,
              kwds...)
    frame(x, y, true; kwds...)
    _curve(x, y, color)
end

"""

```julia
curve!(args...; kwds...)
```

draws a curve in the current frame preserving existing contents.  Arguments are
the same as for [`curve`](@ref).

Possible keywords: `color`, `fig`.

- `color` specifies the color of the lines.
- `fig` specifies the figure where to draw.

"""
function curve!(x::AbstractVector, y::AbstractVector;
               fig::FigureOption = nothing,
               color::ColorOption = nothing)
    select_figure(fig)
    check(x, y)
    _curve(x, y, color)
end

# Helper for `curve` and `curve!` methods.
function _curve(x::AbstractVector, y::AbstractVector,
                color::ColorOption)
    pgsci(get_color(color))
    pgline(x, y)
end

"""

```julia
hist([x,] y; kwds...)
```

draws a histogram of `y` against `x`, that is a stair-style curve.

"""
function hist(x::AbstractVector, y::AbstractVector;
              center::Bool = true,
              color::ColorOption = nothing,
              kwds...)
    frame(x, y, true; kwds...)
    _hist(x, y, center, color)
end

function hist!(x::AbstractVector, y::AbstractVector;
               fig::FigureOption = nothing,
               center::Bool = true,
               color::ColorOption = nothing)
    select_figure(fig)
    check(x, y)
    _hist(x, y, center, color)
end

# Helper for `hist` and `hist!` methods.
function _hist(x, y, center::Bool, color::ColorOption)
    pgsci(get_color(color))
    pgbin(x, y, center)
end

"""

```julia
scatter([x,] y, sym; kwds...)
```

draws points with symbol `sym` at coordinates `(x,y)`.

"""
function scatter(x::AbstractVector, y::AbstractVector,
                 sym::Markers;
                 color::ColorOption = nothing,
                 kwds...)
    frame(x, y, true; kwds...)
    _scatter(x, y, sym, color)
end

function scatter!(x::AbstractVector, y::AbstractVector,
                  sym::Markers;
                  fig::FigureOption = nothing,
                  color::ColorOption = nothing)
    select_figure(fig)
    check(x, y)
    _scatter(x, y, sym, color)
end

# Helper for `scatter` and `scatter!` methods.
function _scatter(x::AbstractVector, y::AbstractVector,
                  sym::Union{Integer,AbstractVector{<:Integer}},
                  color::ColorOption)
    pgsci(get_color(color))
    pgpnts(x, y, sym)
end

# Extend curve-like methods.
for func in (:curve, :curve!, :hist, :hist!)
    @eval begin

        $func(y::AbstractVector; kwds...) =
            $func(axes(y,1), y; kwds...)

        $func(x::AbstractVector, fy; kwds...) =
            $func(x, PGFloat[fy(x[i]) for i in eachindex(x)]; kwds...)

        $func(fx, y::AbstractVector; kwds...) =
            $func(PGFloat[fx(y[i]) for i in eachindex(y)], y; kwds...)

        $func(t::AbstractVector, fx, fy; kwds...) =
            $func(PGFloat[fx(t[i]) for i in eachindex(t)],
                  PGFloat[fy(t[i]) for i in eachindex(t)]; kwds...)
    end
end

for func in (:scatter, :scatter!)
    @eval begin

        $func(y::AbstractVector, sym::Markers; kwds...) =
            $func(axes(y,1), y, sym; kwds...)

        $func(x::AbstractVector, fy, sym::Markers; kwds...) =
            $func(x, PGFloat[fy(x[i]) for i in eachindex(x)], sym; kwds...)

        $func(fx, y::AbstractVector, sym::Markers; kwds...) =
            $func(PGFloat[fx(y[i]) for i in eachindex(y)], y, sym; kwds...)

        $func(t::AbstractVector, fx, fy, sym::Markers; kwds...) =
            $func(PGFloat[fx(t[i]) for i in eachindex(t)],
                  PGFloat[fy(t[i]) for i in eachindex(t)], sym; kwds...)
    end
end

"""

```julia
heatmap(A; kwds...)
```

draws a *heat-map*, that is a map whose cells have pseudo-colors computed
from the values in the 2-dimensional array `A`.  By default, the first and
second dimensions of `A` are assumed to correspond to the abscissa and
ordinate axes respectively and the world coordinates are assumed equal to
the indices of `A` (see below for other possibilities).

Keywords (other keywords are passed to [`frame`](@ref)):

- `vmin` and `vmax` specify the range of values to consider.  It is possible to
  have `vmin > vmax` to reverse the order of the colors.  If unspecified, the
  extreme values of `A` are considered.

- `badcolor` is the color to attribute to bad values (like NaN, *not-a-number*)
  in `A`.

- `cbar` specifies where to put a color-bar: `cbar = 'L'` `'R'`, `'T'` or `'B'`
  to draw a color-bar on the Left, Right, Top or Bottom side of the frame;
  `cbar = nothing` to draw no color-bar.  By default, `cbar = 'R'`.  Keyword
  `zlabel` can be used to specify a label for the color-bar.  For more options,
  use `cbar = nothing` and call [`colorbar`](@ref) directly after drawing the
  heat-map (the range of values are automatically saved by `heatmap` and need
  not be specified to `colorbar` in this case).

- `just` specifies whether both axis should have the same scale or not.  By
  default, `just = true`.

To add another heat-map to an existing plot, call [`heatmap!`](@ref) instead.

To change the default coordiante mapping, the world coordinates `x1`, `x2`,
`y1` and `y2` of the corners of the rectangular region covered by `A` may
be specified:

```julia
heatmap(A, x1, x2, y1, y2; kwds...)
```

The rectangular region covered by `A` has its edges aligned with the world
coordinate axes.  To draw a heat-map of `A` in a region of a different
shape, call:

```julia
heatmap(A, tr; kwds...)
```

with `tr` an affine coordinate transform to map index `(i,j)` to world
coordinates `(x,y)` as follows:

```
x = tr[1] + tr[2]*i + tr[3]*j
y = tr[4] + tr[5]*i + tr[6]*j
```

The coordinate transform `tr` is a 6-element vector or a 2×3 matrix.  Hence
identity is given by: `tr = [0,1,0, 0,0,1]`.

"""
function heatmap(A::AbstractMatrix{T}, args...;
                 fig::FigureOption = nothing,
                 just::Bool = true,
                 cbar::Union{Nothing,Char} = 'R',
                 zlabel::AbstractString = "",
                 cmap::Union{Nothing,AbstractString} = nothing,
                 vmin::Union{Nothing,Real} = nothing,
                 vmax::Union{Nothing,Real} = nothing,
                 badcolor::ColorOption = nothing,
                 kwds...) where {T<:Real}
    # Get size of region covered by the map.
    box = get_extent(Box{PGFloat}, A, args...)

    # Select figure, advance page, set default viewport and set world.
    select_figure(fig)
    newpage()
    setworld(box, just)

    # Select colormap.
    cmap === nothing || palette(cmap)

    # Draw map.
    _heatmap(A, args...; vmin=vmin, vmax=vmax, badcolor=badcolor)

    # Draw the frame.
    frame(box, false; just = just, kwds...)

    # Draw optional color bar.
    cbar === nothing || colorbar(; side = cbar, label = zlabel)
end

function _heatmap(A::AbstractMatrix{<:Real}; kwds...)
    _heatmap(pgarray(PGFloat, A), DEFAULT_XFORM; kwds...)
end

function _heatmap(A::AbstractMatrix{<:Real},
                  tr::AbstractArray{<:Real}; kwds...)
    _heatmap(pgarray(PGFloat, A), get_own_xform(tr); kwds...)
end

function _heatmap(A::DenseMatrix{PGFloat}, tr::DenseArray{PGFloat};
                  vmin::Union{Nothing,Real} = nothing,
                  vmax::Union{Nothing,Real} = nothing,
                  badcolor::ColorOption = nothing)
    bad, lo, hi = rangeofvalues(PGFloat, A, vmin, vmax)
    bad && error("cannot plot image with bad value, call with (x1,x2,y1,y2)")
    memorize_zrange(lo, hi)
    pgimag(A, lo, hi, tr)
end

function _heatmap(A::AbstractMatrix{<:Real},
                  x1::Real, x2::Real, y1::Real, y2::Real; kwds...)
    _heatmap(A, PGFloat(x1), PGFloat(x2), PGFloat(y1), PGFloat(y2); kwds...)
end

function _heatmap(A::AbstractMatrix{<:Integer},
                  x1::PGFloat, x2::PGFloat, y1::PGFloat, y2::PGFloat;
                  vmin::Union{Nothing,Real} = nothing,
                  vmax::Union{Nothing,Real} = nothing,
                  badcolor::ColorOption = BAD_COLOR)
    x1, x2 = adjust_interval(x1, x2)
    y1, y2 = adjust_interval(y1, y2)
    # FIXME: `pgpixl` is much slower than `pgimag` so avoid calling `pgpixl`.
    tr = simple_xform!(get_own_xform(), A, x1, x2, y1, y2)
    bad, lo, hi = rangeofvalues(PGFloat, A, vmin, vmax)
    memorize_zrange(lo, hi)
    @assert bad == false
    pgimag(A, lo, hi, tr)
end

function _heatmap(A::AbstractMatrix{<:AbstractFloat},
                  x1::PGFloat, x2::PGFloat, y1::PGFloat, y2::PGFloat;
                  vmin::Union{Nothing,Real} = nothing,
                  vmax::Union{Nothing,Real} = nothing,
                  badcolor::ColorOption = BAD_COLOR)
    x1, x2 = adjust_interval(x1, x2)
    y1, y2 = adjust_interval(y1, y2)
    select_figure(fig)
    # FIXME: `pgpixl` is much slower than `pgimag` so only call `pgpixl` if
    #        there are bad values.
    bad, lo, hi = rangeofvalues(eltype(A), A, vmin, vmax)
    memorize_zrange(lo, hi)
    if bad
        scaling = LinearScaling(badcolor=badcolor)
        pgpixl(rescale(scaling, A, lo, hi), x1, x2, y1, y2)
    else
        tr = simple_xform!(get_own_xform(), A, x1, x2, y1, y2)
        pgimag(A, lo, hi, tr)
    end
end

"""

```julia
heatmap!(A, tr=[0,1,0, 0,0,1]; kwds...)
```

draws a *heat-map* in the current frame preserving existing contents.
Arguments are the same as for [`heatmap`](@ref).  The only supported keywords
are `vmin`, `vmax` and `fig`.

"""
function heatmap!(A::AbstractMatrix{<:Real}, args...;
                  fig::FigureOption = nothing, kwds...)
    select_figure(fig)
    _heatmap(A, args...; kwds...)
end

# The following are to remember image settings for the color-bar.
const last_bg = Ref{PGFloat}(0)
const last_fg = Ref{PGFloat}(1)

memorize_zrange(bg::Real, fg::Real) =
    memorize_zrange(PGFloat(bg), PGFloat(fg))

function memorize_zrange(bg::PGFloat, fg::PGFloat)
    if bg == fg
        # Quick fix.
        bg -= tiny(bg)
        fg += tiny(fg)
    end
    last_bg[] = bg
    last_fg[] = fg
    return (bg, fg)
end

# Like get_vrange but fix equal endpoints and remember values
# for the color-bar.
function get_zrange(A::AbstractMatrix{PGFloat}, vmin, vmax)
    bg, fg = get_vrange(A, vmin, vmax)
    if bg == fg
        # Quick fix.
        bg -= tiny(bg)
        fg += tiny(fg)
    end
    last_bg[] = bg
    last_fg[] = fg
    return (bg, fg)
end

tiny(val::AbstractFloat) =
    (del = oftype(val, 0.001);
     val == zero(val) ? del : abs(val)*del)

adjust_interval(x1, x2) =
    (x1 == x2 ? (x1 - tiny(x1), x2 + tiny(x2)) : (x1, x2))

"""

```julia
colorbar(vmin, vmax)
```

draws a color-bar for values in the range `vmin` to `vmax`.  It is possible to
have `vmin > vmax` to reverse the order of the colors.  If unspecified, the
last values set by [`heatmap`](@ref) or [`heatmap!`](@ref) are considered.

Keywords:

- `side` specifies where to draw the color-bar: `side = 'L'` `'R'`, `'T'` or
  `'B'` to draw a color-bar on the Left, Right, Top or Bottom side of the
  frame. The default is to draw the color-bar on the right.

- `label` specifies a label for the color-bar.

- `spacing` specifies a distance breween the color-bar and the frame border in
  units of the current character size.  The value may be negative to draw the
  color-bar inside the frame.

- `width` specifies the thickness of the color-bar in units of the current
  character size.

"""
function colorbar(vmin::Real = last_bg[],
                  vmax::Real = last_fg[];
                  side::Char = 'R',
                  label::AbstractString = "",
                  spacing::Real = PGFloat(2),
                  width::Real = PGFloat(5))
    pgwedg(side*'I', spacing, width, vmin, vmax, label)
end

"""

```julia
mouse([mode = :none, [x0, y0]]; kwds...) -> (x, y, c)
```

yields cursor position `(x,y)` and character `c` typed by the user.  The
position is returned in world coordinates. If the device has no cursor or if
some other error occurs, the value `Char(0)` (ASCII NUL character) is returned
for `c`.  For an X-window device, clicking the 1st mouse button is reported as
an `A`, clicking the 2nd mouse button is reported as a `D` and clicking other
mouse buttons is reported as an `X`.

Argument `mode` can be used to specify the feedback for the user:

- If `mode = :none` (the default) no specific feedback is applied.

- If `mode = :line`, a straight line is drawn joining the anchor point and the
  cursor position.

- If `mode = :box`, a hollow rectangle is extended as the cursor is moved, with
  one vertex at the anchor point `(x0,y0)` and the opposite vertex at the
  current cursor position; the edges of the rectangle are horizontal and
  vertical.

- If `mode = :vrange`, two horizontal lines are extended across the width of
  the display, one drawn at ordinate `y0` of the anchor point and the other
  through the moving cursor position. This could be used to select a Y-axis
  range when one end of the range is known.

- If `mode = :hrange`, two vertical lines are extended over the height of the
  display, one drawn at abscissa `x0` of the anchor point and the other through
  the moving cursor position. This could be used to select an X-axis range when
  one end of the range is known.

- If `mode = :hline`, a horizontal line is extended through the cursor position
  over the width of the display. This could be used to select an X-axis value
  such as the start of an X-axis range. The anchor point is ignored.  :vline

- If `mode = :vline`, a vertical line is extended through the cursor position
  over the height of the display. This could be used to select a Y-axis value
  such as the start of a Y-axis range. The anchor point is ignored.  :hline

- If `mode = :cross`, a cross-hair, centered on the cursor, is extended over
  the width and height of the display. The anchor point is ignored.

Argument `(x0,y0)` specifies the position of the anchor point in world
coordinates.  The anchor point may be omitted if `mode` is not `:line`, `:box`,
`:vrange` or `:hrange`.

Keyword `fig` can be used to specify the figure to consider.  The default
is the current one.

Keyword `color` can be used to specify the color of the line(s) drawn for the
feedback.

The returned cursor position can be a `Point` if the first positional argument
is `Point` or if the anchor is specified as a `Point`.  For instance:

```julia
mouse(Point) -> Point(x, y), c
mouse(Point(x0,y0)) -> Point(x, y), c
```

"""
function mouse(mode::Symbol = :none;
               fig::FigureOption = nothing,
               color::ColorOption = nothing)
    md = (mode == :none  ?  PGInt(0) :
          mode == :hline ?  PGInt(5) :
          mode == :vline ?  PGInt(6) :
          mode == :cross ?  PGInt(7) :
          throw_invalid_mouse_mode(mode))
    select_figure(fig)
    mode == :none || pgsci(get_color(color))
    # Start with the (unused) anchor point at the center of the "window".
    xmin, xmax, ymin, ymax = pgqwin()
    pgband(md, false, PGFloat((xmin + xmax)/2), PGFloat((ymin + ymax)/2))
end

function mouse(mode::Symbol, x0::Real, y0::Real;
               fig::FigureOption = nothing,
               color::ColorOption = nothing)
    md = (mode == :none   ?  PGInt(0) :
          mode == :line   ?  PGInt(1) :
          mode == :box    ?  PGInt(2) :
          mode == :vrange ?  PGInt(3) :
          mode == :hrange ?  PGInt(4) :
          mode == :hline  ?  PGInt(5) :
          mode == :vline  ?  PGInt(6) :
          mode == :cross  ?  PGInt(7) :
          throw_invalid_mouse_mode(mode))
    select_figure(fig)
    mode == :none || pgsci(get_color(color))
    pgband(md, false, PGFloat(x0), PGFloat(y0))
end

function mouse(::Type{Point}, args...; kwds...)
    x, y, c = mouse(args...; kwds...)
    return Point(x, y), c
end

function mouse(P::Point; kwds...)
    x, y, c = mouse(P.x, P.y; kwds...)
    return Point(x, y), c
end

@noinline throw_invalid_mouse_mode(val) =
    throw(ArgumentError(val == :line || val == :box ||
                        val == :vrange || val == :hrange ?
                        string("invalid mouse mode ", val,
                               " or specify anchor point") :
                        string("invalid mouse mode ", val)))

"""

```julia
get_xform(arg) -> tr
```

yields coordinate transform `tr` directly usable by PGPlot for a variety of
arguments: `nothing` to get the default *identity* transform, a 2×3 matrix, a
6-element vector, a 6-tuple or an `AffineTransform` from the `TwoDimensional`
package.

"""
get_xform(::Nothing) = DEFAULT_XFORM
get_xform(A::AffineTransform) = PGFloat[A.x, A.xx, A.xy, A.y, A.yx, A.yy]
get_xform(A::NTuple{6,Real}) = PGFloat[A...]
get_xform(A::AbstractVector) = (@assert size(A) == (6,); pgarray(PGFloat, A))
get_xform(A::AbstractMatrix) = (@assert size(A) == (2,3); pgarray(PGFloat, A))

"""

```julia
splat_xform(T, tr) -> Ax,Axx,Axy, Ay,Ayx,Ayy
```

converts the coordinate transform `tr` into a 6-tuple of coefficients of
type `T`.

"""
function splat_xform(::Type{T},
                     tr::AbstractVector{<:Real}) where {T<:AbstractFloat}
    axes(tr,1) == 1:6 || throw_bad_xform_size()
    return (T(tr[1]), T(tr[2]), T(tr[3]),
            T(tr[4]), T(tr[5]), T(tr[6]))
end

function splat_xform(::Type{T},
                     tr::AbstractMatrix{<:Real}) where {T<:AbstractFloat}
    axes(tr) == (1:2, 1:3) || throw_bad_xform_size()
    return (T(tr[1,1]), T(tr[1,2]), T(tr[1,3]),
            T(tr[2,1]), T(tr[2,2]), T(tr[2,3]))
end

"""

```julia
get_extent(T, I, J, tr) -> x1, x2, y1, y2
```

yields the extent of world coordinates after applying coordinates transform
`tr` to the range of indices `I = imin:imax` and `J = jmin:jmax`.  Argument
`T` is the floating-point type of the result.  It can be replace by a
`Box{T}` type to get the resilt as a box instead of a 4-tuplke.

Ranges of indices can be specified by their endpoints:

```julia
get_extent(i1, i2, j1, j2, tr) -> x1, x2, y1, y2
```

In this latter case, the order of `i1` and `i2 (and of `j1` and `j2`) can be
reversed.

Other possibilities:

```julia
get_extent(T, A, tr)
get_extent(T, I, J, tr)
get_extent(T, i1, i2, j1, j2, tr)
```

or

```julia
get_extent(T, A, x1, x2, y1, y2)
get_extent(T, A, box)
```

"""
get_extent(::Type{Box{T}}, args...) where {T<:AbstractFloat} =
    Box{T}(get_extent(T, args...)...)

get_extent(::Type{T}, A::AbstractMatrix) where {T<:AbstractFloat} =
    get_extent(T, axes(A)...)

function get_extent(::Type{T},
                    A::AbstractMatrix,
                    B::Box) where {T<:AbstractFloat}
    return get_extent(T, A, B.x1, B.x2, B.y1, B.y2)
end

function get_extent(::Type{T},
                    A::AbstractMatrix,
                    x1::Real, x2::Real,
                    y1::Real, y2::Real) where {T<:AbstractFloat}
    return get_extent(T, A, T(x1), T(x2), T(y1), T(y2))
end

function get_extent(::Type{T},
                    A::AbstractMatrix,
                    x1::T, x2::T,
                    y1::T, y2::T) where {T<:AbstractFloat}
    return (adjust_interval(x1, x2)..., adjust_interval(y1, y2)...)
end

function get_extent(::Type{T},
                    A::AbstractMatrix,
                    tr::AbstractArray) where {T<:AbstractFloat}
    return get_extent(T, axes(A)..., tr)
end

function get_extent(::Type{T},
                    i1::Integer, i2::Integer,
                    j1::Integer, j2::Integer) where {T<:AbstractFloat}
    return (_get_extent(T, i1, i2)..., _get_extent(T, j1, j2)...)
end

function get_extent(::Type{T},
                    i1::Integer, i2::Integer,
                    j1::Integer, j2::Integer,
                    tr::AbstractArray{<:Real}) where {T<:AbstractFloat}
    _get_extent(T, _get_extent(T, i1, i2)..., _get_extent(T, j1, j2)..., tr)
end

function get_extent(::Type{T},
                    I::AbstractRange{<:Integer},
                    J::AbstractRange{<:Integer}) where {T<:AbstractFloat}
    return (_get_extent(T, I)..., _get_extent(T, J)...)
end

function get_extent(::Type{T},
                    I::AbstractRange{<:Integer},
                    J::AbstractRange{<:Integer},
                    tr::AbstractArray{<:Real}) where {T<:AbstractFloat}
    _get_extent(T, _get_extent(T, I)..., _get_extent(T, J)..., tr)
end

function _get_extent(::Type{T},
                     i1::Integer, i2::Integer) where {T<:AbstractFloat}
    s = (i1 ≤ i2 ? one(T) : -one(T))/T(2)
    return (T(i1) - s, T(I2) + s)
end

function _get_extent(::Type{T},
                     I::AbstractRange{<:Integer}) where {T<:AbstractFloat}
    s = T(step(I))/T(2)
    return (T(first(I)) - s, T(last(I)) + s)
end

function _get_extent(::Type{T},
                     x1::T, x2::T,
                     y1::T, y2::T,
                     tr::AbstractArray{<:Real}) where {T<:AbstractFloat}
    Ax,Axx,Axy, Ay,Ayx,Ayy = splat_xform(T, tr)
    x1p = Ax + Axx*x1 + Axy*y1
    y1p = Ay + Ayx*x1 + Ayy*y1
    x2p = Ax + Axx*x2 + Axy*y1
    y2p = Ay + Ayx*x2 + Ayy*y1
    x3p = Ax + Axx*x1 + Axy*y2
    y3p = Ay + Ayx*x1 + Ayy*y2
    x4p = Ax + Axx*x2 + Axy*y2
    y4p = Ay + Ayx*x2 + Ayy*y2
    xmin = min(x1p, x2p, x3p, x4p)
    xmax = max(x1p, x2p, x3p, x4p)
    ymin = min(y1p, y2p, y3p, y4p)
    ymax = max(y1p, y2p, y3p, y4p)
    return (xmin, xmax, ymin, ymax)
end

"""

```julia
simple_xform!(tr, A, x1p, x2p, y1p, y2p) -> tr
```

overwrites `tr` with the coefficients of the coordinate transform which maps
the indices of the 2-dimensional array `A` to the rectangular region defined by
`x1p`, `x2p`, `y1p` and `y2p`.

Another possibility is to specify the ranges of indices `I` and `J`:

```julia
simple_xform!(tr, I, J, x1p, x2p, y1p, y2p) -> tr
```

Yet another possibility is to specify the limits `x1`, `x2`, `y1` and `y2` of
the rectangular region to map to the output region:

```julia
simple_xform!(tr, x1, x2, y1, y2, x1p, x2p, y1p, y2p) -> tr
```

"""
function simple_xform!(tr::AbstractArray{T},
                       A::AbstractMatrix,
                       x1p::Real, x2p::Real,
                       y1p::Real, y2p::Real) where {T<:AbstractFloat}
    simple_xform!(tr, axes(A)..., x1p, x2p, y1p, y2p)
end

function simple_xform!(tr::AbstractArray{T},
                       I::AbstractRange{<:Integer},
                       J::AbstractRange{<:Integer},
                       x1p::Real, x2p::Real,
                       y1p::Real, y2p::Real) where {T<:AbstractFloat}
    simple_xform!(tr, I, J, T(x1p), T(x2p), T(y1p), T(y2p))
end

function simple_xform!(tr::AbstractArray{T},
                       I::AbstractRange{<:Integer},
                       J::AbstractRange{<:Integer},
                       x1p::T, x2p::T,
                       y1p::T, y2p::T) where {T<:AbstractFloat}
    di = T(step(I)/2)
    x1 = T(first(I)) - di
    x2 = T(last(I)) + di
    dj = T(step(J)/2)
    y1 = T(first(J)) - dj
    y2 = T(last(J)) + dj
    simple_xform!(tr, x1, x2, y1, y2, x1p, x2p, y1p, y2p)
end

function simple_xform!(tr::DenseArray{T},
                       x1::Real, x2::Real,
                       y1::Real, y2::Real,
                       x1p::Real, x2p::Real,
                       y1p::Real, y2p::Real) where {T<:AbstractFloat}
    simple_xform!(tr, T(x1), T(x2), T(y1), T(y2),
                  T(x1p), T(x2p), T(y1p), T(y2p))
end

function simple_xform!(tr::DenseArray{T},
                       x1::T, x2::T,
                       y1::T, y2::T,
                       x1p::T, x2p::T,
                       y1p::T, y2p::T) where {T<:AbstractFloat}
    @assert length(tr) == 6
    ax = (x2p - x1p)/(x2 - x1)
    bx = (x1p*x2 - x2p*x1)/(x2 - x1)
    @assert isfinite(ax) && isfinite(bx)
    ay = (y2p - y1p)/(y2 - y1)
    by = (y1p*y2 - y2p*y1)/(y2 - y1)
    @assert isfinite(ay) && isfinite(by)
    tr[1] = bx
    tr[2] = ax
    tr[3] = zero(T)
    tr[4] = by
    tr[5] = zero(T)
    tr[6] = ay
    return tr
end

"""

```julia
get_vrange([T=eltype(A),] A, lo, hi) -> vmin, vmax
```

yields the range of values in array `A` given optional low and high clipping
levels `lo` and `hi` which both can be either a value or `nothing` to indicate
that the corrsponding extreme (minimum or maximum) value in `A` should be
taken.  The result is a 2-tuple of values of type `T` which is, by default,
that of the elements of `A`.

"""
get_vrange(A::AbstractArray{T}, lo, hi) where {T} =
    get_vrange(T, A, lo, hi)
get_vrange(::Type{T}, A::AbstractArray, lo::Nothing, hi::Nothing) where {T} =
    extrema(A)
get_vrange(::Type{T}, A::AbstractArray, lo::Real, hi::Nothing) where {T} =
    (T(lo), maximum(A))
get_vrange(::Type{T}, A::AbstractArray, lo::Nothing, hi::Real) where {T} =
    (minimum(A), T(hi))
get_vrange(::Type{T}, A::AbstractArray, lo::Real, hi::Real) where {T} =
    (T(lo), T(hi))

function get_vrange(A::AbstractArray{T},
                    opt::Union{Nothing,Tuple},
                    I::NTuple{2,Integer}) where {T}
    get_vrange(T, A, opt, I)
end

function get_vrange(::Type{T}, A::AbstractArray, opt::Nothing,
                    I::NTuple{2,Integer}) where {T}
    get_vrange(T, A, nothing, nothing)
end

function get_vrange(::Type{T}, A::AbstractArray, opt::Tuple,
                    I::NTuple{2,Integer}) where {T}
    get_vrange(T, A, opt[I[1]], opt[I[2]])
end

nicerange(x::Tuple{Real,Real}) = nicerange(x[1], x[2])
nicerange(x1::Real, x2::Real) = nicerange(PGFloat(x1), PGFloat(x2))
function nicerange(x1::PGFloat, x2::PGFloat)
    @assert isfinite(x1) && isfinite(x2)
    if x1 != x2
        return pgrnge(x1, x2)
    elseif x1 == 0
        return (PGFloat(-0.001), PGFloat(+0.001))
    else
        d = PGFloat(10^(round(log10(abs(x1))) - 1))
        return (-d, +d)
    end
end

get_label(::Nothing) = ""
get_label(val::AbstractString) = val
get_label(val::Symbol) = string(val)

get_axis_style(::Nothing) = PGInt(1) # default axis style
get_axis_style(val::Integer) =
    (-2 ≤ val ≤ 2 || val == 10 || val == 20 || val == 30 ? PGInt(val) :
     throw_bad_axis_style(val))
get_axis_style(val::Symbol) =
    # Draw no box, axes or labels:
    (val == :none ? PGInt(-2) :
     # Draw box only:
     val == :box ? PGInt(-1) :
     # Draw box and label it with coordinates:
     val == :label ? PGInt(0) :
     # Same as AXIS=0, but also draw the coordinate axes (X=0, Y=0):
     val == :axis ? PGInt(1) :
     # Same as AXIS=1, but also draw grid lines at major increments of the
     # coordinates:
     val == :grid ? PGInt(2) :
     # Draw box and label X-axis logarithmically:
     val == :loglin ? PGInt(10) :
     # Draw box and label Y-axis logarithmically:
     val == :linlog ? PGInt(20) :
     # Draw box and label both axes logarithmically:
     val == :loglog ? PGInt(30) :
     throw_bad_axis_style(val))

@noinline throw_bad_axis_style(val) =
    throw(ArgumentError(string("invalid axis style: ", val)))

get_color(::Nothing) = FOREGROUND_COLOR_INDEX
get_color(val, def) = get_color(val)
get_color(::Nothing, def) = get_color(def)
get_color(ci::Integer) = PGInt(ci)
get_color(val::UInt32) = get_color(RGB24(val))
get_color(val::Symbol) = get_color(String(val))
get_color(val::AbstractString) = get_color(parse(RGB{PGFloat}, val))
get_color(val::RGBVec) = get_color(convert(RGB{PGFloat}, val))
function get_color(val::Colorant)
    pgscr(USER_COLOR_INDEX, val)
    return USER_COLOR_INDEX
end

get_text_height(val::Real) =
    (@assert isfinite(val) && val > 0; PGFloat(val))

get_font(val::AbstractString) =
  (val == "normal" ? NORMAL_FONT :
   val == "roman"  ? ROMAN_FONT  :
   val == "italic" ? ITALIC_FONT :
   val == "script" ? SCRIPT_FONT : throw_bad_text_font(val))
get_font(val::Symbol) =
  (val == :normal ? NORMAL_FONT :
   val == :roman  ? ROMAN_FONT  :
   val == :italic ? ITALIC_FONT :
   val == :script ? SCRIPT_FONT : throw_bad_text_font(val))
get_font(val::Integer) =
  (val == 1 ? NORMAL_FONT  :
   val == 2 ? ROMAN_FONT   :
   val == 3 ? ITALIC_FONT  :
   val == 4 ? SCRIPT_FONT  : throw_bad_text_font(val))

@noinline throw_bad_text_font(val) =
    throw(ArgumentError(string("invalid text font: ", val)))

"""

```julia
extremefinitevalues(A) -> (bad, amin, amax)
```

yields the extreme finite values in array `A` and a boolean indicating
whether `A` has bad (i.e., non-finite values).  The returned values `amin`
and `amax` are such that `amin ≤ amax` unless there are no finite values in
`A` in which case `(true,typemax(T),typemin(T))` is returned, with `T` the
type of the elements of `A`.

"""
function extremefinitevalues(A::AbstractArray{T}) where {T<:Real}
    amin = typemax(T)
    amax = typemin(T)
    flag = true
    @inbounds @simd for i in eachindex(A)
        val = A[i]
        flag &= isfinite(val)
        amin = ifelse(isfinite(val) & (val < amin), val, amin)
        amax = ifelse(isfinite(val) & (val > amax), val, amax)
    end
    return (!flag, amin, amax)
end

function extremefinitevalues(A::AbstractArray{T}) where {T<:Integer}
    amin = typemax(T)
    amax = typemin(T)
    @inbounds @simd for i in eachindex(A)
        val = A[i]
        amin = ifelse(val < amin, val, amin)
        amax = ifelse(val > amax, val, amax)
    end
    return (false, amin, amax)
end

"""

```julia
minimumfinitevalue(A) -> (bad, amin)
```

yields the smallest finite value in array `A` and a boolean indicating
whether `A` has bad (i.e., non-finite values).  If the elements of `A` are
of a floating-point type `T`, then `amin` is such that `amin < Inf` unless
there are no finite values in `A` in which case `(true,typemax(T))`, that
is `(true,+T(Inf))`, is returned.

"""
function minimumfinitevalue(A::AbstractArray{T}) where {T<:Real}
    amin = typemax(T)
    flag = true
    @inbounds @simd for i in eachindex(A)
        val = A[i]
        flag &= isfinite(val)
        amin = ifelse(isfinite(val) & (val < amin), val, amin)
    end
    return (!flag, amin)
end

function minimumfinitevalue(A::AbstractArray{T}) where {T<:Integer}
    amin = typemax(T)
    @inbounds @simd for i in eachindex(A)
        val = A[i]
        amin = ifelse(val < amin, val, amin)
    end
    return (false, amin)
end

"""

```julia
maximumfinitevalue(A) -> (bad, amax)
```

yields the largest finite value in array `A` and a boolean indicating
whether `A` has bad (i.e., non-finite values).  If the elements of `A` are
of a floating-point type `T`, then `amax` is such that `amax >-Inf` unless
there are no finite values in `A` in which case `(true,typemin(T))`, that
is `(true,-T(Inf))`, is returned.

"""
function maximumfinitevalue(A::AbstractArray{T}) where {T<:Real}
    amax = typemin(T)
    flag = true
    @inbounds @simd for i in eachindex(A)
        val = A[i]
        flag &= isfinite(val)
        amax = ifelse(isfinite(val) & (val > amax), val, amax)
    end
    return (!flag, amax)
end

function maximumfinitevalue(A::AbstractArray{T}) where {T<:Integer}
    amax = typemin(T)
    @inbounds @simd for i in eachindex(A)
        val = A[i]
        amax = ifelse(val > amax, val, amax)
    end
    return (false, amax)
end

"""

```julia
rangeofvalues(T, A, vmin, vmax) -> (bad, amin, amax)
```

"""
function rangeofvalues(::Type{T}, A::AbstractArray,
                       ::Nothing, ::Nothing) where {T}
    bad, amin, amax = extremefinitevalues(A)
    return (bad, T(amin), T(amax))
end

function rangeofvalues(::Type{T}, A::AbstractArray,
                       vmin::Real, ::Nothing) where {T}
    if isfinite(vmin)
        bad, vmax = maximumfinitevalue(A)
        return (bad, T(vmin), T(vmax))
    else
        isnan(vmin) && throw(ArgumentError("vmin must not be a `NaN`"))
        bad, amin, amax = extremefinitevalues(A)
        if vmin > 0
            # Exchange min. and max.
            amin, amax = amax, amin
        end
        return (bad, T(amin), T(amax))
    end
end

function rangeofvalues(::Type{T}, A::AbstractArray,
                       ::Nothing, vmax::Real) where {T}
    if isfinite(vmax)
        bad, vmin = minimumfinitevalue(A)
        return (bad, T(vmin), T(vmax))
    else
        isnan(vmax) && throw(ArgumentError("vmax must not be a `NaN`"))
        bad, amin, amax = extremefinitevalues(A)
        if vmax < 0
            # Exchange min. and max.
            amin, amax = amax, amin
        end
        return (bad, T(amin), T(amax))
    end
end

function rangeofvalues(::Type{T}, A::AbstractArray,
                       vmin::Real, vmax::Real) where {T}
    choice = (isfinite(vmin) ? 1 : 0) | (isfinite(vmax) ? 2 : 0)
    if choice == 3
        return (have_bad_values(A), T(vmin), T(vmax))
    elseif choice == 2
        # vmin is not finite
        isnan(vmin) && throw(ArgumentError("vmin must not be a `NaN`"))
        if vmin ≤ vmax
            bad, amin = minimumfinitevalue(A)
            return (bad, T(amin), T(vmax))
        else
            bad, amax = maximumfinitevalue(A)
            return (bad, T(amax), T(vmax))
        end
    elseif choice == 1
        # vmax is not finite
        isnan(vmax) && throw(ArgumentError("vmax must not be a `NaN`"))
        if vmin ≤ vmax
            bad, amax = maximumfinitevalue(A)
            return (bad, T(vmin), T(amax))
        else
            bad, amin = minimumfinitevalue(A)
            return (bad, T(vmin), T(amin))
        end
    else
        isnan(vmin) && throw(ArgumentError("vmin must not be a `NaN`"))
        isnan(vmax) && throw(ArgumentError("vmax must not be a `NaN`"))
        bad, amin, amax = extremefinitevalues(A)
        if vmax < vmin
            # Exchange min. and max.
            amin, amax = amax, amin
        end
        return (bad, T(amin), T(amax))
    end
end

"""

```julia
rescale(scl, A, vmin, vmax)
```

yields an array of element type `T` (must be an integer) which maps values of
array `A` from the range `vmin:vmax` to `first(scl):last(scl)`.  It is possible
to have `vmin > vmax` or `first(scl) > last(scl)` to reverse the scaling.

If `A` is of floating point type, `bad(scl)` is the value given to bad numbers
(i.e., *NaN*) and, if `vmin` and `vmax` are not both finite, the result is
filled with `bad(scl)`.

See [`rescale!`](@ref) for providing the destination array and
[`rangeofvalues`](@ref) as a mean to determine `vmin` and `vmax`.

"""
function rescale(scl::AbstractScaling{Ti},
                 A::AbstractArray{<:Integer,N},
                 vmin::Real, vmax::Real) where {Ti<:Integer,N}
    Tf = float(promote_type(typeof(vmin), typeof(vmax)))
    rescale!(Array{Ti,N}(undef, size(A)), scl, A, Tf(vmin), Tf(vmax))
end

function rescale(scl::AbstractScaling{Ti},
                 A::AbstractArray{Tf,N},
                 vmin::Real, vmax::Real) where {Ti<:Integer,
                                                Tf<:AbstractFloat,N}
    rescale!(Array{Ti,N}(undef, size(A)), scl, A, Tf(vmin), Tf(vmax))
end

"""

```julia
rescale!(dst, scl, A, vmin, vmax) -> dst
```

overwrites destination array `dst` with the result or rescaling the values in
array `A` from the range `vmin:vmax` and according to `scl`.  See
[`rescale`](@ref) for details.

"""
function rescale!(dst::DenseArray{Ti,N},
                  scl::AbstractScaling{Ti},
                  A::AbstractArray{<:Integer,N},
                  vmin::Real, vmax::Real) where {Ti<:Integer,N}
    Tf = float(promote_type(typeof(vmin), typeof(vmax)))
    rescale!(dst, scl, A, Tf(vmin), Tf(vmax))
end

function rescale!(dst::DenseArray{Ti,N},
                  scl::AbstractScaling{Ti},
                  A::AbstractArray{Tf,N},
                  vmin::Real, vmax::Real) where {Ti<:Integer,
                                                 Tf<:AbstractFloat,N}
    rescale!(dst, scl, A, Tf(vmin), Tf(vmax))
end

function rescale!(dst::DenseArray{Ti,N},
                  scl::AbstractScaling{Ti},
                  A::AbstractArray{<:Integer,N},
                  vmin::Tf, vmax::Tf) where {Ti<:Integer,
                                             Tf<:AbstractFloat,N}
    throw_not_implemented(typeof(scl))
end

function rescale!(dst::DenseArray{Ti,N},
                  scl::AbstractScaling{Ti},
                  A::AbstractArray{Tf,N},
                  vmin::Tf, vmax::Tf) where {Ti<:Integer,
                                             Tf<:AbstractFloat,N}
    throw_not_implemented(typeof(scl))
end

@noinline throw_not_implemented(::Type{T}) where {T<:AbstractScaling} =
    error("scaling of type $T is not implemented")

function rescale!(dst::DenseArray{Ti,N},
                  scl::LinearScaling{Ti},
                  A::AbstractArray{<:Integer,N},
                  vmin::Tf, vmax::Tf) where {Ti<:Integer,
                                             Tf<:AbstractFloat,N}
    cmin, cmax = first(scl), last(scl)
    if vmin == vmax
        cmid = ((cmin + cmax) ÷ 2)
        @inbounds for i in all_indices(dst, A)
            val = A[i]
            dst[i] = (val > vmax ? cmax :
                      val < vmin ? cmin : cmid)
        end
    else
        a = (Tf(cmax) - Tf(cmin))/(vmax - vmin)
        b = Tf(cmin) - a*vmin
        @inbounds for i in all_indices(dst, A)
            val = Tf(A[i])
            dst[i] = (val ≥ vmax ? cmax :
                      val ≤ vmin ? cmin : round(Ti, a*val + b))
        end
    end
    return dst
end

function rescale!(dst::DenseArray{Ti,N},
                  scl::LinearScaling{Ti},
                  A::AbstractArray{Tf,N},
                  vmin::Tf, vmax::Tf) where {Ti<:Integer,
                                             Tf<:AbstractFloat,N}
    cmin, cmax, cbad = first(scl), last(scl), bad(scl)
    if isfinite(vmin) && isfinite(vmin)
        if vmin == vmax
            cmid = ((cmin + cmax) ÷ 2)
            @inbounds for i in all_indices(dst, A)
                val = A[i]
                dst[i] = (isnan(val) ? cbad :
                          val > vmax ? cmax :
                          val < vmin ? cmin : cmid)
            end
        else
            a = (Tf(cmax) - Tf(cmin))/(vmax - vmin)
            b = Tf(cmin) - a*vmin
            @inbounds for i in all_indices(dst, A)
                val = A[i]
                dst[i] = (isnan(val) ? cbad :
                          val ≥ vmax ? cmax :
                          val ≤ vmin ? cmin : round(Ti, a*val + b))
            end
        end
    else
        # There are no finite values in `A`.
        fill!(dst, cbad)
    end
    return dst
end

end # module


text:
################################################################################
###################### DEFINE MASTER PARAMETERS FOR THIS FILE ##################
################################################################################

using Mimi
using Distributions
using CSVFiles
using DataFrames
using CSV
using Random
using StatsBase
using Statistics

# specify model settings
function set_globalbools()
    global use_variability = true

    # set random seed to have similar variability development in the base and the marginal model.
    # set variability seed.
    if use_variability
        global varseed = rand(1:1000000000000)
    end

    global use_linear = false
    global use_logburke = false
    global use_logpopulation = false
    global use_logwherepossible = true
end

set_globalbools()

include("main_model_annualGrowth.jl")

# define the output directory
dir_output = joinpath(@__DIR__, "../../output/")

# define number of Monte Carlo runs
if !@isdefined samplesize
    samplesize = 50000
end

# define the seed
masterseed = 22081994

# define the pulse size
scc_pulse_size = 75000.

################################################################################
################### MC MODEL RUNS    ###########################################
################################################################################

# PART II: triangular rho (growth effects) distributions
df_sccMC = DataFrame(permafr=false, seaice=false, ge_string="-999", scen="-999",
                                  convergence=-999.,  bound=-999., eqwshare=-999.,
                                  civvalue=-999., pulse=-999.,
                                  mean=-999., median=-999., min=-999., max=-999., perc25=-999.,
                                  perc75=-999., sd=-999., varcoeff=-999.,
                                  perc05=-999., perc95=-999., perc10=-999., perc90=-999.)

# get the SCC for three different growth effects distributions and scenarios
for jj_scen in ["RCP4.5 & SSP2", "RCP2.6 & SSP1", "RCP8.5 & SSP5", "1.5 degC Target"]
    for jj_gestring in ["EMPIRICAL"]
        for jj_permafr in [true]
            for jj_seaice in [true]
                for jj_civvalue in [1.]
                    for jj_cbabs in [740.65]
                        for jj_eqwshare in [0.99]
                            for jj_convergence in [1.]
                                for jj_pulse in [scc_pulse_size]

                                    # set the ge_string parameters
                                    if jj_gestring == "EMPIRICAL"
                                        ge_string_min = 0.
                                        ge_string_mode = 0.
                                        ge_string_max = 0.0001
                                        ge_use_empirical = 1.
                                    end

                                    # print out the parameters to track progress
                                    print(string("Scen_", jj_scen, " Ge-mode_", jj_gestring,
                                                 " Permafr_", jj_permafr, " Seaice_", jj_seaice, " Civvalue_", jj_civvalue,
                                                 " cbabs_", jj_cbabs, " eqwshare_", jj_eqwshare,
                                                 " convergence", jj_convergence, "pulse_", jj_pulse))

                                    # define the output for the Monte Carlo files
                                    dir_MCoutput = string(dir_output, "mc_diGE/ge", jj_gestring,
                                                                    "_scen", jj_scen,
                                                                    "_pf", jj_permafr,
                                                                    "_se", jj_seaice,
                                                                    "_co", jj_convergence, "_bd", jj_cbabs, "_eq", jj_eqwshare,
                                                                    "_ci", jj_civvalue,
                                                                    "_p", jj_pulse,
                                                                     "/")

                                    # calculate the stochastic mean SCC
                                    Random.seed!(masterseed)
                                    global scc_mcs_object = get_scc_mcs(samplesize, 2020, dir_MCoutput,
                                                                        scenario=jj_scen,
                                                                        pulse_size=jj_pulse,
                                                                        use_permafrost=jj_permafr,
                                                                        use_seaice=jj_seaice,
                                                                        use_page09damages=false,
                                                                        ge_minimum=ge_string_min,
                                                                        ge_maximum=ge_string_max,
                                                                        ge_mode=ge_string_mode,
                                                                        ge_use_empirical=ge_use_empirical,
                                                                        civvalue_multiplier=jj_civvalue,
                                                                        use_convergence=jj_convergence,
                                                                        cbabs=jj_cbabs,
                                                                        eqwbound=jj_eqwshare)

                                    # write out the full distribution
                                    writedlm(string(dir_output, "SCC_MCS_scen", jj_scen, "_per", jj_permafr, "_sea", jj_seaice,
                                                            "_ge", jj_gestring,
                                                            "_conv", jj_convergence, "_bou", jj_cbabs, "_eqw", jj_eqwshare,
                                                            "_civ", jj_civvalue,
                                                            "_pul", jj_pulse,  ".csv"),
                                                            scc_mcs_object, ",")

                                    # push the results into the data frame
                                    push!(df_sccMC, [jj_permafr, jj_seaice, jj_gestring,  jj_scen,
                                                    jj_convergence,  jj_cbabs, jj_eqwshare,
                                                    jj_civvalue, jj_pulse,
                                                    mean(scc_mcs_object[:, 1]),
                                                    median(scc_mcs_object[:, 1]),
                                                    minimum(scc_mcs_object[:, 1]),
                                                    maximum(scc_mcs_object[:, 1]),
                                                    StatsBase.percentile(scc_mcs_object[:, 1], 25),
                                                    StatsBase.percentile(scc_mcs_object[:, 1], 75),
                                                    Statistics.std(scc_mcs_object[:, 1]),
                                                    Statistics.std(scc_mcs_object[:, 1]) / mean(scc_mcs_object[:, 1]),
                                                    StatsBase.percentile(scc_mcs_object[:, 1], 5),
                                                    StatsBase.percentile(scc_mcs_object[:, 1], 95),
                                                    StatsBase.percentile(scc_mcs_object[:, 1], 10),
                                                    StatsBase.percentile(scc_mcs_object[:, 1], 90)])

                                    # clean out the scc_mcs_object
                                    global scc_mcs_object = nothing
                                    ge_string_min = nothing
                                    ge_string_max = nothing
                                    ge_string_mode = nothing

                                end
                            end
                        end
                    end
                end
            end
        end
    end
end

# remove the first placeholder row
df_sccMC = df_sccMC[df_sccMC[:scen] .!= "-999", :]

# export the results
CSV.write(string(dir_output, "MimiPageGrowthEffectsResultsMonteCarlo.csv"), df_sccMC)

ext: jl


text:
module TestExamplesDGMulti2D

using Test
using Trixi

include("test_trixi.jl")

# pathof(Trixi) returns /path/to/Trixi/src/Trixi.jl, dirname gives the parent directory
EXAMPLES_DIR = joinpath(pathof(Trixi) |> dirname |> dirname, "examples", "dgmulti_2d")

# Start with a clean environment: remove Trixi output directory if it exists
outdir = "out"
isdir(outdir) && rm(outdir, recursive=true)

@testset "DGMulti 2D" begin
  @trixi_testset "elixir_euler_weakform.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_weakform.jl"),
      cells_per_dimension = (4, 4),
      l2 = [0.0013536930300254945, 0.0014315603442106193, 0.001431560344211359, 0.0047393341007602625],
      linf = [0.001514260921466004, 0.0020623991944839215, 0.002062399194485476, 0.004897700392503701],
    )
  end

  @trixi_testset "elixir_euler_weakform.jl (SBP)" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_weakform.jl"),
      cells_per_dimension = (4, 4),
      approximation_type = SBP(),
      l2 = [0.0074706882014934735, 0.005306220583603261, 0.005306220583613591, 0.014724842607716771],
      linf = [0.021563604940952885, 0.01359397832530762, 0.013593978324845324, 0.03270995869587523]
    )
  end

  @trixi_testset "elixir_euler_weakform.jl (Quadrilateral elements)" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_weakform.jl"),
      cells_per_dimension = (4, 4),
      element_type = Quad(),
      l2 = [0.00031892254415307093, 0.00033637562986771894, 0.0003363756298680649, 0.0011100259064243145],
      linf = [0.001073298211445639, 0.0013568139808282087, 0.0013568139808290969, 0.0032249020004324613]
    )
  end

  @trixi_testset "elixir_euler_weakform.jl (EC) " begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_weakform.jl"),
      cells_per_dimension = (4, 4),
      volume_integral = VolumeIntegralFluxDifferencing(flux_ranocha),
      surface_integral = SurfaceIntegralWeakForm(flux_ranocha),
      l2 = [0.007801417730672109, 0.00708583561714128, 0.0070858356171393, 0.015217574294198809],
      linf = [0.011572828457858897, 0.013965298735070686, 0.01396529873508534, 0.04227683691807904],
    )
  end

  @trixi_testset "elixir_euler_weakform.jl (SBP, EC)" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_weakform.jl"),
      cells_per_dimension = (4, 4),
      volume_integral = VolumeIntegralFluxDifferencing(flux_ranocha),
      surface_integral = SurfaceIntegralWeakForm(flux_ranocha),
      approximation_type = SBP(),
      l2 = [0.01280067571168776, 0.010607599608273302, 0.010607599608239775, 0.026408338014056548],
      linf = [0.037983023185674814, 0.05321027922533417, 0.05321027922608157, 0.13392025411844033],
    )
  end

  @trixi_testset "elixir_euler_weakform.jl (Quadrilateral elements, SBP, EC)" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_weakform.jl"),
      cells_per_dimension = (4, 4),
      element_type = Quad(),
      volume_integral = VolumeIntegralFluxDifferencing(flux_ranocha),
      surface_integral = SurfaceIntegralWeakForm(flux_ranocha),
      approximation_type = SBP(),
      l2 = [0.0029373718090697975, 0.0030629360605489465, 0.003062936060545615, 0.0068486089344859755],
      linf = [0.01360165305316885, 0.01267402847925303, 0.012674028479251254, 0.02210545278615017],
    )
  end

  @trixi_testset "elixir_euler_weakform.jl (convergence)" begin
    mean_convergence = convergence_test(@__MODULE__, joinpath(EXAMPLES_DIR, "elixir_euler_weakform.jl"), 2)
    @test isapprox(mean_convergence[:l2], [4.243843382379403, 4.128314378833922, 4.128314378397532, 4.081366752807379], rtol=0.05)
  end

  @trixi_testset "elixir_euler_weakform_periodic.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_weakform_periodic.jl"),
      l2 = [0.0014986508075708323, 0.001528523420746786, 0.0015285234207473158, 0.004846505183839211],
      linf = [0.0015062108658376872, 0.0019373508504645365, 0.0019373508504538783, 0.004742686826709086]
    )
  end

  @trixi_testset "elixir_euler_triangulate_pkg_mesh.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_triangulate_pkg_mesh.jl"),
      l2 = [4.664661209491976e-6, 3.7033509525940745e-6, 4.794877426562555e-6, 1.2682723101532175e-5],
      linf = [2.5099852761334418e-5, 2.2683684021362893e-5, 2.6180448559287584e-5, 5.5752932611508044e-5]
    )
  end

  @trixi_testset "elixir_euler_kelvin_helmholtz_instability.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_kelvin_helmholtz_instability.jl"),
      cells_per_dimension = (32, 32), tspan = (0.0, 0.2),
      l2 = [0.11140378947116614, 0.06598161188703612, 0.10448953167839563, 0.16023209181809595],
      linf = [0.24033843177853664, 0.1659992245272325, 0.1235468309508845, 0.26911424973147735],
    )
  end

  @trixi_testset "elixir_euler_kelvin_helmholtz_instability.jl (Quadrilateral elements, GaussSBP)" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_kelvin_helmholtz_instability.jl"),
      cells_per_dimension = (32, 32), element_type = Quad(), approximation_type=GaussSBP(), tspan = (0.0, 0.2),
      l2 = [0.11141270656347146, 0.06598888014584121, 0.1044902203749932, 0.16023037364774995],
      linf = [0.2414760062126462, 0.1662111846065654, 0.12344140473946856, 0.26978428189564774]
    )
  end

  @trixi_testset "elixir_euler_rayleigh_taylor_instability.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_rayleigh_taylor_instability.jl"),
      cells_per_dimension = (8, 8), tspan = (0.0, 0.2),
      l2 = [0.03548329484912729, 0.002591414376082683, 0.006916327792623457, 0.016235069002818153],
      linf = [0.4783963902824797, 0.022527207050681054, 0.040307056293369226, 0.0852365428206836],
    )
  end

  @trixi_testset "elixir_euler_brown_minion_vortex.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_brown_minion_vortex.jl"),
      cells_per_dimension = 4, tspan = (0.0, 0.1),
      l2 = [0.0066800016110776066, 0.021516763475855016, 0.01069652423536524, 0.15052841129693573],
      linf = [0.01544756362800248, 0.09517304772476806, 0.021957154972646383, 0.33773439650806303],
    )
  end

  @trixi_testset "elixir_mhd_weak_blast_wave.jl (Quad)" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_weak_blast_wave.jl"),
      cells_per_dimension = 4,
      l2 = [0.1562707966203802, 0.19692319035938727, 0.19692319035938863, 0.10641395363893119,
            0.7221962824695998, 0.07678102773723876, 0.07678102773723913, 0.11142445367787217,
            0.006627119909705745],
      linf = [0.16447597822733662, 0.244157345789029, 0.24415734578903472, 0.11982440036793476,
              0.7450328339751362, 0.06357382685763713, 0.0635738268576378, 0.1058830287485999,
              0.005740591170062146],
    )
  end

  @trixi_testset "elixir_mhd_weak_blast_wave.jl (Tri)" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_weak_blast_wave.jl"),
      cells_per_dimension = 4, element_type = Tri(),
      l2 = [0.1348987236501758, 0.1588650593363661, 0.15886505933636905, 0.08571028568614296,
            0.6031732736338957, 0.06286520146660214, 0.06286520146660948, 0.0894610083051161,
            0.003693188791960107],
      linf = [0.16290247390873458, 0.2256891306641319, 0.2256891306641336, 0.09476017042552534,
              0.6906308908961734, 0.05349939593012487, 0.05349939593013042, 0.08830587480616725,
              0.0029551359803035027]
    )
  end

  @trixi_testset "elixir_mhd_weak_blast_wave_SBP.jl (Quad)" begin
    # These setups do not pass CI reliably, see
    # https://github.com/trixi-framework/Trixi.jl/pull/880 and
    # https://github.com/trixi-framework/Trixi.jl/issues/881
    @test_skip @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_weak_blast_wave_SBP.jl"),
      cells_per_dimension = 4,
      l2 = [0.15825983698241494, 0.19897219694837923, 0.19784182473275247, 0.10482833997417325,
            0.7310752391255246, 0.07374056714564853, 0.07371172293240634, 0.10782032253431281,
            0.004921676235111545],
      linf = [0.1765644464978685, 0.2627803272865769, 0.26358136695848144, 0.12347681727447984,
              0.7733289736898254, 0.06695360844467957, 0.06650382120802623, 0.10885097000919097,
              0.007212567638078835],
    )
  end

  @trixi_testset "elixir_mhd_weak_blast_wave_SBP.jl (Tri)" begin
    # These setups do not pass CI reliably, see
    # https://github.com/trixi-framework/Trixi.jl/pull/880 and
    # https://github.com/trixi-framework/Trixi.jl/issues/881
    @test_skip @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_weak_blast_wave_SBP.jl"),
      cells_per_dimension = 4, element_type=Tri(), tspan = (0.0, 0.2),
      l2 = [0.13825044764021147, 0.15472815448314997, 0.1549093274293255, 0.053103596213755405,
            0.7246162776815603, 0.07730777596615901, 0.07733438386480523, 0.109893463921706,
            0.00617678167062838],
      linf = [0.22701306227317952, 0.2905255794821543, 0.2912409425436937, 0.08051361477962096,
              1.0842974228656006, 0.07866000368926784, 0.0786646354518149, 0.1614896380292925,
              0.010358210347485542],
    )
  end

  @trixi_testset "elixir_shallowwater_source_terms.jl (Quad, SBP)" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_source_terms.jl"),
      cells_per_dimension = 8, element_type = Quad(), approximation_type = SBP(),
      l2 = [0.0028731817391463733, 0.033473047741338476, 0.04873655357923285, 2.7341648767587304e-15],
      linf = [0.010385010095182778, 0.08750628939565086, 0.12088392994348407, 9.325873406851315e-15]
    )
  end

  @trixi_testset "elixir_shallowwater_source_terms.jl (Tri, SBP)" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_source_terms.jl"),
      cells_per_dimension = 8, element_type = Tri(), approximation_type = SBP(),
      l2 = [0.005912374812011487, 0.0993653660027428, 0.1670914777143732, 3.2948191109718012e-15],
      linf = [0.02076003852980346, 0.29169601664914424, 0.5674183379872275, 1.1546319456101628e-14]
    )
  end

  @trixi_testset "elixir_shallowwater_source_terms.jl (Tri, Polynomial)" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_source_terms.jl"),
      cells_per_dimension = 8, element_type = Tri(), approximation_type = Polynomial(),
      # The last l2, linf error are the L2 projection error in approximating `b`, so they are not
      # zero for general non-collocated quadrature rules (e.g., for `element_type=Tri()`, `polydeg > 2`).
      l2 = [0.0011751205240013974, 0.021530712606619058, 0.022675456966150455, 1.813056761616414e-5],
      linf = [0.001888045014140971, 0.05466838692127718, 0.06345885709961152, 3.3989933098554914e-5]
    )
  end

  @trixi_testset "elixir_shallowwater_source_terms.jl (Quad, Polynomial)" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_source_terms.jl"),
      cells_per_dimension = 8, element_type = Quad(), approximation_type = Polynomial(),
      # The last l2, linf error are the L2 projection error in approximating `b`. However, this is zero
      # for `Quad()` elements with `Polynomial()` approximations because the quadrature rule defaults to
      # a `(polydeg + 1)`-point Gauss quadrature rule in each coordinate (in general, StartUpDG.jl defaults
      # to the quadrature rule with the fewest number of points which exactly integrates the mass matrix).
      l2 = [0.00010550686466016327, 0.005212211092466598, 0.005529801974869844, 3.062184075681645e-15],
      linf = [0.000259995400729629, 0.0072236204211630906, 0.010364675200833062, 1.021405182655144e-14]
    )
  end

end

# Clean up afterwards: delete Trixi output directory
@test_nowarn isdir(outdir) && rm(outdir, recursive=true)

end # module


text:
function __init__()
    # set up REPL mode
    isdefined(Base, :active_repl) &&
        isinteractive() && typeof(Base.active_repl) != REPL.BasicREPL &&
            !DouglassPrompt.repl_inited(Base.active_repl) && DouglassPrompt.repl_init(Base.active_repl)
end


text:
"""
This function creates the MarketClearingProblem struct for CEM.
"""
function create_cem_mkt_clr_problem(investor_dir::String,
                                    market_names::Vector{Symbol},
                                    expected_portfolio::Vector{<: Project{<: BuildPhase}},
                                    zones::Vector{String},
                                    lines::Vector{ZonalLine},
                                    peak_load::Float64,
                                    rep_hour_weight::Vector{Float64},
                                    average_capital_cost_multiplier::Float64,
                                    scenario::Scenario,
                                    iteration_year::Int64,
                                    yearly_horizon::Int64)

    num_invperiods = yearly_horizon

    #Get updated load growth belief
    load_growth = AxisArrays.AxisArray(zeros(length(zones), num_invperiods), zones, collect(1:num_invperiods))
    for (idx, zone) in enumerate(zones)
        load_growth[idx, :] = get_parameter_values(scenario)[iteration_year]["load_$(zone)", iteration_year:iteration_year + num_invperiods - 1]
    end

    reserve_up_market_bool = false
    reserve_down_market_bool = false
    capacity_market_bool = false
    rec_market_bool = false

    # Find if the investor participates in the following markets:
    if in(:ReserveUp, market_names)
        reserve_up_market_bool = true
    end

    if in(:ReserveDown, market_names)
        reserve_down_market_bool = true
    end

    if in(:Capacity, market_names)
        capacity_market_bool = true
    end

    if in(:REC, market_names)
        rec_market_bool = true
    end

    # Gather markets data-------------------------------------------------------------------------------
    load_data = read_data(joinpath(investor_dir, "timeseries_data_files", "Load", "load_$(iteration_year - 1).csv"))
    reserve_up_demand_data = read_data(joinpath(investor_dir, "timeseries_data_files", "Reserves", "reserve_up_$(iteration_year - 1).csv"))
    reserve_down_demand_data = read_data(joinpath(investor_dir, "timeseries_data_files", "Reserves", "reserve_down_$(iteration_year - 1).csv"))

    num_hours = DataFrames.nrow(load_data)

    energy_mkt_params = read_data(joinpath(investor_dir, "markets_data", "energy_mkt_param.csv"))
    price_cap_energy = AxisArrays.AxisArray(energy_mkt_params.price_cap * 1.0, zones)

    reserve_up_mkt_params = read_data(joinpath(investor_dir, "markets_data", "reserve_up_mkt_param.csv"))
    reserve_up_num_points = AxisArrays.AxisArray(zeros(Int64,length(zones)), zones)
    reserve_up_break_points = AxisArrays.AxisArray([Float64[] for z in zones], zones)
    reserve_up_price_points = AxisArrays.AxisArray([Float64[] for z in zones], zones)

    reserve_down_mkt_params = read_data(joinpath(investor_dir, "markets_data", "reserve_down_mkt_param.csv"))
    price_cap_reservedown = AxisArrays.AxisArray(reserve_down_mkt_params.price_cap * 1.0, zones)
    perc_req_reserve_down = AxisArrays.AxisArray(reserve_down_mkt_params.perc_req, zones)

    zonal_load = AxisArrays.AxisArray(zeros(length(zones), num_hours), zones, (1:num_hours))
    zonal_reserve_up = AxisArrays.AxisArray(zeros(length(zones), num_hours), zones, (1:num_hours))
    zonal_reserve_down = AxisArrays.AxisArray(zeros(length(zones), num_hours), zones, (1:num_hours))

    for (zone_num, zone) in enumerate(zones)
        idx = findfirst(x -> x == "$(zone)", reserve_up_mkt_params.zones)
        if reserve_up_mkt_params.ORDC[idx]
            reserve_up_num_points[zone] = reserve_up_mkt_params.ORDC_points[idx]
            for point in 1:reserve_up_num_points[zone]
                push!(reserve_up_break_points[zone], reserve_up_mkt_params[idx, Symbol("ORDC_x$(point)")] * reserve_up_market_bool)
                push!(reserve_up_price_points[zone], reserve_up_mkt_params[idx, Symbol("ORDC_y$(point)")] * reserve_up_market_bool)
            end
        else
            reserve_up_break_points[zone] = [0.0, 1.0, 1.0, 10.0] * reserve_up_market_bool
            reserve_up_price_points[zone] = [reserve_up_mkt_params.price_cap[idx], reserve_up_mkt_params.price_cap[idx], 0.0, 0.0] * reserve_up_market_bool
        end

        price_cap_reservedown[zone] = reserve_down_mkt_params.price_cap[idx]

        for h in 1:num_hours
            zonal_load[zone, h] = load_data[:, Symbol(zone_num)][h]
            zonal_reserve_up[zone, h] = reserve_up_demand_data[:, Symbol(zone_num)][h] * reserve_up_market_bool
            zonal_reserve_down[zone, h] = reserve_down_demand_data[:, Symbol(zone_num)][h] * reserve_down_market_bool
        end
    end

    energy_annual_increment = AxisArrays.AxisArray(ones(length(zones), num_invperiods), zones, collect(1:num_invperiods))
    reserveup_annual_increment = AxisArrays.AxisArray(ones(length(zones), num_invperiods), zones, collect(1:num_invperiods))
    reservedown_annual_increment = AxisArrays.AxisArray(ones(length(zones), num_invperiods), zones, collect(1:num_invperiods))

    for z in zones
        for p in 1:num_invperiods
            for i in 1:p
                energy_annual_increment[z, p] = energy_annual_increment[z, p] * (1 + load_growth[z, i])
                reserveup_annual_increment[z, p] = energy_annual_increment[z, p] * (1 + load_growth[z, i])
                reservedown_annual_increment[z, p] = energy_annual_increment[z, p] * (1 + load_growth[z, i])
            end
        end
    end

    capacity_mkt_param_file = joinpath(investor_dir, "markets_data", "capacity_mkt_param.csv")
    capacity_annual_increment = load_growth

    REC_mkt_params = read_data(joinpath(investor_dir, "markets_data", "REC_mkt_param.csv"))
    price_cap_rec = REC_mkt_params.price_cap[1]
    rec_req = REC_mkt_params.rec_req[1] * rec_market_bool
    rec_annual_increment = REC_mkt_params.annual_increment[1] * rec_market_bool

    capacity_markets = Vector{CapacityMarket}(undef, num_invperiods)
    energy_markets = Vector{EnergyMarket}(undef, num_invperiods)
    reserve_up_markets = Vector{ReserveUpMarket}(undef, num_invperiods)
    reserve_down_markets = Vector{ReserveDownMarket}(undef, num_invperiods)

    rec_markets = Vector{RECMarket}(undef, num_invperiods)

    average_capacity_growth = Statistics.mean(capacity_annual_increment)

    # Create market products data for the horizon
    for p in 1:num_invperiods
        system_peak_load = (1 + average_capacity_growth) ^ (p) * peak_load
        capacity_markets[p] = create_capacity_demand_curve(capacity_mkt_param_file, system_peak_load, capacity_market_bool)

        energy_markets[p] = EnergyMarket(AxisArrays.AxisArray(zonal_load .* energy_annual_increment[:, p], zones, (1:num_hours)),
                                     price_cap_energy)
        reserve_up_markets[p] = ReserveUpMarket(reserve_up_break_points,
                                                reserve_up_price_points,
                                                AxisArrays.AxisArray(zonal_reserve_up .* reserveup_annual_increment[:, p], zones, (1:num_hours)))
        reserve_down_markets[p] = ReserveDownMarket(AxisArrays.AxisArray(zonal_reserve_down .* reservedown_annual_increment[:, p], zones, (1:num_hours)),
                                                  price_cap_reservedown)
        rec_markets[p] = RECMarket(min(rec_req + rec_annual_increment * (p + iteration_year - 1), 1),
                             price_cap_rec)
    end

    max_peak_loads = AxisArrays.AxisArray([maximum([maximum(market.demand[z, :]) for market in energy_markets]) for z in zones], zones)

    markets = MarketCollection.(capacity_markets,
                                energy_markets,
                                reserve_up_markets,
                                reserve_down_markets,
                                rec_markets)
    #-----------------------------------------------------------------------------------------------------------------------------
    availability_df = read_data(joinpath(investor_dir, "timeseries_data_files", "Availability", "DAY_AHEAD_availability.csv"))

    invested_portfolio = find_active_invested_projects(expected_portfolio)
    option_portfolio = find_option_projects(expected_portfolio)

    cem_projects = MarketProject[]

    for project in invested_portfolio
        tech = get_tech(project)
        zone = get_zone(tech)
        cem_project = create_market_project(project,
                                          price_cap_energy[zone],
                                          maximum(reserve_up_break_points[zone]),
                                          price_cap_reservedown[zone],
                                          max_peak_loads,
                                          iteration_year,
                                          num_hours,
                                          num_invperiods,
                                          availability_df)

        if !isnothing(cem_project)
            push!(cem_projects, cem_project)
        end
    end

    aggregated_options = MarketProject[]
    for option in option_portfolio
        tech = get_tech(option)
        zone = get_zone(tech)
        similar_option = filter(p -> (get_type(tech) == p.tech_type && zone == p.zone), aggregated_options)

        if length(similar_option) < 1
            aggregated_option = create_market_project(option,
                                          price_cap_energy[zone],
                                          maximum(reserve_up_break_points[zone]),
                                          price_cap_reservedown[zone],
                                          max_peak_loads,
                                          iteration_year,
                                          num_hours,
                                          num_invperiods,
                                          availability_df)
            push!(aggregated_options, aggregated_option)
        else
            investment_cost = get_investment_cost(get_finance_data(option))[iteration_year:iteration_year + num_invperiods - 1]
            for (i, cost) in enumerate(similar_option[1].expansion_cost)
                if investment_cost[i] > cost
                    similar_option[1].expansion_cost[i] = investment_cost[i]
                end
            end

            discount_rate = get_discount_rate(get_finance_data(option))
            if discount_rate > similar_option[1].discount_rate
                similar_option[1].discount_rate = discount_rate
            end

            owner = get_ownedby(get_finance_data(option))
            if !in(owner, similar_option[1].ownedby)
               similar_option[1].base_cost_units = similar_option[1].base_cost_units + 1
               push!(similar_option[1].ownedby, owner)
            end
        end
    end

    append!(cem_projects, aggregated_options)

    system = MarketClearingProblem(zones, lines, average_capital_cost_multiplier, markets, cem_projects, rep_hour_weight)

    return system
end


text:
# This is part of Calendars.jl. See the copyright note there.
# ======================= CalendarUtils =====================

# Planetary week
# As defined by Reingold & Dershowitz
const sunday    = 0
const monday    = 1
const tuesday   = 2
const wednesday = 3
const thursday  = 4
const friday    = 5
const saturday  = 6
const weeklen   = 7

# Map weekday specifiers to weekday names.
WeekDays = Dict{Int, String}(
    0 => "Sunday    ",
    1 => "Monday    ",
    2 => "Tuesday   ",
    3 => "Wednesday ",
    4 => "Thursday  ",
    5 => "Friday    ",
    6 => "Saturday  ",
# Well, we keep this to flag potential bugs.
    7 => "Doomsday  "
)

const january   = 1
const february  = 2
const march     = 3
const april     = 4
const may       = 5
const june      = 6
const july      = 7
const august    = 8
const september = 9
const october   = 10
const november  = 11
const december  = 12

const MonthLength = (31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31)

# Calendar specifiers, Int64-enumeration of calendar names
XX = :0  # Unknown
EC = :1  # European Calendar
CE = :2  # Common Era
JD = :3  # Julian (Roman Calendar)
AM = :4  # Anno Mundi
AH = :5  # Anno Hegirae
ID = :6  # ISO Date
EN = :7  # Euro Number
JN = :8  # Julian Number
DN = :9  # FixDay Number

# Map calendar specifiers to calendar names
CalendarNames = Dict{Int64, String}(
    EC => "European ",
    CE => "Common   ",
    JD => "Julian   ",
    AM => "Hebrew   ",
    AH => "Islamic  ",
    ID => "IsoDate  ",
    EN => "EuroNum  ",
    JN => "JulianNum",
    DN => "DayNumber",
    XX => "INVALID  "
)

CalendarName(cal) = CalendarNames[cal]

# Map calendar specifiers to calendar names
SymbolToString = Dict{Int64, String}(
    EC => "EC",
    CE => "CE",
    JD => "JD",
    AM => "AM",
    AH => "AH",
    ID => "ID",
    EN => "EN",
    JN => "JN",
    DN => "DN",
    XX => "00"
)

const DPart = Int64

# Return symbolic name representing a calendar.
# StringToSymbol = CalendarName
function StringToSymbol(calendar::String, warn=false)::DPart
    (calendar == "European"   || calendar == "EC")  && return EC
    (calendar == "Common"     || calendar == "CE")  && return CE
    (calendar == "CurrentEpoch")                    && return CE
    (calendar == "Gregorian"   )                    && return CE
    (calendar == "Julian"     || calendar == "JD")  && return JD
    (calendar == "Roman"      || calendar == "AD")  && return JD
    (calendar == "Hebrew"     || calendar == "AM")  && return AM
    (calendar == "Jewish"      )                    && return AM
    (calendar == "Islamic"    || calendar == "AH")  && return AH
    (calendar == "Hijri"       )                    && return AH
    (calendar == "IsoDate"    || calendar == "ID")  && return ID
    (calendar == "ISODate"     )                    && return ID
    (calendar == "EuroNum"    || calendar == "EN")  && return EN
    (calendar == "JulianNum"  || calendar == "JN")  && return JN
    (calendar == "FixNum"     || calendar == "RD")  && return DN
    (calendar == "DayNum"     || calendar == "DN")  && return DN
    warn && @warn("Unknown calendar: $calendar")
    return XX
end

const InvalidDayNumber = DPart(0)
const InvalidDayOfYear = DPart(0)
const InvalidDuration = DPart(-1)
const InvalidDate = (XX, DPart(0), DPart(0), DPart(0))
const InvalidDateString = "0000-00-00"

"""

A CDate (short for Calender Date) is defined as

```julia
CDate = Tuple{DPart, DPart, DPart, DPart}
```

    A tuple `date` of type `CDate` is unpacked by convention as
    (calendar, year, month, day), where `calendar` is one of

```julia
"European"
"Common"
"Julian"
"Hebrew"
"Islamic"
"IsoDate"
```

    Alternatively the acronyms

```julia
EC, CE, JD, AM, AH, and ID
```

    can be used. "Common" or CE denotes the proleptic Gregorian calendar
    and `DPart` (date part) is a typename for Int64.

```julia
CDateStr(cd::CDate)
```

The function CDateStr converts a date of type CDate to a string representation, where the
numeric part follows the recommendation of ISO 8601 and is prefixed by one of the acronyms
for the calendar names indicated above.

Examples for `CDates` and their string representation given by `CDateStr`:

```julia
("European", 2022,  1,  6)  -> "EC-2022-01-06"
("Common",   2022,  1, 19)  -> "CE-2022-01-19"
("Julian",   2022,  1,  6)  -> "JD-2022-01-06"
("Hebrew",   5782, 11, 17)  -> "AM-5782-11-17"
("Islamic",  1443,  6, 15)  -> "AH-1443-06-15"
("IsoDate",  2022,  3,  3)  -> "ID-2022-03-03"
```

The components of a `CDate` can be accessed by applying the functions
```julia
Calendar, Year, Month, Day and Date
```
to a calendar date.

The function CalendarName(CC) returns the calendar name from a calendar specifier CC.
For example both CalendarName(ID) and CalendarName(6) return the string "IsoDate".
"""
const CDate = NTuple{4, DPart}
const YMD = NTuple{3, DPart}
const DateTable = NTuple{8, CDate}

Calendar(date::CDate) = date[1]
Year(date::CDate) = date[2]
Month(date::CDate) = date[3]
Day(date::CDate) = date[4]
Date(date::CDate) = (date[2], date[3], date[4])

Calendar(d::Tuple{String, DPart, DPart, DPart}) = StringToSymbol(d[1])
Year(d::Tuple{String, DPart, DPart, DPart}) = d[2]
Month(d::Tuple{String, DPart, DPart, DPart}) = d[3]
Day(d::Tuple{String, DPart, DPart, DPart}) = d[4]
Date(d::Tuple{String, DPart, DPart, DPart}) = (d[2], d[3], d[4])

# TODO Write a 'get-' interface for the DateTable, like:
GetFixNum(dt::DateTable) = dt[7][4] - 2

# Return the name of the weekday for the given fix day number.
function WeekDay(dn::DPart)::String
    day = rem(dn, weeklen)
    return WeekDays[day]
end
WeekDay(dt::DateTable) = WeekDay(GetFixNum(dt))

# CDate stringified, no validation
function DateStr(year::DPart, month::DPart, day::DPart)
    y = lpad(year,  4, "0")
    m = lpad(month, 2, "0")
    d = lpad(day,   2, "0")
    return "$y-$m-$d"
end
DateStr(d::CDate) = DateStr(d[2], d[3], d[4])

function CDateStr(cd::CDate)
    cal, year, month, day = cd
    if cal == DN || cal == JN || cal == EN
        s = lpad(day, 7, "0")
        return SymbolToString[cal]*"#$s"
    end
        s = DateStr(year, month, day)
    return SymbolToString[cal]*"-$s"
end

CDateStr(cal::String, d::YMD) =
CDateStr((StringToSymbol(cal), d[1], d[2], d[3]))

CDateStr(cal::DPart, d::YMD) =
CDateStr((cal, d[1], d[2], d[3]))

CDateStr(cal::String, day::DPart) =
CDateStr((StringToSymbol(cal), DPart(0), DPart(0), day))

CDateStr(cal::DPart, day::DPart) =
CDateStr((cal, DPart(0), DPart(0), day))

function CDateStr(cal::DPart, year::DPart, month::DPart)
    y = lpad(year,  4, "0")
    m = lpad(month, 2, "0")
    return SymbolToString[cal] * "-$y-$m"
end

CDateStr(cal::String, year::DPart, month::DPart) =
CDateStr(StringToSymbol(cal), year, month)

CDateStr(cal::String, year::DPart, month::DPart, day::DPart) =
CDateStr((StringToSymbol(cal), year, month, day))

"""

PrintDateTable(D::DateTable, io=stdout)

Print a `DateTable` (which normally is build by `CalendarDates`) to `io`.
"""
function PrintDateTable(D::DateTable, io=stdout)
    for d in D
        if d[1] != XX  # Don't print invalid dates
            cname = CalendarNames[d[1]]
            println(io, cname, " ", CDateStr(d))
        end
    end
end

"""

```julia
PrintDateLine(date::CDate, io=stdout)
```

Given a `date` print a line with all representations of this date to `io`. The parts of the date can also be given separately.
For example:

```julia
julia> PrintDateLine(EC, 2022, 2, 2)
```
| CE-2022-02-02 | EC-2022-02-02 | JD-2022-01-20 | AM-5782-12-01 | AH-1443-06-29 | ID-2022-05-03 | EN#0738190 | JN#2459612 |
"""
function PrintDateLine(date::CDate, io=stdout)
    dn = DayNumberFromDate(date)
    println(io,  "| ",
        CDateStr(DateFromDayNumber(CE, dn)), " | ",
        CDateStr(DateFromDayNumber(EC, dn)), " | ",
        CDateStr(DateFromDayNumber(JD, dn)), " | ",
        CDateStr(DateFromDayNumber(AM, dn)), " | ",
        CDateStr(DateFromDayNumber(AH, dn)), " | ",
        CDateStr(DateFromDayNumber(ID, dn)), " | ",
        CDateStr(DateFromDayNumber(EN, dn)), " | ",
        CDateStr(DateFromDayNumber(JN, dn)), " |"
    )
end

PrintDateLine(cal, year, month, day, io=stdout) = PrintDateLine((cal, year, month, day), io)

header = "|    Common     |   European    |    Julian     |    Hebrew     |    Islamic    |    IsoDate    |  EuroNum   | JulianNum  |"
markup = "|    :---:      |     :---:     |    :---:      |     :---:     |    :---:      |     :---:     |   :---:    |   :---:    |"

"""

```julia
PrintEuropeanMonth(year::DPart, month::DPart, io=stdout)
```

Print a calendar for the given European month as a markdown table to `io`.

For example:
```julia
julia> PrintEuropeanMonth(2022, 2)
```
"""
function PrintEuropeanMonth(year::DPart, month::DPart, io=stdout)
    datestr = CDateStr(EC, year, month)
    println(io, "# Calendars of the month ", datestr)
    println(io, header)
    println(io, markup)
    dn = DayNumberFromDate((JD, year, month, 1))
    if dn < GregorysBreak
        for day in 1:LastDayOfMonthJulian(year, month)
            PrintDateLine((JD, year, month, day), io)
        end
    else
    for day in 1:LastDayOfMonthGregorian(year, month)
            PrintDateLine((CE, year, month, day), io)
        end
    end
end

datadir = realpath(joinpath(dirname(@__FILE__)))

"""

```julia
SaveEuropeanMonth(year::DPart, month::DPart, dirname:String)
```

Save the calendar of the given European month as a markdown table to a file in the directory
`dirname`. If no directory is given the file is written to the execution directory.

For example:
```julia
julia> SaveEuropeanMonth(2022, 2)
```
"""
function SaveEuropeanMonth(year::DPart, month::DPart, dir=datadir)
    datestr = CDateStr(EC, year, month)
    filename = joinpath(dir, "Calendar$datestr.md")
    f = open(filename, "w") do f
        PrintEuropeanMonth(year, month, f)
    end
    @info("Saved European month $datestr to: $filename")
end

function PrintWeekLine(date::CDate, io=stdout)
    dn = DayNumberFromDate(date)
    println(io,  "| ",
        DateStr(DateFromDayNumber(EC, dn)), " | ",
        DateStr(DateFromDayNumber(AM, dn)), " | ",
        DateStr(DateFromDayNumber(AH, dn)), " | "
    )
end

wheader = "|  Weekday  |  European  |   Hebrew   |  Islamic   |"
wmarkup = "|   :---:   |    :---:   |   :---:    |   :---:    |"

"""

```julia
PrintIsoWeek(year::DPart, week::DPart, io=stdout)
```

Print a calendar for the given Iso week as a markdown table to `io`.

For example:
```julia
julia> PrintIsoWeek(2022, 2)
```
"""
function PrintIsoWeek(year::DPart, week::DPart, io=stdout)
    println(io, "# Days of week ", week, " of the year ", year, ".")
    println(io, wheader)
    println(io, wmarkup)
    for day in 1:7
        ind = day != 7 ? day : 0
        print("| ", WeekDays[ind])
        PrintWeekLine((ID, year, week, day), io)
    end
end

NextDay(date::CDate)     = DateFromDayNumber(date[1],  1 + DayNumberFromDate(date))
PreviousDay(date::CDate) = DateFromDayNumber(date[1], -1 + DayNumberFromDate(date))

Warning(d) = string(d) * " is not a valid date!"


text:
module DoubletSource

using LinearAlgebra
using Base.Iterators
using StaticArrays
using Statistics

using ..AeroMDAO: AbstractPanel2D, Panel2D, WakePanel2D, Point2D, collocation_point, p1, p2, transform_panel, affine_2D, panel_length, panel_angle, panel_tangent, panel_normal, panel_dist, rotation, inverse_rotation, midpair_map, pressure_coefficient, wake_panel, wake_panels, panel_points, panel_vector

## Doublet-source Dirichlet boundary condition
#===========================================================================#

source_potential(str, x, z, x1, x2) = str / 4π * ((x - x1) * log((x - x1)^2 + z^2) - (x - x2) * log((x - x2)^2 + z^2) + 2z * (atan(z, x - x2) - atan(z, x - x1)))

source_velocity(str, x, z, x1, x2) = SVector(str / (4π) * log(((x - x1)^2 + z^2) / ((x - x2)^2 + z^2)), doublet_potential(str, x, z, x1, x2))

doublet_potential(str, x, z, x1, x2) = str / (2π) * (atan(z, x - x1) - atan(z, x - x2))

doublet_velocity(str, x, z, x1, x2) = SVector(str / (2π) * - (z / ((x - x1)^2 + z^2) - z / ((x - x2)^2 + z^2) ), str / (2π) * ( (x - x1) / ((x - x1)^2 + z^2) - (x - x2) / ((x - x2)^2 + z^2)))

## Matrix helpers
#===========================================================================#

function doublet_influence(panel_j :: AbstractPanel2D, panel_i :: AbstractPanel2D)
    xp, yp = transform_panel(panel_j, panel_i)
    ifelse(panel_i == panel_j, 0.5, doublet_potential(1., xp, yp, 0., panel_length(panel_j)))
end

function source_influence(panel_j :: AbstractPanel2D, panel_i :: AbstractPanel2D)
    xp, yp = transform_panel(panel_j, panel_i)
    source_potential(1., xp, yp, 0., panel_length(panel_j))
end

boundary_condition(panel_j :: AbstractPanel2D, panel_i :: AbstractPanel2D, u) = -source_influence(panel_j, panel_i) * dot(u, panel_normal(panel_j))

## Aerodynamic coefficients
#===========================================================================#

panel_velocity(dφ, dr, u, α) = dφ / dr + dot(u, α)

lift_coefficient(cp, dist_colpoints, panel_angle) = - cp * dist_colpoints * cos(panel_angle)

lift_coefficient(wake_strength, speed) = 2. * wake_strength / speed

"""
    aerodynamic_coefficients(vels, Δrs, panel_angles, speed, α)

Compute the lift, moment, and pressure coefficients given associated arrays of edge speeds, adjacent collocation point distances, panel angles, the freestream speed, and angle of attack ``\\alpha``.
"""
function eval_coefficients(vels, Δrs, xjs, panel_angles, speed, α)
    cps   = @. pressure_coefficient(speed, vels)
    cls   = @. lift_coefficient(cps, Δrs, panel_angles)
    cms   = @. -cls * xjs * cos(α)

    cls, cms, cps
end

## Matrix assembly
#===========================================================================#

include("matrix_func.jl")

export solve_problem

function solve_problem(panels, u, α, sources :: Bool, wake_length)
    speed 			= norm(u)
    xs	 			= getindex.(panel_points(panels)[2:end-1], 1)

    # Blunt trailing edge tests
    te_panel 		= Panel2D((p2 ∘ last)(panels), (p1 ∘ first)(panels))
    r_te 			= panel_vector(te_panel)
    φ_TE  			= dot(u, r_te)

    # Solve for doublet strengths
    φs 				= solve_strengths(panels, u, α, r_te, sources; bound = wake_length)

    # Evaluate inviscid edge velocities
    u_es, Δrs 		= tangential_velocities(panels, φs, u, sources)

    # Compute coefficients
    cls, cms, cps 	= eval_coefficients(u_es, Δrs, xs, panel_angle.(panels[2:end]), speed, α)

    # Evaluate lift coefficient from wake doublet strength
    cl_wake 		= lift_coefficient(φs[end] - φs[1] + φ_TE, speed)

    cls, cms, cps, cl_wake
end

function solve_problem(panels, u, α, num_wake :: Integer, wake_length)
    speed 			= norm(u)
    xs	 			= getindex.(panel_points(panels)[2:end-1], 1)
    wakes 			= wake_panels(panels, wake_length, num_wake)
    φs				= solve_strengths(panels, u, α, wakes; bound = wake_length)
    u_es, Δrs 		= tangential_velocities(panels, φs[1:end-1], u, false)
    cls, cms, cps 	= evaluate_coefficients(u_es, Δrs, xs, panel_angle.(panels[2:end]), speed, α)
    cl_wake 		= lift_coefficient(last(φs), speed)

    cls, cms, cps, cl_wake
end

end

max_forks_repo_name: GodotMisogi/AeroMDAO