text:
# # Using MLJ to classifiy the MNIST image dataset

using Pkg
const DIR = @__DIR__
Pkg.activate(DIR)
Pkg.instantiate()

# **Julia version** is assumed to be ^1.6

using MLJ
using Flux
import MLJFlux
import MLJIteration # for `skip`

MLJ.color_off()

using Plots
pyplot(size=(600, 300*(sqrt(5)-1)));

# ## Basic training

# Downloading the MNIST image dataset:

import MLDatasets: MNIST

ENV["DATADEPS_ALWAYS_ACCEPT"] = true
images, labels = MNIST.traindata();

# In MLJ, integers cannot be used for encoding categorical data, so we
# must force the labels to have the `Multiclass` [scientific
# type](https://juliaai.github.io/ScientificTypes.jl/dev/). For
# more on this, see [Working with Categorical
# Data](https://alan-turing-institute.github.io/MLJ.jl/dev/working_with_categorical_data/).

labels = coerce(labels, Multiclass);
images = coerce(images, GrayImage);

# Checking scientific types:

@assert scitype(images) <: AbstractVector{<:Image}
@assert scitype(labels) <: AbstractVector{<:Finite}

# Looks good.

# For general instructions on coercing image data, see [Type coercion
# for image
# data](https://alan-turing-institute.github.io/ScientificTypes.jl/dev/#Type-coercion-for-image-data-1)

images[1]

# We start by defining a suitable `Builder` object. This is a recipe
# for building the neural network. Our builder will work for images of
# any (constant) size, whether they be color or black and white (ie,
# single or multi-channel).  The architecture always consists of six
# alternating convolution and max-pool layers, and a final dense
# layer; the filter size and the number of channels after each
# convolution layer is customisable.

import MLJFlux
struct MyConvBuilder
    filter_size::Int
    channels1::Int
    channels2::Int
    channels3::Int
end

make2d(x::AbstractArray) = reshape(x, :, size(x)[end])

function MLJFlux.build(b::MyConvBuilder, rng, n_in, n_out, n_channels)
    k, c1, c2, c3 = b.filter_size, b.channels1, b.channels2, b.channels3
    mod(k, 2) == 1 || error("`filter_size` must be odd. ")
    p = div(k - 1, 2) # padding to preserve image size
    init = Flux.glorot_uniform(rng)
    front = Chain(
        Conv((k, k), n_channels => c1, pad=(p, p), relu, init=init),
        MaxPool((2, 2)),
        Conv((k, k), c1 => c2, pad=(p, p), relu, init=init),
        MaxPool((2, 2)),
        Conv((k, k), c2 => c3, pad=(p, p), relu, init=init),
        MaxPool((2 ,2)),
        make2d)
    d = Flux.outputsize(front, (n_in..., n_channels, 1)) |> first
    return Chain(front, Dense(d, n_out, init=init))
end

# **Note.** There is no final `softmax` here, as this is applied by
# default in all MLJFLux classifiers. Customisation of this behaviour
# is controlled using using the `finaliser` hyperparameter of the
# classifier.

# We now define the MLJ model. If you have a GPU, substitute
# `acceleration=CUDALibs()` below:

ImageClassifier = @load ImageClassifier
clf = ImageClassifier(builder=MyConvBuilder(3, 16, 32, 32),
                      batch_size=50,
                      epochs=10,
                      rng=123)

# You can add Flux options `optimiser=...` and `loss=...` here. At
# present, `loss` must be a Flux-compatible loss, not an MLJ
# measure. To run on a GPU, set `acceleration=CUDALib()`.

# Binding the model with data in an MLJ machine:
mach = machine(clf, images, labels);

# Training for 10 epochs on the first 500 images:

fit!(mach, rows=1:500, verbosity=2);

# Inspecting:

report(mach)

#-

chain = fitted_params(mach)

#-

Flux.params(chain)[2]

#-

# Adding 20 more epochs:

clf.epochs = clf.epochs + 20
fit!(mach, rows=1:500);

# Computing an out-of-sample estimate of the loss:

predicted_labels = predict(mach, rows=501:1000);
cross_entropy(predicted_labels, labels[501:1000]) |> mean

# Or, in one line:

evaluate!(mach,
          resampling=Holdout(fraction_train=0.5),
          measure=cross_entropy,
          rows=1:1000,
          verbosity=0)

# ## Wrapping the MLJFlux model with iteration controls

# Any iterative MLJFlux model can be wrapped in *iteration controls*,
# as we demonstrate next. For more on MLJ's `IteratedModel` wrapper,
# see the [MLJ
# documentation](https://alan-turing-institute.github.io/MLJ.jl/dev/controlling_iterative_models/).

# The "self-iterating" classifier, called `iterated_clf` below, is for
# iterating the image classifier defined above until one of the
# following stopping criterion apply:

# - `Patience(3)`: 3 consecutive increases in the loss
# - `InvalidValue()`: an out-of-sample loss, or a training loss, is `NaN`, `Inf`, or `-Inf`
# - `TimeLimit(t=5/60)`: training time has exceeded 5 minutes

# These checks (and other controls) will be applied every two epochs
# (because of the `Step(2)` control). Additionally, training a
# machine bound to `iterated_clf` will:

# - save a snapshot of the machine every three control cycles (every six epochs)
# - record traces of the out-of-sample loss and training losses for plotting
# - record mean value traces of each Flux parameter for plotting

# For a complete list of controls, see [this
# table](https://alan-turing-institute.github.io/MLJ.jl/dev/controlling_iterative_models/#Controls-provided).

# ### Wrapping the classifier

# Some helpers

make2d(x::AbstractArray) = reshape(x, :, size(x)[end])
make1d(x::AbstractArray) = reshape(x, length(x));

# To extract Flux params from an MLJFlux machine

parameters(mach) = make1d.(Flux.params(fitted_params(mach)));

# To store the traces:

losses = []
training_losses = []
parameter_means = Float32[];
epochs = []

# To update the traces:

update_loss(loss) = push!(losses, loss)
update_training_loss(losses) = push!(training_losses, losses[end])
update_means(mach) = append!(parameter_means, mean.(parameters(mach)));
update_epochs(epoch) = push!(epochs, epoch)

# The controls to apply:

save_control =
    MLJIteration.skip(Save(joinpath(DIR, "mnist.jlso")), predicate=3)

controls=[Step(2),
          Patience(3),
          InvalidValue(),
          TimeLimit(5/60),
          save_control,
          WithLossDo(),
          WithLossDo(update_loss),
          WithTrainingLossesDo(update_training_loss),
          Callback(update_means),
          WithIterationsDo(update_epochs)
];

# The "self-iterating" classifier:

iterated_clf = IteratedModel(model=clf,
                       controls=controls,
                       resampling=Holdout(fraction_train=0.7),
                       measure=log_loss)

# ### Binding the wrapped model to data:

mach = machine(iterated_clf, images, labels);

# ### Training

fit!(mach, rows=1:500);

# ### Comparison of the training and out-of-sample losses:

plot(epochs, losses,
     xlab = "epoch",
     ylab = "root squared error",
     label="out-of-sample")
plot!(epochs, training_losses, label="training")

savefig(joinpath(DIR, "loss.png"))

# ### Evolution of weights

n_epochs =  length(losses)
n_parameters = div(length(parameter_means), n_epochs)
parameter_means2 = reshape(copy(parameter_means), n_parameters, n_epochs)'
plot(epochs, parameter_means2,
     title="Flux parameter mean weights",
     xlab = "epoch")

# **Note.** The the higher the number, the deeper the chain parameter.

savefig(joinpath(DIR, "weights.png"))

# ### Retrieving a snapshot for a prediction:

mach2 = machine(joinpath(DIR, "mnist3.jlso"))
predict_mode(mach2, images[501:503])

# ### Restarting training

# Mutating `iterated_clf.controls` or `clf.epochs` (which is otherwise
# ignored) will allow you to restart training from where it left off.

iterated_clf.controls[2] = Patience(4)
fit!(mach, rows=1:500)

plot(epochs, losses,
     xlab = "epoch",
     ylab = "root squared error",
     label="out-of-sample")
plot!(epochs, training_losses, label="training")

using Literate #src
Literate.markdown(@__FILE__, @__DIR__, execute=false) #src
Literate.notebook(@__FILE__, @__DIR__, execute=false) #src


text:
module BenchUniformLoops

using BenchmarkTools
using SyncBarriers

function work()
    t1 = time_ns() + 100
    while t1 > time_ns()
    end
end

#=
const SINKS = Float64[]
const SINKS_STRIDE = 32
function __init__()
    resize!(SINKS, SINKS_STRIDE * (Threads.nthreads() + 1))
    fill!(SINKS, -1)
end

function work(n = 10)
    i = SINKS_STRIDE * Threads.threadid()
    x = SINKS[i]
    for _ in 1:n
        x = sin(3 * x)
    end
    SINKS[i] = x
end
=#

function spawn_and_wait(ncycles, ntasks)
    tasks = Vector{Task}(undef, ntasks)
    for _ in 1:ncycles
        empty!(tasks)
        for _ in 1:ntasks
            t = Threads.@spawn work()
            push!(tasks, t)
        end
        foreach(wait, tasks)
    end
end

# bad idea?
function spawn_and_barrier(ncycles, ntasks)
    tasks = Vector{Task}(undef, ntasks)
    barrier = SyncBarriers.CentralizedBarrier(ntasks + 1)
    for k in 1:ncycles
        s = mod(k, 2^7) == 0
        s && empty!(tasks)
        for i in 1:ntasks
            t = Threads.@spawn begin
                work()
                SyncBarriers.cycle!(barrier[i])
            end
            s && push!(tasks, t)
        end
        SyncBarriers.cycle!(barrier[ntasks+1])
        s && foreach(wait, tasks)
    end
end

# => works well with dissemination
function barrier_with_static(factory, ncycles, ntasks, spin)
    barrier = factory(ntasks)
    Threads.@threads :static for i in 1:ntasks
        for k in 1:ncycles
            work()
            SyncBarriers.cycle!(barrier[i], spin)
        end
    end
end

# => works well with centralized
function barrier_with_spawn(factory, ncycles, ntasks, spin)
    barrier = factory(ntasks)
    tasks = empty!(Vector{Task}(undef, ntasks))
    for i in 1:ntasks
        t = Threads.@spawn begin
            for k in 1:ncycles
                work()
                SyncBarriers.cycle!(barrier[i], spin)
            end
        end
        push!(tasks, t)
    end
    foreach(wait, tasks)
end

function setup(;
    ncycles::Integer = 2^10,
    ntasks::Integer = Threads.nthreads(),
    spin::Integer = 10_000,
    nbranches_list::AbstractVector{<:Integer} = Int[],
)
    @debug "BenchUniformLoops: ncycles=$ncycles ntasks=$ntasks spin=$spin"

    barrier_list = []
    push_barrier!(name, value) = push!(barrier_list, (name = name, value = value))
    push_barrier!("dissemination", SyncBarriers.DisseminationBarrier),
    push_barrier!("centralized", SyncBarriers.CentralizedBarrier),
    if isempty(nbranches_list)
        if ntasks > 8
            push_barrier!("tree-8", SyncBarriers.TreeBarrier{8})
            push_barrier!("flat-tree-8", SyncBarriers.FlatTreeBarrier{8})
            push_barrier!("static-tree-8-8", SyncBarriers.StaticTreeBarrier{8,8})
        elseif ntasks > 4
            push_barrier!("tree-4", SyncBarriers.TreeBarrier{4})
            push_barrier!("flat-tree-4", SyncBarriers.FlatTreeBarrier{4})
            push_barrier!("static-tree-4-4", SyncBarriers.StaticTreeBarrier{4,4})
        elseif ntasks > 2
            push_barrier!("tree-2", SyncBarriers.TreeBarrier{2})
            push_barrier!("flat-tree-2", SyncBarriers.FlatTreeBarrier{2})
            push_barrier!("static-tree-2-2", SyncBarriers.StaticTreeBarrier{2,2})
        end
    else
        for n in nbranches_list
            push_barrier!("tree-$n", SyncBarriers.TreeBarrier{n})
            push_barrier!("flat-tree-$n", SyncBarriers.FlatTreeBarrier{n})
            push_barrier!("static-tree-$n-$n", SyncBarriers.StaticTreeBarrier{n,n})
        end
    end

    spin_list = [
        (name = "nospin", value = nothing),
        (name = "spin", value = spin),
        # ...
    ]
    loop_list = [
        (name = "static", value = barrier_with_static),
        (name = "spawn", value = barrier_with_spawn),
        # ...
    ]

    suite = BenchmarkGroup()
    suite["wait"] = @benchmarkable spawn_and_wait($ncycles, $ntasks)
    # suite["spawn_and_barrier"] = @benchmarkable spawn_and_barrier($ncycles, $ntasks)
    for barrier in barrier_list
        s1 = suite[string(barrier.name)] = BenchmarkGroup()
        for spinarg in spin_list
            s2 = s1[string(spinarg.name)] = BenchmarkGroup()
            for loop in loop_list
                s2[string(loop.name)] = @benchmarkable(
                    loop($(barrier.value), $ncycles, $ntasks, $(spinarg.value)),
                    setup = begin
                        # Workaround the error " function argument and static parameter
                        # names must be distinct":
                        loop = $(loop.value)
                    end
                )
            end
        end
    end
    return suite
end

function clear() end

end  # module

ext: jl


text:
# Autogenerated wrapper script for SCS_GPU_jll for x86_64-linux-gnu
export libscsgpuindir

