text:
using MacroTools, Reactive

attribute_doc(name, func, indent) = ""
function attribute_doc(name, func::Symbol, indent)
    io = IOBuffer()
    f = getfield(Makie, func)
    println(io, "Attribute `$name`, conversion function [`$func`](@ref)")
    String(take!(io))
end

"""
Creates the expression that fetches the default for an attribute, converts it to a node
and inserts it into the kw_arg dictionary
"""
function convert_expr(var, cfunc, args, mainfunc, fargs, dictsym)
    scene_sym, kw_sym = fargs
    var_sym = QuoteNode(var)
    tmpsym = gensym("tmp")
    if length(args) == 1 && args[1] == var
        quote
            $tmpsym = find_default($scene_sym, $kw_sym, $(QuoteNode(mainfunc)), $var_sym)
            $(esc(var)) = to_node($tmpsym, x-> ($(esc(cfunc))($scene_sym, x)))
            $dictsym[$var_sym] = $(esc(var))
        end
    else # case when a secondary convert function gets called, which builds up on other attributes
        args = esc.(args)
        quote
            $(esc(var)) = $(esc(cfunc))($scene_sym, $(args...))
            $dictsym[$var_sym] = $(esc(var))
        end
    end
end

function process_body_element(elem, fargs, mainfunc, dictsym, kwarg_keys, indent = 1)
    docs = []; symbols = []
    if isa(elem, Expr) && elem.head == :kw
        var, call = elem.args
        @capture(call,
            f_(args__) |
            args__::f_
        ) || error("Use a call or type assert on right hand side")

        expr = convert_expr(var, f, args, mainfunc, fargs, dictsym)
        return expr, [var], [attribute_doc(var, f, indent)]
    end
    if isa(elem, Expr) && elem.head == :block
        syms = Symbol[]
        docs = []
        result = Expr(:block)
        for arg in elem.args
            expr, _syms, _docs = process_body_element(arg, fargs, mainfunc, dictsym, kwarg_keys, indent)
            push!(result.args, expr)
            append!(syms, _syms)
            append!(docs, _docs)
        end
        return result, syms, docs
    end
    # If is documented
    if isa(elem, Expr) && elem.head == :macrocall &&
            length(elem.args) == 3 && elem.args[1].head == :core &&
            elem.args[1].args[1] == Symbol("@doc")
        push!(docs, elem.args[2])
        # the expression that is documented
        return process_body_element(elem.args[3], fargs, mainfunc, dictsym, kwarg_keys, indent)
    end
    # exclusive blocks
    if isa(elem, Expr) && elem.head == :call && elem.args[1] == :xor
        xor_expr = Expr(:block)
        first_expr = Expr(:block)
        current_expr = xor_expr
        docio = IOBuffer()
        println(docio, "## Exclusive Attribute sets:\n---\n")
        for (i, arg) in enumerate(elem.args[2:end])
            expression, docs, condition = if isa(arg, Expr) && arg.head == :if
                condition_syms = arg.args[1]
                if isa(condition_syms, Expr) && condition_syms.head == :tuple
                    map!(QuoteNode, condition_syms.args, condition_syms.args)
                    condition = :(all(x-> x in $kwarg_keys, $condition_syms))
                    expression, syms, docs = process_body_element(arg.args[2], fargs, mainfunc, dictsym, kwarg_keys, indent + 1)
                    expression, docs, condition
                else
                    error("Needs if (sym1, sym2, sym3), found if $(condition_syms)")
                end
            else
                expression, syms, docs = process_body_element(arg, fargs, mainfunc, dictsym, kwarg_keys, indent + 1)
                condition = :(any(x-> x in $kwarg_keys, $syms))
                expression, docs, condition
            end
            println(docio, join(docs, "\n"))
            if i == length(elem.args[2:end])
                println(docio, "## end\n---\n")
            else
                println(docio, "## or\n---\n")
            end
            if i == 1
                first_expr = expression
                # since a xor block defaults to first definition, we don't actually
                # need to check for it's args and create an if - since it will land
                # in the last else block anyways
                continue
            end
            else_expr = Expr(:block)
            ifelse = Expr(:if, condition, expression, else_expr)
            push!(current_expr.args, ifelse)
            current_expr = else_expr # now we need to insert into else

        end
        # defaults to first block in the last else block
        current_expr.args = first_expr.args
        return xor_expr, Symbol[], [String(take!(docio))]
    end
    found = @capture(elem,
        (var_ = f_(args__)) |
        (var_ = args__::f_)
    )
    result = if found
        push!(symbols, var)
        push!(docs, attribute_doc(var, f, indent))
        convert_expr(var, f, args, mainfunc, fargs, dictsym)
    else
        elem
    end
    return result, symbols, docs
end