using OpenBLAS32_jll
JLLWrappers.@generate_wrapper_header("SCS_GPU")
JLLWrappers.@declare_library_product(libscsgpuindir, "libscsgpuindir.so")
function __init__()
    JLLWrappers.@generate_init_header(OpenBLAS32_jll)
    JLLWrappers.@init_library_product(
        libscsgpuindir,
        "lib/libscsgpuindir.so",
        nothing,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using TermWin

tshow( TermWin, title="TermWin" )


text:
#=
In England the currency is made up of pound, £, and pence, p, and there are eight coins in general circulation:

1p, 2p, 5p, 10p, 20p, 50p, £1 (100p) and £2 (200p).
It is possible to make £2 in the following way:

1×£1 + 1×50p + 2×20p + 1×5p + 1×2p + 3×1p
How many different ways can £2 be made using any number of coins?
=#
coins = [1,2,5,10,20,50,100,200]
# Number of ways of making change for n cents
# using only the first k types of coins
function coin(n,k)
  if k < 1 || n < 0
    return 0
  elseif n == 0
    return 1
  end
  coin(n, k-1)+coin(n-coins[k], k)
end
@time println(coin(200,8))

max_stars_repo_name: daniel-beard/JuliaProjectEuler

max_issues_repo_name: daniel-beard/JuliaProjectEuler

max_forks_repo_name: daniel-beard/JuliaProjectEuler


text:
# This file contains the core functions of WignerSymbols and CG-coefficient.

"""
    HalfInt = Union{Integer, Rational}
Angular momentum quantum numbers may be half integers and integers. With `HalfInt` type, you can use both `2` and `3//2` as parameters.
But the parameter like `5//3`, who's denominator is not `2` while gives out error.
"""
const HalfInt = Union{Integer,Rational}

"""
basic CG coefficient calculation function
"""
function _dCG(dj1::BigInt, dj2::BigInt, dj3::BigInt, dm1::BigInt, dm2::BigInt, dm3::BigInt)
    check_jm(dj1, dm1) || error(jm_mismatching_msg(dj1, dm1))
    check_jm(dj2, dm2) || error(jm_mismatching_msg(dj2, dm2))
    check_jm(dj3, dm3) || error(jm_mismatching_msg(dj3, dm3))
    check_couple(dj1, dj2, dj3) || error(miss_couple_msg(dj1, dj2, dj3))
    dm1 + dm2 == dm3 || return zero(SqrtRational)
    J::BigInt = div(dj1 + dj2 + dj3, 2)
    Jm1::BigInt = J - dj1
    Jm2::BigInt = J - dj2
    Jm3::BigInt = J - dj3
    j1mm1::BigInt = div(dj1 - dm1, 2)
    j2mm2::BigInt = div(dj2 - dm2, 2)
    j3mm3::BigInt = div(dj3 - dm3, 2)
    j2pm2::BigInt = div(dj2 + dm2, 2)
    A = SqrtRational(
        binomial(dj1, Jm2) * binomial(dj2, Jm3) // (
            binomial(J + 1, Jm3) * binomial(dj1, j1mm1) * binomial(dj2, j2mm2) * binomial(dj3, j3mm3)
        )
    )
    B::BigInt = zero(BigInt)
    low::BigInt = max(zero(BigInt), j1mm1 - Jm2, j2pm2 - Jm1)
    high::BigInt = min(Jm3, j1mm1, j2pm2)
    for z in low:high
        B = -B + binomial(Jm3, z) * binomial(Jm2, j1mm1 - z) * binomial(Jm1, j2pm2 - z)
    end
    return A * (iphase(high) * B)
end

"""
use CG coefficient to calculate 3j-symbol
"""
function _d3j(dj1::BigInt, dj2::BigInt, dj3::BigInt, dm1::BigInt, dm2::BigInt, dm3::BigInt)
    iphase(dj1 + div(dj3 + dm3, 2)) * SqrtRational(1 // (dj3 + 1)) * _dCG(dj1, dj2, dj3, -dm1, -dm2, dm3)
end

"""
basic 6j-symbol calculation funciton
"""
function _d6j(dj1::BigInt, dj2::BigInt, dj3::BigInt, dj4::BigInt, dj5::BigInt, dj6::BigInt)
    check_couple(dj1, dj2, dj3) || error(miss_couple_msg(dj1, dj2, dj3))
    check_couple(dj1, dj5, dj6) || error(miss_couple_msg(dj1, dj5, dj6))
    check_couple(dj4, dj2, dj6) || error(miss_couple_msg(dj4, dj2, dj6))
    check_couple(dj4, dj5, dj3) || error(miss_couple_msg(dj4, dj5, dj3))
    j123::BigInt = div(dj1 + dj2 + dj3, 2)
    j156::BigInt = div(dj1 + dj5 + dj6, 2)
    j426::BigInt = div(dj4 + dj2 + dj6, 2)
    j453::BigInt = div(dj4 + dj5 + dj3, 2)
    jpm123::BigInt = div(dj1 + dj2 - dj3, 2)
    jpm132::BigInt = div(dj1 + dj3 - dj2, 2)
    jpm231::BigInt = div(dj2 + dj3 - dj1, 2)
    jpm156::BigInt = div(dj1 + dj5 - dj6, 2)
    jpm426::BigInt = div(dj4 + dj2 - dj6, 2)
    jpm453::BigInt = div(dj4 + dj5 - dj3, 2)
    A = SqrtRational(
        binomial(j123 + 1, dj1 + 1) * binomial(dj1, jpm123) // (
            binomial(j156 + 1, dj1 + 1) * binomial(dj1, jpm156) *
            binomial(j453 + 1, dj4 + 1) * binomial(dj4, jpm453) *
            binomial(j426 + 1, dj4 + 1) * binomial(dj4, jpm426)
        )
    )
    A = A * Rational{BigInt}(1 // (dj4 + 1))
    B::BigInt = zero(BigInt)
    low::BigInt = max(j123, j453, j426, j156)
    high::BigInt = min(jpm123 + j453, jpm132 + j426, jpm231 + j156)
    for n = low:high
        B = -B + binomial(n + 1, j123 + 1) * binomial(jpm123, n - j453) *
                binomial(jpm132, n - j426) * binomial(jpm231, n - j156)
    end
    return A * (iphase(high) * B)
end

"""
basic 9j-symbol calculation funciton
"""
function _d9j(dj1::BigInt, dj2::BigInt, dj3::BigInt,
              dj4::BigInt, dj5::BigInt, dj6::BigInt,
              dj7::BigInt, dj8::BigInt, dj9::BigInt)
    check_couple(dj1, dj2, dj3) || error(miss_couple_msg(dj1, dj2, dj3))
    check_couple(dj4, dj5, dj6) || error(miss_couple_msg(dj4, dj5, dj6))
    check_couple(dj7, dj8, dj9) || error(miss_couple_msg(dj7, dj8, dj9))
    check_couple(dj1, dj4, dj7) || error(miss_couple_msg(dj1, dj4, dj7))
    check_couple(dj2, dj5, dj8) || error(miss_couple_msg(dj2, dj5, dj8))
    check_couple(dj3, dj6, dj9) || error(miss_couple_msg(dj3, dj6, dj9))
    j123::BigInt = div(dj1 + dj2 + dj3, 2)
    j456::BigInt = div(dj4 + dj5 + dj6, 2)
    j789::BigInt = div(dj7 + dj8 + dj9, 2)
    j147::BigInt = div(dj1 + dj4 + dj7, 2)
    j258::BigInt = div(dj2 + dj5 + dj8, 2)
    j369::BigInt = div(dj3 + dj6 + dj9, 2)
    pm123::BigInt = div(dj1 + dj2 - dj3, 2)
    pm132::BigInt = div(dj1 + dj3 - dj2, 2)
    pm231::BigInt = div(dj2 + dj3 - dj1, 2)
    pm456::BigInt = div(dj4 + dj5 - dj6, 2)
    pm465::BigInt = div(dj4 + dj6 - dj5, 2)
    pm564::BigInt = div(dj5 + dj6 - dj4, 2)
    pm789::BigInt = div(dj7 + dj8 - dj9, 2)
    pm798::BigInt = div(dj7 + dj9 - dj8, 2)
    pm897::BigInt = div(dj8 + dj9 - dj7, 2)
    P0_nu::BigInt = binomial(j123 + 1, dj1 + 1) * binomial(dj1, pm123) *
                    binomial(j456 + 1, dj5 + 1) * binomial(dj5, pm456) *
                    binomial(j789 + 1, dj9 + 1) * binomial(dj9, pm798)
    P0_de::BigInt = binomial(j147 + 1, dj1 + 1) * binomial(dj1, div(dj1 + dj4 - dj7, 2)) *
                    binomial(j258 + 1, dj5 + 1) * binomial(dj5, div(dj2 + dj5 - dj8, 2)) *
                    binomial(j369 + 1, dj9 + 1) * binomial(dj9, div(dj3 + dj9 - dj6, 2))
    P0 = SqrtRational(P0_nu // P0_de)
    dtl::BigInt = max(abs(dj2 - dj6), abs(dj4 - dj8), abs(dj1 - dj9))
    dth::BigInt = min(dj2 + dj6, dj4 + dj8, dj1 + dj9)
    PABC::Rational{BigInt} = zero(Rational{BigInt})
    for dt::BigInt = dtl:2:dth
        j19t::BigInt = div(dj1 + dj9 + dt, 2)
        j26t::BigInt = div(dj2 + dj6 + dt, 2)
        j48t::BigInt = div(dj4 + dj8 + dt, 2)
        Pt_de::BigInt = binomial(j19t + 1, dt + 1) * binomial(dt, div(dj1 + dt - dj9, 2)) *
                        binomial(j26t + 1, dt + 1) * binomial(dt, div(dj2 + dt - dj6, 2)) *
                        binomial(j48t + 1, dt + 1) * binomial(dt, div(dj4 + dt - dj8, 2))
        Pt_de *= (dt + 1)^2
        xl::BigInt = max(j123, j369, j26t, j19t)
        xh::BigInt = min(pm123 + j369, pm132 + j26t, pm231 + j19t)
        At::BigInt = zero(BigInt)
        for x = xl:xh
            At = -At + binomial(x + 1, j123 + 1) * binomial(pm123, x - j369) * binomial(pm132, x - j26t) * binomial(pm231, x - j19t)
        end
        yl::BigInt = max(j456, j26t, j258, j48t)
        yh::BigInt = min(pm456 + j26t, pm465 + j258, pm564 + j48t)
        Bt::BigInt = zero(BigInt)
        for y = yl:yh
            Bt = -Bt + binomial(y + 1, j456 + 1) * binomial(pm456, y - j26t) * binomial(pm465, y - j258) * binomial(pm564, y - j48t)
        end
        zl::BigInt = max(j789, j19t, j48t, j147)
        zh::BigInt = min(pm789 + j19t, pm798 + j48t, pm897 + j147)
        Ct::BigInt = zero(BigInt)
        for z = zl:zh
            Ct = -Ct + binomial(z + 1, j789 + 1) * binomial(pm789, z - j19t) * binomial(pm798, z - j48t) * binomial(pm897, z - j147)
        end
        PABC += (iphase(xh + yh + zh) * At * Bt * Ct) // Pt_de
    end
    return P0 * (iphase(dth) * PABC)
end

"""
    dCG(dj1::Integer, dj2::Integer, dj3::Integer, dm1::Integer, dm2::Integer, dm3::Integer)
CG coefficient function with double angular monentum number parameters, so that the parameters can be integer.
You can calculate `dCG(1, 1, 2, 1, 1, 2)` to calculate the real `CG(1//2, 1//2, 1, 1/2, 1//2, 1)`
"""
dCG(dj1::Integer, dj2::Integer, dj3::Integer, dm1::Integer, dm2::Integer, dm3::Integer) = _dCG(BigInt.((dj1, dj2, dj3, dm1, dm2, dm3))...)

"""
    d3j(dj1::Integer, dj2::Integer, dj3::Integer, dm1::Integer, dm2::Integer, dm3::Integer)
3j-symbol function with double angular monentum number parameters, so that the parameters can be integer.
"""
d3j(dj1::Integer, dj2::Integer, dj3::Integer, dm1::Integer, dm2::Integer, dm3::Integer) = _d3j(BigInt.((dj1, dj2, dj3, dm1, dm2, dm3))...)

"""
    d6j(dj1::Integer, dj2::Integer, dj3::Integer, dj4::Integer, dj5::Integer, dj6::Integer)
6j-symbol function with double angular monentum number parameters, so that the parameters can be integer.
"""
d6j(dj1::Integer, dj2::Integer, dj3::Integer, dj4::Integer, dj5::Integer, dj6::Integer) = _d6j(BigInt.((dj1, dj2, dj3, dj4, dj5, dj6))...)

"""
    d9j(dj1::Integer, dj2::Integer, dj3::Integer,
        dj4::Integer, dj5::Integer, dj6::Integer,
        dj7::Integer, dj8::Integer, dj9::Integer)
9j-symbol function with double angular monentum number parameters, so that the parameters can be integer.
"""
d9j(dj1::Integer, dj2::Integer, dj3::Integer,
    dj4::Integer, dj5::Integer, dj6::Integer,
    dj7::Integer, dj8::Integer, dj9::Integer) = _d9j(BigInt.((dj1, dj2, dj3, dj4, dj5, dj6, dj7, dj8, dj9))...)

@doc raw"""
    CG(j1::HalfInt, j2::HalfInt, j3::HalfInt, m1::HalfInt, m2::HalfInt, m3::HalfInt)
CG coefficient ``\langle j_1m_1 j_2m_2 | j_3m_3 \rangle``
"""
CG(j1::HalfInt, j2::HalfInt, j3::HalfInt, m1::HalfInt, m2::HalfInt, m3::HalfInt) = _dCG(BigInt.((2j1, 2j2, 2j3, 2m1, 2m2, 2m3))...)

@doc raw"""
    threeJ(j1::HalfInt, j2::HalfInt, j3::HalfInt, m1::HalfInt, m2::HalfInt, m3::HalfInt)
Wigner 3j-symbol
```math
\begin{pmatrix}
j_1 & j_2 & j_3 \\
m_1 & m_2 & m_3
\end{pmatrix}
```
"""
threeJ(j1::HalfInt, j2::HalfInt, j3::HalfInt, m1::HalfInt, m2::HalfInt, m3::HalfInt) = _d3j(BigInt.((2j1, 2j2, 2j3, 2m1, 2m2, 2m3))...)

@doc raw"""
    sixJ(j1::HalfInt, j2::HalfInt, j3::HalfInt, j4::HalfInt, j5::HalfInt, j6::HalfInt)
Wigner 6j-symbol
```math
\begin{Bmatrix}
j_1 & j_2 & j_3 \\
j_4 & j_5 & j_6
\end{Bmatrix}
```
"""
sixJ(j1::HalfInt, j2::HalfInt, j3::HalfInt, j4::HalfInt, j5::HalfInt, j6::HalfInt) = _d6j(BigInt.((2j1, 2j2, 2j3, 2j4, 2j5, 2j6))...)

@doc raw"""
    nineJ(j1::HalfInt, j2::HalfInt, j3::HalfInt,
          j4::HalfInt, j5::HalfInt, j6::HalfInt,
          j7::HalfInt, j8::HalfInt, j9::HalfInt)
Wigner 9j-symbol
```math
\begin{Bmatrix}
j_1 & j_2 & j_3 \\
j_4 & j_5 & j_6 \\
j_7 & j_8 & j_9
\end{Bmatrix}
```
"""
nineJ(j1::HalfInt, j2::HalfInt, j3::HalfInt,
      j4::HalfInt, j5::HalfInt, j6::HalfInt,
      j7::HalfInt, j8::HalfInt, j9::HalfInt) = _d9j(BigInt.((2j1, 2j2, 2j3, 2j4, 2j5, 2j6, 2j7, 2j8, 2j9))...)


text:
# ---
# title: 643. Maximum Average Subarray I
# id: problem643
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories: Array
# link: <https://leetcode.com/problems/maximum-average-subarray-i/description/>
# hidden: true
# ---
#
# Given an array consisting of `n` integers, find the contiguous subarray of
# given length `k` that has the maximum average value. And you need to output
# the maximum average value.
#
# **Example 1:**
#
#
#
#     Input: [1,12,-5,-6,50,3], k = 4
#     Output: 12.75
#     Explanation: Maximum average is (12-5-6+50)/4 = 51/4 = 12.75
#
#
#
#
# **Note:**
#
#   1. 1 <= `k` <= `n` <= 30,000.
#   2. Elements of the given array will be in the range [-10,000, 10,000].
#
#
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
cd(@__DIR__); include("setups/grid7x3.jl");

## ROT2
D = natural_eigdist(𝚽, 𝛌, Q; α = 0.1, distance = :ROT)
E = transform(fit(MDS, D, maxoutdim=2, distances=true))
E[1, :] .*= sign(E[1, 2]); E[2, :] .*= -sign(E[1, 2])
plt = grid7x3_mds_heatmaps(E, 𝚽; backend = :pyplot)
savefig(plt, "../figs/Grid7x3_MDS_ROT2_alpha01.png")
display(plt)


text:
struct DimsSelector{N}
    dims::NTuple{N, Int}
end
dims(args...) = DimsSelector(args)

function (d::DimsSelector)(c::CartesianIndex{N}) where N
    t = ntuple(N) do n
        return n in d.dims ? c[n] : 1
    end
    return CartesianIndex(t)
end

compute_label(data, name::StringLike) = string(name)
compute_label(data, name::Integer) = string(columnnames(data)[name])
compute_label(data, name::DimsSelector) = ""

struct NameTransformationLabel
    name::Any
    transformation::Any
    label::String
end

function NameTransformationLabel(name, transformation, label::Symbol)
    return NameTransformationLabel(name, transformation, string(label))
end

function NameTransformationLabel(data, x::Union{StringLike, Integer, DimsSelector})
    return NameTransformationLabel(x, identity, compute_label(data, x))
end

function NameTransformationLabel(data, x::Pair{<:Any, <:StringLike})
    name, label = x
    return NameTransformationLabel(name, identity, label)
end

function NameTransformationLabel(data, x::Pair{<:Any, <:Any})
    name, transformation = x
    label = compute_label(data, name)
    return NameTransformationLabel(name, transformation, label)
end

function NameTransformationLabel(data, x::Pair{<:Any, <:Pair})
    name, transformation_label = x
    transformation, label = transformation_label
    return NameTransformationLabel(name, transformation, label)
end

function apply_context(data, axs, names::ArrayLike)
    return map(name -> apply_context(data, axs, name), names)
end

apply_context(data, axs, name::StringLike) = getcolumn(data, Symbol(name))

function apply_context(data, axs, idx::Integer)
    name = columnnames(data)[idx]
    return getcolumn(data, name)
end

function apply_context(data, axs::NTuple{N, Any}, d::DimsSelector) where N
    sz = ntuple(N) do n
        return n in d.dims ? length(axs[n]) : 1
    end
    return reshape(CartesianIndices(sz), 1, sz...)
end

struct Labeled{T}
    label::AbstractString
    value::T
end

Labeled(x) = Labeled(getlabel(x), getvalue(x))

getlabel(x::Labeled) = x.label
getvalue(x::Labeled) = x.value

getlabel(x) = ""
getvalue(x) = x

function process_data(data, mappings′)
    mappings = map(mappings′) do x
        return map(Base.Fix1(NameTransformationLabel, data), maybewrap(x))
    end
    axs = Broadcast.combine_axes(mappings.positional..., values(mappings.named)...)
    labeledarrays = map(mappings) do ntls
        names = map(ntl -> ntl.name, ntls)
        transformations = map(ntl -> ntl.transformation, ntls)
        labels = map(ntl -> ntl.label, ntls)
        res = map(transformations, names) do transformation, name
            cols = apply_context(data, axs, maybewrap(name))
            map(transformation, cols...)
        end
        return Labeled(join(unique(labels), ' '), unnest(res))
    end
    return Entry(Any, labeledarrays, Dict{Symbol, Any}())
end

process_transformations(layers::Layers) = map(process_transformations, layers)

function process_transformations(layer::Layer)
    init = process_data(layer.data, layer.mappings)
    res = foldl(process_transformations, layer.transformations; init)
    return res isa Entry ? splitapply(res) : res
end

process_transformations(v::AbstractArray{Entry}, f) = map(f, v)

process_transformations(le::Entry, f) = f(le)


text:
"""
	contour(cmd0::String="", arg1=nothing; kwargs...)

Reads a table data and produces a raw contour plot by triangulation.

Full option list at [`contour`]($(GMTdoc)contour.html)

Parameters
----------

- $(GMT.opt_J)
- **A** | **annot** :: [Type => Str | Number]       ``Arg = [-|[+]annot_int][labelinfo]``

    *annot_int* is annotation interval in data units; it is ignored if contour levels are given in a file.
    ($(GMTdoc)contour.html#a)
- $(GMT.opt_B)
- **C** | **cont** | **contour** | **contours** | **levels** :: [Type => Str | Number | GMTcpt]  ``Arg = [+]cont_int``

    Contours to be drawn may be specified in one of three possible ways.
    ($(GMTdoc)contour.html#c)
- **D** | **dump** :: [Type => Str]

    Dump contours as data line segments; no plotting takes place.
    ($(GMTdoc)contour.html#d)
- **E** | **index** :: [Type => Str | Mx3 array]

    Give name of file with network information. Each record must contain triplets of node
    numbers for a triangle.
    ($(GMTdoc)contour.html#e)
- **G** | **labels** :: [Type => Str]

    Controls the placement of labels along the quoted lines.
    ($(GMTdoc)contour.html#g)
- **I** | **fill** | **colorize** :: [Type => Bool]

    Color the triangles using the color scale provided via **C**.
    ($(GMTdoc)contour.html#i)
- $(GMT.opt_Jz)
- **L** | **mesh** :: [Type => Str | Number]

    Draw the underlying triangular mesh using the specified pen attributes (if not provided, use default pen)
    ($(GMTdoc)contour.html#l)
- **N** | **no_clip** :: [Type => Bool]

    Do NOT clip contours or image at the boundaries [Default will clip to fit inside region].
    ($(GMTdoc)contour.html#n)
- $(GMT.opt_P)
- **Q** | **cut** :: [Type => Str | Number]         ``Arg = [cut[unit]][+z]]``

    Do not draw contours with less than cut number of points.
    ($(GMTdoc)contour.html#q)
- **S** | **skip** :: [Type => Str | []]            ``Arg = [p|t]``

    Skip all input xyz points that fall outside the region.
    ($(GMTdoc)contour.html#s)
- **T** | **ticks** :: [Type => Str]                 ``Arg = [+|-][+a][+dgap[/length]][+l[labels]]``

    Draw tick marks pointing in the downward direction every *gap* along the innermost closed contours.
    ($(GMTdoc)contour.html#t)
- $(GMT.opt_R)
- $(GMT.opt_U)
- $(GMT.opt_V)
- **W** | **pen** :: [Type => Str | Number]

    Sets the attributes for the particular line.
    ($(GMTdoc)contour.html#w)
- $(GMT.opt_X)
- $(GMT.opt_Y)
- **Z** | **scale** :: [Type => Str]

    Use to subtract shift from the data and multiply the results by factor before contouring starts.
    ($(GMTdoc)contour.html#z)
- $(GMT.opt_bi)
- $(GMT.opt_bo)
- $(GMT.opt_d)
- $(GMT.opt_di)
- $(GMT.opt_do)
- $(GMT.opt_e)
- $(GMT.opt_h)
- $(GMT.opt_i)
- $(GMT.opt_p)
- $(GMT.opt_t)
- $(GMT.opt_swap_xy)
- $(GMT.opt_savefig)
"""
function contour(cmd0::String="", arg1=nothing; first=true, kwargs...)

    gmt_proggy = (IamModern[1]) ? "contour "  : "pscontour "
	length(kwargs) == 0 && return monolitic(gmt_proggy, cmd0, arg1)

	d, K, O = init_module(first, kwargs...)		# Also checks if the user wants ONLY the HELP mode
	dict_auto_add!(d)			# The ternary module may send options via another channel

	cmd, opt_B, opt_J, opt_R = parse_BJR(d, "", "", O, " -JX" * split(def_fig_size, '/')[1] * "/0")
	cmd, = parse_common_opts(d, cmd, [:UVXY :bo :c :d :do :e :p :t :params], first)
	cmd  = parse_these_opts(cmd, d, [[:D :dump], [:I :fill :colorize], [:N :no_clip], [:Q :cut], [:S :skip]])
	cmd *= add_opt_pen(d, [:L :mesh], "L", true)     # TRUE to also seek (lw,lc,ls)
	cmd, opt_W = parse_contour_AGTW(d, cmd)

	# If file name sent in, read it and compute a tight -R if this was not provided
	arg2 = nothing;		arg3 = nothing
	cmd, arg1, opt_R, info = read_data(d, cmd0, cmd, arg1, opt_R, false, true)
	if (occursin(" -I", cmd) || occursin("+c", opt_W))			# Only try to load cpt if -I was set
		cmd, N_used_, arg1, arg2, arg3 = get_cpt_set_R(d, "", cmd, opt_R, (arg1 === nothing ? 1 : 0), arg1, arg2, arg3, "pscontour")
	end
	N_used = (arg1 !== nothing) + (arg2 !== nothing) + (arg3 !== nothing)

	if (!occursin(" -C", cmd))			# Otherwise ignore an eventual :cont because we already have it
		cmd, args, n, = add_opt(d, cmd, "C", [:C :cont :contour :contours :levels], :data, Array{Any,1}([arg1, arg2, arg3]), (x="",))
		if (n > 0)
			for k = 3:-1:1
				(args[k] === nothing) && continue
				if (isa(args[k], Array{<:Number}))
					cmd *= arg2str(args[k], ',')
					if (length(args[k]) == 1)  cmd *= ","  end		# A single contour needs to end with a ","
					break
				elseif (isa(args[k], GMTcpt))
					arg1, arg2, arg3 = args[:]
					break
				end
			end
		end
	end

	if ((val = find_in_dict(d, [:E :index])[1]) !== nothing)
		cmd *= " -E"
		if (isa(val, Matrix{<:Real}) || isa(val, GMTdataset) || isa(val, Vector{<:GMTdataset}))   # Now need to find the free slot where to store the indices array
			(N_used == 0) ? arg1 = val : (N_used == 1 ? arg2 = val : arg3 = val)
		else
			cmd *= arg2str(val)
		end
	end

	if (!occursin(" -W", cmd) && !occursin(" -I", cmd) && !occursin(" -D", cmd))  cmd *= " -W"  end	# Use default pen

	if (occursin("-I", cmd) && !occursin("-C", cmd))
		r = round_wesn([info[1].data[5], info[1].data[6], info[1].data[5], info[1].data[6]])
		info[1].data[5], info[1].data[6] = r[1], r[2]
		opt_T = (isempty(current_cpt[1])) ? @sprintf(" -T%.14g/%.14g/11+n", info[1].data[5], info[1].data[6]) : ""
		if (N_used <= 1)
			cmd, arg1, arg2, = add_opt_cpt(d, cmd, [:C], 'C', N_used, arg1, arg2, true, true, opt_T, true)
		else
			cmd, arg2, arg3, = add_opt_cpt(d, cmd, [:C], 'C', N_used, arg2, arg3, true, true, opt_T, true)
		end
	end

#	if (occursin(" -I", cmd) && (!isa(arg1, GMTcpt) && !isa(arg2, GMTcpt) && !isa(arg3, GMTcpt)))
#		# If arg to a -C ends with .cpt, accept it
#		if ((ind = findfirst("-C", cmd)) !== nothing && !endswith(split(cmd[ind[2]+1:end])[1], ".cpt"))
#			error("fill option rquires passing a CPT")
#		end
#	end

	_cmd = finish_PS_nested(d, [gmt_proggy * cmd])
	return finish_PS_module(d, _cmd, "-D", K, O, true, arg1, arg2, arg3)
end

# ---------------------------------------------------------------------------------------------------
contour!(cmd0::String="", arg1=nothing; kw...) = contour(cmd0, arg1; first=false, kw...)
contour(arg1, cmd0::String=""; kw...) = contour(cmd0, arg1; first=true, kw...)
contour!(arg1, cmd0::String=""; kw...) = contour(cmd0, arg1; first=false, kw...)

# ---------------------------------------------------------------------------------------------------
const pscontour  = contour
const pscontour! = contour!


text:
#=
It is possible to show that the square root of two can be expressed as an infinite continued fraction.

√ 2 = 1 + 1/(2 + 1/(2 + 1/(2 + ... ))) = 1.414213...

By expanding this for the first four iterations, we get:

1 + 1/2 = 3/2 = 1.5
1 + 1/(2 + 1/2) = 7/5 = 1.4
1 + 1/(2 + 1/(2 + 1/2)) = 17/12 = 1.41666...
1 + 1/(2 + 1/(2 + 1/(2 + 1/2))) = 41/29 = 1.41379...

The next three expansions are 99/70, 239/169, and 577/408, but the eighth expansion, 1393/985, is the first example where the number of digits in the numerator exceeds the number of digits in the denominator.

In the first one-thousand expansions, how many fractions contain a numerator with more digits than denominator?
=#
function calc()
  p = BigInt(1)
  q = BigInt(1)
  count = 1
  total = 0
  while count <= 1000
    np = p + 2q
    nq = p + q
    p,q = np,nq
    if length(string(p)) > length(string(q))
      total = total + 1
    end
    count = count + 1
  end
  total
end
@time println(calc())


text:
module ntest
using  ..nonlineq

function test_broyden()

    fcn(x) = [ (x[1] + 3) * (x[2]^2 - 7) + 18,  sin( x[2] * exp(x[1]) - 1) ]

    x0 = [-0.5, 1.4]

    (x_sol, iter) = nonlineq.broyden(fcn, x0)

    println("Solution at: $(x_sol)")
    println("Iter: $(iter)" )

end

end # ntest

max_stars_repo_name: HomoModelicus/julia

max_issues_repo_name: HomoModelicus/julia

max_forks_repo_name: HomoModelicus/julia


text:
using GeometricBase.Config
using GeometricIntegrators.Integrators
using GeometricIntegrators.Integrators: CacheType, Parameters
using Test

struct MyParameters{DT,TT} <: Parameters{DT,TT} end
struct MyIntegrator{DT,TT} <: Integrator{DT,TT} end

params = MyParameters{Float64,Float64}()
int    = MyIntegrator{Float64,Float64}()

@test_throws ErrorException IntegratorCache(params)
@test_throws ErrorException IntegratorCache(int)
@test_throws ErrorException CacheType(Float32, params)


text:
# TODO: Make this fast
function expl2(a::AbstractSet, b::AbstractSet)
	d = sqrt(length(union(a, b)) - length(intersect(a, b)))
	return exp(-d)
end

expl2(a::AbstractVector, b::AbstractVector) = exp(-euclidean(a, b))


text:
"""
Computes the extinction coefficient \$\\alpha_\\nu\$ for hydrogen-related
bound-free and free-free transitions.

Includes:
* Neutral Hydrogen bound-free and free-free.
* H\$^-\$ bound-free and free-free.
* H\$_2^-\$ free-free.
* H\$_2^+\$ free-free.
* Rayleigh scattering by molecular H\$_2\$
"""

#=----------------------------------------------------------------------------
                            Catch-all functions
----------------------------------------------------------------------------=#
"""
    hminus_ff(
        λ::Unitful.Length,
        temperature::Unitful.Temperature,
        h_neutral_density::NumberDensity,
        electron_density::NumberDensity;
        recipe::String="stilley"
    )

Compute free-free extinction from H minus ion. Recipe can be one of:

- `stilley` (default): Interpolates table from
  [Stilley & Callaway (1970)](https://ui.adsabs.harvard.edu/abs/1970ApJ...160..245S/abstract),
  which is valid for λ up to 9113 nm.

- `john`: Follows
  [John (1988)](https://ui.adsabs.harvard.edu/abs/1988A%26A...193..189J/abstract),
  which is valid beyond 9113 nm but may not be good below 364.5 nm.
"""
function hminus_ff(
    λ::Unitful.Length,
    temperature::Unitful.Temperature,
    h_neutral_density::NumberDensity,
    electron_density::NumberDensity;
    recipe::String="stilley"
)
    if recipe == "stilley"
        hminus_ff_stilley(λ, temperature, h_neutral_density, electron_density)
    elseif recipe == "john"
        hminus_ff_john(λ, temperature, h_neutral_density, electron_density)
    else
        throw("NotImplemented recipe $recipe")
    end
end

"""
    hminus_bf(
        λ::Unitful.Length,
        temperature::Unitful.Temperature,
        h_neutral_density::NumberDensity,
        electron_density::NumberDensity;
        recipe::String="geltman"
    )

Compute bound-free extinction from H minus ion. Recipe can be one of:

- `geltman` (default): Uses recipe from
  [Geltman (1962)](https://ui.adsabs.harvard.edu/abs/1962ApJ...136..935G/abstract)

- `john`: Follows
  [John (1988)](https://ui.adsabs.harvard.edu/abs/1988A%26A...193..189J/abstract),
  which is valid beyond 9113 nm but may not be good below 364.5 nm.

- `wbr`: Follows
  [Wishart (1979)](https://ui.adsabs.harvard.edu/abs/1979MNRAS.187P..59W) for λ > 175 nm,
  and [Broad and Reinhardt (1976)](https://ui.adsabs.harvard.edu/abs/1976PhRvA..14.2159B)
  for λ <= 164 nm.
"""
function hminus_bf(
    λ::Unitful.Length,
    temperature::Unitful.Temperature,
    h_neutral_density::NumberDensity,
    electron_density::NumberDensity;
    recipe::String="geltman"
)
    if recipe == "geltman"
        hminus_bf_geltman(λ, temperature, h_neutral_density, electron_density)
    elseif recipe == "john"
        hminus_bf_john(λ, temperature, h_neutral_density, electron_density)
    elseif recipe == "wbr"
        hminus_bf_wbr(λ, temperature, h_neutral_density, electron_density)
    else
        throw("NotImplemented recipe $recipe")
    end
end

#=----------------------------------------------------------------------------
                   Recipes from Karzas and Latter / Kurucz
----------------------------------------------------------------------------=#
#= Tabulated Gaunt factors from Kurucz (1970, SAO Special Report no. 309), page 77,
 which is a fit to  figures 3-5 of Karzas and Latter (1961, ApJ Suppl 6, 167)

- Note from Mats Carlsson in RH:
    There is extrapolation outside the range of lg(gamma2) (T outside [1570,1.57e8K])
    or outside range of lg(U) (for lambda=3.2mm for T>49 kK). This should be OK for
    reasonable extrapolation distances (and better than setting to constant end value
    or zero). Interpolation tested against table in Gustafsson (1973) with results
    within 1%.
=#
const kurucz_ff_table = [5.53  5.49  5.46  5.43  5.40  5.25  5.00  4.69  4.48  4.16  3.85
                         4.91  4.87  4.84  4.80  4.77  4.63  4.40  4.13  3.87  3.52  3.27
                         4.29  4.25  4.22  4.18  4.15  4.02  3.80  3.57  3.27  2.98  2.70
                         3.64  3.61  3.59  3.56  3.54  3.41  3.22  2.97  2.70  2.45  2.20
                         3.00  2.98  2.97  2.95  2.94  2.81  2.65  2.44  2.21  2.01  1.81
                         2.41  2.41  2.41  2.41  2.41  2.32  2.19  2.02  1.84  1.67  1.50
                         1.87  1.89  1.91  1.93  1.95  1.90  1.80  1.68  1.52  1.41  1.30
                         1.33  1.39  1.44  1.49  1.55  1.56  1.51  1.42  1.33  1.25  1.17
                         0.90  0.95  1.00  1.08  1.17  1.30  1.32  1.30  1.20  1.15  1.11
                         0.55  0.58  0.62  0.70  0.85  1.01  1.15  1.18  1.15  1.11  1.08
                         0.33  0.36  0.39  0.46  0.59  0.76  0.97  1.09  1.13  1.10  1.08
                         0.19  0.21  0.24  0.28  0.38  0.53  0.76  0.96  1.08  1.09  1.09]
# table_y is log10(hν / kT)
const kurucz_ff_table_y = -4:0.5:1.5
# table_x is log10(3.28805e15 * (Z^2 h) / (kT))
const kurucz_ff_table_x = -3:0.5:2.0
const h_k = h / k_B
const hc_k = h * c_0 / k_B
const hminusχ = 0.754u"eV"
const saha_const = h^2 / (2 * π * m_e * k_B)  # Constant for Saha equation
const αff_const = (4 / (3 * h * c_0) * (e^2 / (4 * π * ε_0))^3 *
                   sqrt(2 * π / (3 * m_e^3 * k_B))) |> u"K^(1/2) * m^5 / s^3"
const αbf_const = (4 * e^2 / (3 * π * sqrt(3) * ε_0 * m_e * c_0^2 * R_∞)) |> u"m^2"
const kurucz_ff_interp = LinearInterpolation((kurucz_ff_table_y, kurucz_ff_table_x),
                                             kurucz_ff_table, extrapolation_bc=Line())

"""
    gaunt_ff(ν::Unitful.Frequency, temperature::Unitful.Temperature, charge::Int)
    gaunt_ff(λ::Unitful.Length, temperature::Unitful.Temperature, charge::Int)

Compute Gaunt factor for free-free based on [Karzas and Latter (1961, ApJ Suppl 6, 167)]
(https://ui.adsabs.harvard.edu/abs/1961ApJS....6..167K/abstract)
fit in
[Kurucz (1970, SAO Special Report no. 309), page 77](https://ui.adsabs.harvard.edu/abs/1970SAOSR.309.....K/abstract)
"""
function gaunt_ff(ν::Unitful.Frequency, temperature::Unitful.Temperature, charge::Int)
    lookup_y = log10(h_k * ν / temperature)
    lookup_x = log10(3.28805e15u"Hz" * charge^2 * h_k / temperature)
    return kurucz_ff_interp(lookup_y, lookup_x)::Float64
end

function gaunt_ff(λ::Unitful.Length, temperature::Unitful.Temperature, charge::Int)
    return gaunt_ff(c_0 / λ, temperature, charge)
end

#=----------------------------------------------------------------------------
                            Recipes from Seaton
----------------------------------------------------------------------------=#
"""
    gaunt_bf(charge::Int, n_eff::Number, λ::Unitful.Length)::Float64

Compute bound-free Gaunt factor for a given charge, effective principal
quantum number and wavelength λ. Taken from RH. Formula from
[Seaton (1960), Rep. Prog. Phys. 23, 313](https://ui.adsabs.harvard.edu/abs/1960RPPh...23..313S/abstract),
page 316.
"""
function gaunt_bf(λ::Unitful.Length, charge::Real, n_eff::Real)::Float64
    x = ustrip(1 / (λ * R_∞ * charge^2) |> u"m/m")
    x3 = x^(1/3)
    nsqx = 1 / (n_eff^2 * x)
    g_bf = 1 + 0.1728 * x3 * (1 - 2 * nsqx) - 0.0496 * x3^2 * (1 - (1 - nsqx) * 0.66666667 * nsqx)
    @assert g_bf >= 0 "gaunt_bf negative, calculation will not be reliable"
    return g_bf
end

function gaunt_bf(ν::Unitful.Frequency, charge::Real, n_eff::Real)::Float64
    return gaunt_bf(c_0 / ν, charge, n_eff)
end

"""
    n_eff(energy_upper::Unitful.Energy, energy_lower::Unitful.Energy, Z::Integer)

Compute the effective principal quantum number for a given energy difference
and atomic number `Z`.
"""
function n_eff(energy_upper::Unitful.Energy, energy_lower::Unitful.Energy, Z::Integer)
    return Z * sqrt(Ryh / (energy_upper - energy_lower))
end

#=----------------------------------------------------------------------------
                            Recipes from Stilley
----------------------------------------------------------------------------=#
const stilley_ff_λ = [0.0, 303.8, 455.6, 506.3, 569.5, 650.9, 759.4, 911.3,
                      1013.0, 1139.0, 1302.0, 1519.0, 1823.0, 2278.0, 3038.0,
                      4556.0, 9113.0]  # in nm
const stilley_ff_t = 5040.0 ./ [0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 1.1, 1.2,
                                1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9, 2.0]  # in K
const stilley_ff_table =
   [[0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00  #=    0.0 nm
=#   0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00];
    [3.44e-02 4.18e-02 4.91e-02 5.65e-02 6.39e-02 7.13e-02 7.87e-02 8.62e-02  #=  303.8 nm
=#   9.36e-02 1.01e-01 1.08e-01 1.16e-01 1.23e-01 1.30e-01 1.38e-01 1.45e-01];
    [7.80e-02 9.41e-02 1.10e-01 1.25e-01 1.40e-01 1.56e-01 1.71e-01 1.86e-01  #=  455.6 nm
=#   2.01e-01 2.16e-01 2.31e-01 2.45e-01 2.60e-01 2.75e-01 2.89e-01 3.03e-01];
    [9.59e-02 1.16e-01 1.35e-01 1.53e-01 1.72e-01 1.90e-01 2.08e-01 2.25e-01  #=  506.3 nm
=#   2.43e-01 2.61e-01 2.78e-01 2.96e-01 3.13e-01 3.30e-01 3.47e-01 3.64e-01];
    [1.21e-01 1.45e-01 1.69e-01 1.92e-01 2.14e-01 2.36e-01 2.58e-01 2.80e-01  #=  569.5 nm
=#   3.01e-01 3.22e-01 3.43e-01 3.64e-01 3.85e-01 4.06e-01 4.26e-01 4.46e-01];
    [1.56e-01 1.88e-01 2.18e-01 2.47e-01 2.76e-01 3.03e-01 3.31e-01 3.57e-01  #=  650.9 nm
=#   3.84e-01 4.10e-01 4.36e-01 4.62e-01 4.87e-01 5.12e-01 5.37e-01 5.62e-01];
    [2.10e-01 2.53e-01 2.93e-01 3.32e-01 3.69e-01 4.06e-01 4.41e-01 4.75e-01  #=  759.4 nm
=#   5.09e-01 5.43e-01 5.76e-01 6.08e-01 6.40e-01 6.72e-01 7.03e-01 7.34e-01];
    [2.98e-01 3.59e-01 4.16e-01 4.70e-01 5.22e-01 5.73e-01 6.21e-01 6.68e-01  #=  911.3 nm
=#   7.15e-01 7.60e-01 8.04e-01 8.47e-01 8.90e-01 9.32e-01 9.73e-01 1.01e+00];
    [3.65e-01 4.39e-01 5.09e-01 5.75e-01 6.39e-01 7.00e-01 7.58e-01 8.15e-01  #= 1013.0 nm
=#   8.71e-01 9.25e-01 9.77e-01 1.03e+00 1.08e+00 1.13e+00 1.18e+00 1.23e+00];
    [4.58e-01 5.50e-01 6.37e-01 7.21e-01 8.00e-01 8.76e-01 9.49e-01 1.02e+00  #= 1139.0 nm
=#   1.09e+00 1.15e+00 1.22e+00 1.28e+00 1.34e+00 1.40e+00 1.46e+00 1.52e+00];
    [5.92e-01 7.11e-01 8.24e-01 9.31e-01 1.03e+00 1.13e+00 1.23e+00 1.32e+00  #= 1302.0 nm
=#   1.40e+00 1.49e+00 1.57e+00 1.65e+00 1.73e+00 1.80e+00 1.88e+00 1.95e+00];
    [7.98e-01 9.58e-01 1.11e+00 1.25e+00 1.39e+00 1.52e+00 1.65e+00 1.77e+00  #= 1519.0 nm
=#   1.89e+00 2.00e+00 2.11e+00 2.21e+00 2.32e+00 2.42e+00 2.51e+00 2.61e+00];
    [1.14e+00 1.36e+00 1.58e+00 1.78e+00 1.98e+00 2.17e+00 2.34e+00 2.52e+00  #= 1823.0 nm
=#   2.68e+00 2.84e+00 3.00e+00 3.15e+00 3.29e+00 3.43e+00 3.57e+00 3.70e+00];
    [1.77e+00 2.11e+00 2.44e+00 2.75e+00 3.05e+00 3.34e+00 3.62e+00 3.89e+00  #= 2278.0 nm
=#   4.14e+00 4.39e+00 4.63e+00 4.86e+00 5.08e+00 5.30e+00 5.51e+00 5.71e+00];
    [3.10e+00 3.71e+00 4.29e+00 4.84e+00 5.37e+00 5.87e+00 6.36e+00 6.83e+00  #= 3038.0 nm
=#   7.28e+00 7.72e+00 8.14e+00 8.55e+00 8.95e+00 9.33e+00 9.71e+00 1.01e+01];
    [6.92e+00 8.27e+00 9.56e+00 1.08e+01 1.19e+01 1.31e+01 1.42e+01 1.52e+01  #= 4556.0 nm
=#   1.62e+01 1.72e+01 1.82e+01 1.91e+01 2.00e+01 2.09e+01 2.17e+01 2.25e+01];
    [2.75e+01 3.29e+01 3.80e+01 4.28e+01 4.75e+01 5.19e+01 5.62e+01 6.04e+01  #= 9133.0 nm
=#   6.45e+01 6.84e+01 7.23e+01 7.60e+01 7.97e+01 8.32e+01 8.67e+01 9.01e+01]]
const stilley_ff_interp = LinearInterpolation((stilley_ff_λ, stilley_ff_t[end:-1:1]),
                             stilley_ff_table[:, end:-1:1], extrapolation_bc=Line())

"""
    hminus_ff_stilley(
        λ::Unitful.Length,
        temperature::Unitful.Temperature,
        h_neutral_density::NumberDensity,
        electron_density::NumberDensity
    )

Compute free-free extinction from H minus ion, for a given wavelength, temperature,
density of neutral hydrogen atoms `h_neutral_density`, and electron density.
Interpolates table from
[Stilley & Callaway (1970)](https://ui.adsabs.harvard.edu/abs/1970ApJ...160..245S/abstract),
page 255, which is valid for λ up to 9113 nm.
"""
function hminus_ff_stilley(
    λ::Unitful.Length,
    temperature::Unitful.Temperature,
    h_neutral_density::NumberDensity,
    electron_density::NumberDensity
)
    electron_pressure = electron_density * k_B * temperature
    λi = ustrip(λ |> u"nm")   # convert to units of table
    temp = ustrip(temperature |> u"K")
    kappa = stilley_ff_interp(λi, temp)::Float64 * 1e-29u"m^4/N"
    return h_neutral_density * electron_pressure * kappa |> u"m^-1"
end

#=----------------------------------------------------------------------------
                            Recipes from Geltman
----------------------------------------------------------------------------=#
const geltman_bf_λ = [   0.0,   50.0,  100.0,  150.0,  200.0,  250.0,  300.0,
                       350.0,  400.0,  450.0,  500.0,  550.0,  600.0,  650.0,
                       700.0,  750.0,  800.0,  850.0,  900.0,  950.0, 1000.0,
                      1050.0, 1100.0, 1150.0, 1200.0, 1250.0, 1300.0, 1350.0,
                      1400.0, 1450.0, 1500.0, 1550.0, 1600.0, 1641.9]  # in nm
const geltman_bf_σ = [0.00, 0.15, 0.33, 0.57, 0.85, 1.17, 1.52, 1.89, 2.23,
                      2.55, 2.84, 3.11, 3.35, 3.56, 3.71, 3.83, 3.92, 3.95,
                      3.93, 3.85, 3.73, 3.58, 3.38, 3.14, 2.85, 2.54, 2.20,
                      1.83, 1.46, 1.06, 0.71, 0.40, 0.17, 0.0]  # in 1e-21 m^2

const geltman_bf_interp = LinearInterpolation(geltman_bf_λ, geltman_bf_σ, extrapolation_bc=0)

"""
    hminus_bf_geltman(
        λ::Unitful.Length,
        temperature::Unitful.Temperature,
        h_minus_density::NumberDensity
    )
    hminus_bf_geltman(
        λ::Unitful.Length,
        temperature::Unitful.Temperature,
        h_neutral_density::NumberDensity,
        electron_density::NumberDensity
    )

Compute extinction from H minus ion, from input H minus populations. Uses recipe from
[Geltman (1962)](https://ui.adsabs.harvard.edu/abs/1962ApJ...136..935G/abstract)
"""
function hminus_bf_geltman(
    λ::Unitful.Length,
    temperature::Unitful.Temperature,
    h_minus_density::NumberDensity
)
    λi = ustrip(λ |> u"nm")   # convert to units of table
    κ = geltman_bf_interp(λi)::Float64 * 1e-21u"m^2"
    stimulated_emission = exp(-hc_k / (λ * temperature))
    return h_minus_density * (1 - stimulated_emission) * κ
end

function hminus_bf_geltman(
    λ::Unitful.Length,
    temperature::Unitful.Temperature,
    h_neutral_density::NumberDensity,
    electron_density::NumberDensity
)
    h_minus_density = calc_hminus_density(h_neutral_density, temperature, electron_density)
    return hminus_bf_geltman(λ, temperature, h_minus_density)
end

"""
    calc_hminus_density(
        h_neutral_density::NumberDensity,
        temperature::Unitful.Temperature,
        electron_density::NumberDensity
    )

Compute H minus populations based on electron density, temperature, and
density of neutral hydrogen atoms.
"""
function calc_hminus_density(
    h_neutral_density::NumberDensity,
    temperature::Unitful.Temperature,
    electron_density::NumberDensity
)
    tmp = (ustrip(saha_const) / ustrip(temperature |> u"K"))^(3/2) * u"m^3"
    ϕ = tmp * exp(hminusχ / (k_B * temperature)) / 4
    return h_neutral_density * electron_density * ϕ
end

#=----------------------------------------------------------------------------
                            Recipes from John
----------------------------------------------------------------------------=#
"""
    hminus_ff_john(
        λ::Unitful.Length,
        temperature::Unitful.Temperature,
        h_neutral_density::NumberDensity,
        electron_density::NumberDensity
    )

Compute free-free extinction from H minus ion. Uses recipe from
[John (1988)](https://ui.adsabs.harvard.edu/abs/1988A%26A...193..189J/abstract),
which is valid beyond 9113 nm but may not be good below 364.5 nm.
"""
function hminus_ff_john(
    λ::Unitful.Length,
    temperature::Unitful.Temperature,
    h_neutral_density::NumberDensity,
    electron_density::NumberDensity
)
    λμ = ustrip(λ |> u"μm")
    if λμ > 0.3645
        table =
          [    0.0000     0.0000      0.0000      0.0000     0.0000    0.0000
            2483.3460   285.8270  -2054.2910   2827.7760 -1341.5370  208.9520
           -3449.8890 -1158.3820   8746.5230 -11485.6320  5303.6090 -812.9390
            2200.0400  2427.7190 -13651.1050  16755.5240 -7510.4940 1132.7380
            -696.2710 -1841.4000   8624.9700 -10051.5300  4400.0670 -655.0200
              88.2830   444.5170  -1863.8640   2095.2880  -901.7880  132.9850]
    else
        table =
          [  518.1021  -734.8666   1021.1775   -479.0721    93.1373   -6.4285
             473.2636  1443.4137  -1977.3395    922.3575  -178.9275   12.3600
            -482.2089  -737.1616   1096.8827   -521.1341   101.7963   -7.0571
             115.5291   169.6374   -245.6490    114.2430   -21.9972    1.5097
               0.0000     0.0000      0.0000      0.0000     0.0000    0.0000
               0.0000     0.0000      0.0000      0.0000     0.0000    0.0000]
    end
    sqrtθ = sqrt((5040.0u"K" / temperature) |> u"K/K")
    λinv = 1.0 / λμ
    κ = 0.0
    for i in 1:6
        κ += sqrtθ^(1 + i) * (λμ^2 * table[i, 1] + table[i, 2] +
                                λinv * (table[i, 3] + λinv * (table[i, 4] +
                                λinv * (table[i, 5] + λinv * table[i, 6]))))
    end
    # NOTE: in RH temperature from electron pressure is set to 5040 K!
    # Also, RH uses sometimes nHminus = nH * pe, other times atmos.nHmin...
    return κ * 1e-32u"m^4/N" * h_neutral_density * electron_density * k_B * temperature
end

"""
    hminus_bf_john(
        λ::Unitful.Length,
        temperature::Unitful.Temperature,
        h_neutral_density::NumberDensity,
        electron_density::NumberDensity
    )

Compute extinction from H minus ion. Uses recipe from
[John (1988)](https://ui.adsabs.harvard.edu/abs/1988A%26A...193..189J/abstract),
which is valid for  0.125 <= λ (μm) <= 1.6419. Seems to already include stimulated
emission.
"""
function hminus_bf_john(
    λ::Unitful.Length,
    temperature::Unitful.Temperature,
    h_neutral_density::NumberDensity,
    electron_density::NumberDensity
)
    table = [152.519, 49.534, -118.858, 92.536, -34.194, 4.982]
    λμ = ustrip(λ |> u"μm")
    λ0 = ustrip(1.6419u"μm")
    λ1 = ustrip(0.125u"μm")  # edge wavelength when approximation for fλ no longer valid
    temp = ustrip(temperature |> u"K")
    λidiff = max(0.0, 1.0 / λμ - 1.0 / λ0)  # cases beyond λ0 set to zero
    σλ = 1e-18 * λμ^3 * λidiff^1.5
    fλ = 0.0
    if λμ < λ1
        λidiff = 1.0 / λ1 - 1.0 / λ0
    end
    for n in 1:6
        fλ += table[n] * λidiff^((n-1)/2)
    end
    σλ *= fλ
    α = h_k * c_0
    κ = 0.750 * sqrt(temp)^-5 * exp(α / (λ0*u"μm" * temperature)) *
       (1 - exp(-α / (λ * temperature))) * σλ * 1e-3u"m^4/N"
    return κ * h_neutral_density * electron_density * k_B * temperature
end

#=----------------------------------------------------------------------------
        Recipes from Wishart (1979) and Broad and Reinhardt (1976)
----------------------------------------------------------------------------=#
const wbr_λ = [     18, 19.6, 21.4, 23.6, 26.4, 29.8, 34.3, 40.4, 49.1, 62.6,  121, 139,
                   164,  175,  200,  225,  250,  275,  300,  325,  350,  375,  400, 425,
                   450,  475,  500,  525,  550,  575,  600,  625,  650,  675,  700, 725,
                   750,  775,  800,  825,  850,  875,  900,  925,  950,  975, 1000, 1025,
                  1050, 1075, 1100, 1125, 1150, 1175, 1200, 1225, 1250, 1275, 1300, 1325,
                  1350, 1375, 1400, 1425, 1450, 1475, 1500, 1525, 1550, 1575, 1600, 1610,
                  1620, 1630]   # in nm
const wbr_σ = [0.067, 0.088, 0.117, 0.155, 0.206, 0.283, 0.414, 0.703,  1.24,  2.33,
                5.43,  5.91,  7.29, 7.918, 9.453, 11.08, 12.75, 14.46, 16.19, 17.92,
               19.65, 21.35, 23.02, 24.65, 26.24, 27.77, 29.23, 30.62, 31.94, 33.17,
               34.32, 35.37, 36.32, 37.17, 37.91, 38.54, 39.07, 39.48, 39.77, 39.95,
               40.01, 39.95, 39.77, 39.48, 39.06, 38.53, 37.89, 37.13, 36.25, 35.28,
               34.19, 33.01, 31.72, 30.34, 28.87, 27.33, 25.71, 24.02, 22.26, 20.46,
               18.62, 16.74, 14.85, 12.95, 11.07, 9.211, 7.407, 5.677, 4.052, 2.575,
               1.302, 0.8697, 0.4974, 0.1989]  # in 1e-22 m^2
const wbr_bf_interp = LinearInterpolation(wbr_λ, wbr_σ, extrapolation_bc=0)

"""
    hminus_bf_wbr(
        λ::Unitful.Length,
        temperature::Unitful.Temperature,
        h_minus_density::NumberDensity
    )
    hminus_bf_wbr(
        λ::Unitful.Length,
        temperature::Unitful.Temperature,
        h_neutral_density::NumberDensity,
        electron_density::NumberDensity
    )

Compute extinction from H minus ion, from input H minus populations. Uses recipe from
[Wishart (1979)](https://ui.adsabs.harvard.edu/abs/1979MNRAS.187P..59W) for λ down to 175 nm,
and recipe from [Broad and Reinhardt (1976)](https://ui.adsabs.harvard.edu/abs/1976PhRvA..14.2159B)
for λ=164 nm and below, following the recommendation from Mathisen (1984, MSC thesis).
"""
function hminus_bf_wbr(
    λ::Unitful.Length,
    temperature::Unitful.Temperature,
    h_minus_density::NumberDensity
)
    λi = ustrip(λ |> u"nm")   # convert to units of table
    κ = wbr_bf_interp(λi)::Float64 * 1e-22u"m^2"
    stimulated_emission = exp(-hc_k / (λ * temperature))
    return h_minus_density * (1 - stimulated_emission) * κ
end

function hminus_bf_wbr(
    λ::Unitful.Length,
    temperature::Unitful.Temperature,
    h_neutral_density::NumberDensity,
    electron_density::NumberDensity
)
    h_minus_density = calc_hminus_density(h_neutral_density, temperature, electron_density)
    return hminus_bf_wbr(λ, temperature, h_minus_density)
end

#=----------------------------------------------------------------------------
                            Recipes from Mihalas
----------------------------------------------------------------------------=#
"""
    hydrogenic_ff(
        ν::Unitful.Frequency,
        temperature::Unitful.Temperature,
        electron_density::NumberDensity,
        species_density::NumberDensity,
        charge::Int
    )

Compute free-free extinction for a hydrogen-like species. Following
Mihalas (1978) p. 101 and
[Rutten's IART](https://www.uio.no/studier/emner/matnat/astro/AST4310/h20/pensumliste/iart.pdf)
p 68. For the hydrogen case, `ion_density` is the proton density (H II).
"""
function hydrogenic_ff(
    ν::Unitful.Frequency,
    temperature::Unitful.Temperature,
    electron_density::NumberDensity,
    species_density::NumberDensity,
    charge::Int
)
    ν = ν |> u"s^-1"
    stimulated_emission = exp(-h_k * ν / temperature)
    return (αff_const * charge^2 / sqrt(temperature) * ν^-3 * electron_density *
       species_density * (1 - stimulated_emission) * gaunt_ff(ν, temperature, charge))
end

"""
    hydrogenic_bf(
        ν::Unitful.Frequency,
        ν_edge::Unitful.Frequency,
        temperature::Unitful.Temperature,
        species_density::NumberDensity,
        charge::Real,
        n_eff::AbstractFloat
    )

Compute bound-free extinction for a hydrogen-like species. Following
Mihalas (1978) p. 99 and
[Rutten's IART](https://www.uio.no/studier/emner/matnat/astro/AST4310/h20/pensumliste/iart.pdf)
p 70. Using simplified constant and expression for threshold cross section.
"""
function hydrogenic_bf(
    ν::Unitful.Frequency,
    ν_edge::Unitful.Frequency,
    temperature::Unitful.Temperature,
    species_density::NumberDensity,
    charge::Real,
    n_eff::AbstractFloat
)
    if ν < ν_edge
        return 0 * αbf_const * species_density
    else
        ν3_ratio = (ν_edge / ν)^3
        stimulated_emission = exp(-h_k * ν / temperature)
        return (αbf_const * charge^4 * ν3_ratio * species_density * n_eff *
                (1 - stimulated_emission) * gaunt_bf(ν, charge, n_eff))
    end
end

#=----------------------------------------------------------------------------
                            Recipes from Bell (1980)
----------------------------------------------------------------------------=#
const bell_ff_λ = [   0.0,  350.5,  414.2,   506.3,   569.6,  650.9,  759.4,
                    911.3, 1139.1, 1518.8,  1822.6,  2278.3, 3037.7, 3645.2,
                   4556.5, 6075.3, 9113.0, 11391.3, 15188.3]   # in nm
const bell_ff_t = 5040.0 ./ [0.5, 0.8, 1.0, 1.2, 1.6, 2.0, 2.8, 3.6]  # in K
const bell_ff_κ =
    [0.00e+00  0.00e+00  0.00e+00  0.00e+00  0.00e+00  0.00e+00  0.00e+00  0.00e+00
     4.17e-02  6.10e-02  7.34e-02  8.59e-02  1.11e-01  1.37e-01  1.87e-01  2.40e-01
     5.84e-02  8.43e-02  1.01e-01  1.17e-01  1.49e-01  1.82e-01  2.49e-01  3.16e-01
     8.70e-02  1.24e-01  1.46e-01  1.67e-01  2.10e-01  2.53e-01  3.39e-01  4.27e-01
     1.10e-01  1.54e-01  1.80e-01  2.06e-01  2.55e-01  3.05e-01  4.06e-01  5.07e-01
     1.43e-01  1.98e-01  2.30e-01  2.59e-01  3.17e-01  3.75e-01  4.92e-01  6.09e-01
     1.92e-01  2.64e-01  3.03e-01  3.39e-01  4.08e-01  4.76e-01  6.13e-01  7.51e-01
     2.73e-01  3.71e-01  4.22e-01  4.67e-01  5.52e-01  6.33e-01  7.97e-01  9.63e-01
     4.20e-01  5.64e-01  6.35e-01  6.97e-01  8.06e-01  9.09e-01  1.11e+00  1.32e+00
     7.36e-01  9.75e-01  1.09e+00  1.18e+00  1.34e+00  1.48e+00  1.74e+00  2.01e+00
     1.05e+00  1.39e+00  1.54e+00  1.66e+00  1.87e+00  2.04e+00  2.36e+00  2.68e+00
     1.63e+00  2.14e+00  2.36e+00  2.55e+00  2.84e+00  3.07e+00  3.49e+00  3.90e+00
     2.89e+00  3.76e+00  4.14e+00  4.44e+00  4.91e+00  5.28e+00  5.90e+00  6.44e+00
     4.15e+00  5.38e+00  5.92e+00  6.35e+00  6.99e+00  7.50e+00  8.32e+00  9.02e+00
     6.47e+00  8.37e+00  9.20e+00  9.84e+00  1.08e+01  1.16e+01  1.28e+01  1.38e+01
     1.15e+01  1.48e+01  1.63e+01  1.74e+01  1.91e+01  2.04e+01  2.24e+01  2.40e+01
     2.58e+01  3.33e+01  3.65e+01  3.90e+01  4.27e+01  4.54e+01  4.98e+01  5.33e+01
     4.03e+01  5.20e+01  5.70e+01  6.08e+01  6.65e+01  7.08e+01  7.76e+01  8.30e+01
     7.16e+01  9.23e+01  1.01e+02  1.08e+02  1.18e+02  1.26e+02  1.38e+02  1.47e+02]
# Linear extrapolation in λ, flat extrapolation in θ
const bell_ff_interp = LinearInterpolation((bell_ff_λ, bell_ff_t[end:-1:1]),
                   bell_ff_κ[:, end:-1:1], extrapolation_bc=(Line(), Flat()))

"""
    h2minus_ff(
        λ::Unitful.Length,
        temperature::Unitful.Temperature,
        h2_density::NumberDensity,
        electron_density::NumberDensity
    )

Compute extinction from H2^- molecule, according to recipe from
[Bell (1980)](https://ui.adsabs.harvard.edu/abs/1980JPhB...13.1859B/abstract),
page 1863. Stimulated emission is included.
"""
function h2minus_ff(
    λ::Unitful.Length,
    temperature::Unitful.Temperature,
    h2_density::NumberDensity,
    electron_density::NumberDensity
)
    electron_pressure = electron_density * k_B * temperature
    λi = ustrip(λ |> u"nm")   # convert to units of table
    temp = ustrip(temperature |> u"K")
    κ = bell_ff_interp(λi, temp)::Float64 * 1e-29u"m^4/N"
    return h2_density * electron_pressure * κ |> u"m^-1"
end

#=----------------------------------------------------------------------------
                            Recipes from Bates (1952)
----------------------------------------------------------------------------=#
const bates_λ = 1f7 ./ [   500,   1000,   1500,   2000,   2500,   3000,
                          3500,   4000,   5000,   6000,   7000,   8000,
                          9000, 10_000, 12_000, 14_000, 16_000, 18_000,
                        20_000, 22_000, 24_000, 26_000,  Inf32]  # in nm
const bates_t = [2.5f+03, 3.0f+03, 3.5f+03, 4.0f+03, 5.0f+03,
                    6.0f+03, 7.0f+03, 8.0f+03, 1.0f+04, 1.2f+04]  # in K
# H2plus bf + ff extinction coefficient in 1e-49 m^5
const bates_κ = convert(Array{Float32, 2},
                [  1.14  0.94  0.80  0.69  0.55  0.46  0.39  0.34  0.27  0.226
                   1.61  1.32  1.12  0.97  0.77  0.63  0.54  0.47  0.38  0.31
                   1.97  1.60  1.35  1.17  0.92  0.76  0.64  0.56  0.44  0.37
                   2.28  1.83  1.53  1.32  1.03  0.85  0.72  0.63  0.50  0.41
                   2.56  2.04  1.69  1.45  1.13  0.92  0.78  0.68  0.54  0.44
                   2.84  2.23  1.84  1.56  1.21  0.99  0.83  0.72  0.57  0.47
                   3.1   2.42  1.97  1.67  1.28  1.04  0.88  0.76  0.60  0.49
                   3.4   2.60  2.10  1.77  1.35  1.09  0.92  0.79  0.62  0.51
                   4.0   2.98  2.36  1.96  1.47  1.17  0.98  0.84  0.66  0.54
                   4.8   3.4   2.63  2.15  1.57  1.25  1.04  0.89  0.69  0.57
                   5.6   3.9   2.91  2.33  1.67  1.31  1.08  0.92  0.71  0.58
                   6.7   4.4   3.2   2.53  1.77  1.37  1.12  0.95  0.73  0.59
                   7.9   5.0   3.5   2.74  1.87  1.43  1.16  0.97  0.74  0.60
                   9.3   5.6   3.9   2.95  1.97  1.48  1.19  0.99  0.75  0.61
                  13.0   7.2   4.7   3.4   2.18  1.58  1.25  1.03  0.77  0.62
                  18.1   9.3   5.7   4.0   2.40  1.69  1.30  1.06  0.78  0.62
                  25.1  11.9   7.0   4.7   2.64  1.80  1.36  1.09  0.78  0.62
                  35.   15.2   8.4   5.4   2.91  1.91  1.41  1.11  0.79  0.61
                  47.   19.3  10.2   6.3   3.2   2.03  1.46  1.14  0.79  0.61
                  64.   24.3  12.2   7.3   3.5   2.16  1.52  1.16  0.79  0.60
                  86.   31.   14.6   8.4   3.8   2.29  1.57  1.18  0.79  0.59
                 114.   38.   17.3   9.6   4.2   2.42  1.63  1.21  0.79  0.58
                   0.    0.    0.    0.    0.    0.    0.    0.    0.    0.  ])
const bates_bf_fraction = convert(Array{Float32, 2},
        [0.059  0.046  0.037  0.031  0.022  0.017  0.014  0.011  0.008  0.006
         0.135  0.107  0.087  0.072  0.053  0.041  0.033  0.027  0.020  0.015
         0.214  0.171  0.141  0.118  0.088  0.069  0.056  0.046  0.034  0.026
         0.291  0.236  0.196  0.166  0.125  0.098  0.080  0.067  0.049  0.038
         0.363  0.298  0.250  0.214  0.162  0.129  0.105  0.088  0.065  0.050
         0.430  0.357  0.303  0.260  0.200  0.159  0.131  0.110  0.082  0.064
         0.490  0.413  0.353  0.305  0.237  0.190  0.157  0.132  0.099  0.077
         0.546  0.464  0.400  0.349  0.273  0.221  0.183  0.155  0.116  0.091
         0.640  0.556  0.486  0.429  0.342  0.280  0.234  0.200  0.151  0.120
         0.715  0.632  0.561  0.501  0.406  0.336  0.284  0.243  0.186  0.148
         0.775  0.696  0.625  0.564  0.464  0.388  0.331  0.285  0.220  0.176
         0.822  0.748  0.680  0.619  0.517  0.437  0.375  0.326  0.254  0.204
         0.859  0.792  0.727  0.667  0.564  0.482  0.417  0.364  0.286  0.232
         0.888  0.827  0.767  0.709  0.607  0.524  0.456  0.400  0.317  0.258
         0.929  0.881  0.829  0.777  0.680  0.597  0.526  0.467  0.375  0.309
         0.954  0.917  0.874  0.829  0.739  0.658  0.586  0.525  0.428  0.356
         0.970  0.942  0.906  0.867  0.786  0.708  0.638  0.577  0.476  0.400
         0.980  0.959  0.930  0.897  0.824  0.751  0.683  0.622  0.519  0.440
         0.987  0.970  0.947  0.919  0.854  0.786  0.721  0.661  0.558  0.476
         0.991  0.978  0.960  0.936  0.879  0.816  0.754  0.695  0.593  0.510
         0.994  0.984  0.969  0.949  0.898  0.841  0.782  0.726  0.625  0.541
         0.996  0.988  0.976  0.959  0.914  0.862  0.806  0.752  0.653  0.569
         0.     0.     0.     0.     0.     0.     0.     0.     0.     0.   ])
const bates_ff_fraction = 1 .- bates_bf_fraction
const bates_ff_κ = bates_κ .* bates_ff_fraction
const bates_bf_κ = bates_κ .* bates_bf_fraction
const bates_ff_interp = LinearInterpolation((bates_λ[end:-1:1], bates_t),
                          bates_ff_κ[end:-1:1, :], extrapolation_bc=(Line(), Flat()))
const bates_bf_interp = LinearInterpolation((bates_λ[end:-1:1], bates_t),
                          bates_bf_κ[end:-1:1, :], extrapolation_bc=(Line(), Flat()))

"""
    h2plus_ff(λ::Unitful.Length, temperature::Unitful.Temperature,
              h_neutral_density::NumberDensity, proton_density::NumberDensity)

Compute free-free extinction from H2plus molecule, according to recipe from
[Bates (1952)](https://ui.adsabs.harvard.edu/abs/1952MNRAS.112...40B/abstract),
page 43.
"""
function h2plus_ff(
    λ::Unitful.Length,
    temperature::Unitful.Temperature,
    h_neutral_density::NumberDensity,
    proton_density::NumberDensity
)
    λi = convert(Float32, ustrip(λ |> u"nm"))   # convert to units of table
    temp = convert(Float32, ustrip(temperature |> u"K"))
    κ = bates_ff_interp(λi, temp)::Float32 * u"m^5"
    # Table in 1e-49 m^5, splitting in two to prevent overflow in Float32 inputs
    return κ * (h_neutral_density * 1f-25) * (proton_density * 1f-24)
end

"""
    h2plus_bf(λ::Unitful.Length, temperature::Unitful.Temperature,
              h_neutral_density::NumberDensity, proton_density::NumberDensity)

Compute bound-free extinction from H2plus molecule, according to recipe from
[Bates (1952)](https://ui.adsabs.harvard.edu/abs/1952MNRAS.112...40B/abstract),
page 43.
"""
function h2plus_bf(
    λ::Unitful.Length,
    temperature::Unitful.Temperature,
    h_neutral_density::NumberDensity,
    proton_density::NumberDensity
)
    λi = convert(Float32, ustrip(λ |> u"nm"))  # convert to units of table
    temp = convert(Float32, ustrip(temperature |> u"K"))
    κ = bates_bf_interp(λi, temp)::Float32 * u"m^5"
    # Table in 1e-49 m^5, splitting in two to prevent overflow in Float32 inputs
    return κ * (h_neutral_density * 1f-25) * (proton_density * 1f-24)
end

#=----------------------------------------------------------------------------
                  Recipes from Victor & Dalgarno (1969)
----------------------------------------------------------------------------=#
const victor_h2_λ = [121.57, 130.00, 140.00, 150.00, 160.00, 170.00, 185.46,
                     186.27, 193.58, 199.05, 230.29, 237.91, 253.56, 275.36,
                     296.81, 334.24, 404.77, 407.90, 435.96, 546.23, 632.80]  # in nm
const victor_h2_σ = [2.35E-28, 1.22E-28, 6.80E-29, 4.24E-29, 2.84E-29, 2.00E-29,
                     1.25E-29, 1.22E-29, 1.00E-29, 8.70E-30, 4.29E-30, 3.68E-30,
                     2.75E-30, 1.89E-30, 1.36E-30, 8.11E-31, 3.60E-31, 3.48E-31,
                     2.64E-31, 1.04E-31, 5.69E-32]  # in m2
const victor_h2_interp = LinearInterpolation(victor_h2_λ, victor_h2_σ, extrapolation_bc=Line())

"""
    function rayleigh_h2(λ::Unitful.Length, h2_density::NumberDensity)

Compute extinction from Rayleigh scattering from H2 molecules. Uses recipe from
[Victor and Dalgarno (1969)](https://aip.scitation.org/doi/pdf/10.1063/1.1671412),
J. Chem. Phys. 50, 2535, page 2538 for λ <= 632.80, and the recipe from
[Tarafdar & Vardya (1973)](https://ui.adsabs.harvard.edu/abs/1973MNRAS.163..261T/abstract)
page 272, for λ > 632.80.
"""
function rayleigh_h2(λ::Unitful.Length, h2_density::NumberDensity)
    λi = ustrip(λ |> u"nm")
    if λi >= victor_h2_λ[1]
        if λi <= victor_h2_λ[end]
            σ_h2 = victor_h2_interp(λi)::Float64 * u"m^2"
        else
            λ2 =1 / λi^2
            # Tarafdar coeffs converted from λ^4, λ^6, λ^8 in Å to nm and cm^2 to m^2:
            σ_h2 = (8.779e-21 + (1.323e-16 + 2.245e-12 * λ2) * λ2) * λ2^2 * u"m^2"
        end
    else
        σ_h2 = 0.0u"m^2"
    end
    return σ_h2 * h2_density
end

#=----------------------------------------------------------------------------
                  Recipes from Dalgarno (1962)
----------------------------------------------------------------------------=#
"""
    function rayleigh_h(λ::Unitful.Length, h_ground_density::NumberDensity)

Compute extinction from Rayleigh scattering from H atoms. `h_ground_density`
is the number density of neutral hydrogen atoms in the ground state.
Uses recipe from Dalgarno (1962), Geophysics Corp. of America, Technical Report No. 62-28-A
(unavailable), which is accurate to 1% for λ > 125.0 nm.
"""
function rayleigh_h(λ::Unitful.Length, h_ground_density::NumberDensity)
    λi = ustrip(λ |> u"Å")
    if λi >= 1215.7
        λ2 = 1 / λi^2
        # First coefficient has conversion from Mbarn to m^2. From RH:
        σ_h = (5.81e-17 * λ2^2 * (1 + 2.452e6 * λ2 +  4.801e12 * λ2^2)) * u"m^2"
    else
        σ_h = 0.0u"m^2"
    end
    return σ_h * h_ground_density
end


text:
export body!, content!, loadcss!, loadjs!, load!, importhtml!

content!(o, sel, html::AbstractString; fade = true) =
  fade ?
    @js_(o, Blink.fill($sel, $html)) :
    @js_ o document.querySelector($sel).innerHTML = $html

content!(o, sel, html; fade = true) =
  content!(o, sel, stringmime(MIME"text/html"(), html), fade = fade)

body!(w, html; fade = true) = content!(w, "body", html, fade = fade)

function loadcss!(w, url)
  @js_ w begin
    @var link = document.createElement("link")
    link.type = "text/css"
    link.rel = "stylesheet"
    link.href = $url
    document.head.appendChild(link)
  end
end

function importhtml!(w, url; async=false)
  if async
    @js_ w begin
      @var link = document.createElement("link")
      link.rel = "import"
      link.href = $url
      document.head.appendChild(link)
    end
  else
    @js w begin
      @new Promise(function (resolve, reject)
          @var link = document.createElement("link")
          link.rel = "import"
          link.href = $url
          link.onload = (e) -> resolve(true)
          link.onerror = (e) -> reject(false)
          document.head.appendChild(link)
      end)
    end
  end
end

function loadjs!(w, url)
  @js w @new Promise(function (resolve, reject)
    @var script = document.createElement("script")
    script.src = $url
    script.onload = resolve
    script.onerror = (e) -> reject(
                               Dict("name"=>"JSLoadError",
                                    "message"=>"failed to load " + this.src)
                                    )
    document.head.appendChild(script)
  end)
end

isurl(f) = ismatch(r"^https?://", f)

function load!(w, file)
  if !isurl(file)
    resource(file)
    file = basename(file)
  end
  ext = Mux.extension(file)
  if ext == "js"
    loadjs!(w, file)
  elseif ext == "css"
    loadcss!(w, file)
  elseif ext == "html"
    importhtml!(w, file)
  else
    error("Blink: Unsupported file type")
  end
end


text:
import NNlib: DenseConvDims, PoolDims

## generic convolution

function convNd(x, w; stride=1, padding=0, dilation=1)
    # TODO: unflipweights
    cdims = DenseConvDims(x, w; stride=stride, padding=padding, dilation=dilation)
    return NNlib.conv(x, w, cdims)
end

function rrule(::typeof(convNd), x, w; stride=1, padding=0, dilation=1)
    cdims = DenseConvDims(x, w; stride=stride, padding=padding, dilation=dilation)
    y = NNlib.conv(x, w, cdims)
    function convNd_pullback(dy)
        # TODO: unflipweights
        dy = unthunk(dy)
        return NoTangent(), NNlib.∇conv_data(dy, w, cdims), NNlib.∇conv_filter(x, dy, cdims)
    end
    return y, convNd_pullback
end

## specific convolution functions

const conv1d = convNd
const conv2d = convNd
const conv3d = convNd

## maxpool2d

function maxpool2d(x, kernel_size; stride=kernel_size, padding=0, dilation=1)
    pdims = PoolDims(x, kernel_size; stride=stride, padding=padding, dilation=dilation)
    return NNlib.maxpool(x, pdims)
end

function rrule(::typeof(maxpool2d), x, kernel_size; stride=kernel_size, padding=0, dilation=1)
    pdims = PoolDims(x, kernel_size; stride=stride, padding=padding, dilation=dilation)
    y = NNlib.maxpool(x, pdims)
    function maxpool_pullback(dy)
        dy = unthunk(dy)
        return NoTangent(), NNlib.∇maxpool(dy, y, x, pdims), NoTangent()
    end
    return y, maxpool_pullback
end


text:
# # Leaf Level Photosynthesis
#
# Here, we will go one step further and look at the entire leaf-level response to different environmental conditions, including the impact of leaf diffusive conductance here. See next tutorial for how stomatal conductance responds to the environment.

# ## Leaf Diffusive Conductance and Stomatal Conductance
# Before, we focussed mainly on the demand-driven constraints through Rubisco and RuBP regeneration. Leaf diffusive conductance is highly important as it drives the suppy-side of photosynthesis and interacts with the energy balance and leaf temperature as latent heat fluxes are a major factor in surface cooling.
#
# Before, we have derived net rates of photosynthesis $A_n$, which have to be matched with the supply side through CO₂ diffusion:
#
# $$A_n^{diff} = g_{leaf,C}(C_a-C_c) = g_{leaf,C}\frac{P_a-P_c}{P_{atm}}$$
#
# which can be separated into diffusion from the air to the leaf surface with a boundary layer conductance $g_{b,C}$, diffusion from the surface to the interstitial air-space with stomatal conductance $g_{s,C}$ and diffusion from the  interstitial air-space to the chloroplast with mesophyll conductance $g_{m,C}$ (the letter C stands for CO₂ here, as the diffusion constants vary with species, e.g. H$_2$O stomatal conductance is a factor 1.6 higher than those for CO₂):
#
# $$A_n^{diff} = g_{b,C}(C_a-C_s) = g_{s,C}(C_s-C_i) = g_{m,C}(C_i-C_c) = g_{leaf,C}(C_a-C_c)$$
#
# $$g_{leaf,C} = \left(\frac{1}{g_{b,C}}+\frac{1}{g_{s,C}}+\frac{1}{g_{m,C}}\right)^{-1}$$
#
# The importance of diffusive conductance and its interplay with photosynthesis here is that the supply and demand rates determine $C_c$. A reduction in $C_c$ reduces the demand-limited rates while it increase the diffusion rates (at a given $g$). Most models run internal so-called A-$C_c$ iterations to ensure both rates are balanced and in steady-state ($\partial C_c/\partial t=0$). We implement this option but also opt to run stomatal conductance prognostically, as stomata don't open and close instantanously but have a response time of around 15 minutes.
#
# Below, we show simple examples of how the environmental conditions and the leaf diffusive conductance impact leaf photosynthesis.
#

## load packages
using PyPlot
using Land.Photosynthesis

## set the default floating type
FT     = Float32;

## define default photosynthesis parameter sets
c3_set = C3CLM(FT);
c4_set = C4CLM(FT);

## define leaf photosynthetic parameters and environmental conditions
leaf_3 = LeafPhotoContainer{FT}(APAR=1200, Vcmax25=60, Jmax25=120, Vpmax25=80, Rd25=1);
leaf_4 = LeafPhotoContainer{FT}(APAR=1200, Vcmax25=60, Jmax25=120, Vpmax25=80, Rd25=1);
envir  = EnvironmentConditions{FT}();
#----------------------------------------------------------------------------

# ## $A$ ~ $P_c$

## set the internal CO₂ to 30 Pa and temperature to 298.15 K and PAR to 1500
leaf_3.p_i=30; leaf_3.T=298.15; leaf_3.APAR=1500;
leaf_4.p_i=30; leaf_4.T=298.15; leaf_4.APAR=1500;

## set an array of internal CO₂
p_array = collect(FT, 5:200);
ac_3 = similar(p_array); aj_3 = similar(p_array); ap_3 = similar(p_array); ag_3 = similar(p_array); an_3 = similar(p_array);
ac_4 = similar(p_array); aj_4 = similar(p_array); ap_4 = similar(p_array); ag_4 = similar(p_array); an_4 = similar(p_array);

## update the temperature and PAR dependent photosynthesis
Photosynthesis.photosynthesis_TD!(c3_set, leaf_3, envir);
Photosynthesis.photosynthesis_TD!(c4_set, leaf_4, envir);

## update the CO₂ dependent photosynthesis
for i in eachindex(p_array);
    p_i = p_array[i];
    photosynthesis_PD!(c3_set, leaf_3, p_i);
    photosynthesis_PD!(c4_set, leaf_4, p_i);

    ac_3[i], aj_3[i], ap_3[i], ag_3[i], an_3[i] = leaf_3.Ac, leaf_3.Aj, leaf_3.Ap, leaf_3.Ag, leaf_3.An;
    ac_4[i], aj_4[i], ap_4[i], ag_4[i], an_4[i] = leaf_4.Ac, leaf_4.Aj, leaf_4.Ap, leaf_4.Ag, leaf_4.An;
end

## plot the CO₂ responses for C3 and C4 photosynthesis
figure(figsize=(10,6), dpi=100);
tight_layout(true);
subplot(1,2,1);
plot(p_array, ac_3, "r:", lw=2, label="Ac");
plot(p_array, aj_3, "b:", lw=2, label="Aj");
plot(p_array, ap_3, "g:", lw=2, label="Ap");
plot(p_array, ag_3, "k-", lw=1, label="Ag");
plot(p_array, an_3, "k:", lw=1, label="An");
xlabel("Leaf Internal CO₂ (Pa)"            , fontsize=16);
ylabel("Photosynthetic Rate (μmol m⁻² s⁻¹)", fontsize=16);
legend(loc="lower right")
subplot(1,2,2);
plot(p_array, ac_4, "r:", lw=2, label="Ac");
plot(p_array, aj_4, "b:", lw=2, label="Aj");
plot(p_array, ap_4, "g:", lw=2, label="Ap");
plot(p_array, ag_4, "k-", lw=1, label="Ag");
plot(p_array, an_4, "k:", lw=1, label="An");
xlabel("Leaf Internal CO₂ (Pa)"            , fontsize=16);
legend(loc="lower right")
gcf()
#----------------------------------------------------------------------------

# ## $A$ ~ $PAR$

## set the internal CO₂ to 30 Pa and temperature to 298.15 K and PAR to 1500
leaf_3.p_i=30; leaf_3.T=298.15; leaf_3.APAR=1500;
leaf_4.p_i=30; leaf_4.T=298.15; leaf_4.APAR=1500;

## set an array of APAR
par_array = collect(FT, 10:10:2000);
ac_3 = similar(par_array); aj_3 = similar(par_array); ap_3 = similar(par_array); ag_3 = similar(par_array); an_3 = similar(par_array);
ac_4 = similar(par_array); aj_4 = similar(par_array); ap_4 = similar(par_array); ag_4 = similar(par_array); an_4 = similar(par_array);

for i in eachindex(par_array);
    leaf_3.APAR = par_array[i];
    leaf_4.APAR = par_array[i];

    ## update the temperature dependent photosynthesis
    Photosynthesis.photosynthesis_TD!(c3_set, leaf_3, envir);
    Photosynthesis.photosynthesis_TD!(c4_set, leaf_4, envir);

    ## update the CO₂ dependent photosynthesis
    photosynthesis_PD!(c3_set, leaf_3, leaf_3.p_i);
    photosynthesis_PD!(c4_set, leaf_4, leaf_4.p_i);

    ac_3[i], aj_3[i], ap_3[i], ag_3[i], an_3[i] = leaf_3.Ac, leaf_3.Aj, leaf_3.Ap, leaf_3.Ag, leaf_3.An;
    ac_4[i], aj_4[i], ap_4[i], ag_4[i], an_4[i] = leaf_4.Ac, leaf_4.Aj, leaf_4.Ap, leaf_4.Ag, leaf_4.An;
end

## plot the CO₂ responses for C3 and C4 photosynthesis
figure(figsize=(10,6), dpi=100);
tight_layout(true);
subplot(1,2,1);
plot(par_array, ac_3, "r:", lw=2, label="Ac");
plot(par_array, aj_3, "b:", lw=2, label="Aj");
plot(par_array, ap_3, "g:", lw=2, label="Ap");
plot(par_array, ag_3, "k-", lw=1, label="Ag");
plot(par_array, an_3, "k:", lw=1, label="An");
xlabel("Absorbed PAR (μmol m⁻² s⁻¹)"       , fontsize=16);
ylabel("Photosynthetic Rate (μmol m⁻² s⁻¹)", fontsize=16);
legend(loc="lower right")
subplot(1,2,2);
plot(par_array, ac_4, "r:", lw=2, label="Ac");
plot(par_array, aj_4, "b:", lw=2, label="Aj");
plot(par_array, ap_4, "g:", lw=2, label="Ap");
plot(par_array, ag_4, "k-", lw=1, label="Ag");
plot(par_array, an_4, "k:", lw=1, label="An");
xlabel("Absorbed PAR (μmol m⁻² s⁻¹)"       , fontsize=16);
legend(loc="lower right")
gcf()
#----------------------------------------------------------------------------

# ## $A$ ~ $T$

## set the internal CO₂ to 30 Pa and temperature to 298.15 K and PAR to 1500
leaf_3.p_i=30; leaf_3.T=298.15; leaf_3.APAR=1500;
leaf_4.p_i=30; leaf_4.T=298.15; leaf_4.APAR=1500;

## set an array of APAR
t_array = collect(FT, 10:50) .+ FT(273.15);
ac_3 = similar(t_array); aj_3 = similar(t_array); ap_3 = similar(t_array); ag_3 = similar(t_array); an_3 = similar(t_array);
ac_4 = similar(t_array); aj_4 = similar(t_array); ap_4 = similar(t_array); ag_4 = similar(t_array); an_4 = similar(t_array);

for i in eachindex(t_array);
    leaf_3.T = t_array[i];
    leaf_4.T = t_array[i];

    ## update the temperature dependent photosynthesis
    Photosynthesis.photosynthesis_TD!(c3_set, leaf_3, envir);
    Photosynthesis.photosynthesis_TD!(c4_set, leaf_4, envir);

    ## update the CO₂ dependent photosynthesis
    photosynthesis_PD!(c3_set, leaf_3, leaf_3.p_i);
    photosynthesis_PD!(c4_set, leaf_4, leaf_4.p_i);

    ac_3[i], aj_3[i], ap_3[i], ag_3[i], an_3[i] = leaf_3.Ac, leaf_3.Aj, leaf_3.Ap, leaf_3.Ag, leaf_3.An;
    ac_4[i], aj_4[i], ap_4[i], ag_4[i], an_4[i] = leaf_4.Ac, leaf_4.Aj, leaf_4.Ap, leaf_4.Ag, leaf_4.An;
end

## plot the CO₂ responses for C3 and C4 photosynthesis
figure(figsize=(10,6), dpi=100);
tight_layout(true);
subplot(1,2,1);
plot(t_array, ac_3, "r:", lw=2, label="Ac");
plot(t_array, aj_3, "b:", lw=2, label="Aj");
plot(t_array, ap_3, "g:", lw=2, label="Ap");
plot(t_array, ag_3, "k-", lw=1, label="Ag");
plot(t_array, an_3, "k:", lw=1, label="An");
xlabel("Temperature (K)"                   , fontsize=16);
ylabel("Photosynthetic Rate (μmol m⁻² s⁻¹)", fontsize=16);
legend(loc="upper left")
subplot(1,2,2);
plot(t_array, ac_4, "r:", lw=2, label="Ac");
plot(t_array, aj_4, "b:", lw=2, label="Aj");
plot(t_array, ap_4, "g:", lw=2, label="Ap");
plot(t_array, ag_4, "k-", lw=1, label="Ag");
plot(t_array, an_4, "k:", lw=1, label="An");
xlabel("Temperature (K)"                   , fontsize=16);
legend(loc="upper left")
gcf()
#----------------------------------------------------------------------------

max_stars_repo_name: bbuman/Land

max_issues_repo_name: bbuman/Land

max_forks_repo_name: bbuman/Land


text:
using QuantumCumulants
using Test

@testset "mixed-order" begin

# Hilbertspace
hf = FockSpace(:cavity)
ha = NLevelSpace(:atom, 2)
h = hf⊗ha

# Operators
a = Destroy(h,:a,1)
σ(i,j) = Transition(h,:σ,i,j)

@test isequal(average(a'*a, [2,1]), average(a'*a))
@test isequal(average(a'*σ(1,2), [2,1]; mix_choice=minimum), average(a')*average(σ(1,2)))

he = meanfield([a'*a,σ(2,2)],a'*a + σ(2,2) + a'*σ(1,2) + a*σ(2,1))
he_avg1 = cumulant_expansion(he,2)
he_avg2 = cumulant_expansion(he,[2,1])
he_avg3 = cumulant_expansion(he,[2,1];mix_choice=minimum)

@test isequal(he_avg1.equations, he_avg2.equations)
@test !isequal(he_avg1.equations, he_avg3.equations)

he = meanfield(a'*σ(1,2), a'*a + σ(2,2) + a*σ(2,1))

# N-atom laser
# Parameters
N = 2 #number of atoms
κ, g, Γ23, Γ13, Γ12, Ω, Δc, Δ3 = cnumbers("κ g Γ_{23} Γ_{13} Γ_{12} Ω Δ_c Δ_3")

# Hilbertspace
hf = FockSpace(:cavity)
ha = ⊗([NLevelSpace(Symbol(:atom,i),3) for i=1:N]...)
h = hf ⊗ ha

# Operators
a = Destroy(h,:a)
σ(i,j,k) = Transition(h,Symbol("σ_{$k}"),i,j,k+1)

# Hamiltonian
H = -Δc*a'a + sum(g*(a'*σ(1,2,i) + a*σ(2,1,i)) for i=1:N) + sum(Ω*(σ(3,1,i) + σ(1,3,i)) for i=1:N) - sum(Δ3*σ(3,3,i) for i=1:N)

# Jumps
J = [a;[σ(1,2,i) for i=1:N];[σ(1,3,i) for i=1:N];[σ(2,3,i) for i=1:N]]

# Rates
rates = [κ;[Γ12 for i=1:N];[Γ13 for i=1:N];[Γ23 for i=1:N]]

# list of operators
ops = [a'a, σ(2,2,1), σ(3,3,1)]

he = meanfield(ops,H,J; rates=rates)
he_avg_ = cumulant_expansion(he,[2,1,1]) #second order average

he_avg = complete(he_avg_;order=[2,1,1])
@test isempty(find_missing(he_avg))
@test isempty(findall(x -> (aon = acts_on(x); 2 in aon && 3 in aon), he_avg.states))

he_avg_ = cumulant_expansion(he,[2,1,1];mix_choice=minimum)
he_avg = complete(he_avg_;order=[2,1,1],mix_choice=minimum)
@test isempty(find_missing(he_avg))
@test isempty(findall(x -> (aon = acts_on(x); 2 in aon && 3 in aon), he_avg.states))
@test isempty(findall(x -> (aon = acts_on(x); 1 in aon && 3 in aon), he_avg.states))
@test isempty(findall(x -> (aon = acts_on(x); 1 in aon && 2 in aon), he_avg.states))

end # testset


text:
type mycell
  color   #颜色
  index   #对于堆栈中的坐标
end

type maze
  arr # 迷宫主体
  # 长宽
  m
  n
  # 初始点、结束点
  st
  ed
  # 路径
  path
  # 栈指针
  p
end

max_stars_repo_name: RandomMaze/makeMaze

max_issues_repo_name: RandomMaze/makeMaze

max_forks_repo_name: RandomMaze/makeMaze


text:
import Pkg
psg = "C:/Users/dusty/Documents/PlaymatProjects/PlaymatSimulator"
gd = "C:/Users/dusty/Documents/PlaymatProjects/PlaymatGames/Mage"
Pkg.activate(psg)

using Colors
using Serialization
using PlaymatSimulator

"""
Elder Dragon Highlander aka Commander
EDH is a variant of an MtG game that uses a 100 card singleton deck with 1 or more cards designated as the "Commander" and a 40 point lifepool. The Commander should start on the battlefield and be visible to all players.
"""

USER_SETTINGS = deserialize("C:/Users/dusty/Documents/PlaymatProjects/PlaymatSimulator/tmp/user_selection.jls")
const GAME_NAME = USER_SETTINGS[:GAME_NAME]
const DECK_NAME = USER_SETTINGS[:DECK_NAME]

game_include("$gd/Base/notebooks/Base.jl")

game_include("$gd/Base/notebooks/game_settings.jl")

game_include("$gd/Base/notebooks/game_rules.jl")

game_include("$gd/$GAME_NAME/game_state.jl")

merge!(gs, USER_SETTINGS)
deck = deserialize("$gd/$GAME_NAME/decks/$DECK_NAME/$DECK_NAME.jls")
gs[:deck] = deck

add_texts!(gs)

reset_stage!(gs)

# required GameOne variables
const SCREEN_HEIGHT = gs[:SCREEN_HEIGHT]
const SCREEN_WIDTH = gs[:SCREEN_WIDTH]

play_music(gd * "/Base/music/chad_crouch_-_moonrise.mp3", 1)  # play_music(name, loops=-1)

PlaymatSimulator.GameOne.SDL2.ShowCursor(Int32(0))  # hides system mouse cursor