"""
    `@default function name(args...) end`

Macro that allows for concise default creations.
From this expression:

    ```julia
        @default function sprites(backend, scene, kw_args)
            attribute = convert_function(attribute)
            attribute2 = convert_function2(attribute, attribute2)
        end
    ```

It creates a function `sprites_default(scene, kw_args)::Dict{Symbol, Any}`

Which will first look in kw_args for `:attribute`, if found it will call `convert_function(kw_args[:attribute])`
and insert it in the returned attribute dictionary.
If it's not found in kw_args, it will search in scene.theme.sprites for `:attribute` and if not found there it will
search one level higher (scene.theme).
This can be manually achieved by calling: `find_default(scene, kw_args, func, :attribute)`.

The same will be done for attribute2, which also demonstrate that you can reference previously defined attributes and that you can
use as many inputs to convert_function as you want.
You can optionally define a doc string for an attribute like this:
    ```julia
        Attribute 1 is great
        attribute = convert_fun(attribute)
    ```
If you don't define any doc string, it will default to the doc string of the convert function.

For attributes that don't need a complex convert function, you can simply use a type
assert:
    ```Julia
        attribute = to_float(attribute)
        # will become
        attribute = convert(Float32, find_default(scene, kw_args, sprites, :attribute))
    ```

# xor blocks - exlusive sets of attributes
You can define blocks of exclusive attributes like so:
    ```Julia
    xor(
        begin
            color = to_color(color)
        end,
        begin
            colormap = to_colormap(colormap)
            intensity = to_intensity(intensity)
            colornorm = to_colornorm(colornorm, intensity)
        end
    )
```
This gets desugared to:

    ```julia
    if any(x-> x in keys(kw_args), (:color, :intensity, :colornorm))
        colormap = to_colormap(colormap)
        intensity = to_intensity(intensity)
        colornorm = to_colornorm(colornorm, intensity)
    else
        color = to_color(color)
    end
    ```
    So the first block becomes the default block if no key is in the kw_args.
    You can also explicitely define what keys kw_args needs to contain to get selected:
    ```julia
    xor(
        begin
            color = to_color(color)
        end,
        if (colormap, intensity) # now kw_args needs to contain attributes colormap && intensity
            colormap = to_colormap(colormap)
            intensity = to_intensity(intensity)
            colornorm = to_colornorm(colornorm, intensity)
        end
    )
    ```
"""
macro default(func)
    @capture(func,
        function mainfunc_(fargs__)
            body__
        end
    ) || error("Please use @default with a function declaration, e.g.
        @default function myfunc(args...)
            ...
        end
    ")

    length(fargs) == 2 || error("Function should have 2 arguments, namely scene and kw_args. Found: $fargs")
    docs = []
    result = []
    dictsym = gensym(:attributes)
    kwarg_keys = gensym(:keys)
    for elem in body
        expr, syms, _docs = process_body_element(elem, fargs, mainfunc, dictsym, kwarg_keys)
        push!(result, expr)
        append!(docs, _docs)
    end
    io = IOBuffer()
    for elem in docs
        # the `    ` before Attribute get it parsed as a code block, which destroys cross references
        println(io, replace(elem, "    Attribute", "Attribute"))
        println(io, "\n____________________\n")
    end
    docstr = String(take!(io))
    expr = quote
        """
        $($(docstr))
        """
        function $(esc(Symbol("$(mainfunc)_defaults")))($(fargs...))
            $dictsym = Dict{Symbol, Any}()
            $kwarg_keys = keys($(fargs[2]))
            $(result...)
            merge!($(fargs[2]), $dictsym)
            return $(fargs[2])
        end
    end
    expr
end

nice_dump(x, intent = 0) = (print("    "^intent); show(x); println())
function nice_dump(x::Expr, intent = 0)
    Base.is_linenumber(x) && return
    println("    "^intent, "head: ", x.head)
    intent += 1
    for elem in x.args
        nice_dump(elem, intent)
    end
end

"""
Billboard attribute to always have a primitive face the camera.
Can be used for rotation.
"""
immutable Billboard end


text:
using Pkg

# Activate environment at current directory
Pkg.activate(dirname(@__FILE__))

# Add Spine Julia registry
pkg"registry add https://github.com/Spine-project/SpineJuliaRegistry.git"

# Download and install all required packages
Pkg.instantiate()

# Set Python executable to current and re-build PyCall if necessary
ENV["PYTHON"] = Sys.which("python")
try
    Pkg.build("PyCall")
catch e
    if !isa(e, Pkg.Types.PkgError)
        throw(e)
    end
end

# Precompile everything
Pkg.API.precompile()

ext: jl


text:
# Tests of the AbstractModel objects in FiniteVolumes.jl

using Test
using StaticArrays
using FiniteVolumes
using ForwardDiff, LinearAlgebra

@testset "Fluxes" begin

@testset "Anonynous models" begin
    @testset "(1 var, 1D)" begin
        flux(α, n) = 0.5*α.^2 * n
        @test flux(3.0, 1.0) == 4.5
        @test FiniteVolumes.jacobian(flux, 2.0, 1.0) == 2.0
    end

    @testset "(1 var, 2D)" begin
        flux(α, n) = n' * [α, 1.0] * α
        @test flux(1.0, SVector(1.0, 0.0)) == 1.0
        @test flux(1.0, SVector(0.0, 1.0)) == 1.0
        @test FiniteVolumes.jacobian(flux, 1.0, SVector(1.0, 0.0)) == 2.0
        @test FiniteVolumes.jacobian(flux, 1.0, SVector(0.0, 1.0)) == 1.0
    end

    @testset "(2 var, 1D)" begin
        flux(u, n) = [0.5*u[1].^2 + u[2], u[2]] .* n
        @test flux([1.0, 2.0], 1.0) == [2.5, 2.0]
        @test flux(SVector(1.0, 2.0), 1.0) == SVector(2.5, 2.0)
        @test FiniteVolumes.jacobian(flux, [1.0, 1.0], 1.0) == [1.0 1.0; 0.0 1.0]
    end

    @testset "(2 var, 2D)" begin
        flux(u, n) = n' * [[u[1], u[2]], [0.0, 0.0]]
        @test flux([1.0, 2.0], [1.0, 0.0]) == [1.0, 2.0]
        @test flux([1.0, 2.0], [0.0, 1.0]) == [0.0, 0.0]
        @test FiniteVolumes.jacobian(flux, [1.0, 2.0], [1.0, 0.0]) == [1.0 0.0; 0.0 1.0]
        @test FiniteVolumes.jacobian(flux, [1.0, 2.0], [0.0, 1.0]) == [0.0 0.0; 0.0 0.0]
    end
end

@testset "LinearAdvection" begin
    mesh = CartesianMesh(10)
    w0 = [sin(2π*x) for x in cell_centers(mesh)]
    i_face = (FiniteVolumes.CartesianMeshes.Half(11),)

    flux = LinearAdvectionFlux(1.0)
    @test (Upwind())(flux, mesh, w0, i_face) == w0[5]
    # @btime (Upwind())($f, $mesh, $w, $i_face)  # ~4ns

    flux = LinearAdvectionFlux(-1.0)
    @test (Upwind())(flux, mesh, w0, i_face) == -w0[6]
end

@testset "Shallow water" begin
    @testset "1D" begin
        f = ShallowWater()
        @test f(SVector(1.0, 0.0), 1.0) == SVector(0.0, 9.81/2)
        @test f(SVector(1.0, 1.0), 1.0) == SVector(1.0, 1 + 9.81/2)
        @test f(SVector(1.0, 1.0), -1.0) == -SVector(-1.0, 1 + 9.81/2)

        for h in [1.0, 2.0], u in [-1.0, 0.0, 1.0]
            @test all(FiniteVolumes.eigvals(f, SVector(h, h*u), 1.0) .≈ SVector(u-sqrt(h*f.g), u+sqrt(h*f.g)))
        end

        f2(v, n) = n * SVector(v[2], v[2]^2/v[1] + v[1]^2*9.81/2)
        mesh = CartesianMesh(2)
        w = [SVector(1.0, 0.0), SVector(2.0, 0.0)]
        i_face = (FiniteVolumes.CartesianMeshes.Half(3),)
        @test (Upwind())(f, mesh, w, i_face) ≈ (Upwind())(f2, mesh, w, i_face)
    end

    @testset "2D" begin
        flux = ShallowWater(9.8)

        for v in [SVector(1.0, 0.0, 0.0), SVector(1.0, 1.0, 1.0)], n in [SVector(1.0, 0.0), SVector(0.0, 1.0), sqrt(2)/2*SVector(1.0, 1.0)]
            J1 = ForwardDiff.jacobian(v -> flux(v, n), v)
            J2 = FiniteVolumes.jacobian(flux, v, n)
            @test J1 == J2
            λ1 = LinearAlgebra.eigvals(Array(J1))
            λ2 = FiniteVolumes.eigvals(flux, v, n)
            @test λ1 ≈ λ2
            eg1 = LinearAlgebra.eigen(Array(J1)).vectors
            eg2 = FiniteVolumes.eigen(flux, v, n)[2]
            for (c1, c2) in zip(eachcol(eg1), eachcol(eg2))
                @test maximum(abs.(c2)) .* abs.(c1) ≈ maximum(abs.(c1)) .* abs.(c2)
            end
        end
    end
end

end


text:
### A Pluto.jl notebook ###
# v0.12.20

using Markdown
using InteractiveUtils

# ╔═╡ 372d3cf2-6173-11eb-356e-23c959c3fd89
# edit the code below to set your name and UGent username

student = (name = "Sam Janssen", email = "Sam.Janssen@UGent.be");

# press the ▶ button in the bottom right of this cell to run your edits
# or use Shift+Enter

# you might need to wait until all other cells in this notebook have completed running.
# scroll down the page to see what's up

# ╔═╡ eb0428ac-5d8c-11eb-09a3-2b3cfc77f3f4
begin
	using DSJulia;
	tracker = ProgressTracker(student.name, student.email);
	md"""
	Submission by: **_$(student.name)_**
	"""
end

# ╔═╡ 70be3952-5d8c-11eb-1509-b3f7077d57e0
md"""
# Composite and parametric types

The previous notebook showed that it is easy to extend the type system with our own types in the hierarchy. The interesting thing is to have the concrete 'leafs' on this tree contain data that can be used in the functions.

## Composite types

Composite types, sometimes call records, structs (matlab) or objects (python), can store several values in its *fields*.

When defining a new composite type, we can choose them to be mutable or immutable:
- mutable types are defined using `mutable struct ... end`, they allow the fields to be changed after the object is created;
- immutable types are defined similarly using `struct ... end`, after creating the object its fields cannot be changed (immutable by default).
Mutable types are a bit more flexible, though might be less safe and are more difficult to work with. As the compiler
knows everything in advance, it often optimizes better for immutable types. The choice depends on the application, though
generally immutable types are the better choice!

As an example, let us define an agent type for an ecological individual-based model (IBM). We create the abstract type `Agent` for which we
can then specify several children types.
"""

# ╔═╡ acd7de0c-5d8c-11eb-120a-8b79f2b8eb3b
abstract type Agent end

# ╔═╡ af8c6460-5d8c-11eb-3ba8-c16e8855e992
md"The concrete types in such an IBM might represent an animal type you want to model, for example preys and predators. Making a concrete type of
a prey animal, we want each to have an unique identifier (represented by an integer) and a position. As we expect the agent to move, hence changing
its position when our simulation runs, we choose a mutable type."

# ╔═╡ b96a33b8-5d8c-11eb-01de-439f53cdc355
mutable struct Prey <: Agent
    id::Int
    pos
end

# ╔═╡ bcacf89e-5d8c-11eb-0077-e5761b8855a3
md"Notice the type annotation for `id`, which we choose to always represent by an integer."

# ╔═╡ c7a077ba-5d8c-11eb-08bb-07c343ea8ab1
md"After defining a composite type, a constructor for that type is immediately available. A constructor is a function to instantiate a type."

# ╔═╡ cd0d8636-5d8c-11eb-19f9-4da4550d306f
deer = Prey(1, (0.5, 1.9))

# ╔═╡ ce25e25c-5d8c-11eb-2e8e-b5b1e7350d70
md"You can always check which field names are available."

# ╔═╡ d672cc72-5d8c-11eb-2c06-0341181e3a3d
fieldnames(Prey)

# ╔═╡ dac11770-5d8c-11eb-1058-2d043e172931
md"The fields can be accessed easily:"

# ╔═╡ 61c43794-6174-11eb-1545-2db114b929e4
deer.pos

# ╔═╡ e7b9023a-5d8c-11eb-1387-cfa7c41ab6ca
fyi(md"This is just syntactic sugar for `getfield`, e.g. `getfield(deer, :id)`")

# ╔═╡ fb6e62d4-5d8c-11eb-34f3-bf3df7cd4cb3
md"Similarly, a predator type can be defined. In addition to an `id` and `position`, which each agent has, they also have a size, determining its mobility."

# ╔═╡ 01fe6f9a-5d8d-11eb-0519-03aefcd587bb
mutable struct Predator <: Agent
    id::Int
    pos
    size::Float64
end

# ╔═╡ 05ae1a2a-5d8d-11eb-0ed8-496f48194232
wolf = Predator(2, (0.0, 0.0), 40.0)  # 40 kg wolf

# ╔═╡ 18d04a36-5d8d-11eb-1986-693eaad5d5be
md"Using the `.` syntax for accessing the fields is very \"object-oriented\" and not  very juliaesque! We should define custom getter functions
for the user to access the relevant fields. We could define `id` and `position` methods to get the respective
fields for the two agents. However, since these fields should be defined for every `Agent` type, we can just create
these for the Agent type!"

# ╔═╡ 1df0419c-5d8d-11eb-3444-d5741bf10d32
id(agent::Agent) = agent.id

# ╔═╡ 234dd2e4-5d8d-11eb-366c-c15474da799f
position(agent::Agent) = agent.pos

# ╔═╡ 23f7e8c4-5d8d-11eb-04af-b5e268c13a14
md"Here, we could theoretically have ommited the type annotation in the function. Then the function would accept objects of the non-agent type and likely yield an error because they don't have the `id` or `pos` field. Now, these functions will return a `MethodError` when given a non-`Agent` input."

# ╔═╡ 330d2090-5d8d-11eb-1003-a52a078514b2
md"A slightly more interesting example is by extending `size`."

# ╔═╡ 375bf8d8-5d8d-11eb-1a49-69905d38effe
Base.size(agent::Predator) = agent.size

# ╔═╡ 39d6d5c4-5d8d-11eb-0e07-11d891ff87a3
size(wolf)

# ╔═╡ 3f5740f6-5d8d-11eb-19c3-ddedf6003e53
md"Here, we had to import `size` because we are extending a function from the `Base` library to work with a new type (doing something vastly different than its original function)."

# ╔═╡ 47049bf0-5d8d-11eb-18c2-733287b69420
md"Similarly, we can program the interaction behaviour between the agents."

# ╔═╡ 50f17586-5d8d-11eb-0eec-579467b787d0
begin
	interact(agent1::Agent, agent2::Agent) = nothing
	interact(agent1::Predator, agent2::Prey) = "eat"
	interact(agent1::Prey, agent2::Predator) = "run"
end

# ╔═╡ 585995ae-5d8d-11eb-256f-bd8e9eb52063
md"We have chosen the default behaviour that two Agents of unspecified types do not interact at all, this will now be the case when a prey meets other prey, a predator an other predator or a new third type comes into the equation."

# ╔═╡ 62e49c94-5d8d-11eb-39ac-f30febf282ff
fyi(md"Since in these simple examples, the `interact` methods do not use their arguments, merely perform type checking, we could have written this as `interact(::Agent,::Agent) = ...` etc.")

# ╔═╡ e5c61ace-654a-11eb-05be-79153b6222b6
QuestionBlock(
	title=md"**Assignment: the wildlife**",
	description = md"""
	Test a few random wild life encounters.
	"""
)

# ╔═╡ 5c141870-654b-11eb-26f3-4b58b0b2e25d
mutable struct Neutral <: Agent end

# ╔═╡ 4d2d4e6c-654b-11eb-02c5-698af6d43b37
fly = Neutral()

# ╔═╡ 6f69292e-654b-11eb-012e-317c7b2f33ea
wolf

# ╔═╡ 78f6ea1c-654b-11eb-217e-6da9345f2fe3
deer

# ╔═╡ 7ea3de2a-654b-11eb-18c2-0fa369696cb9
fauna = [fly, wolf, deer]

# ╔═╡ 3422d2fc-654b-11eb-37c8-93a84229ab21
rand(fauna)

# ╔═╡ 340538c8-654b-11eb-3733-272fb9a8f106
#... interact

# ╔═╡ b1525c16-654b-11eb-0428-43008354b402

# ╔═╡ 6e8548b4-5d8d-11eb-3fcc-45cb005e5c5e
md"""
## Parametric types

Sometimes we want more flexiblility in defining types. Think of designing a new type of matrix. Here you would like to work them for all
numeric datatypes, `Int`, `Int8`, `Float6`, `Rational`, in addition to new datatypes that might not even be defined yet! To this end, we use
*parametric types*, types that **depend** on another type.

For example, consider a 2-dimensional coordinate:
(x,y)
"""

# ╔═╡ cf6dea8c-5d8d-11eb-3f54-4d947305f5e5
begin
	struct Point{T}
		x::T
		y::T
	end

# PASTE YOUR CONSTRUCTORS HERE!

end

# ╔═╡ d305007e-5d8d-11eb-2505-3347d7d1a561
md"Here, each coordinate of the type `Point` has two attributes, `x` and `y`, of the same type. The specific type of `Point` can vary."

# ╔═╡ d84355d4-5d8d-11eb-2c2c-21daf0364c21
p = Point(1.0, 2.0)

# ╔═╡ db651312-5d8d-11eb-09f1-bb693144d8f7
p_int = Point(1, 2)

# ╔═╡ 8b7688ce-5d8e-11eb-0100-dd2712247ff6
p_string = Point("hello", "goodbye")  # might not make sense but currently allowed!

# ╔═╡ 9c0f7916-5d8e-11eb-15cd-11badf9f9c01
Point([1 2; 3 4], [2 3; 9 4])  # maybe relevant to quantum physics?

# ╔═╡ dd5706ec-5d8d-11eb-29df-75170b547cf8
md"Note that"

# ╔═╡ e799c73e-5d8d-11eb-3df4-13dbef13a342
p isa Point

# ╔═╡ eb563ede-5d8d-11eb-1e91-55bd3b36f4eb
p isa Point{Float64}  # more specific

# ╔═╡ fc2d5c92-5d8d-11eb-08c3-29d0478f1884
p isa Point{Int}  # obviously not true

# ╔═╡ 072f2382-5d8e-11eb-377a-d39d0821559e
p isa Point{Real}  # unexpectedly not true!

# ╔═╡ 17022b06-5d8e-11eb-391b-7194962a2c18
md"The observations above is very important! Even though `Float16 <: Real`, this does not hold for the corresponding parametric types."

# ╔═╡ 393e2f9e-5d8e-11eb-11d7-511c32ce6e48
Point(1, 2.0)  # should error initially, but will be fixed later in the notebook!

# ╔═╡ 4b262a72-5d8e-11eb-02bb-5fb2bd07f280
md"Parametric types can be used in dispatch. For example, if we want to compute the norm of a Point, this would only make sense if Point is a number."

# ╔═╡ 605779e6-5d8e-11eb-3e08-c7420ef76aba
norm(p::Point{T} where {T<:Number}) = sqrt(p.x^2 + p.y^2)

# ╔═╡ 6646eb5c-5d8e-11eb-1e01-f3011c4230de
norm(p)

# ╔═╡ 6960f8d2-5d8e-11eb-0215-2de7b54e3081
norm(p_int)  # dispatch creates a method for this type

# ╔═╡ 8321136a-5d8e-11eb-0da2-e750b81d7ce9
md"""
## Constructors

Constructors are functions that create new objects.

### Outer constructors

We have already seen that when creating a new `struct`, this immediately initiates the constructor (e.g., `Point(1.0, 2.0)`). These can also be made explicitly:

```julia
Point(x::T, y::T) where {T<:Real} = Point{T}(x,y)
```

In a normal file, we could run this, but Pluto does not allow us to have constructors at different places. Copy-paste it to the cell where we defined the `Point` type.

"""

# ╔═╡ e246757e-5d8e-11eb-313c-61a63246cf9a
Point(1, 2.0)

# ╔═╡ ecc14f96-5d8f-11eb-2b18-cb23fad26b6e
Point(1.2)

# ╔═╡ 9d4b9eee-5d8f-11eb-07bb-557415e4ac4a
md"""

The line above yields an error, because there is no constructor method when two inputs are of a different type. The following constructors will resolve such cases satisfactorily.

```julia
Point(x::Real, y::Real) = Point(promote(x, y)...)
```

Add it to the definition and see the change!

We can write other constructors just like functions. For example, support that when we provide a single `x`, we want to create a point (x, y):

```julia
Point(x) = Point(x, x)
```
"""

# ╔═╡ 142666de-5d90-11eb-3231-efe53ddc9b0d
md"""

### Inner constructors

The above examples show *outer constructors*. These are defined outside the structure. We can also use *inner constructors*, which are declared within the definition of the type. These make use of the keyword `new`. For example, let us define an ordered pair.
"""

# ╔═╡ 27fcaede-5d90-11eb-1cea-91fcc4b6b0fe
struct OrderedPair
  x
  y
  function OrderedPair(x, y)
    if x < y
      new(x, y)
    else
      new(y, x)
    end
  end
end

# ╔═╡ 2a224fde-5d90-11eb-1c46-3fd248350914
OrderedPair(18, 23)

# ╔═╡ 2c937298-5d90-11eb-06e6-ab70b9d3701e
OrderedPair(8, 2)

# ╔═╡ 2efe030e-5d90-11eb-38ae-092222d3a8d4
fyi(md"For parametric types, the `new` keyword should be type annotated. So, for in the `Point` example one would use `new{T}(x,y)`.")

# ╔═╡ 46df4eb0-5d90-11eb-1fdf-f34a7bcb7191
md"""

## Example: iterators

We can extend Julia by making use of establised interfaces, such as for iterators like `1:0.1:10`. For example, suppose we want to iterate over the first $n$ squares of natural numbers.
"""

# ╔═╡ a4a1cb18-5d90-11eb-08ee-8570368a056b
struct Squares
    count::Int
end

# ╔═╡ aca1c930-5d90-11eb-29d9-954e097bbe3b
md"To make this an iterator, we just have to extend the `iterate` function of `Base`. This is what is needed for Julia to treat this as an iterator because all functions fall back to this."

# ╔═╡ a9502b64-5d90-11eb-144c-3d7ce0949e67
Base.iterate(S::Squares, state=1) = state > S.count ? nothing : (state*state, state+1)

# ╔═╡ da6cc5b8-5d90-11eb-07ff-db8bdb504054
@terminal for i in Squares(7)
    println(i)
end

# ╔═╡ e3759d4c-5d90-11eb-0bea-bb4247623ec2
25 ∈ Squares(10)

# ╔═╡ 07998440-5d91-11eb-1a65-8de428eac89c
sum(Squares(18093))

# ╔═╡ e11b0b10-6621-11eb-0bdb-f3719cc92a20
@elapsed sum(Squares(18093))

# ╔═╡ 0b8980aa-5d91-11eb-385b-71568ec0e325
md"Let's give the compiler some additional information!"

# ╔═╡ 192d9fd4-5d91-11eb-1cb9-c706aad03480
Base.eltype(::Type{Squares}) = Int

# ╔═╡ 1fa68c72-5d91-11eb-1102-c18460e92ee6
Base.length(S::Squares) = S.count

# ╔═╡ 26a70358-5d91-11eb-0241-699b5dc4783f
md"Now this works:"

# ╔═╡ 2270e790-5d91-11eb-20e5-29905f232734
collect(Squares(4))

# ╔═╡ 2cbe2370-5d91-11eb-130a-5da52d1a62c8
md"I remember there is a fancy formula to compute the sum of squared natural numbers. Uncomment it to activate!"

# ╔═╡ 49f1d98c-5d91-11eb-1657-f320e9fcdc0e
#Base.sum(S::Squares) = (n = S.count; return n*(n+1)*(2n+1)÷6)

# ╔═╡ 4cb68744-5d91-11eb-2b3e-e7df55888c93
sum(Squares(18093))  # much faster now!

# ╔═╡ e99af5c0-6621-11eb-058b-45c3719930d0
@elapsed sum(Squares(18093))

# ╔═╡ 579e3828-5d91-11eb-1d33-b94628d61fc0
md"""

## Example: custom matrices

Similarly, we can make our very own matrix types. Consider the Strang matrix, a [tridiagonal matrix](https://en.wikipedia.org/wiki/Tridiagonal_matrix) with 2 on the main diagonal, and -1 on the first diagonal below and above the main diagonal.
"""

# ╔═╡ e9a99a00-5d91-11eb-2c50-8be452cab83f
struct Strang <: AbstractMatrix{Int}
    n::Int
end

# ╔═╡ ec62c35c-5d91-11eb-3773-b9385f312f7f
Base.size(S::Strang) = (S.n, S.n)

# ╔═╡ efb0b460-5d91-11eb-2534-496df689dc60
Base.getindex(S::Strang, i, j) = i==j ? 2 : (abs(i - j) == 1 ?  -1 : 0)

# ╔═╡ f3c3114c-5d91-11eb-1d37-6d97ea6d267f
S = Strang(1000)  # holy cow! Looks just like a real matrix!

# ╔═╡ 04dcda58-5d92-11eb-10ba-396947081338
sum(S)  # works, but slow...

# ╔═╡ fbdb2958-6621-11eb-3cb6-a9bdeea3bdb7
@elapsed sum(S)

# ╔═╡ 066901bc-6622-11eb-23ad-ed869adad36d
md"Uncomment the line below."

# ╔═╡ 0f878dea-5d92-11eb-0000-b7484532ee70
#Base.sum(S::Strang) = 2

# ╔═╡ 11630c02-5d92-11eb-1746-4dabf327fbbe
sum(S)

# ╔═╡ 046ce4f8-6622-11eb-3c4f-7b6bf21fb77b
@elapsed sum(S)

# ╔═╡ 1e65cb9c-5d92-11eb-3526-332169917fd9
v = randn(1000)

# ╔═╡ 201f59ee-5d92-11eb-33ae-51904d249dd4
S * v  # works, but slow

# ╔═╡ 310d8966-6176-11eb-1b23-9d942f197fbd
md"Luckily there is a trick to computing the product of a tridiagonal matrix and a vector:"

# ╔═╡ 276e9af4-5d92-11eb-1399-993570859698
function Base.:*(S::Strang, v::Vector)
    n = length(v)
    @assert size(S, 2) == n
    x = similar(v)
    for i in 1:n
        x[i] = v[i]
        i > 1 && (x[i] += v[i-1])
        i < n && (x[i] += v[i+1])
    end
    return x
end

# ╔═╡ 300a8428-5d92-11eb-188b-05d00df4f6a7
S * v  # fast (linear time in v)

# ╔═╡ cf146bf4-6177-11eb-1eaa-c35efde57b3e
md"## Exercises"

# ╔═╡ 3ae60e88-5d94-11eb-0c50-1d74ea104758
struct WizCur
	# complete me!
end

# ╔═╡ 48301af2-5d94-11eb-0019-7737667c9cea
galleons(money::WizCur) = missing

# ╔═╡ 4ea80eda-5d94-11eb-3882-21a41d2d65f8
sickles(money::WizCur) = missing

# ╔═╡ 5af60d90-5d94-11eb-2ee4-b7bfc2caf53b
knuts(money::WizCur) = missing

# ╔═╡ 5f7c75ac-5d94-11eb-137a-7914cd009821
#=
function Base.show(io::IO, money::WizCur)
    print(io, "I am printed for WIZCUR, make me say something meaningful")
end
=#

# ╔═╡ 678ca64a-5d94-11eb-2b85-0b706526e35b
Base.isless(m1::WizCur, m2::WizCur) = missing

# ╔═╡ 8afbd434-5d94-11eb-366b-d3d719189ef7
Base.isgreater(m1::WizCur, m2::WizCur) = missing

# ╔═╡ 9226b2d8-5d94-11eb-2bea-491eb7dc1da7
Base.isequal(m1::WizCur, m2::WizCur) = missing

# ╔═╡ 95146d46-5d94-11eb-22aa-c1a544e0d784
Base.:+(m1::WizCur, m2::WizCur) = missing

# ╔═╡ 9eab40be-5d94-11eb-0c59-21f5824fb812
money_ron = missing

# ╔═╡ a137e0f8-5d94-11eb-2209-73acad549307
money_harry = missing

# ╔═╡ a79ba114-5d94-11eb-16ae-9906c6cdf54f
dungbomb_fund = money_ron + money_harry

# ╔═╡ dd0b7c7a-6177-11eb-2cb7-8b194a75d776
begin

	q_wc = Question(
		validators = @safe[
			Solutions.WizCur(732, 19, 0) |> Solutions.galleons ==
				WizCur(732, 19, 0) |> galleons,
			Solutions.WizCur(732, 19, 0) |> Solutions.knuts ==
				WizCur(732, 19, 0) |> knuts,
			Solutions.WizCur(732, 19, 0) |> Solutions.sickles ==
					WizCur(732, 19, 0) |> sickles,
			galleons(dungbomb_fund) == Solutions.galleons(Solutions.dungbomb_fund),
			sickles(dungbomb_fund) == Solutions.sickles(Solutions.dungbomb_fund),
			knuts(dungbomb_fund) == Solutions.knuts(Solutions.dungbomb_fund)
		]
	)

	qb_wc = QuestionBlock(
		title=md"**Exercise: wizarding currency**",
		description = md"""

The British Wizarding World uses Galleons, Sickles, and Knuts as a currency. There are 17 Sickles in a Galleon, and 29 Knuts in a Sickle, meaning there are 493 Knuts to a Galleon. We will make a structure `WizCur` to represent wizarding currency. This structure has three integer-valued fields: `galleons`, `sickles`, and `knuts`. The constructor should always create tidy representations, meaning that, for example, if the number of knuts is 29 or more, it just adds an appropriate number of sickles such that the number knuts is less than 29 (it's magical money). The same applies to the sickles, which can also never exceed 17.

Overload `Base.show` such that Julia prints your currency as, for example, `7G, 2S, 9K`.

Also, overload the function `+` to add two instances of `WizCur` and the `>` and `<` operators to compare two instances of wizarding currency.

The piggy bank with Ron's life savings contains 19 Sickles and 732 Knuts. Harry has 3 Galleons, 1 Sickle, and 7 Knuts pocket change. Who has the most money? How many do they have together?
		""",
		hints = [
			hint(md"You might find `%` and `div` useful here.")
		],
		questions = [q_wc]
	)
		validate(qb_wc, tracker)

end

# ╔═╡ 392228e2-617d-11eb-09a5-c9e5649356eb

# ╔═╡ d448a2e0-5d92-11eb-18a6-9ff817992154
begin
	struct Vandermonde{T,VT} <: AbstractMatrix{T}
		α::VT
		m::Int
		Vandermonde(α::AbstractVector{T}, m) where {T} = missing
	end

	# take length of α as a default value of m
	Vandermonde(α::Vector{<:Number}) = missing
end

# ╔═╡ 902f4dfe-617d-11eb-2957-71130adca3ae
begin
	q_vm = Question(
		validators = @safe[

			Solutions.Vandermonde(Solutions.α, length(Solutions.α)) ==
				Vandermonde(Solutions.α, length(Solutions.α)),

			Solutions.Vandermonde(Solutions.α) == Vandermonde(Solutions.α),

			Solutions.Vandermonde(Solutions.α)[1] == Vandermonde(Solutions.α)[1],

			size(Solutions.Vandermonde(Solutions.α)) ==
				size(Vandermonde(Solutions.α))
		]
	)

	qb_vm = QuestionBlock(
		title=md"**Exercise: Vandermonde matrix**",
		description = md"""
		The [Vandermonde matrix](https://en.wikipedia.org/wiki/Vandermonde_matrix) can be obtained from a vector by taking the powers from 0 till $m-1$.

		$${\displaystyle V={\begin{bmatrix}1&\alpha _{1}&\alpha _{1}^{2}&\dots &\alpha _{1}^{n-1}\\1&\alpha _{2}&\alpha _{2}^{2}&\dots &\alpha _{2}^{n-1}\\1&\alpha _{3}&\alpha _{3}^{2}&\dots &\alpha _{3}^{n-1}\\\vdots &\vdots &\vdots &\ddots &\vdots \\1&\alpha _{m}&\alpha _{m}^{2}&\dots &\alpha _{m}^{n-1}\end{bmatrix}},}$$

		$$V_{i,j} = \alpha_i^{j-1}$$

		Complete the implementation to store and process this matrix.
		""",
		questions = [q_vm]
	)
		validate(qb_vm, tracker)

end

# ╔═╡ bd91a60e-5d93-11eb-09d4-830ca69439bf
Base.size(V::Vandermonde) = missing

# ╔═╡ c2ecfec8-5d93-11eb-2640-07bc07f3da98
Base.getindex(V::Vandermonde, i, j) = missing

# ╔═╡ cb3e91cc-5d93-11eb-020c-d73c10131755
α = [1, 2, 3, 4]

# ╔═╡ d107c75e-5d93-11eb-0e6f-097b1291e460
V = Vandermonde(α, 4)

# ╔═╡ 7f02b0a0-617f-11eb-1263-91423840def3

# ╔═╡ f01448f0-617d-11eb-1829-0fcfe19b3115
begin
	qb_dvm = QuestionBlock(
		title=md"**Exercise: determinant of a Vandermonde matrix**",
		description = md"""
		The determinant of a Vandermonde matrix is easy to compute:

		$${\displaystyle \det(V)=\prod _{1\leq i<j\leq n}(x_{j}-x_{i}).}$$

		Overload this for the Vandermonde matrix!

		""",
		hints = [
			hint(md"`prod`")
		]
	)
end

# ╔═╡ d2a076ea-5d93-11eb-216e-f5c37d330b40
import LinearAlgebra

# ╔═╡ dc945902-5d93-11eb-1121-a7ae99c5862e
#LinearAlgebra.det(V::Vandermonde) = missing

# ╔═╡ Cell order:
# ╠═372d3cf2-6173-11eb-356e-23c959c3fd89
# ╟─eb0428ac-5d8c-11eb-09a3-2b3cfc77f3f4
# ╟─70be3952-5d8c-11eb-1509-b3f7077d57e0
# ╠═acd7de0c-5d8c-11eb-120a-8b79f2b8eb3b
# ╟─af8c6460-5d8c-11eb-3ba8-c16e8855e992
# ╠═b96a33b8-5d8c-11eb-01de-439f53cdc355
# ╟─bcacf89e-5d8c-11eb-0077-e5761b8855a3
# ╟─c7a077ba-5d8c-11eb-08bb-07c343ea8ab1
# ╠═cd0d8636-5d8c-11eb-19f9-4da4550d306f
# ╟─ce25e25c-5d8c-11eb-2e8e-b5b1e7350d70
# ╠═d672cc72-5d8c-11eb-2c06-0341181e3a3d
# ╟─dac11770-5d8c-11eb-1058-2d043e172931
# ╠═61c43794-6174-11eb-1545-2db114b929e4
# ╟─e7b9023a-5d8c-11eb-1387-cfa7c41ab6ca
# ╟─fb6e62d4-5d8c-11eb-34f3-bf3df7cd4cb3
# ╠═01fe6f9a-5d8d-11eb-0519-03aefcd587bb
# ╠═05ae1a2a-5d8d-11eb-0ed8-496f48194232
# ╟─18d04a36-5d8d-11eb-1986-693eaad5d5be
# ╠═1df0419c-5d8d-11eb-3444-d5741bf10d32
# ╠═234dd2e4-5d8d-11eb-366c-c15474da799f
# ╟─23f7e8c4-5d8d-11eb-04af-b5e268c13a14
# ╟─330d2090-5d8d-11eb-1003-a52a078514b2
# ╠═375bf8d8-5d8d-11eb-1a49-69905d38effe
# ╠═39d6d5c4-5d8d-11eb-0e07-11d891ff87a3
# ╟─3f5740f6-5d8d-11eb-19c3-ddedf6003e53
# ╟─47049bf0-5d8d-11eb-18c2-733287b69420
# ╠═50f17586-5d8d-11eb-0eec-579467b787d0
# ╟─585995ae-5d8d-11eb-256f-bd8e9eb52063
# ╟─62e49c94-5d8d-11eb-39ac-f30febf282ff
# ╟─e5c61ace-654a-11eb-05be-79153b6222b6
# ╠═5c141870-654b-11eb-26f3-4b58b0b2e25d
# ╠═4d2d4e6c-654b-11eb-02c5-698af6d43b37
# ╠═6f69292e-654b-11eb-012e-317c7b2f33ea
# ╠═78f6ea1c-654b-11eb-217e-6da9345f2fe3
# ╠═7ea3de2a-654b-11eb-18c2-0fa369696cb9
# ╠═3422d2fc-654b-11eb-37c8-93a84229ab21
# ╠═340538c8-654b-11eb-3733-272fb9a8f106
# ╠═b1525c16-654b-11eb-0428-43008354b402
# ╟─6e8548b4-5d8d-11eb-3fcc-45cb005e5c5e
# ╠═cf6dea8c-5d8d-11eb-3f54-4d947305f5e5
# ╟─d305007e-5d8d-11eb-2505-3347d7d1a561
# ╠═d84355d4-5d8d-11eb-2c2c-21daf0364c21
# ╠═db651312-5d8d-11eb-09f1-bb693144d8f7
# ╠═8b7688ce-5d8e-11eb-0100-dd2712247ff6
# ╠═9c0f7916-5d8e-11eb-15cd-11badf9f9c01
# ╟─dd5706ec-5d8d-11eb-29df-75170b547cf8
# ╠═e799c73e-5d8d-11eb-3df4-13dbef13a342
# ╠═eb563ede-5d8d-11eb-1e91-55bd3b36f4eb
# ╠═fc2d5c92-5d8d-11eb-08c3-29d0478f1884
# ╠═072f2382-5d8e-11eb-377a-d39d0821559e
# ╟─17022b06-5d8e-11eb-391b-7194962a2c18
# ╠═393e2f9e-5d8e-11eb-11d7-511c32ce6e48
# ╟─4b262a72-5d8e-11eb-02bb-5fb2bd07f280
# ╠═605779e6-5d8e-11eb-3e08-c7420ef76aba
# ╠═6646eb5c-5d8e-11eb-1e01-f3011c4230de
# ╠═6960f8d2-5d8e-11eb-0215-2de7b54e3081
# ╟─8321136a-5d8e-11eb-0da2-e750b81d7ce9
# ╠═e246757e-5d8e-11eb-313c-61a63246cf9a
# ╠═ecc14f96-5d8f-11eb-2b18-cb23fad26b6e
# ╟─9d4b9eee-5d8f-11eb-07bb-557415e4ac4a
# ╟─142666de-5d90-11eb-3231-efe53ddc9b0d
# ╠═27fcaede-5d90-11eb-1cea-91fcc4b6b0fe
# ╠═2a224fde-5d90-11eb-1c46-3fd248350914
# ╠═2c937298-5d90-11eb-06e6-ab70b9d3701e
# ╟─2efe030e-5d90-11eb-38ae-092222d3a8d4
# ╟─46df4eb0-5d90-11eb-1fdf-f34a7bcb7191
# ╠═a4a1cb18-5d90-11eb-08ee-8570368a056b
# ╟─aca1c930-5d90-11eb-29d9-954e097bbe3b
# ╠═a9502b64-5d90-11eb-144c-3d7ce0949e67
# ╠═da6cc5b8-5d90-11eb-07ff-db8bdb504054
# ╠═e3759d4c-5d90-11eb-0bea-bb4247623ec2
# ╠═07998440-5d91-11eb-1a65-8de428eac89c
# ╠═e11b0b10-6621-11eb-0bdb-f3719cc92a20
# ╟─0b8980aa-5d91-11eb-385b-71568ec0e325
# ╠═192d9fd4-5d91-11eb-1cb9-c706aad03480
# ╠═1fa68c72-5d91-11eb-1102-c18460e92ee6
# ╟─26a70358-5d91-11eb-0241-699b5dc4783f
# ╠═2270e790-5d91-11eb-20e5-29905f232734
# ╟─2cbe2370-5d91-11eb-130a-5da52d1a62c8
# ╠═49f1d98c-5d91-11eb-1657-f320e9fcdc0e
# ╠═4cb68744-5d91-11eb-2b3e-e7df55888c93
# ╠═e99af5c0-6621-11eb-058b-45c3719930d0
# ╟─579e3828-5d91-11eb-1d33-b94628d61fc0
# ╠═e9a99a00-5d91-11eb-2c50-8be452cab83f
# ╠═ec62c35c-5d91-11eb-3773-b9385f312f7f
# ╠═efb0b460-5d91-11eb-2534-496df689dc60
# ╠═f3c3114c-5d91-11eb-1d37-6d97ea6d267f
# ╠═04dcda58-5d92-11eb-10ba-396947081338
# ╠═fbdb2958-6621-11eb-3cb6-a9bdeea3bdb7
# ╠═066901bc-6622-11eb-23ad-ed869adad36d
# ╠═0f878dea-5d92-11eb-0000-b7484532ee70
# ╠═11630c02-5d92-11eb-1746-4dabf327fbbe
# ╠═046ce4f8-6622-11eb-3c4f-7b6bf21fb77b
# ╠═1e65cb9c-5d92-11eb-3526-332169917fd9
# ╠═201f59ee-5d92-11eb-33ae-51904d249dd4
# ╟─310d8966-6176-11eb-1b23-9d942f197fbd
# ╠═276e9af4-5d92-11eb-1399-993570859698
# ╠═300a8428-5d92-11eb-188b-05d00df4f6a7
# ╟─cf146bf4-6177-11eb-1eaa-c35efde57b3e
# ╟─dd0b7c7a-6177-11eb-2cb7-8b194a75d776
# ╠═3ae60e88-5d94-11eb-0c50-1d74ea104758
# ╠═48301af2-5d94-11eb-0019-7737667c9cea
# ╠═4ea80eda-5d94-11eb-3882-21a41d2d65f8
# ╠═5af60d90-5d94-11eb-2ee4-b7bfc2caf53b
# ╠═5f7c75ac-5d94-11eb-137a-7914cd009821
# ╠═678ca64a-5d94-11eb-2b85-0b706526e35b
# ╠═8afbd434-5d94-11eb-366b-d3d719189ef7
# ╠═9226b2d8-5d94-11eb-2bea-491eb7dc1da7
# ╠═95146d46-5d94-11eb-22aa-c1a544e0d784
# ╠═9eab40be-5d94-11eb-0c59-21f5824fb812
# ╠═a137e0f8-5d94-11eb-2209-73acad549307
# ╠═a79ba114-5d94-11eb-16ae-9906c6cdf54f
# ╟─392228e2-617d-11eb-09a5-c9e5649356eb
# ╟─902f4dfe-617d-11eb-2957-71130adca3ae
# ╠═d448a2e0-5d92-11eb-18a6-9ff817992154
# ╠═bd91a60e-5d93-11eb-09d4-830ca69439bf
# ╠═c2ecfec8-5d93-11eb-2640-07bc07f3da98
# ╠═cb3e91cc-5d93-11eb-020c-d73c10131755
# ╠═d107c75e-5d93-11eb-0e6f-097b1291e460
# ╟─7f02b0a0-617f-11eb-1263-91423840def3
# ╟─f01448f0-617d-11eb-1829-0fcfe19b3115
# ╠═d2a076ea-5d93-11eb-216e-f5c37d330b40
# ╠═dc945902-5d93-11eb-1121-a7ae99c5862e


text:
export arrow_out, arrow_in, arrow_in_out, winding_plot,
    count_coils, winding_mmf, mmf_sum, mmf_plot, winding_mmf_plot,
    slot_label, mmf_label

function arrow_out(x, y; r = 0.2, color = "black", linewidth = 1, ri = 0.05)
    phi = collect(0:pi/36:2*pi)
    fill(x .+ r .* cos.(phi), y .+ r .* sin.(phi),
        color = "white", linewidth = linewidth)
    plot(x .+ r .* cos.(phi), y .+ r .* sin.(phi),
        color = color, linewidth = linewidth)
    fill(x .+ ri .* cos.(phi), y .+ ri .* sin.(phi),
        color = color, linewidth = linewidth)
end

function arrow_in(x, y; r = 0.2, color = "black", linewidth = 1)
    phi = collect(0:pi/36:2*pi)
    fill(x .+ r .* cos.(phi), y .+ r .* sin.(phi),
        color = "white", linewidth = linewidth)
    plot(x .+ r .* cos.(phi), y .+ r .* sin.(phi),
        color = color, linewidth = linewidth)
    plot([x - r / sqrt(2), x + r/sqrt(2)], [y - r / sqrt(2), y + r/sqrt(2)],
        color = color, linewidth = linewidth)
    plot([x - r / sqrt(2), x + r/sqrt(2)], [y + r / sqrt(2), y - r/sqrt(2)],
        color = color, linewidth = linewidth)
end

function arrow_in_out(x, y; mag = 1, r = 0.2, color = "black", linewidth = 1,
    ri = 0.05)
    # mag > 0 plot arrow_out; regular size only with mag = 1
    # mag < 0 plot arrow_in; regular size only with mag = -1
    # mag = 0  does not plot an arrow
    # Size of arrow is scaled with mag
    if mag > 0
        arrow_out(x, y, r = r * mag, color = color, linewidth = linewidth,
            ri = ri)
    elseif mag < 0
        arrow_in(x, y,r = -r * mag, color = color, linewidth = linewidth)
    end
end

function winding_plot(w; r = 0.2,
    color = ["#E6550D", "#8CA252", "#6B6ECF"],
    slot = 1, yoke = 1, fillcolor = "#D9D9D9")
    # Determine the number of layers per slot
    Nl = size(w,1)
    Ns = size(w,2)
    for y in 1:Nl
        for x in 1:Ns
            # plot arrows
            arrow_in_out(x, Nl - y + 1, mag = sign(w[y,x]), color = color[abs(w[y,x])])
        end
    end
    # Create slot structure
    if abs(slot) == 1
        xcontour = cat([0.5, 0.5],
            repeat(collect(0.75:0.5:Ns + 0.4), inner=[2]),
            [Ns + 0.5, Ns + 0.5], dims = 1)
        if slot == 1
            ycontour = cat([0.5 - yoke],
                repeat([Nl+0.5, Nl+0.5, 0.5, 0.5], outer=[Ns]),
                [Nl+0.5, Nl+0.5, 0.5 - yoke], dims = 1)
        else
            ycontour = cat([Nl + yoke + 0.5],
                repeat([0.5, 0.5, Nl+0.5, Nl+0.5], outer=[Ns]),
                [0.5, 0.5, Nl + yoke + 0.5], dims = 1)
        end
    end
    # Plot slots
    fill(xcontour, ycontour, color = fillcolor)
    # Limit horiziontal range
    xlim(0,Ns+1)
    # Make a circle to appear as circle
    axis("equal")
    # Add tickst at the slot locations
    xticks(collect(1:Ns),fill("",Ns))
    # Do not use vertical ticks
    yticks([])
    # gca().spines["left"].set_visible(false)
    # gca().spines["top"].set_visible(false)
    # gca().spines["right"].set_visible(false)
    # gca().spines["bottom"].set_visible(false)
end

function count_coils(w; positive = true)
    # Number of phase
    m = max(maximum(w),minimum(-w))
    # Count positive and negative coil sides
    coils_pos = zeros(Int,m)
    coils_neg = zeros(Int,m)
    # Zero matrix for comparison
    coils_zer = zeros(Int,m)
    # Loop counting
    for k in 1:size(w,1)
        for j in 1:m
            coils_pos[j] += count(w[k,:] .== +j)
            coils_neg[j] += count(w[k,:] .== -j)
        end
    end

    # In case there are more positive than negative coil sides, a warning is
    # caused
    if coils_pos - coils_neg != coils_zer
        @warn "count_coils: Coils sides are not balanced"
        for j in 1:m
            if  coils_pos[j] != coils_neg[j]
                println("Phase ", j, " has ", coils_pos[j], " positive coil sides")
                println("Phase ", j, " has ", coils_neg[j], " negative coil sides")
            end
        end
    end
    # Return number of coils
    if positive
        return coils_pos
    else
        return coils_neg
    end
end

@doc raw"""
# Function call

`mmfwinding(w, i)`

# Description

Calculates the magnetomotive force distribution based on the winding layout
and the actual currents of polyphase winding

# Input variable

`w[:,:]` Winding array with equal number of windings per slot, where

- `size(w,1)` = number of layers, e.g. 1 or 2
- `size(w,2)` = number of slots

`i[:]` Actual phase currents for, e.g.,  the three phases, scaled to 1

# Output variables

`mmf[:,:]` Magnetomotive force, where

- `size(mmf,1)` = number phases
- `size(mmm,1)` = number of slots
"""
function winding_mmf(w, i)
    # Number of slots
    Ns = size(w, 2);
    # Number of layers
    Nl = size(w, 1)
    # Determine number of coils
    coils = count_coils(w)
    # Determine number of phases
    m = size(coils,1)

    # Create relative MMF curves
    wp = zeros(m, Ns);
    for l = 1:Nl
        for k = 1:Ns
            phase = abs(w[l, k])
            orientation = sign(w[l, k])
            # Do not consider invalid phase indices
            if phase > 0
                wp[phase, k] += orientation
            end
        end
    end
    # Initialize mmf
    mmf = zeros(size(wp))
    # Make mean value equal to zero
    for k = 1:m
        # Create cumulative sum
        mmf[k,:] = Base.cumsum(wp[k,:])
        mmf[k,:] = mmf[k,:] .- Statistics.mean(mmf[k,:]);
    end
    # Scale with I and N
    for k = 1:m
        mmf[k,:] = mmf[k,:] * i[k];
    end
    return mmf
end

function mmf_sum(mmf)
    return sum(mmf, dims = 1)
end

function slot_label(Ns; start = 1, inc = 1)
    ticks_label = fill("", Ns)
    for k in start:Ns
        if mod(k + start - 2, inc) == 0
            ticks_label[k] = string(k)
        end
    end
    xticks(collect(1:Ns), ticks_label)
end

function mmf_label(mmf1_max = 0; mmf_max = 0, inc = 0.5, showmmf = true,
    mmf1_label_neg = L"$-\hat V_{m,1}$", mmf1_label_pos = L"$+\hat V_{m,1}$",
    mmf_label_neg = L"$-\hat V_{m}$", mmf_label_pos = L"$+\hat V_{m}$")

    # Pick the greater of the two maxima
    mmf_limit = max(mmf1_max, mmf_max)

    # Create vector of positive ticks
    ticks_pos = collect(+inc:+inc:+mmf_limit)
    # Create vector of negative ticks
    ticks_neg = collect(-mmf_limit:+inc:-inc)
    # Create zero tick
    ticks_zer = [0]
    # Concatenate ticks
    ticks = vcat(ticks_neg,ticks_zer,ticks_pos)
    # Add only zero as label
    ticks_label = vcat(fill("",length(ticks_neg)),
        L"$0$", fill("",length(ticks_pos)))
    # Add MMF labels only if showmmf = true
    if showmmf
        # Add ±mmf1_max if not zero
        if mmf1_max != 0 && mmf1_max != mmf_max
            ticks = vcat(ticks,[-mmf1_max],[+mmf1_max])
            ticks_label = vcat(ticks_label,
                [mmf1_label_neg],[mmf1_label_pos])
        end
        # Add ±mmf_max if not zero
        if mmf_max != 0
            ticks = vcat(ticks,[-mmf_max],[+mmf_max])
            ticks_label = vcat(ticks_label,
                [mmf_label_neg],[mmf_label_pos])
        end
    end
    # Plot ticks
    yticks(ticks, ticks_label)
end

function mmf_plot(mmf; index = collect(1:size(mmf,1)),
    color = ["#E6550D", "#8CA252", "#6B6ECF", "black"],
    showsum = true, label = ["1", "2", "3", L"$\Sigma$"],
    showlegend = true, loc = "best",
    linestyle=[lineStyle1,lineStyle2,lineStyle3,lineStyle4],
    linewidth=[lineWidth1,lineWidth3,lineWidth2,lineWidth4])

    # Number of phases
    m = size(mmf,1)
    # Number of slots
    Ns = size(mmf,2)
    # Array of slot indices
    slots = collect(1:Ns)
    # Plot step curve
    for k in index
        step([0.5; slots; Ns + 0.5], [ mmf[k,end]; mmf[k,end]; mmf[k,:]],
            color = color[k], linestyle = linestyle[k],
            linewidth = linewidth[k], label = label[k]);
    end
    # Local maximum and minimum of the largest phase MMF
    mmf1_max = maximum(maximum(mmf,dims=2))
    mmf1_min = minimum(minimum(mmf,dims=2))

    if showsum
        # Determine sum of all MMF curves
        mmfsum = mmf_sum(mmf)
        # Determine the maximum of the total MMF
        mmf_max = max(maximum(maximum(mmf,dims=2)),maximum(mmfsum))
        mmf_min = min(minimum(minimum(mmf,dims=2)),minimum(mmfsum))
        step([0.5; slots; Ns + 0.5], [mmfsum[end]; mmfsum[end]; mmfsum'],
            linewidth = linewidth[m+1], color = color[m+1] ,
            label = label[m+1])
    else
        mmf_max = mmf1_max
        mmf_min = mmf1_min
    end
    # Limit horizontal range
    xlim(0,Ns+1)
    # Add tickst at the slot locations
    xticks(collect(1:Ns),fill("",Ns))
    # Do not use vertical ticks

    yt = collect(mmf_min:0.5:mmf_max)
    ylim(mmf_min-0.2,mmf_max+0.2)
    Nyt = size(yt,1)
    yticks(yt,fill("",Nyt))
    grid(true)
    legend(loc =loc, fontsize = legendFontSize);
    return (mmf1_max,mmf_max)
end

function winding_mmf_plot(w, i; r = 0.2,
    color = ["#E6550D", "#8CA252", "#6B6ECF", "black"],
    slot = 1, yoke = 1, fillcolor = "#D9D9D9",
    showsum = true, label = ["1", "2", "3", L"$\Sigma$"],
    linestyle=[lineStyle1,lineStyle3,lineStyle2,lineStyle4],
    linewidth=[lineWidth1,lineWidth3,lineWidth2,lineWidth4],
    showlegend = true, loc = "best",
    index = collect(1:size(winding_mmf(w,i),1)),
    showslot = true, start = 1, inc = 1, showmmf = true,
    mmf1_label_neg = L"$-\hat V_{m,1}$", mmf1_label_pos = L"$+\hat V_{m,1}$",
    mmf_label_neg = L"$-\hat V_{m}$", mmf_label_pos = L"$+\hat V_{m}$")

    # Plot winding layout
    # https://matplotlib.org/3.1.1/tutorials/intermediate/tight_layout_guide.html#sphx-glr-tutorials-intermediate-tight-layout-guide-py
    subplot2grid((3, 1), (0, 0))
    winding_plot(w, r = 0.2,
        color = color,
        slot = slot, yoke = yoke, fillcolor = fillcolor)
    # Determine exact extension of plot area
    ax = gca()

    # Plot MMF
    subplot2grid((3, 1), (1, 0), rowspan = 2)
    mmf = winding_mmf(w,i)
    # Plot MMF distribution
    (mmf1_max, mmf_max) = mmf_plot(mmf, index = collect(1:size(mmf,1)),
        color = color,
        showsum = showsum, label = label,
        linestyle = linestyle,
        linewidth = linewidth,
        showlegend = showlegend, loc = loc)
    if showslot
        slot_label(size(mmf,2), start = start, inc = inc)
    end
    # Re-adjust by doing everything twice
    mmf_label(mmf1_max, mmf_max = mmf_max, showmmf = showmmf,
        mmf1_label_neg = mmf1_label_neg, mmf1_label_pos = mmf1_label_pos,
        mmf_label_neg = mmf_label_neg, mmf_label_pos = mmf_label_pos)
    tight_layout()
    xlim(ax.axis()[1],ax.axis()[2])
    tight_layout()
    xlim(ax.axis()[1],ax.axis()[2])
end


text:
"""
`module  JAC.PairAnnihilation1Photon`  ... a submodel of JAC that contains all methods for computing positron-bound-electron pair annihilation
                                           (PEPA) with single-photon emission cross sections and rates; e^+ + |i(N)> --> |f(N-1)> + photon;
                                           it is using JAC, JAC.ManyElectron, JAC.Radial.
"""
module PairAnnihilation1Photon

    using JAC, JAC.ManyElectron, JAC.Radial
    global JAC_counter = 0

    """
    `struct  PairAnnihilation1Photon.Settings`  ... defines a type for the details and parameters in computing positron-bound-electron pair
             annihilation (PEPA) with single-photon emission cross sections and rates; e^+ + |i(N)> --> |f(N-1)> + photon.

        + multipoles              ::Array{EmMultipole}           ... Specifies the multipoles of the radiation field that are to be included.
        + gauges                  ::Array{UseGauge}              ... Specifies the gauges to be included into the computations.
        + positronEnergies        ::Array{Float64,1}             ... List of positron energies.
        + printBeforeComputation  ::Bool                         ... True, if all energies and lines are printed before their evaluation.
        + selectLines             ::Bool                         ... True, if lines are selected individually for the computations.
        + selectedLines           ::Array{Tuple{Int64,Int64},1}  ... List of lines, given by tupels (inital-level, final-level).
    """
    struct Settings
        multipoles                ::Array{EmMultipole}
        gauges                    ::Array{UseGauge}
        positronEnergies          ::Array{Float64,1}
        printBeforeComputation    ::Bool
        selectLines               ::Bool
        selectedLines             ::Array{Tuple{Int64,Int64},1}
    end

    """
    `JAC.PairAnnihilation1Photon.Settings()`  ... constructor for the default values of pair-annihilation photon line computations
    """
    function Settings()
        Settings(EmMultipole[], UseGauge[], Float64[], false, false, Tuple{Int64,Int64}[])
    end

     # `Base.show(io::IO, settings::PairAnnihilation1Photon.Settings)`
     #		... prepares a proper printout of the variable settings::PairAnnihilation1Photon.Settings.
    function Base.show(io::IO, settings::PairAnnihilation1Photon.Settings)
        println(io, "multipoles:               $(settings.multipoles)  ")
        println(io, "gauges:                   $(settings.gauges)  ")
        println(io, "positronEnergies:         $(settings.positronEnergies)  ")
        println(io, "printBeforeComputation:   $(settings.printBeforeComputation)  ")
        println(io, "selectLines:              $(settings.selectLines)  ")
        println(io, "selectedLines:            $(settings.selectedLines)  ")
    end

    """
    `struct  PairAnnihilation1Photon.Channel`  ... defines a type for a positron-bound-electron pair annihilation (PEPA) with single-photon
             emission channel that specifies all quantum numbers, phases and amplitudes.

        + multipole      ::EmMultipole          ... Multipole of the photon absorption.
        + gauge          ::EmGauge              ... Gauge for dealing with the (coupled) radiation field.
        + kappa          ::Int64                ... partial-wave of the incoming free positron
        + symmetry       ::LevelSymmetry        ... total angular momentum and parity of the scattering state
        + phase          ::Float64              ... phase of the partial wave
        + amplitude      ::Complex{Float64}     ... PairAnnihilation1PhotonChannel amplitude associated with the given channel.
   """
    struct  Channel
        multipole        ::EmMultipole
        gauge            ::EmGauge
        kappa            ::Int64
        symmetry         ::LevelSymmetry
        phase            ::Float64
        amplitude        ::Complex{Float64}
    end

    """
    `struct  PairAnnihilation1Photon.Line`  ... defines a type for a positron-bound-electron pair-annihilation (photon) line that may
                                                include the definition of channels.

        + initialLevel   ::Level                 ... initial-(state) level
        + finalLevel     ::Level                 ... final-(state) level
        + positronEnergy ::Float64                ... Energy of the (incoming free) positron.
        + photonEnergy   ::Float64                ... Energy of the emitted photon.
        + crossSection   ::EmProperty             ... Cross section for this pair-annihilation (photon) line.
        + hasChannels    ::Bool                   ... Determines whether the individual (sub-) channels are defined in terms of their
                                                      free-positron energy, kappa, multipole, etc., or not.
        + channels       ::Array{PairAnnihilation1Photon.Channel,1}  ... List of PairAnnihilation1Photon.Channels of this line.
    """
    struct  Line
        initialLevel     ::Level
        finalLevel       ::Level
        positronEnergy   ::Float64
        photonEnergy     ::Float64
        crossSection     ::EmProperty
        hasChannels      ::Bool
        channels         ::Array{PairAnnihilation1Photon.Channel,1}
    end

    """
    `JAC.PairAnnihilation1Photon.Line()`  ... 'empty' constructor for a pair-annihilation (photon) line between a specified initial and
                                              final level.
    """
    function Line()
        Line(Level(), Level(), 0., 0., EmProperty(0., 0.), false, PairAnnihilation1Photon[] )
    end

    # `Base.show(io::IO, line::PairAnnihilation1Photon.Line)`  ... prepares a proper printout of the variable line::PairAnnihilation1Photon.Line.
    function Base.show(io::IO, line::PairAnnihilation1Photon.Line)
        println(io, "initialLevel:      $(line.initialLevel)  ")
        println(io, "finalLevel:        $(line.finalLevel)  ")
        println(io, "positronEnergy:    $(line.positronEnergy)  ")
        println(io, "photonEnergy:      $(line.photonEnergy)  ")
        println(io, "crossSection:      $(line.crossSection)  ")
        println(io, "hasChannels:       $(line.hasChannels)  ")
        println(io, "channels:          $(line.channels)  ")
    end

    """
    `JAC.PairAnnihilation1Photon.computeAmplitudesProperties(line::PairAnnihilation1Photon.Line, grid::Radial.Grid,
                                                             settings::PairAnnihilation1Photon.Settings)`  ... to compute all amplitudes and
         properties of the given line; a line::PairAnnihilation1Photon.Line is returned for which the amplitudes and properties are now evaluated.
    """
    function  computeAmplitudesProperties(line::PairAnnihilation1Photon.Line, grid::Radial.Grid, settings::PairAnnihilation1Photon.Settings)
        global JAC_counter
        newChannels = PairAnnihilation1Photon.Channel[]
        for channel in line.channels
            # Generate a continuum orbital
            JAC_counter = JAC_counter + 1
            if   JAC_counter < 20   println("PairAnnihilation1Photon.computeAmplitudes..-aa: warning ... no cont. orbital is generated.")  end
            phase  = 0.
            # Define a proper continuum basis from the initialLevel.basis and the continuum orbital
            JAC_counter = JAC_counter + 1
            if   JAC_counter < 20   println("PairAnnihilation1Photon.computeAmplitudes..-ab: warning ... no coninuum basis is yet generated.") end
            # Compute the transition matrix for the continuum and the initial-state basis
            JAC_counter = JAC_counter + 1
            if   JAC_counter < 20   println("PairAnnihilation1Photon.computeAmplitudes..-ac: warning ... no transition matrix is computed.") end
            # matrix    = JAC.PhotoIonization.computeMatrix(channel.multipole, channel.gauge, line.omega, line.finalLevel.basis,
            #                                               line.initialLevel.basis, grid, settings)
            # amplitude = line.finalLevel.mc * matrix * line.initialLevel.mc
            amplitude = 1.0
            push!( newChannels, PairAnnihilation1Photon.Channel( channel.multipole, channel.gauge, channel.kappa, channel.symmetry,
                                                                 phase, amplitude) )
        end
        # Calculate the photonrate and angular beta if requested
        JAC_counter = JAC_counter + 1
        if   JAC_counter < 20   println("PairAnnihilation1Photon.computeAmplitudesProperties-ba: warning ... cs set to -1.") end
        crossSection = EmProperty(-1., -1.)
        line = PairAnnihilation1Photon.Line( line.initialLevel, line.finalLevel, line.positronEnergy, line.photonEnergy,
                                             crossSection, true, newChannels)
        return( line )
    end

    """
    `JAC.PairAnnihilation1Photon.computeLines(finalMultiplet::Multiplet, initialMultiplet::Multiplet, grid::Radial.Grid,
                                              settings::PairAnnihilation1Photon.Settings; output::Bool=true)`  ... to compute the
         pair-annihilation single-photon emission amplitudes and all properties as requested by the given settings. A list of
         lines::Array{PairAnnihilation1Photon.Lines} is returned.
    """
    function  computeLines(finalMultiplet::Multiplet, initialMultiplet::Multiplet, grid::Radial.Grid, settings::PairAnnihilation1Photon.Settings;
                           output::Bool=true)
        lines = JAC.PairAnnihilation1Photon.determineLines(finalMultiplet, initialMultiplet, settings)
        # Display all selected lines before the computations start
        if  settings.printBeforeComputation    JAC.PairAnnihilation1Photon.displayLines(lines)    end
        # Calculate all amplitudes and requested properties
        newLines = PairAnnihilation1Photon.Line[]
        for  line in lines
            newLine = JAC.PairAnnihilation1Photon.computeAmplitudesProperties(line, grid, settings)
            push!( newLines, newLine)
        end
        # Print all results to screen
        JAC.PairAnnihilation1Photon.displayResults(lines)
        #
        if    output    return( lines )
        else            return( nothing )
        end
    end

    """
    `JAC.PairAnnihilation1Photon.determineChannels(finalLevel::Level, initialLevel::Level, settings::PairAnnihilation1Photon.Settings)`
         ... to determine a list of pair-annihilation single-photon emission Channel for a transitions from the initial to final level and by
         taking into account the particular settings of for this computation; an Array{PairAnnihilation1Photon.Channel,1} is returned.
    """
    function determineChannels(finalLevel::Level, initialLevel::Level, settings::PairAnnihilation1Photon.Settings)
        channels = PairAnnihilation1Photon.Channel[];
        symi = LevelSymmetry(initialLevel.J, initialLevel.parity);    symf = LevelSymmetry(finalLevel.J, finalLevel.parity)
        for  mp in settings.multipoles
            for  gauge in settings.gauges
                symList = JAC.AngularMomentum.allowedMultipoleSymmetries(symi, mp)
                ##x println("mp = $mp   symi = $symi   symList = $symList")
                for  symt in symList
                    kappaList = JAC.AngularMomentum.allowedKappaSymmetries(symt, symf)
                    for  kappa in kappaList
                        # Include further restrictions if appropriate
                        if     string(mp)[1] == 'E'  &&   gauge == Basics.UseCoulomb
                            push!(channels, PairAnnihilation1Photon.Channel(mp, JAC.Coulomb,   kappa, symt, 0., Complex(0.)) )
                        elseif string(mp)[1] == 'E'  &&   gauge == Basics.UseBabushkin
                            push!(channels, PairAnnihilation1Photon.Channel(mp, JAC.Babushkin, kappa, symt, 0., Complex(0.)) )
                        elseif string(mp)[1] == 'M'
                            push!(channels, PairAnnihilation1Photon.Channel(mp, JAC.Magnetic,  kappa, symt, 0., Complex(0.)) )
                        end
                    end
                end
            end
        end
        return( channels )
    end

    """
    `JAC.PairAnnihilation1Photon.determineLines(finalMultiplet::Multiplet, initialMultiplet::Multiplet,
                                                settings::PairAnnihilation1Photon.Settings)` ... to determine a list of
         PairAnnihilation1Photon.Line's for transitions between levels from the initial- and final-state multiplets, and  by taking into account
         the particular selections and settings for this computation; an Array{PairAnnihilation1Photon.Line,1} is returned. Apart from the
         level specification, all physical properties are set to zero during the initialization process.
    """
    function  determineLines(finalMultiplet::Multiplet, initialMultiplet::Multiplet, settings::PairAnnihilation1Photon.Settings)
        if    settings.selectLines    selectLines   = true;   selectedLines = Basics.determine("selected lines", settings.selectedLines)
        else                          selectLines   = false
        end

        lines = PairAnnihilation1Photon.Line[]
        for  i = 1:length(initialMultiplet.levels)
            for  f = 1:length(finalMultiplet.levels)
                if  selectLines  &&  !(haskey(selectedLines, (i,f)) )    continue   end
                for  ep in settings.positronEnergies
                    omega = ep - (finalMultiplet.levels[f].energy - initialMultiplet.levels[i].energy) + 2* Defaults.getDefaults("mc^2")
                    if  omega < 0    continue   end

                    channels = JAC.PairAnnihilation1Photon.determineChannels(finalMultiplet.levels[f], initialMultiplet.levels[i], settings)
                    push!( lines, PairAnnihilation1Photon.Line(initialMultiplet.levels[i], finalMultiplet.levels[f], ep, omega,
                                                               EmProperty(0., 0.), true, channels) )
                end
            end
        end
        return( lines )
    end

    """
    `JAC.PairAnnihilation1Photon.displayLines(lines::Array{PairAnnihilation1Photon.Line,1})`  ... to display a list of lines and channels that
         have been selected due to the prior settings. A neat table of all selected transitions and energies is printed but nothing is returned
         otherwise.
    """
    function  displayLines(lines::Array{PairAnnihilation1Photon.Line,1})
        println(" ")
        println("  Selected pair-annihilation single-photon emission lines:")
        println(" ")
        println("  ", JAC.TableStrings.hLine(175))
        sa = "  ";   sb = "  "
        sa = sa * JAC.TableStrings.center(18, "i-level-f"; na=2);                                sb = sb * JAC.TableStrings.hBlank(20)
        sa = sa * JAC.TableStrings.center(18, "i--J^P--f"; na=4);                                sb = sb * JAC.TableStrings.hBlank(22)
        sa = sa * JAC.TableStrings.center(38, "Energies  " * JAC.TableStrings.inUnits("energy"); na=4);
        sb = sb * JAC.TableStrings.center(38, " i--f       positron      omega"; na=4)
        sa = sa * JAC.TableStrings.flushleft(57, "List of multipoles, gauges, kappas and total symmetries"; na=4)
        sb = sb * JAC.TableStrings.flushleft(57, "partial (multipole, gauge, total J^P)                  "; na=4)
        println(sa);    println(sb);    println("  ", JAC.TableStrings.hLine(175))
        #
        for  line in lines
            sa  = "  ";    isym = LevelSymmetry( line.initialLevel.J, line.initialLevel.parity)
                           fsym = LevelSymmetry( line.finalLevel.J,   line.finalLevel.parity)
            sa = sa * JAC.TableStrings.center(18, JAC.TableStrings.levels_if(line.initialLevel.index, line.finalLevel.index); na=2)
            sa = sa * JAC.TableStrings.center(18, JAC.TableStrings.symmetries_if(isym, fsym); na=4)
            energy = line.initialLevel.energy - line.finalLevel.energy
            sa = sa * @sprintf("%.5e", Defaults.convertUnits("energy: from atomic", energy))              * "  "
            sa = sa * @sprintf("%.5e", Defaults.convertUnits("energy: from atomic", line.positronEnergy)) * "  "
            sa = sa * @sprintf("%.5e", Defaults.convertUnits("energy: from atomic", line.photonEnergy))   * "    "
            kappaMultipoleSymmetryList = Tuple{Int64,EmMultipole,EmGauge,LevelSymmetry}[]
            for  i in 1:length(line.channels)
                push!( kappaMultipoleSymmetryList, (line.channels[i].kappa, line.channels[i].multipole, line.channels[i].gauge,
                                                    line.channels[i].symmetry) )
            end
            wa = JAC.TableStrings.kappaMultipoleSymmetryTupels(85, kappaMultipoleSymmetryList)
            sb = sa * wa[1];    println( sb )
            for  i = 2:length(wa)
                sb = JAC.TableStrings.hBlank( length(sa) ) * wa[i];    println( sb )
            end
        end
        println("  ", JAC.TableStrings.hLine(175))
        #
        return( nothing )
    end

    """
    `JAC.PairAnnihilation1Photon.displayResults(lines::Array{PairAnnihilation1Photon.Line,1})`  ... to list all results, energies, cross sections,
         etc. of the selected lines. A neat table is printed but nothing is returned otherwise.
    """
    function  displayResults(lines::Array{PairAnnihilation1Photon.Line,1})
        println(" ")
        println("  Pair-annihilation single-photon emission cross sections:")
        println(" ")
        println("  ", JAC.TableStrings.hLine(128))
        sa = "  ";   sb = "  "
        sa = sa * JAC.TableStrings.center(18, "i-level-f"; na=2);                                sb = sb * JAC.TableStrings.hBlank(20)
        sa = sa * JAC.TableStrings.center(18, "i--J^P--f"; na=4);                                sb = sb * JAC.TableStrings.hBlank(22)
        sa = sa * JAC.TableStrings.center(38, "Energies  " * JAC.TableStrings.inUnits("energy"); na=4);
        sb = sb * JAC.TableStrings.center(38, " i--f       positron      omega"; na=4)
        sa = sa * JAC.TableStrings.center(10, "Multipoles"; na=1);                              sb = sb * JAC.TableStrings.hBlank(13)
        sa = sa * JAC.TableStrings.center(30, "Cou -- Cross section -- Bab"; na=3)
        sb = sb * JAC.TableStrings.center(30, JAC.TableStrings.inUnits("cross section") * "          " *
                                              JAC.TableStrings.inUnits("cross section"); na=3)
        println(sa);    println(sb);    println("  ", JAC.TableStrings.hLine(128))
        #
        for  line in lines
            sa  = "  ";    isym = LevelSymmetry( line.initialLevel.J, line.initialLevel.parity)
                           fsym = LevelSymmetry( line.finalLevel.J,   line.finalLevel.parity)
            sa = sa * JAC.TableStrings.center(18, JAC.TableStrings.levels_if(line.initialLevel.index, line.finalLevel.index); na=2)
            sa = sa * JAC.TableStrings.center(18, JAC.TableStrings.symmetries_if(isym, fsym); na=4)
            energy = line.initialLevel.energy - line.finalLevel.energy
            sa = sa * @sprintf("%.5e", Defaults.convertUnits("energy: from atomic", energy))              * "  "
            sa = sa * @sprintf("%.5e", Defaults.convertUnits("energy: from atomic", line.positronEnergy)) * "  "
            sa = sa * @sprintf("%.5e", Defaults.convertUnits("energy: from atomic", line.photonEnergy))   * "     "
            multipoles = EmMultipole[]
            for  ch in line.channels
                multipoles = push!( multipoles, ch.multipole)
            end
            multipoles = unique(multipoles);   mpString = JAC.TableStrings.multipoleList(multipoles) * "          "
            sa = sa * JAC.TableStrings.flushleft(11, mpString[1:10];  na=2)
            sa = sa * @sprintf("%.6e", Defaults.convertUnits("cross section: from atomic", line.crossSection.Coulomb))     * "    "
            sa = sa * @sprintf("%.6e", Defaults.convertUnits("cross section: from atomic", line.crossSection.Babushkin))   * "    "
            println(sa)
        end
        println("  ", JAC.TableStrings.hLine(128))
        #
        return( nothing )
    end

end # module


text:
export AWG, CoilGeometry, WindingLayers

import Base: ==, ≈
const Length = Unitful.Length

struct AWG{T<:Real, L<:Unitful.Length, R<:ResistanceLength, Q<:Unitful.Quantity}
    gauge::T
    d_bare::L
    d_insul::L
    Ω::R
    ft_lb::Q
end

struct CoilGeometry{T<:Unitful.Length}
    id::T
    od::T
    len::T
end
CoilGeometry{T}(x::CoilGeometry) where {T<:Length} = CoilGeometry{T}(x.id,x.od,x.len)
Base.length(::CoilGeometry)=1
==(a::CoilGeometry, b::CoilGeometry) = a.id == b.id && a.od == b.od && a.len == b.len
≈(a::CoilGeometry, b::CoilGeometry) = a.id ≈ b.id && a.od ≈ b.od && a.len ≈ b.len
Base.:convert(::Type{T}, x::CoilGeometry) where {T<:CoilGeometry} = T(x)
Base.:promote_rule(::Type{CoilGeometry{T}},::Type{CoilGeometry{S}}) where {T,S} = CoilGeometry{promote_type(T,S)}

function CoilGeometry(id::Length,od::Length,len::Length, gap::Length=zero(id))
    CoilGeometry(promote(id+2gap, od-2gap, len-2gap)...)
end
function CoilGeometry(coil::CoilGeometry;id=coil.id,od=coil.od,len=coil.len)
    CoilGeometry(id,od,len)
end

struct WindingLayers
    layers::Int
    odd_layer::Int
    even_layer::Int
    odd_turns::Int
    even_turns::Int
    turns::Int
end

function WindingLayers(layers::T, odd_layer::T, even_layer::T) where {T<:Integer}
    div,rem = divrem(layers,2)
    odd_turns = odd_layer * (div+rem)
    even_turns = even_layer * div
    WindingLayers(layers, odd_layer, even_layer, odd_turns, even_turns, odd_turns+even_turns)
end


text:
const _GLOBALS_EXT = ".global.jls"

function globfile(Proj::Module, gid::Symbol, gids::Symbol...)
    gids = string.([gid, gids...])
    datdir(Proj, join(gids, "."), _GLOBALS_EXT)
end

function lglob(Proj::Module, gid::Symbol, gids::Symbol...)
    fname = globfile(Proj, gid, gids...)
    !isfile(fname) && error("global is missing, glob file '", basename(fname), "' not found")
    return ldat(fname; verbose = false)
end

function lglob(Proj::Module, gidv::Vector)
    map(gidv) do arg
        (arg isa Symbol) ? lglob(Proj, arg) : lglob(Proj, arg...)
    end
end

macro lglob(Proj::Symbol, ex, exs...)
    _lglob_macro_ex(Proj, ex, exs...)
end

macro lglob(ex, exs...)
    error("You must specify a project module Ex: `Proj gid`")
end

## ------------------------------------------------------------
function sglob(Proj::Module, dat, gid::Symbol, gids::Symbol...)
    fname = globfile(Proj, gid, gids...)
    sdat(Proj, dat, fname; verbose = false)
end

function sglob(Proj::Module; gids...)
    for (gid, dat) in gids
        sglob(Proj, dat, gid)
    end
end

sglob(f::Function, Proj::Module, gid::Symbol, gids::Symbol...) = sglob(Proj, f(), gid, gids...)

macro sglob(Proj::Symbol, ex, exs...)
    _sglob_macro_ex(Proj, ex, exs...)
end

macro sglob(ex, exs...)
    error("You must specify a project module Ex: `Proj gid=dat`")
end

## ------------------------------------------------------------
delglob(Proj::Module, gid::Symbol, gids::Symbol...) = rm(globfile(Proj, gid, gids...); force = true)


text:
## --- Load observed cooling against which to fit

    # McDonough bulk silicate earth
    mcdbse = importdataset("data/mcdbse.csv",',')

    # # Average cooling curve
    # obs = importdataset("data/CoolingAverage.csv",',')
    # TrelObs_time = obs["Age"]/1000
    # TrelObs = obs["Trel"]
    # TrelObs_sigma = obs["Trel_sigma"]
    # TrelObs_sigma[1] = 0.5;

    # Keller and Schoene cooling curve
    obs = importdataset("data/KellerSchoene.psv",'|')
    TrelObs_time = obs["Age"]/1000
    TrelObs = obs["T_p"] .- obs["T_p"][1]
    TrelObs_sigma = obs["T_p_sigma"]
    TrelObs_sigma[1] = 0.5;

## -- Constant parameters

    mutable struct Parameters
        R::Float64 # Universal gas constant (J/Mol/K)
        g::Float64 # Gravitational acceleration (m/s^2)
        s_yr::Float64 # Seconds per year
        J_MeV::Float64 # Joules per megaelectron-volt
        mol::Float64 # Avogadro's number (Atoms per mole)
        lambda40K::Float64   # K-40 decay constant (1/years)
        lambda87Rb::Float64  # Rb-87 decay constant (1/years)
        lambda235U::Float64  # U-235 decay constant (1/years)
        lambda238U::Float64  # U-338 decay constant (1/years)
        lambda232Th::Float64 # Th-232 decay constant (1/years)

        Mm::Float64 # Mass of mantle (kg)
        Mc::Float64 # Mass of core (kg)
        Re::Float64 # Radius of Earth (m)
        Rp::Float64 # Radius of "plume feeding area"
        Cp_m::Float64  # Specific heat of mantle (J/kg/C)
        Cp_c::Float64  # Specific heat of core (J/kg/C)
        rho_m::Float64 # Density of mantle (kg/m^3)
        rho_c::Float64 # Density of core (kg/m^3)

        Qm_initial::Float64 # Initial mantle heat flux (W)
        Qc_initial::Float64 # Initial core heat flux (W)
        Tm_initial::Float64 # Initial mantle temperature (K)
        Tc_initial::Float64 # Initial core temperature (K)
        Qm_initial_sigma::Float64
        Qc_initial_sigma::Float64
        Tm_initial_sigma::Float64
        Tc_initial_sigma::Float64

        Qm_now::Float64 # Present mantle heat flux (W)
        Qc_now::Float64 # Present core heat flux (W)
        Tm_now::Float64 # Present mantle temperature (K)
        Tc_now::Float64 # Present core temperature (K)
        Qm_now_sigma::Float64
        Qc_now_sigma::Float64
        Tm_now_sigma::Float64
        Tc_now_sigma::Float64

        Qrcc_now::Float64   # Present-day crustal radiogenic heat production (W)
        Qrcc_now_sigma::Float64
        Qrm_now::Float64   # Present-day mantle radiogenic heat production (W)
        Qrm_now_sigma::Float64

        alpha::Float64 # Thermal expansion coefficient
        kappa::Float64 # Thermal diffusivity (m^2/s)

        Trl::Float64 # Reference temperature, CMB (C)
        mur_l::Float64 # Reference viscosity, CMB (Pa s)
        nc::Float64 # Core exponent
        Kc::Float64 # Condutivity of the core (W/m/C)

        Tru::Float64 # Reference temperature, upper mantle (C)
        mur_u::Float64 # Reference viscosity, upper mantle (Pa s)
        Ea::Float64 # Activation energy, upper mantle (J/mol)
        Ha::Float64 # Activation enthalpy, CMB (J/mol)
        Km::Float64 # Conductivity of the mantle (W/m/C)
        Xm::Float64 # Mantle heat capacity correction for average temperature
        Xc::Float64 # Core heat capacity correction for average temperature

        eta_L::Float64 # Lithospheric viscosity (Pa s)
        Rc::Float64 # Lithspheric plate bending radius of curvature (m)

        Q_addtl::Float64 # Constant additional heat flux to base of mantle (for sensitivity testing)
    end

    Parameters() = Parameters(8.3145, 9.8, 3.1556926e7, 1.602176565e-13, 6.022e23, 5.552729156131902e-10, 1.4203835667211993e-11, 9.845840632953767e-10, 1.551358953804712e-10, 4.933431890106372e-11, 4.0e24, 2.0e24, 6.371e6, 1.0e6, 1000.0, 800.0, 3500.0, 12300.0, 4.35e14, 2.58e13, 1800.0, 7000.0, 2.0e14, 1.0e13, 200.0, 2000.0, 3.5e13, 6.0e12, 1573.15, 4673.15, 1.0e13, 1.0e13, 100.0, 1000.0, 7.2e12, 1.0e12, 1.2e13, 3.0e12, 2.5e-5, 8.6e-7, 3773.15, 1.0e19, 0.2, 30.0, 1573.15, 1.0e21, 400000.0, 800000.0, 3.0, 1.2, 1.11, 1.0e23, 450000.0, 0.0)

    # # Define parameters
    # p = Parameters() # Create struct
    # p.R = 8.3145; # Universal gas constant (J/Mol/K)
    # p.g = 9.8; # Gravitational acceleration (m/s^2)
    # p.s_yr = 31556926; # Seconds per year
    # p.J_MeV=1.602176565E-13; # Joules per megaelectron-volt
    # p.mol=6.022E23; # Avogadro's number (Atoms per mole)
    # p.lambda40K  = log(2) / 1.2483E9 # K-40 decay constant (1/years)
    # p.lambda87Rb = log(2) / 4.88E10  # Rb-87 decay constant (1/years)
    # p.lambda235U = log(2) / 7.04E8   # U-235 decay constant (1/years)
    # p.lambda238U = log(2) / 4.468E9  # U-238 decay constant (1/years)
    # p.lambda232Th = log(2) / 1.405E10 # Th-232 decay constant (1/years)
    #
    # p.Mm = 4E24; # Mass of mantle (kg)
    # p.Mc = 2E24; # Mass of core (kg)
    # p.Re = 6371E3; # Radius of Earth (m)
    # p.Rp = 1000E3; # Radius of Davies' "plume feeding area"
    # p.Cp_m = 1000;  # Specific heat of mantle (J/kg/C)
    # p.Cp_c = 800;  # Specific heat of core (J/kg/C)
    # p.rho_m = 3500; # Density of mantle (kg/m^3)
    # p.rho_c = 12300; # Density of core (kg/m^3)
    #
    # p.Qm_initial = 43.5E12 # Initial mantle heat flux (W)
    # p.Qc_initial = 25.8E12 # Initial core heat flux (W)
    # p.Tm_initial = 1800.0 # Initial mantle temperature (K)
    # p.Tc_initial = 7000.0 # Initial core temperature (K).
    # p.Qm_initial_sigma = 20E12
    # p.Qc_initial_sigma = 10E12
    # p.Tm_initial_sigma = 200
    # p.Tc_initial_sigma = 2000
    #
    # p.Qm_now = 35E12 # Present mantle heat flux (W)
    # p.Qc_now = 6E12 # Present core heat flux (W)
    # p.Tm_now = 1300+273.15 # Present mantle temperature (K)
    # p.Tc_now = 4400+273.15 # Present core temperature (K). I/O boundary 6230 +/- 500 according to 10.1126/science.1233514
    # p.Qm_now_sigma = 10E12
    # p.Qc_now_sigma = 10E12
    # p.Tm_now_sigma = 100.0
    # p.Tc_now_sigma = 1000.0
    #
    # p.Qrcc_now       =  7.2e12 # Crustal radiogenic heat production today (W) (Jaupart average 7.6+/-1.9)
    # p.Qrcc_now_sigma =  1.0e12 # uncertainty (W)
    # p.Qrm_now        = 12.0e12 # Mantle radiogenic heat production today(W) (Jaupart average: 10.2+/-3)
    # p.Qrm_now_sigma  =  3.0e12 # uncertainty (W)
    #
    # p.alpha = 2.5E-5; # Thermal expansion coefficient
    # p.kappa = 0.86E-6; # Thermal diffusivity (m^2/s)
    #
    # p.Trl = 3500+273.15; # Reference temperature, CMB (C)
    # p.mur_l = 1E19; # Reference viscosity, CMB (Pa s)
    # p.nc = 1/5; # Core exponent
    # p.Kc = 30; # Condutivity of the core (W/m/C)
    #
    # p.Tru = 1300+273.15; # Reference temperature, upper mantle (C)
    # p.mur_u = 1E21; # Reference viscosity, upper mantle (Pa s)
    # p.Ea = 400E3; # Activation energy, upper mantle (J/mol)
    # p.Ha = 800E3; # Activation enthalpy CMB (J/mol)
    # p.Km = 3; # Conductivity of the mantle (W/m/C)
    # p.Xm = 1.2; # Mantle heat capacity correction for average temperature
    # p.Xc = 1.11; # Core heat capacity correction for average temperature
    #
    # p.eta_L = 1E23; # Plate viscsity (Pa s)
    # p.Rc = 4.5E5; # Lithspheric plate bending radius of curvature (m)
    #
    # p.Q_addtl = 0 # Additional (constant) heat flux
    # # Done defining parameters

## --- Calculate past mantle heat production model for a given time vector

    function bse_heat_production(p::Parameters,composition::Dict,timevec)
        # Composition: Bulk Silicate Earth
        K40i  = 0.00012 .* composition["K"]  .* exp.(p.lambda40K  .* timevec*10^6);
        Rb87i = 0.27835 .* composition["Rb"] .* exp.(p.lambda87Rb .* timevec*10^6);
        U235i = 0.00720 .* composition["U"]  .* exp.(p.lambda235U .* timevec*10^6);
        U238i = 0.99274 .* composition["U"]  .* exp.(p.lambda238U .* timevec*10^6);
        Th232i = composition["Th"] .* exp.(p.lambda232Th .* timevec*10^6)

        # Heat production in W/Kg, from Rybach, 1985
        H = 1000 .*K40i./10^6 .* p.lambda40K/p.s_yr  .* (0.6*0.8928 + 1.460*0.1072)*p.J_MeV*p.mol/40 +
            1000 .* Rb87i  ./10^6 .* p.lambda87Rb/p.s_yr .* 1/3*0.283*p.J_MeV*p.mol/87 +
            1000 .* U235i  ./10^6 .* p.lambda235U/p.s_yr .* 45.26*p.J_MeV*p.mol/235 +
            1000 .* U238i  ./10^6 .* p.lambda238U/p.s_yr .* 46.34*p.J_MeV*p.mol/238 +
            1000 .* Th232i ./10^6 .* p.lambda232Th/p.s_yr .* 38.99*p.J_MeV*p.mol/232

        return H
    end

## --- Likelihood and residual functions

    function LL(x_prop,x,sigma)
        return sum( -(x_prop .- x).^2 ./ (2*sigma.^2) .- log.(sqrt.(2*pi*sigma)))
    end

    # Return the signed log-likelihood for a given temperature curve
    function signed_LL(Tm, timevec, TrelAve_time, TrelAve, TrelAve_sigma)
        Trel = linterp1(timevec, Tm .- Tm[1], TrelAve_time)
        signed_LL = sum(sign.(Trel - TrelAve) .*
                (Trel - TrelAve).^2 ./ (2 .* TrelAve_sigma.^2)
                .- log.(sqrt.(2*pi*TrelAve_sigma))
            )
        return signed_LL
    end

    # Return the residual for a given temperature curve
    function resid(Tm, timevec, TrelAve_time, TrelAve, TrelAve_sigma)
        Trel = linterp1(timevec, Tm .- Tm[1], TrelAve_time)
        return sum((Trel-TrelObs).^2 ./ (2 .* TrelObs_sigma.^2))
    end

## --- Model functions (reverse time)

    # Finite-difference cooling model
    function coolingmodel(p::Parameters, Qrm, dt, timevec)

        # Allocate arrays
        Qc = Array{Float64}(undef, size(timevec))
        Qm = Array{Float64}(undef, size(timevec))
        dTm_dt = Array{Float64}(undef, size(timevec))
        dTc_dt = Array{Float64}(undef, size(timevec))
        mu_u = Array{Float64}(undef, size(timevec))
        mu_l = Array{Float64}(undef, size(timevec))
        Tm = Array{Float64}(undef, size(timevec))
        Tc = Array{Float64}(undef, size(timevec))
        dp = Array{Float64}(undef, size(timevec))

        # Consolidated core constants
        a_c = 4pi * p.Rc^2 * p.Kc * ((2 * p.g * p.rho_c * p.alpha) / (3 * p.kappa * p.Rp^2))^p.nc

        # ---  Variable parameters to adjust --- #
        Qm[1] = p.Qm_now # Present mantle heat flux (W)
        Qc[1] = p.Qc_now # Present core heat flux (W)
        Tm[1] = p.Tm_now # Present mantle temperature (K)
        Tc[1] = p.Tc_now # Present core temperature (K)

        mu_u[1] = p.mur_u * exp(p.Ea/p.R*(1/Tm[1]-1/p.Tru)) # Inital mantle viscosity
        mu_l[1] = p.mur_l * exp(p.Ha/p.R*(1/Tc[1]-1/p.Trl)) # Initial lower mantle / CMB viscosity
        epsilon = Tm[1]^2 / (p.Ha/p.R*(Tc[1]-Tm[1])) # Initial CMB boundary layer thickness ratio

        # Intial lithospheric plate thickness (piecewise fit to Korenaga)
        dp[1] = 115E3 + 275*(Tm[1]-(1400+273.15)) +
            (Tm[1]<(1410+273.15)) *
            (0.003649*(Tm[1]-273.15).^2 - 10.4*(Tm[1]-273.15) + 7409.5)

        # Adjustment to fit modern (imposed) core and mantle heat fluxes
        Beta_c = Qc[1] /
            (a_c * epsilon^(4*p.nc) * ((Tc[1]-Tm[1])^(1+p.nc) / mu_l[1]^p.nc))
        Beta_d = Qm[1] /
            (4*pi*p.Re^2*p.Km*Tm[1]*
                (p.g*p.rho_m*p.alpha*Tm[1]*dp[1] /
                    (4*mu_u[1]*p.kappa + (4/3*p.eta_L*p.kappa)*(dp[1]/p.Rc)^3)
                )^(1/2)
            )

        # Initial core and mantle cooling rates
        dTc_dt[1] = (p.Mc.*0 - Qc[1]) /
            (p.Xc * p.Mc * p.Cp_c) * p.s_yr*1E6
        dTm_dt[1] = (Qrm[1] + Qc[1] - Qm[1]) /
            (p.Xm * p.Mm * p.Cp_m) * p.s_yr * 1E6

        for i=2:length(timevec)
            Tm[i] = Tm[i-1] - dTm_dt[i-1]*dt # New mantle temperature
            Tc[i] = Tc[i-1] - dTc_dt[i-1]*dt # New core temperature

            mu_u[i] = p.mur_u * exp(p.Ea/p.R*(1/Tm[i]-1/p.Tru)) # New mantle viscosity
            mu_l[i] = p.mur_l * exp(p.Ha/p.R*(1/Tc[i]-1/p.Trl)) # New lower mantle / CMB viscosity
            epsilon = Tm[i]^2 / (p.Ha/p.R*(Tc[i]-Tm[i])) # New CMB boundary layer thickness ratio

            # Lithospheric plate thickness (piecewise fit to Korenaga)
            dp[i] = 115E3 + 275*(Tm[i]-(1400+273.15)) +
                (Tm[i]<(1410+273.15)) *
                (0.003649*(Tm[i]-273.15).^2 - 10.4*(Tm[i]-273.15) + 7409.5)

            # Core heat flux
            Qc[i] = Beta_c * a_c * abs(epsilon)^(4*p.nc) *
                abs(Tc[i]-Tm[i])^(1+p.nc) / abs(mu_l[i])^p.nc

            # Mantle heat flux
            Qm[i] = Beta_d*4*pi*p.Re^2*p.Km*Tm[i] *
                sqrt(abs(p.g*p.rho_m*p.alpha*Tm[i]*dp[i]/(4*mu_u[i]*p.kappa +
                    (4/3*p.eta_L*p.kappa)*(dp[i]/p.Rc)^3)
                ))

            # Core and mantle cooling rates by energy balance
            dTc_dt[i] = (p.Mc.*0 - Qc[i]) /
                (p.Xc * p.Mc * p.Cp_c) * p.s_yr * 1E6
            dTm_dt[i] = (Qrm[i] + Qc[i] + p.Q_addtl - Qm[i]) /
                (p.Xm * p.Mm * p.Cp_m) * p.s_yr * 1E6
        end

        return (Tm, Tc, Qm, Qc, dp)
    end

    # Run model and calculate its signed log likelihood
    function coolingmodel_LL(p::Parameters, Qrm, dt, timevec, TrelAve_time, TrelAve, TrelAve_sigma)
        Tm = coolingmodel(p, Qrm, dt, timevec)[1]
        return signed_LL(Tm, timevec, TrelAve_time, TrelAve, TrelAve_sigma)
    end

    # Run model and calculate its residual
    function coolingmodel_resid(p::Parameters, Qrm, dt, timevec, TrelAve_time, TrelAve, TrelAve_sigma)
        Tm = coolingmodel(p, Qrm, dt, timevec)[1]
        return resid(Tm, timevec, TrelAve_time, TrelAve, TrelAve_sigma)
    end

    # Find best-fit present-day Qrm with newton's method
    function coolingmodel_Newton_Qrm(p::Parameters, nSteps, Qrm, dt, timevec, TrelObs_time, TrelObs, TrelObs_sigma)
        Qrm_now = p.Qrm_now
        dQrm_now = p.Qrm_now*1E-9
        for i=1:nSteps-1
            ll = coolingmodel_LL(p, Qrm .* (Qrm_now / Qrm[1]), dt, timevec, TrelObs_time, TrelObs, TrelObs_sigma)
            llp = coolingmodel_LL(p, Qrm .* ((Qrm_now+dQrm_now) / Qrm[1]), dt, timevec, TrelObs_time, TrelObs, TrelObs_sigma)
            dll_dQrm = (llp-ll) / dQrm_now
            Qrm_now = Qrm_now - ll/dll_dQrm/2
        end
        return Qrm_now > 0 ? Qrm_now : 0.0
    end

## --- Model functions (forward time)

    # Finite-difference cooling model
    function coolingmodel_forward(p::Parameters, Qrm, dt, timevec)

        # Allocate arrays
        Qc = Array{Float64}(undef, size(timevec))
        Qm = Array{Float64}(undef, size(timevec))
        dTm_dt = Array{Float64}(undef, size(timevec))
        dTc_dt = Array{Float64}(undef, size(timevec))
        mu_u = Array{Float64}(undef, size(timevec))
        mu_l = Array{Float64}(undef, size(timevec))
        Tm = Array{Float64}(undef, size(timevec))
        Tc = Array{Float64}(undef, size(timevec))
        dp = Array{Float64}(undef, size(timevec))

        # Consolidated core constants
        a_c = 4pi * p.Rc^2 * p.Kc * ((2 * p.g * p.rho_c * p.alpha) / (3 * p.kappa * p.Rp^2))^p.nc

        # ---  Variable parameters to adjust --- #
        Qm[end] = p.Qm_initial # Present mantle heat flux (W)
        Qc[end] = p.Qc_initial # Present core heat flux (W)
        Tm[end] = p.Tm_initial # Present mantle temperature (K)
        Tc[end] = p.Tc_initial # Present core temperature (K)

        mu_u[end] = p.mur_u * exp(p.Ea/p.R*(1/Tm[end]-1/p.Tru)) # Inital mantle viscosity
        mu_l[end] = p.mur_l * exp(p.Ha/p.R*(1/Tc[end]-1/p.Trl)) # Initial lower mantle / CMB viscosity
        epsilon = Tm[end]^2 / (p.Ha/p.R*(Tc[end]-Tm[end])) # Initial CMB boundary layer thickness ratio

        # Intial lithospheric plate thickness (piecewise fit to Korenaga)
        dp[end] = 115E3 + 275*(Tm[end]-(1400+273.15)) +
            (Tm[end]<(1410+273.15)) *
            (0.003649*(Tm[end]-273.15).^2 - 10.4*(Tm[end]-273.15) + 7409.5)

        # Adjustment to reconcile modern core and mantle heat fluxes
        dp_now = 115E3 + 275*(p.Tm_now-(1400+273.15)) +
            (p.Tm_now<(1410+273.15)) *
            (0.003649*(p.Tm_now-273.15).^2 - 10.4*(p.Tm_now-273.15) + 7409.5)
            # Present-day lithospheric thickness (piecewise fit to Korenaga)
        mu_u_now = p.mur_u * exp(p.Ea/p.R*(1/p.Tm_now-1/p.Tru)) # Present-day mantle viscosity
        mu_l_now = p.mur_l * exp(p.Ha/p.R*(1/p.Tc_now-1/p.Trl)) # Present-day lower mantle / CMB viscosity
        epsilon_now = p.Tm_now^2 / (p.Ha/p.R*(p.Tc_now-p.Tm_now)) # Present-day CMB boundary layer thickness ratio
        Beta_c = p.Qc_now /
            (a_c * epsilon_now^(4*p.nc) * ((p.Tc_now-p.Tm_now)^(1+p.nc) / mu_l_now^p.nc))
            # Core calibration factor
        Beta_d = p.Qm_now /
            (4*pi*p.Re^2*p.Km*p.Tm_now*
                (p.g*p.rho_m*p.alpha*p.Tm_now*dp_now /
                    (4*mu_u_now*p.kappa + (4/3*p.eta_L*p.kappa)*(dp_now/p.Rc)^3)
                )^(1/2)
            ) # Mantle calibration factor

        # Beta_c = 2.670859522430628
        # Beta_d = 0.08052045415991978

        # Initial core and mantle cooling rates
        dTc_dt[end] = (p.Mc.*0 - Qc[end]) /
            (p.Xc * p.Mc * p.Cp_c) * p.s_yr*1E6
        dTm_dt[end] = (Qrm[end] + Qc[end] - Qm[end]) /
            (p.Xm * p.Mm * p.Cp_m) * p.s_yr * 1E6

        for i=(length(timevec)-1):-1:1
            Tm[i] = Tm[i+1] + dTm_dt[i+1]*dt # New mantle temperature
            Tc[i] = Tc[i+1] + dTc_dt[i+1]*dt # New core temperature

            mu_u[i] = p.mur_u * exp(p.Ea/p.R*(1/Tm[i]-1/p.Tru)) # New mantle viscosity
            mu_l[i] = p.mur_l * exp(p.Ha/p.R*(1/Tc[i]-1/p.Trl)) # New lower mantle / CMB viscosity
            epsilon = Tm[i]^2 / (p.Ha/p.R*(Tc[i]-Tm[i])) # New CMB boundary layer thickness ratio

            # Lithospheric plate thickness (piecewise fit to Korenaga)
            dp[i] = 115E3 + 275*(Tm[i]-(1400+273.15)) +
                (Tm[i]<(1410+273.15)) *
                (0.003649*(Tm[i]-273.15).^2 - 10.4*(Tm[i]-273.15) + 7409.5)

            # Core heat flux
            Qc[i] = Beta_c * a_c * abs(epsilon)^(4*p.nc) *
                abs(Tc[i]-Tm[i])^(1+p.nc) / abs(mu_l[i])^p.nc

            # Mantle heat flux
            Qm[i] = Beta_d*4*pi*p.Re^2*p.Km*Tm[i] *
                sqrt(abs(p.g*p.rho_m*p.alpha*Tm[i]*dp[i]/(4*mu_u[i]*p.kappa +
                    (4/3*p.eta_L*p.kappa)*(dp[i]/p.Rc)^3)
                ))

            # Core and mantle cooling rates by energy balance
            dTc_dt[i] = (p.Mc.*0 - Qc[i]) /
                (p.Xc * p.Mc * p.Cp_c) * p.s_yr * 1E6
            dTm_dt[i] = (Qrm[i] + Qc[i] + p.Q_addtl - Qm[i]) /
                (p.Xm * p.Mm * p.Cp_m) * p.s_yr * 1E6
        end

        return (Tm, Tc, Qm, Qc, dp)
    end

    # Run model and calculate its signed log likelihood
    function coolingmodel_forward_LL(p::Parameters, Qrm, dt, timevec, TrelAve_time, TrelAve, TrelAve_sigma)
        Tm = coolingmodel_forward(p, Qrm, dt, timevec)[1]
        return signed_LL(Tm, timevec, TrelAve_time, TrelAve, TrelAve_sigma)
    end

    # Run model and calculate its residual
    function coolingmodel_forward_resid(p::Parameters, Qrm, dt, timevec, TrelAve_time, TrelAve, TrelAve_sigma)
        Tm = coolingmodel_forward(p, Qrm, dt, timevec)[1]
        return resid(Tm, timevec, TrelAve_time, TrelAve, TrelAve_sigma)
    end

    # Find best-fit present-day Qrm with newton's method
    function coolingmodel_forward_Newton_Qrm(p::Parameters, nSteps, Qrm, dt, timevec, TrelObs_time, TrelObs, TrelObs_sigma)
        Qrm_now = p.Qrm_now
        dQrm_now = p.Qrm_now*1E-9
        for i=1:nSteps-1
            ll = coolingmodel_forward_LL(p, Qrm .* (Qrm_now / Qrm[1]), dt, timevec, TrelObs_time, TrelObs, TrelObs_sigma)
            llp = coolingmodel_forward_LL(p, Qrm .* ((Qrm_now+dQrm_now) / Qrm[1]), dt, timevec, TrelObs_time, TrelObs, TrelObs_sigma)
            dll_dQrm = (llp-ll) / dQrm_now
            Qrm_now = Qrm_now - ll/dll_dQrm/2
        end
        return Qrm_now > 0 ? Qrm_now : 0.0
    end

## --- End of File

max_stars_repo_name: brenhinkeller/SecularCooling

max_issues_repo_name: brenhinkeller/SecularCooling

max_forks_repo_name: brenhinkeller/SecularCooling


text:
using Test
using StackOverflow

qHolder = getrecentquestionsfortag()
@test_nowarn qHolder[1].answer_count

@test_nowarn getrecentbadgerecipients()
@test_nowarn geterrors()

siteinfoholder = getsiteinfo()
@test_nowarn dump(siteinfoholder)

@test_nowarn searchtag("julia")

@test_nowarn StackOverflow.searcherror("404")

aHolder = getanswers()
@test typeof(aHolder[1].question_id) == Int

@test_nowarn StackOverflow.getanswerfromquestion(qHolder[1])
@test_nowarn StackOverflow.getquestionfromanswer(aHolder[1])

@test_nowarn getusers()
@test_nowarn getcomments()
@test_nowarn getsuggested_edits()

@test_nowarn StackOverflow.getmoderators()
@test_nowarn StackOverflow.getopenquestions()
@test_nowarn StackOverflow.getnoanswerquestions()
@test_nowarn StackOverflow.getfeaturedquestions()


text:
module TestDivisionFilters

    using Test
    using NCBITaxonomy

    @test typeof(taxon(virusfilter(), "Ebolavirus")) <: NCBITaxon
    @test typeof(taxon(plantfilter(), "Acer")) <: NCBITaxon

    @test_throws NameHasNoDirectMatch taxon(mammalfilter(false), "Homo")
    @test_throws NameHasNoDirectMatch taxon(mammalfilter(false), "Mus musculus")
    @test_throws NameHasNoDirectMatch taxon(vertebratefilter(false), "Homo")

    @test typeof(taxon(vertebratefilter(true), "Homo")) <: NCBITaxon
    @test typeof(taxon(primatefilter(), "Homo")) <: NCBITaxon
    @test typeof(taxon(vertebratefilter(true), "Homo")) <: NCBITaxon
    @test typeof(taxon(bacteriafilter(), "Pseudomonas")) <: NCBITaxon
    @test typeof(taxon(invertebratefilter(), "Lamellodiscus")) <: NCBITaxon
    @test typeof(taxon(mammalfilter(true), "Mus musculus")) <: NCBITaxon
    @test typeof(taxon(rodentfilter(), "Mus musculus")) <: NCBITaxon

end


text:
using Xpress, Compat.Test

m = Xpress.Model()

@testset "Attributes" begin
    @test Xpress.get_intattr(m, Xpress.XPRS_COLS) == 0
    @test Xpress.get_dblattr(m, Xpress.XPRS_OBJRHS) == 0.0
    @test Xpress.get_strattr(m, Xpress.XPRS_MATRIXNAME) == ""
end


text:
"""
    SubgridScale

Sub-grid scale parameters, shared by multiple models
including the atmosphere, ocean and surface fluxes.
"""
module SubgridScale

export von_karman_const

""" Von Karman constant (dimensionless) """
function von_karman_const end

end # module SubgridScale

max_stars_repo_name: climate-machine/CLIMAParameters

max_issues_repo_name: climate-machine/CLIMAParameters

max_forks_repo_name: climate-machine/CLIMAParameters


text:
"""
    unshortlink(url)

unshorten a short url with the service of http://checkshorturl.com,
 and return the url if it is not a short link
"""
function unshortlink(url)
    cmd = `curl "http://checkshorturl.com/expand.php?u=$url"`
    @show cmd
    html = read(cmd, String)
    m = match(r".*<td .*Long URL<\/td>\s*<td .*><a href=\"([^ ]*)\"", html)
    m === nothing ? url : m.captures[1]
end

isgooglesheet(url) = occursin("docs.google.com/spreadsheets", url)
isgoogledrive(url) = occursin("drive.google.com", url)

function googlesheet_handler(url; format=:csv)
    link, expo = splitdir(url)
    if startswith(expo, "edit") || expo == ""
        url = link * "/export?format=$format"
    elseif startswith(expo, "export")
        url = replace(url, r"format=([a-zA-Z]*)(.*)"=>SubstitutionString("format=$format\\2"))
    end
    url
end

function mybegoogle_download(url, localdir)
    long_url = unshortlink(url)
    if isgooglesheet(long_url)
        long_url = googlesheet_handler(long_url)
    end

    if isgoogledrive(url)
        download_gdrive(long_url, localdir)
    else
        DataDeps.fetch_http(long_url, localdir)
    end
end

function find_gcode(ckj)
    for cookie ∈ ckj
        if match(r"_warning_", cookie.name) !== nothing
            return cookie.value
        end
    end

    nothing
end

function download_gdrive(url, localdir)
    rq = HTTP.request("GET", url; cookies=true)
    ckj = HTTP.CookieRequest.default_cookiejar["drive.google.com"]
    gcode = find_gcode(ckj)
    @assert gcode !== nothing

    rq = HTTP.request("GET", "$url&confirm=$gcode"; cookies=true)
    hcd = HTTP.header(rq, "Content-Disposition")
    m = match(r"filename=\\\"(.*)\\\"", hcd)
    if m === nothing
        filename = "gdrive_downloaded"
    else
        filename = m.captures[]
    end

    filepath = joinpath(localdir, filename)

    open(filepath, "w+") do f
        write(f, rq.body)
    end

    filepath
end


text:
## support routines

function clopts(opts...)
    args = ["", opts...]
    LLVM.API.LLVMParseCommandLineOptions(Int32(length(args)),
        [Base.unsafe_convert(Cstring, arg) for arg in args], C_NULL)
end