text:
"""
    group_generator_basic(generators::Array{T}; prnt=false) -> Set(group), num_all_mult

Calculate all other elements of this closed group by multiplying all generators and storing new elements and then multiplying the generators again with the new elements.
Very efficient with ``O(\\Omega \\cdot N)`` multiplications.
It is the most efficient algorithm for a given set of generators. Only the use of smaller sets of generators can reduce the number of multiplications further.

The `prnt` variable can be used to specify if the programm should output the progressbar. This can slow things down.
"""
function group_generator_basic(generators::Array{T}; prnt = false)::Tuple{Set{T},Int64} where T

    # Sets to store generated elements in different levels of generation
    group_prev_level = Set([])
    group_current_level = Set(copy(generators))
    group_next_level = Set([])

    level = 0
    size_after = 0
    new_elements = 0

    num_all_mult = 0

    while length(group_current_level) > 0
        level += 1
        num_multiplications = 0

        # calculate the number of multiplications in this level
        num_max_mult = length(group_current_level) * (length(group_current_level) + length(group_prev_level))

        if prnt
            prog = Progress(num_max_mult)
        end

        # empty next for new iteration
        group_next_level = Set([])

        # multiplicate current elements with all generators and save new elements to next
        for M in generators
            for N in group_current_level
                num_multiplications += 1
                el = M * N

                #union!(group_next_level, [el]) # slower alternative
                push!(group_next_level, el)

                # space optimization was replaced by switching to the Set data-type.
                #if num_multiplications % 1_000_000 == 0
                    # free up some space every 1 000 000 iterations (decrease RAM usage from 20+ GB to 2-3 GB)
                    #unique!(group_next_level)
                #end

                if prnt && num_multiplications % 1_000 in [0,1]
                    size_next_level = length(group_next_level)
                    update!(prog, num_multiplications; showvalues = [(:i,"$num_multiplications/$num_max_mult"),(:level, level),(:size,size_after),(:new_elements, new_elements),(:size_next_level, "$size_next_level ($(size_next_level/num_multiplications * 100)%)")])
                end
            end
        end

        if prnt print("current → prev...") end

        # current is now done, move to previous
        union!(group_prev_level, group_current_level)

        if prnt println(" done.") end

        if prnt print("next → current...") end

        # move all from next to current that are not in previous
        group_current_level = setdiff(group_next_level, group_prev_level)

        if prnt println(" done.") end

        # some stats to display
        size_after = length(group_prev_level)
        new_elements = length(group_current_level)
        num_all_mult += num_multiplications
        println("level #",level," size ", size_after, " and ", new_elements," new elements (",num_multiplications," multiplications)")
        if prnt
            println("level #",level," size ", size_after, " and ", new_elements," new elements (",num_multiplications," multiplications)")
            # countdown to prevent directly overwriting with progressbar
            print("5..."); sleep(1); print("4..."); sleep(1); print("3..."); sleep(1); print("2..."); sleep(1); println("1..."); sleep(1);
        end
    end

    println("\nEnded after ", level, " iterations. \nThe resulting group has ", size_after, " elements.")

    # returns the final group and the number of multiplications needed
    return group_prev_level, num_all_mult
end


text:
function raindrops(number)
    pling = ""
    plang = ""
    plong = ""

    if number % 3 == 0
        pling = "Pling"
    end
    if number % 5 == 0
        plang = "Plang"
    end
    if number % 7 == 0
        plong = "Plong"
    end

    plingplangplong= pling*plang*plong
    plingplangplong == "" ? string(number) : plingplangplong

end

ext: jl

max_stars_repo_name: Rigante-pl/Exercism-Julia

max_issues_repo_name: Rigante-pl/Exercism-Julia

max_forks_repo_name: Rigante-pl/Exercism-Julia


text:
workspace()
using BoxChopper

box = [0 1; 0 1; 0 1]
r₀ = [0.5, 0.5, 0.5]
nout = [1, -2, 1]

@code_warntype volfrac(box, nout, r₀)


text:
# Generate large types where each field has a default value to which it's
# initialized.

macro varset(name::Symbol, table)
    @assert table.head == :block
    table = table.args

    names = Any[]
    vars = Any[]
    defaults = Any[]
    parameters = Any[]
    parameters_expr = Expr(:parameters)

    for row in table
        if typeof(row) == Expr && row.head == :line
            continue
        end

        if typeof(row) == Symbol
            var = row
            typ = :Any
            default = :nothing
        elseif row.head == :tuple
            @assert 2 <= length(row.args) <= 3
            var = row.args[1]
            typ = row.args[2]
            default = length(row.args) > 2 ? row.args[3] : :nothing
        else
            error("Bad varset systax")
        end

        push!(names, var)
        push!(vars, :($(var)::$(typ)))
        push!(defaults, default)
        push!(parameters, Expr(:kw, var, default))
        parameters_expr = Expr(:parameters, parameters...)
    end

    new_with_defaults = Expr(:call, :new, names...)

    ex =
    quote
        type $(name)
            $(vars...)

            function $(name)()
                new($(defaults...))
            end

            function $(name)($(parameters_expr))
                $(new_with_defaults)
            end

            # shallow copy constructor
            function $(name)(a::$(name))
                b = new()
                for name in fieldnames($(name))
                    setfield!(b, name, getfield(a, name))
                end
                b
            end
        end

        function copy(a::$(name))
            $(name)(a)
        end
    end

    esc(ex)
end


text:
"""
$(TYPEDEF)
"""
struct EnsembleProblem{T,T2,T3,T4,T5} <: AbstractEnsembleProblem
  prob::T
  prob_func::T2
  output_func::T3
  reduction::T4
  u_init::T5
  safetycopy::Bool
end

DEFAULT_PROB_FUNC(prob,i,repeat) = prob
DEFAULT_OUTPUT_FUNC(sol,i) = (sol,false)
DEFAULT_REDUCTION(u,data,I) = append!(u, data), false
EnsembleProblem(prob;
    output_func = DEFAULT_OUTPUT_FUNC,
    prob_func= DEFAULT_PROB_FUNC,
    reduction = DEFAULT_REDUCTION,
    u_init = nothing,
    safetycopy = prob_func !== DEFAULT_PROB_FUNC) =
    EnsembleProblem(prob,prob_func,output_func,reduction,u_init,safetycopy)

EnsembleProblem(;prob,
    output_func = DEFAULT_OUTPUT_FUNC,
    prob_func= DEFAULT_PROB_FUNC,
    reduction = DEFAULT_REDUCTION,
    u_init = nothing, p = nothing, safetycopy = prob_func !== DEFAULT_PROB_FUNC) =
    EnsembleProblem(prob,prob_func,output_func,reduction,u_init,safetycopy)


text:
using Soss

m = @model x begin
    α ~ StudentT(ν = 3.0)
    β ~ Normal()
    y ~ For(1:1000) do j
        Normal(μ = α + β * x[j])
    end
    return y
end

#####################################

julia> prior(m, :y)
@model begin
        β ~ Normal()
        α ~ StudentT(ν = 3.0)
    end

julia> likelihood(m, :y)
@model (α, β, x) begin
        y ~ For(1:1000) do j
                Normal(μ = α + β * x[j])
            end
    end

############################################

julia> markovBlanket(m, :α)
@model (β, x) begin
        α ~ StudentT(ν = 3.0)
        y ~ For(1:1000) do j
                Normal(μ = α + β * x[j])
            end
    end

######################################################

julia> dag = digraph(m)
SimpleDigraph{Symbol} (n=4, m=3)

julia> dag.N
Dict{Symbol, Set{Symbol}} with 4 entries:
  :α => Set([:y])
  :y => Set()
  :β => Set([:y])
  :x => Set([:y])

julia> dag.NN
Dict{Symbol, Set{Symbol}} with 4 entries:
  :α => Set()
  :y => Set([:α, :β, :x])
  :β => Set()
  :x => Set()

############################################################

x = randn(1000);
y = rand(m(;x));

julia> symlogdensity(m(;x) | (;y))
3.4166191036395746 - (0.5(5687.909333430344 + (4008.967498395781α) + (1282.559235689009β) + (1000(α^2)) + (985.752713011282(β^2)) - (173.80770668906146α*β))) - (2.0(log(3.0 + α^2)))

#############################################

julia> dynamicHMC(m(;x) | (;y))
1000-element TupleArray with schema (β = Float64, α = Float64)
(β = -0.835456 ± 0.0315286, α = -2.07694 ± 0.030521)

m1 = @model μ, σ begin
    x ~ Normal(μ,σ)
    y ~ Normal(x,1) |> iid(20)
end

m2 = @model μ, σ begin
    z ~ Normal()
    x = σ * z + μ
    y ~ Normal(x,1) |> iid(20)
end


text:
module MapPlot

using UrbanMaps
using PyPlot

export plot

function PyPlot.plot(map::UrbanMap)
    buildings = map.buildings
    for b in buildings
        pts = b.points
        # extra point for connecting building
        n_p = size(pts,1)
        n_d = size(pts,2)
        new_points = zeros(n_p+1, n_d)
        new_points[1:n_p,:] = pts
        new_points[end,:] = pts[1,:]
        plot(new_points[:,1], new_points[:,2] ,"black",lw=2.0)
    end
end

end # module

max_stars_repo_name: etotheipluspi/ParticleFilters

max_issues_repo_name: etotheipluspi/ParticleFilters

max_forks_repo_name: etotheipluspi/ParticleFilters


text:
module CyclicArrays

 include("utils.jl")
 export CyclicArray,shiftc,stagger,diff

 """
     CyclicArray
 CyclicArray data structure. Available constructors:
 ```
 CyclicArray(data::AbstractArray{T,N}, connections::AbstractArray)
 ```
 To generate a new CyclicArray with the connections of x:
 ```
 CyclicArray(x::CyclicArray)
 ```
 To generate a new CyclicArray with the connections:
 ```
 CyclicArray(connections)
 ```
 To generate a new CyclicArray with the connections of y and data from x:
 ```
 CyclicArray(x::AbstractArray,y::CyclicArray)
 ```
 To generate a new CyclicArray with predifined connections:
 ```
 CyclicArray(x::AbstractArray,str::String)
 ```
 Available options:

  "1D" - one dimensional array, one face

  "1D2F" - one dimensional array, two faces

  "2D" - two dimensional array, one face

  "2DXY" - two dimensional array, one face, switching x-y directions

  "2DFL" - two dimensional array, one face, flip sides

  "3D" - three dimensional array, one face

  "3DXY" - three dimensional array, one face, switching x-y directions

  "cubed" - cubed sphere, three dimensional array, six faces

  "cubed2D" - cubed sphere, two dimensional array, six faces
 """
 struct CyclicArray{T,N} <: AbstractArray{T,N} #inherits from AbstractArray
    data::AbstractArray{T,N}
    connections::Array{T,4}
    function CyclicArray(data::AbstractArray{T,N},connections) where {T,N}
      #creates the type only with a vector x
      if ((length(size(data)))<(size(connections)[2]))&(~isempty(data))
        error(string("number of array dimensions should be larger or equal to the spatial dimensions (",size(connections)[2],")"))
      else
        return new{T,N}(data,connections)
      end
    end
 end

 Base.show(io::IO, ::MIME"text/plain", A::CyclicArray{T,1}) where{T} =
           print(io, length(A),"-element CyclicArray{$T,1}:\n   ", A.data)

 Base.size(A::CyclicArray) = size(A.data)
 CyclicArray(x::CyclicArray)=CyclicArray([],x.connections)
 CyclicArray(x::AbstractArray,y::CyclicArray)=CyclicArray(x,y.connections)
 CyclicArray(connections::Array{T,4}) where T =CyclicArray([],connections)

 function MakeCyclicArray(x::AbstractArray,str::String)
  if str=="1D"
    nfaces=1;
    faces=zeros(nfaces,1,2,4);
    faces[1,1,1,:]=[1,1,2,0];
    faces[1,1,2,:]=[1,1,1,0];
  elseif str=="1D2F"
      nfaces=2;
    faces=zeros(nfaces,1,2,4);
    faces[1,1,1,:]=[2,1,2,0];
    faces[1,1,2,:]=[2,1,1,0];
    faces[2,1,1,:]=[1,1,2,0];
    faces[2,1,2,:]=[1,1,1,0];
  elseif str=="2D"
    nfaces=1;
    faces=zeros(nfaces,2,2,4);
    faces[1,1,1,:]=[1,1,2,0];
    faces[1,1,2,:]=[1,1,1,0];
    faces[1,2,1,:]=[1,2,2,0];
    faces[1,2,2,:]=[1,2,1,0];
  elseif str=="2DXY"
    nfaces=1;
    faces=zeros(nfaces,2,2,4);
    faces[1,1,1,:]=[1,2,2,0];
    faces[1,1,2,:]=[1,2,1,0];
    faces[1,2,1,:]=[1,1,2,0];
    faces[1,2,2,:]=[1,1,1,0];
  elseif str=="2DFL"
    nfaces=1;
    faces=zeros(nfaces,2,2,4);
    faces[1,1,1,:]=[1,1,2,1];
    faces[1,1,2,:]=[1,1,1,1];
    faces[1,2,1,:]=[1,2,2,1];
    faces[1,2,2,:]=[1,2,1,1];
  elseif str=="3D"
    nfaces=1;
    faces=zeros(nfaces,3,2,4);
    faces[1,1,1,:]=[1,1,2,0];
    faces[1,1,2,:]=[1,1,1,0];
    faces[1,2,1,:]=[1,2,2,0];
    faces[1,2,2,:]=[1,2,1,0];
    faces[1,3,1,:].=-1;
    faces[1,3,2,:].=-1;
  elseif str=="3DXY"
    nfaces=1;
    faces=zeros(nfaces,3,2,4);
    faces[1,1,1,:]=[1,2,2,0];
    faces[1,1,2,:]=[1,2,1,0];
    faces[1,2,1,:]=[1,1,2,0];
    faces[1,2,2,:]=[1,1,1,0];
    faces[1,3,1,:].=-1;
    faces[1,3,2,:].=-1;
  elseif str=="cubed"
    nfaces=6
    faces=zeros(Int64,nfaces,3,2,4);

    faces[1,1,1,:]=[4,1,2,0]
    faces[1,1,2,:]=[2,1,1,0]
    faces[1,2,1,:]=[6,2,2,0]
    faces[1,2,2,:]=[5,2,1,0]
    faces[1,3,1,:].=-1
    faces[1,3,2,:].=-1

    faces[2,1,1,:]=[1,1,2,0]
    faces[2,1,2,:]=[3,1,1,0]
    faces[2,2,1,:]=[6,1,2,1]
    faces[2,2,2,:]=[5,1,2,0]
    faces[2,3,1,:].=-1
    faces[2,3,2,:].=-1

    faces[3,1,1,:]=[2,1,2,0]
    faces[3,1,2,:]=[4,1,1,0]
    faces[3,2,1,:]=[6,2,1,1]
    faces[3,2,2,:]=[5,2,2,1]
    faces[3,3,1,:].=-1
    faces[3,3,2,:].=-1

    faces[4,1,1,:]=[3,1,2,0]
    faces[4,1,2,:]=[1,1,1,0]
    faces[4,2,1,:]=[6,1,1,0]
    faces[4,2,2,:]=[5,1,1,1]
    faces[4,3,1,:].=-1
    faces[4,3,2,:].=-1

    faces[5,1,1,:]=[4,2,2,1]
    faces[5,1,2,:]=[2,2,2,0]
    faces[5,2,1,:]=[1,2,2,0]
    faces[5,2,2,:]=[3,2,2,0]
    faces[5,3,1,:].=-1
    faces[5,3,2,:].=-1

    faces[6,1,1,:]=[4,2,1,0]
    faces[6,1,2,:]=[2,2,1,1]
    faces[6,2,1,:]=[3,2,1,1]
    faces[6,2,2,:]=[1,2,1,0]
    faces[6,3,1,:].=-1
    faces[6,3,2,:].=-1
  elseif str=="cubed2D"
    nfaces=6
    faces=zeros(Int64,nfaces,2,2,4);

    faces[1,1,1,:]=[4,1,2,0]
    faces[1,1,2,:]=[2,1,1,0]
    faces[1,2,1,:]=[6,2,2,0]
    faces[1,2,2,:]=[5,2,1,0]

    faces[2,1,1,:]=[1,1,2,0]
    faces[2,1,2,:]=[3,1,1,0]
    faces[2,2,1,:]=[6,1,2,1]
    faces[2,2,2,:]=[5,1,2,0]

    faces[3,1,1,:]=[2,1,2,0]
    faces[3,1,2,:]=[4,1,1,0]
    faces[3,2,1,:]=[6,2,1,1]
    faces[3,2,2,:]=[5,2,2,1]

    faces[4,1,1,:]=[3,1,2,0]
    faces[4,1,2,:]=[1,1,1,0]
    faces[4,2,1,:]=[6,1,1,0]
    faces[4,2,2,:]=[5,1,1,1]

    faces[5,1,1,:]=[4,2,2,1]
    faces[5,1,2,:]=[2,2,2,0]
    faces[5,2,1,:]=[1,2,2,0]
    faces[5,2,2,:]=[3,2,2,0]

    faces[6,1,1,:]=[4,2,1,0]
    faces[6,1,2,:]=[2,2,1,1]
    faces[6,2,1,:]=[3,2,1,1]
    faces[6,2,2,:]=[1,2,1,0]
  end
  return CyclicArray(x,faces);
 end

 CyclicArray(x::AbstractArray,str::String)=MakeCyclicArray(x,str)
 #CyclicArray(str::String)=MakeCyclicArray([],str)

 Base.view(A::CyclicArray) = Base.view(A.data)
 Base.ndims(A::CyclicArray) = ndims(A.data)
 Base.Dims(A::CyclicArray) = Dims(A.data)

 Base.length(A::CyclicArray)=length(A.data)
 Base.similar(A::CyclicArray)=CyclicArray(Array{eltype(A.data)}(undef, size(A.data)),A.connections)
 Base.similar(A::CyclicArray, S::DataType)=CyclicArray(Array{S}(undef, size(A.data)),A.connections)
 Base.similar(A::CyclicArray, dims::Dims)=CyclicArray(similar(A.data,dims),A.connections)
 # The next line produce error, not clear why
 #Base.similar(A::CyclicArray, S::DataType, dims::Dims)=CyclicArray(Array{S}(undef, dims),A.connections)
 Base.axes(A::CyclicArray) = axes(A.data)

 Base.findall(A::CyclicArray) = findall(A.data)
 Base.findfirst(A::CyclicArray) = findfirst(A.data)
 Base.findlast(A::CyclicArray) = findlast(A.data)
 Base.findprev(A::CyclicArray, I...) = findprev(A.data, I...)
 Base.findnext(A::CyclicArray, I...) = findnext(A.data, I...)
 Base.checkbounds(A::CyclicArray, I...) = nothing
 Base.CartesianIndices(A::CyclicArray) = CartesianIndices(A.data)
 Base.maximum(A::CyclicArray,dims) = maximum(A.data, dims=dims)
 Base.minimum(A::CyclicArray,dims) = minimum(A.data, dims=dims)

 """
    diff

 ```
 diff(A::CyclicArray)
 diff(A::CyclicArray; dims::Integer)
 ```

 Finite difference operator on a vector or a multidimensional array A. In the latter case the dimension to operate on needs to be
 specified with the dims keyword argument.
 """
 function Base.diff(A::CyclicArray; dims=1::Integer)
  I=size(A.data)
  I1=[UnitRange(1:I[i]) for i in 1:length(I)]
  I2=copy(I1)
  I2[dims]=I1[dims].+1
  B=CyclicArray(zeros(size(A.data)),A.connections)
  CI1=CartesianIndices(Tuple(I1))
  CI2=CartesianIndices(Tuple(I2))
  N=length(CI1[1])
  for n in 1:length(CI1)
    I1b=[CI1[n][i]  for i in 1:N]
    I2b=[CI2[n][i]  for i in 1:N]
    B[I1b...]=A[I2b...]-A[I1b...]
  end
  return B
 end

 """
    stagger
 ```
 stagger(A::CyclicArray; dims=1::Integer, frac=0.5::Real)
 ```
 Linear interpolation of the array to intermidiate grid positions
 """
 function stagger(A::CyclicArray; dims=1::Integer, frac=0.5::Real)
  I=size(A.data)
  I1=[UnitRange(1:I[i]) for i in 1:length(I)]
  I2=copy(I1)
  I2[dims]=I1[dims].+Int(sign(frac))
  B=CyclicArray(zeros(size(A.data)),A.connections)
  CI1=CartesianIndices(Tuple(I1))
  CI2=CartesianIndices(Tuple(I2))
  N=length(CI1[1])
  for n in 1:length(CI1)
    I1b=[CI1[n][i]  for i in 1:N]
    I2b=[CI2[n][i]  for i in 1:N]
    B[I1b...]=(1-frac)*A[I1b...]+frac*A[I2b...]
  end
  return B
 end

"""
    shiftc
```
shiftc(A::CyclicArray; dims=1::Integer, shift=1::Real)
```
Shifts array by an integer
"""
 function shiftc(A::CyclicArray; dims=1::Integer, shift=1::Real)
   I=size(A.data)
   I1=[UnitRange(1:I[i]) for i in 1:length(I)]
   I2=copy(I1)
   I2[dims]=I1[dims].-shift
   B=CyclicArray(zeros(size(A.data)),A.connections)
   CI1=CartesianIndices(Tuple(I1))
   CI2=CartesianIndices(Tuple(I2))
   N=length(CI1[1])
   for n in 1:length(CI1)
     I1b=[CI1[n][i]  for i in 1:N]
     I2b=[CI2[n][i]  for i in 1:N]
     B[I1b...]=A[I2b...]
   end
   return B
  end

 Base.:-(A::CyclicArray)=CyclicArray(.-A.data,A.connections)

 Base.:*(A::CyclicArray, B::CyclicArray)=CyclicArray(.*(A.data,B.data),A.connections)
 Base.:*(A::Array, B::CyclicArray)=CyclicArray(.*(A,B.data),B.connections)
 Base.:*(A::CyclicArray, B::Array)=CyclicArray(.*(A.data,B),A.connections)
 Base.:*(A::Number, B::CyclicArray)=CyclicArray(.*(A,B.data),B.connections)
 Base.:*(A::CyclicArray, B::Number)=CyclicArray(.*(A.data,B),A.connections)

 Base.:-(A::CyclicArray, B::CyclicArray)=CyclicArray(.-(A.data,B.data),A.connections)
 Base.:-(A::Number, B::CyclicArray)=CyclicArray(.-(A,B.data),B.connections)
 Base.:-(A::CyclicArray, B::Number)=CyclicArray(.-(A.data,B),A.connections)
 Base.:-(A::Array, B::CyclicArray)=CyclicArray(.-(A,B.data),B.connections)
 Base.:-(A::CyclicArray, B::Array)=CyclicArray(.-(A.data,B),A.connections)

 Base.:+(A::CyclicArray, B::CyclicArray)=CyclicArray(.+(A.data,B.data),A.connections)
 Base.:+(A::Number, B::CyclicArray)=CyclicArray(.+(A,B.data),B.connections)
 Base.:+(A::CyclicArray, B::Number)=CyclicArray(.+(A.data,B),A.connections)
 Base.:+(A::Array, B::CyclicArray)=CyclicArray(.+(A,B.data),B.connections)
 Base.:+(A::CyclicArray, B::Array)=CyclicArray(.+(A.data,B),A.connections)

 Base.:^(A::CyclicArray, B::CyclicArray)=CyclicArray(.^(A.data,B.data),A.connections)
 Base.:^(A::Number, B::CyclicArray)=CyclicArray(.^(A,B.data),B.connections)
 Base.:^(A::CyclicArray, B::Number)=CyclicArray(.^(A.data,B),A.connections)
 Base.:^(A::Array, B::CyclicArray)=CyclicArray(.^(A,B.data),B.connections)
 Base.:^(A::CyclicArray, B::Array)=CyclicArray(.^(A.data,B),A.connections)

 Base.:/(A::CyclicArray, B::CyclicArray)=CyclicArray(./(A.data,B.data),A.connections)
 Base.:/(A::Number, B::CyclicArray)=CyclicArray(./(A,B.data),B.connections)
 Base.:/(A::CyclicArray, B::Number)=CyclicArray(./(A.data,B),A.connections)
 Base.:/(A::Array, B::CyclicArray)=CyclicArray(./(A,B.data),B.connections)
 Base.:/(A::CyclicArray, B::Array)=CyclicArray(./(A.data,B),A.connections)

 Base.:\(A::CyclicArray, B::CyclicArray)=CyclicArray(.\(A.data,B.data),A.connections)
 Base.:\(A::Number, B::CyclicArray)=CyclicArray(.\(A,B.data),B.connections)
 Base.:\(A::CyclicArray, B::Number)=CyclicArray(.\(A.data,B),A.connections)
 Base.:\(A::Array, B::CyclicArray)=CyclicArray(.\(A,B.data),B.connections)
 Base.:\(A::CyclicArray, B::Array)=CyclicArray(.\(A.data,B),A.connections)

 Base.:>=(A::CyclicArray, B::CyclicArray)=CyclicArray(Array{Int}(.>=(A.data,B.data)),A.connections)
 Base.:>=(A::Number, B::CyclicArray)=CyclicArray(Array{Int}(.>=(A,B.data)),B.connections)
 Base.:>=(A::CyclicArray, B::Number)=CyclicArray(Array{Int}(.>=(A.data,B)),A.connections)
 Base.:>=(A::Array, B::CyclicArray)=CyclicArray(Array{Int}(.>=(A,B.data)),B.connections)
 Base.:>=(A::CyclicArray, B::Array)=CyclicArray(Array{Int}(.>=(A.data,B)),A.connections)

 Base.:<=(A::CyclicArray, B::CyclicArray)=CyclicArray(Array{Int}(.<=(A.data,B.data)),A.connections)
 Base.:<=(A::Number, B::CyclicArray)=CyclicArray(Array{Int}(.<=(A,B.data)),B.connections)
 Base.:<=(A::CyclicArray, B::Number)=CyclicArray(Array{Int}(.<=(A.data,B)),A.connections)
 Base.:<=(A::Array, B::CyclicArray)=CyclicArray(Array{Int}(.<=(A,B.data)),B.connections)
 Base.:<=(A::CyclicArray, B::Array)=CyclicArray(Array{Int}(.<=(A.data,B)),A.connections)

 Base.:>(A::CyclicArray, B::CyclicArray)=CyclicArray(Array{Int}(.>(A.data,B.data)),A.connections)
 Base.:>(A::Number, B::CyclicArray)=CyclicArray(Array{Int}(.>(A,B.data)),B.connections)
 Base.:>(A::CyclicArray, B::Number)=CyclicArray(Array{Int}(.>(A.data,B)),A.connections)
 Base.:>(A::Array, B::CyclicArray)=CyclicArray(Array{Int}(.>(A,B.data)),B.connections)
 Base.:>(A::CyclicArray, B::Array)=CyclicArray(Array{Int}(.>(A.data,B)),A.connections)

 Base.:<(A::CyclicArray, B::CyclicArray)=CyclicArray(Array{Int}(.<(A.data,B.data)),A.connections)
 Base.:<(A::Number, B::CyclicArray)=CyclicArray(Array{Int}(.<(A,B.data)),B.connections)
 Base.:<(A::CyclicArray, B::Number)=CyclicArray(Array{Int}(.<(A.data,B)),A.connections)
 Base.:<(A::Array, B::CyclicArray)=CyclicArray(Array{Int}(.<(A,B.data)),B.connections)
 Base.:<(A::CyclicArray, B::Array)=CyclicArray(Array{Int}(.<(A.data,B)),A.connections)

 Base.:(==)(A::CyclicArray, B::CyclicArray)=CyclicArray(Array{Int}(.==(A.data,B.data)),A.connections)
 Base.:(==)(A::Number, B::CyclicArray)=CyclicArray(Array{Int}(.==(A,B.data)),B.connections)
 Base.:(==)(A::CyclicArray, B::Number)=CyclicArray(Array{Int}(.==(A.data,B)),A.connections)
 Base.:(==)(A::Array, B::CyclicArray)=CyclicArray(Array{Int}(.==(A,B.data)),B.connections)
 Base.:(==)(A::CyclicArray, B::Array)=CyclicArray(Array{Int}(.==(A.data,B)),A.connections)

 function CalcIndex(A::CyclicArray, I::Vararg{Int, N}) where N # implements A[I]
  connections=A.connections
  nfaces=size(A.connections)[1]
  nspatial=size(A.connections)[2]
  I1=[i for i in I]
  N0=length(I1)
  S=size(A)
  if nspatial==3
    nx=S[N0];ny=S[N0-1];nz=S[N0-2];
    while (I1[N0]<1 || I1[N0]>nx) || (I1[N0-1]<1 || I1[N0-1]>ny) || (I1[N0-2]<1 || I1[N0-2]>nz)
      I2=I1
      i=I1[N0]; j=I1[N0-1]; k=I1[N0-2]
      if nfaces>1
        f=I1[N0-3]
      else
        f=1
      end
      if I1[N0]<1
        if connections[f,1,1,1]==-1
          return NaN
        else
          f1=connections[f,1,1,1]
          axis=connections[f,1,1,2];
          side=connections[f,1,1,3];
          flip=connections[f,1,1,4];
          I2[N0] = kd(axis,1) * kd(side,1) * (1-i)  +
                   kd(axis,1) * kd(side,2) * (nx+i) +
                   kd(axis,2) * kd(flip,0) *  j     +
                   kd(axis,2) * kd(flip,1) * (nx+1-j)
          I2[N0-1] = kd(axis,1) * kd(flip,0) *  j +
                     kd(axis,1) * kd(flip,1) * (ny+1-j) +
                     kd(axis,2) * kd(side,1) * (1-i) +
                     kd(axis,2) * kd(side,2) * (ny+i)
          I2[N0-2] = k
          if nfaces>1; I2[N0-3]=f1 end
          I1=I2
          i=I1[N0]; j=I1[N0-1]; k=I1[N0-2]
        end
      end
      if I1[N0]>nx
        if connections[f,1,2,1]==-1
          return NaN
        else
          f1=connections[f,1,2,1];
          axis=connections[f,1,2,2];
          side=connections[f,1,2,3];
          flip=connections[f,1,2,4];
          I2[N0] = kd(axis,1) * kd(side,1) * (i-nx)   +
                   kd(axis,1) * kd(side,2) * (nx-(i-nx-1)) +
                   kd(axis,2) * kd(flip,0) *  j     +
                   kd(axis,2) * kd(flip,1) * (nx+1-j)
          I2[N0-1] = kd(axis,1) * kd(flip,0) *  j +
                     kd(axis,1) * kd(flip,1) * (ny+1-j) +
                     kd(axis,2) * kd(side,1) * (i-ny) +
                     kd(axis,2) * kd(side,2) * ((ny+1)-(i-ny))
          I2[N0-2] = k
          if nfaces>1; I2[N0-3]=f1 end
          I1=I2
          i=I1[N0]; j=I1[N0-1]; k=I1[N0-2]
        end
      end

      if I1[N0-1]<1
        if connections[f,2,1,1]==-1
          return NaN
        else
          f1=connections[f,2,1,1]
          axis=connections[f,2,1,2];
          side=connections[f,2,1,3];
          flip=connections[f,2,1,4];
          I2[N0] = kd(axis,1) * kd(side,1) * (1-j) +
                   kd(axis,1) * kd(side,2) *  (nx+j) +
                   kd(axis,2) * kd(flip,0) * i +
                   kd(axis,2) * kd(flip,1) * (nx+1-i)
          I2[N0-1] = kd(axis,1) * kd(flip,0) *  i +
                     kd(axis,1) * kd(flip,1) * (ny+1-i) +
                     kd(axis,2) * kd(side,1) * (1-j) +
                     kd(axis,2) * kd(side,2) * (ny+j)
          I2[N0-2] = k
          if nfaces>1; I2[N0-3]=f1 end
          I1=I2
          i=I1[N0]; j=I1[N0-1]; k=I1[N0-2]
        end
      end
      if I1[N0-1]>ny
        if connections[f,2,2,1]==-1
          return NaN
        else
          f1=connections[f,2,2,1]
          axis=connections[f,2,2,2];
          side=connections[f,2,2,3];
          flip=connections[f,2,2,4];
          I2[N0] = kd(axis,1) * kd(side,1) * (j-ny)   +
                  kd(axis,1) * kd(side,2) * (ny+1-(j-ny)) +
                  kd(axis,2) * kd(flip,0) *  i     +
                  kd(axis,2) * kd(flip,1) * (nx+1-i)
          I2[N0-1] = kd(axis,1) * kd(flip,0) *  i +
                    kd(axis,1) * kd(flip,1) * (ny+1-i) +
                    kd(axis,2) * kd(side,1) * (j-ny) +
                    kd(axis,2) * kd(side,2) * (ny+1-(j-ny))
          I2[N0-2] = k
          if nfaces>1; I2[N0-3]=f1 end
          I1=I2
          i=I1[N0]; j=I1[N0-1]; k=I1[N0-2]
        end
      end

      if I1[N0-2]<1
        if connections[f,3,1,1]==-1
          return NaN
        else
          f1=connections[f,3,1,1]
          axis=connections[f,3,1,2];
          side=connections[f,3,1,3];
          flip=connections[f,3,1,4];
          I2[N0] = i
          I2[N0-1] = k
          I2[N0-2] = kd(axis,3) * kd(side,1) * (k-nz)   +
                    kd(axis,3) * kd(side,2) * (nz+1-(k-nz))
          if nfaces>1; I2[N0-3]=f1 end
          I1=I2
          i=I1[N0]; j=I1[N0-1]; k=I1[N0-2]
        end
      end
      if I1[N0-2]>nz
        if connections[f,3,2,1]==-1
          return NaN
        else
          f1=connections[f,3,2,1]
          axis=connections[f,3,2,2];
          side=connections[f,3,2,3];
          flip=connections[f,3,2,4];
          I2[N0] = i
          I2[N0-1] = j
          I2[N0-2] = kd(axis,3) * kd(side,1) * (k-nz)   +
                    kd(axis,3) * kd(side,2) * (nz+1-(k-nz))
          if nfaces>1; I2[N0-3]=f1 end
          I1=I2
          i=I1[N0]; j=I1[N0-1]; k=I1[N0-2]
        end
      end
    end
  # 2d array
  elseif nspatial==2
    nx=S[N0];ny=S[N0-1]
    while (I1[N0]<1 || I1[N0]>nx) || (I1[N0-1]<1 || I1[N0-1]>ny)
      I2=I1
      i=I1[N0]; j=I1[N0-1]
      if nfaces>1
        f=I1[N0-2]
      else
        f=1
      end
      if I1[N0]<1
        if connections[f,1,1,1]==-1
          return NaN
        else
          f1=connections[f,1,1,1]
          axis=connections[f,1,1,2];
          side=connections[f,1,1,3];
          flip=connections[f,1,1,4];
          I2[N0] = kd(axis,1) * kd(side,1) * (1-i)  +
                   kd(axis,1) * kd(side,2) * (nx+i) +
                   kd(axis,2) * kd(flip,0) *  j     +
                   kd(axis,2) * kd(flip,1) * (nx+1-j)
          I2[N0-1] = kd(axis,1) * kd(flip,0) *  j +
                     kd(axis,1) * kd(flip,1) * (ny+1-j) +
                     kd(axis,2) * kd(side,1) * (1-i) +
                     kd(axis,2) * kd(side,2) * (ny+i)
          if nfaces>1; I2[N0-2]=f1 end
          I1=I2
          i=I1[N0]; j=I1[N0-1]
        end
      end
      if I1[N0]>nx
        if connections[f,1,2,1]==-1
          return NaN
        else
          f1=connections[f,1,2,1];
          axis=connections[f,1,2,2];
          side=connections[f,1,2,3];
          flip=connections[f,1,2,4];
          I2[N0] = kd(axis,1) * kd(side,1) * (i-nx)   +
                   kd(axis,1) * kd(side,2) * (nx-(i-nx-1)) +
                   kd(axis,2) * kd(flip,0) *  j     +
                   kd(axis,2) * kd(flip,1) * (nx+1-j)
          I2[N0-1] = kd(axis,1) * kd(flip,0) *  j +
                     kd(axis,1) * kd(flip,1) * (ny+1-j) +
                     kd(axis,2) * kd(side,1) * (i-ny) +
                     kd(axis,2) * kd(side,2) * ((ny+1)-(i-ny))
          if nfaces>1; I2[N0-2]=f1 end
          I1=I2
          i=I1[N0]; j=I1[N0-1]
        end
      end

      if I1[N0-1]<1
        if connections[f,2,1,1]==-1
          return NaN
        else
          f1=connections[f,2,1,1]
          axis=connections[f,2,1,2];
          side=connections[f,2,1,3];
          flip=connections[f,2,1,4];
          I2[N0] = kd(axis,1) * kd(side,1) * (1-j) +
                   kd(axis,1) * kd(side,2) *  (nx+j) +
                   kd(axis,2) * kd(flip,0) * i +
                   kd(axis,2) * kd(flip,1) * (nx+1-i)
          I2[N0-1] = kd(axis,1) * kd(flip,0) *  i +
                     kd(axis,1) * kd(flip,1) * (ny+1-i) +
                     kd(axis,2) * kd(side,1) * (1-j) +
                     kd(axis,2) * kd(side,2) * (ny+j)
          if nfaces>1; I2[N0-2]=f1 end
          I1=I2
          i=I1[N0]; j=I1[N0-1]
        end
      end
      if I1[N0-1]>ny
        if connections[f,2,2,1]==-1
          return NaN
        else
          f1=connections[f,2,2,1]
          axis=connections[f,2,2,2];
          side=connections[f,2,2,3];
          flip=connections[f,2,2,4];
          I2[N0] = kd(axis,1) * kd(side,1) * (j-ny)   +
                  kd(axis,1) * kd(side,2) * (ny+1-(j-ny)) +
                  kd(axis,2) * kd(flip,0) *  i     +
                  kd(axis,2) * kd(flip,1) * (nx+1-i)
          I2[N0-1] = kd(axis,1) * kd(flip,0) *  i +
                    kd(axis,1) * kd(flip,1) * (ny+1-i) +
                    kd(axis,2) * kd(side,1) * (j-ny) +
                    kd(axis,2) * kd(side,2) * (ny+1-(j-ny))
          if nfaces>1; I2[N0-2]=f1 end
          I1=I2
          i=I1[N0]; j=I1[N0-1]
        end
      end
    end
  # 1d array
  elseif nspatial==1
    nx=size(A)[N0];
    while (I1[N0]<1 || I1[N0]>nx)
      I2=I1
      i=I1[N0]
      if nfaces>1
        f=I1[N0-1]
      else
        f=1
      end
      if I1[N0]<1
        if connections[f,1,1,1]==-1
          return NaN
        else
          f1=connections[f,1,1,1]
          axis=connections[f,1,1,2];
          side=connections[f,1,1,3];
          flip=connections[f,1,1,4];
          I2[N0] = kd(axis,1) * kd(side,1) * (1-i) +
                   kd(axis,1) * kd(side,2) * (nx+i)
          if nfaces>1; I2[N0-1]=f1 end
          I1=I2
          i=I1[N0]
        end
      end
      if I1[N0]>nx
        if connections[f,1,2,1]==-1
          return NaN
        else
          f1=connections[f,1,2,1];
          axis=connections[f,1,2,2];
          side=connections[f,1,2,3];
          flip=connections[f,1,2,4];
          I2[N0] = kd(axis,1) * kd(side,1) * (i-nx)   +
                  kd(axis,1) * kd(side,2) * ((nx+1)-(i-nx))
          if nfaces>1; I2[N0-1]=f1 end
          I1=I2
          i=I1[N0]
        end
      end
    end
  end
  return I1[N0-length(S)+1:N0]
end

function Base.getindex(A::CyclicArray, I::Vararg{Int, N}) where N # implements A[I]
  I1=CalcIndex(A,I...)
  if all(isnan.([I1[i] for i in N]))
    return NaN
  else
    return Base.getindex(A.data,I1...)
  end
end

function Base.setindex!(A::CyclicArray,value,I::Vararg{Int, N}) where N # A[I] = value
    I1=CalcIndex(A,I...)
  return Base.setindex!(A.data,value,I1...)
end

end

max_issues_repo_name: udistr/CircularArrays


text:
module PhysData

import CoolProp
import PhysicalConstants: CODATA2014
import Unitful: ustrip
import CSV
import DelimitedFiles: readdlm
import Polynomials
import Luna: Maths, Utils

include("data/lookup_tables.jl")

"Speed of light"
const c = ustrip(CODATA2014.SpeedOfLightInVacuum)
"Pressure in Pascal at standard conditions (atmospheric pressure)"
const atm = ustrip(CODATA2014.atm)
"Pressure in Pascal of 1 bar"
const bar = 100000
"Boltzmann constant"
const k_B = ustrip(CODATA2014.k_B)
"Permittivity of vacuum"
const ε_0 = ustrip(CODATA2014.ε_0)
"Permeability of vacuum"
const μ_0 = ustrip(CODATA2014.μ_0)
"Electron charge"
const electron = ustrip(CODATA2014.e)
"Electron mass"
const m_e = ustrip(CODATA2014.m_e)
"Ratio of electron charge squared to electron mass (for plasma)"
const e_ratio = electron^2/m_e
"Reduced Planck's constant"
const ħ = ustrip(CODATA2014.ħ)
"Atomic unit of energy"
const au_energy = ħ*c*ustrip(CODATA2014.α)/ustrip(CODATA2014.a_0)
"Atomic unit of time"
const au_time = ħ/au_energy
"Atomic unit of electric field"
const au_Efield = au_energy/(electron*ustrip(CODATA2014.a_0))
"Room temperature in Kelvin (ca 20 deg C)"
const roomtemp = 293.15
"Avogadro constant"
const N_A = ustrip(CODATA2014.N_A)
"Amagat (Loschmidt constant)"
const amg = atm/(k_B*273.15)
"Atomic mass unit"
const m_u = ustrip(CODATA2014.m_u)

const gas = (:Air, :He, :HeJ, :Ne, :Ar, :Kr, :Xe, :N2, :H2, :O2, :CH4, :SF6, :N2O)
const gas_str = Dict(
    :He => "He",
    :HeJ => "He",
    :Ar => "Ar",
    :Ne => "Neon",
    :Kr => "Krypton",
    :Xe => "Xenon",
    :Air => "Air",
    :N2 => "Nitrogen",
    :H2 => "Hydrogen",
    :O2 => "Oxygen",
    :CH4 => "Methane",
    :SF6 => "SulfurHexafluoride",
    :N2O => "NitrousOxide"
)
const glass = (:SiO2, :BK7, :KBr, :CaF2, :BaF2, :Si, :MgF2, :ADPo, :ADPe, :KDPo, :KDPe)
const metal = (:Ag,:Al)

"""
    wlfreq(ωλ)

Change from ω (angular frequency) to λ (wavelength) and vice versa
"""
wlfreq(ωλ) = 2π*c/ωλ

"""
    ΔλΔω(Δλ, λ)

Convert Δλ (wavelength bandwidth) at λ (central wavelength) to Δω (angular frequency bandwidth)
"""
ΔλΔω(Δλ, λ) = (2π*c)*Δλ/λ^2

eV_to_m(eV) = wlfreq(electron*eV/ħ)

"""
    γ_Börzsönyi(B1, C1, B2, C2)

Sellmeier expansion for linear susceptibility from Applied Optics 47, 27, 4856 (2008) at
room temperature and atmospheric pressure
"""
function γ_Börzsönyi(B1, C1, B2, C2)
    return μm -> (B1 * μm^2 / (μm^2 - C1) + B2 * μm^2 / (μm^2 - C2))
end

"""
    γ_JCT(B1, C1, B2, C2, B3, C3)

Adapted Sellmeier expansion for helium made to fit high frequency data
Phys. Rev. A 92, 033821 (2015)
"""
function γ_JCT(B1, C1, B2, C2, B3, C3)
    return μm -> (B1 * μm^2 / (μm^2 - C1)
                  + B2 * μm^2 / (μm^2 - C2)
                  + B3 * μm^2 / (μm^2 - C3))
end

"""
    γ_Peck(B1, C1, B2, C2, dens)

Sellmeier expansion for linear susceptibility from
J. Opt. Soc. Am. 67, 1550 (1977)
"""
function γ_Peck(B1, C1, B2, C2, dens)
    return μm -> @. (((B1 / (C1 - 1/μm^2) + B2 / (C2 - 1/μm^2)) + 1)^2 - 1)/dens
end

"""
    γ_Zhang(A, B, C, dens)

Sellmeier expansion for Oxygen from Applied Optics 50, 35, 6484 (2011)
"""
function γ_Zhang(A, B, C, dens)
    return μm -> ((1 + A + B/(C-1/μm^2))^2 - 1)/dens
end

"""
    γ_QuanfuHe(A, B, C, dens)

Sellmeier expansion for CH4, SF6 and N2O from Atmospheric Chemistry and Physics 2021, 21 (19), 14927–14940.
https://doi.org/10.5194/acp-21-14927-2021.

"""
function γ_QuanfuHe(A, B, C, dens)
    return μm -> ((1 + 1e-8*(A + B/(C - (1e4/μm)^2))))/dens
end

"""
    sellmeier_gas(material::Symbol)

Return function for linear polarisability γ, i.e. susceptibility of a single particle,
calculated from Sellmeier expansions.
"""
function sellmeier_gas(material::Symbol)
    dens = density(material, 1.0, 273.15)
    if material == :He
        B1 = 4977.77e-8
        C1 = 28.54e-6
        B2 = 1856.94e-8
        C2 = 7.76e-3
        return γ_Börzsönyi(B1/dens, C1, B2/dens, C2)
    elseif material == :HeJ
        B1 = 2.16463842e-05
        C1 = -6.80769781e-04
        B2 = 2.10561127e-07
        C2 = 5.13251289e-03
        B3 = 4.75092720e-05
        C3 = 3.18621354e-03
        return γ_JCT(B1/dens, C1, B2/dens, C2, B3/dens, C3)
    elseif material == :Ne
        B1 = 9154.48e-8
        C1 = 656.97e-6
        B2 = 4018.63e-8
        C2 = 5.728e-3
        return γ_Börzsönyi(B1/dens, C1, B2/dens, C2)
    elseif material == :Ar
        B1 = 20332.29e-8
        C1 = 206.12e-6
        B2 = 34458.31e-8
        C2 = 8.066e-3
        return γ_Börzsönyi(B1/dens, C1, B2/dens, C2)
    elseif material == :Kr
        B1 = 26102.88e-8
        C1 = 2.01e-6
        B2 = 56946.82e-8
        C2 = 10.043e-3
        return γ_Börzsönyi(B1/dens, C1, B2/dens, C2)
    elseif material == :Xe
        B1 = 103701.61e-8
        C1 = 12750e-6
        B2 = 31228.61e-8
        C2 = 0.561e-3
        return γ_Börzsönyi(B1/dens, C1, B2/dens, C2)
    elseif material == :Air
        B1 = 14926.44e-8
        C1 = 19.36e-6
        B2 = 41807.57e-8
        C2 = 7.434e-3
        return γ_Börzsönyi(B1/dens, C1, B2/dens, C2)
    elseif material == :N2
        B1 = 39209.95e-8
        C1 = 1146.24e-6
        B2 = 18806.48e-8
        C2 = 13.476e-3
        return γ_Börzsönyi(B1/dens, C1, B2/dens, C2)
    elseif material == :H2
        B1 = 14895.6e-6
        C1 = 180.7
        B2 = 4903.7e-6
        C2 = 92.0
        return γ_Peck(B1, C1, B2, C2, density(material, atm/bar, 273.15))
    elseif material == :O2
        # Applied Optics 50, 35, 6484 (2011)
        A = 1.181494e-4
        B = 9.708931e-3
        C = 75.4
        return γ_Zhang(A, B, C, density(material, atm/bar, roomtemp))
    elseif material == :CH4
        # Atmospheric Chemistry and Physics 2021, 21 (19), 14927–14940.
        A = 3603.09
        B = 4.40362e14
        C = 1.1741e10
        return  γ_QuanfuHe(A, B, C, density(material, atm/bar, 288.15))
    elseif material == :N2O
        # Atmospheric Chemistry and Physics 2021, 21 (19), 14927–14940.
        A = 22095.0
        B = 1.66291e14
        C = 6.75226e9
        return  γ_QuanfuHe(A, B, C, density(material, atm/bar, 288.15))
    elseif material == :SF6
        # Atmospheric Chemistry and Physics 2021, 21 (19), 14927–14940.
        A = 18997.7
        B = 8.27663e14
        C = 1.56833e10
        return  γ_QuanfuHe(A, B, C, density(material, atm/bar, 288.15))
    else
        throw(DomainError(material, "Unknown gas $material"))
    end
end

"""
    sellmeier_glass(material::Symbol)

Sellmeier for glasses. Returns function of wavelength in μm which in turn returns the
refractive index directly
"""
function sellmeier_glass(material::Symbol)
    if material == :SiO2
        #  J. Opt. Soc. Am. 55, 1205-1208 (1965)
        # TODO: Deal with sqrt of negative values better (somehow...)
        return μm -> sqrt(complex(1
             + 0.6961663/(1-(0.0684043/μm)^2)
             + 0.4079426/(1-(0.1162414/μm)^2)
             + 0.8974794/(1-(9.896161/μm)^2)
             ))
    elseif material == :BK7
        # ref index info (SCHOTT catalogue)
        return μm -> sqrt(complex(1
             + 1.03961212/(1-0.00600069867/μm^2)
             + 0.231792344 / (1-0.0200179144/μm^2)
             + 1.01046945/(1-103.560653/μm^2)
             ))
    elseif material == :CaF2
        # Appl. Opt. 41, 5275-5281 (2002)
        return μm -> sqrt(complex(1
             + 0.443749998/(1-0.00178027854/μm^2)
             + 0.444930066/(1-0.00788536061/μm^2)
             + 0.150133991/(1-0.0124119491/μm^2)
             + 8.85319946/(1-2752.28175/μm^2)
             ))
    elseif material == :KBr
        # J. Phys. Chem. Ref. Data 5, 329-528 (1976)
        return μm -> sqrt(complex(1
             + 0.39408
             + 0.79221/(1-(0.146/μm)^2)
             + 0.01981/(1-(0.173/μm)^2)
             + 0.15587/(1-(0.187/μm)^2)
             + 0.17673/(1-(60.61/μm)^2)
             + 2.06217/(1-(87.72/μm)^2)
             ))
    elseif material == :BaF2
        # J. Phys. Chem. Ref. Data 9, 161-289 (1980)
        return μm -> sqrt(complex(1
             + 0.33973
             + 0.81070/(1-(0.10065/μm)^2)
             + 0.19652/(1-(29.87/μm)^2)
             + 4.52469/(1-(53.82/μm)^2)
             ))
    elseif material == :Si
        # J. Opt. Soc. Am., 47, 244-246 (1957)
        return μm -> sqrt(complex(1
             + 10.6684293/(1-(0.301516485/μm)^2)
             + 0.0030434748/(1-(1.13475115/μm)^2)
             + 1.54133408/(1-(1104/μm)^2)
             ))
    elseif material == :MgF2
        return μm -> @. sqrt(complex(1
            + 0.27620
            + 0.60967/(1-(0.08636/μm)^2)
            + 0.0080/(1-(18.0/μm)^2)
            + 2.14973/(1-(25.0/μm)^2)
            ))
    elseif material == :ADPo
        return μm -> @. sqrt(complex(
            2.302842
            + 15.102464*μm^2/(μm^2-400)
            + 0.011125165/(μm^2-0.01325366)
        ))
    elseif material == :ADPe
        return μm -> @. sqrt(complex(
            2.163510
            + 5.919896*μm^2/(μm^2-400)
            + 0.009616676/(μm^2-0.01298912)
        ))
    elseif material == :KDPo
        return μm -> @. sqrt(complex(
            2.259276
            + 13.00522*μm^2/(μm^2-400)
            + 0.01008956/(μm^2-0.0129426)
        ))
    elseif material == :KDPe
        return μm -> @. sqrt(complex(
            2.132668
            + 3.2279924*μm^2/(μm^2-400)
            + 0.008637494/(μm^2-0.0122810)
        ))
    else
        throw(DomainError(material, "Unknown glass $material"))
    end
end

"""
    sellmeier_crystal(material, axis)

Sellmeier for crystals. Returns function of wavelength in μm which in turn returns the
refractive index directly. Possible values for `axis` depend on the type of crystal.
"""
function sellmeier_crystal(material, axis)
    if material == :BBO
        if axis == :o
            return μm -> sqrt(complex(1
                + 0.90291/(1-0.003926/μm^2)
                + 0.83155/(1-0.018786/μm^2)
                + 0.76536/(1-60.01/μm^2)
                ))
        elseif axis == :e
            return μm -> sqrt(complex(1
                + 1.151075/(1-0.007142/μm^2)
                + 0.21803/(1-0.02259/μm^2)
                + 0.656/(1-263/μm^2)
                ))
        else
            throw(DomainError(axis, "Unknown BBO axis $axis"))
        end
    elseif material == :LBO
        # C Chen et al., J Opt. Soc. Am. 6, 616-621 (1989)
        # F. Hanson and D. Dick., Opt. Lett. 16, 205-207 (1991).
        if axis == :x
            return μm -> sqrt(complex(
                2.45768
                + 0.0098877/(μm^2-0.026095)
                - 0.013847*μm^2
            ))
        elseif axis == :y
            return μm -> sqrt(complex(
                2.52500
                + 0.017123/(μm^2+0.0060517)
                - 0.0087838*μm^2
            ))
        elseif axis == :z
            return μm -> sqrt(complex(
                2.58488
                + 0.012737/(μm^2-0.021414)
                - 0.016293*μm^2
            ))
        else
            throw(DomainError(axis, "Unknown LBO axis $axis"))
        end
    else
        throw(DomainError(material, "Unknown crystal $material"))
    end
end

function ref_index_fun_uniax(material; axes=(:o, :e))
    n_o = sellmeier_crystal(material, axes[1])
    n_e = sellmeier_crystal(material, axes[2])
    n(λ, θ) = sqrt(1/((cos(θ)/n_o(λ*1e6))^2+(sin(θ)/n_e(λ*1e6))^2))
    return n
end

"""
    χ1_fun(gas::Symbol)

Get function to return χ1 (linear susceptibility) for gases as a function of
wavelength in SI units, pressure in bar, and temperature in Kelvin.
"""
function χ1_fun(gas::Symbol)
    γ = sellmeier_gas(gas)
    f = let γ=γ, gas=gas
        function χ1(λ, P, T)
            γ(λ*1e6)*density(gas, P, T)
        end
    end
    return f
end

function χ1_fun(gas::Symbol, P, T)
    γ = sellmeier_gas(gas)
    dens = density(gas, P, T)
    return λ -> γ(λ*1e6)*dens
end

"""
    χ1(gas::Symbol, λ, P=1.0, T=roomtemp)

Calculate χ1 at wavelength λ in SI units, pressure P in bar and temperature T in Kelvin.
Gases only.
"""
function χ1(gas::Symbol, λ, P=1.0, T=roomtemp)
    return χ1_fun(gas)(λ, P, T)
end

"""
    ref_index(material, λ, P=1.0, T=roomtemp; lookup=nothing)

Get refractive index for any material at wavelength given in SI units.
"""
function ref_index(material, λ, P=1.0, T=roomtemp; lookup=nothing)
    return ref_index_fun(material, P, T; lookup=lookup)(λ)
end

"""
    ref_index_fun(material::Symbol, P=1.0, T=roomtemp; lookup=nothing)

Get function which returns refractive index.
"""
function ref_index_fun(material::Symbol, P=1.0, T=roomtemp; lookup=nothing)
    if material in gas
        χ1 = χ1_fun(material, P, T)
        return λ -> sqrt(1 + χ1(λ))
    elseif material in glass
        if isnothing(lookup)
            lookup = (material == :SiO2)
        end
        if lookup
            spl = lookup_glass(material)
            return λ -> spl(λ*1e6)
        else
            sell = sellmeier_glass(material)
            return λ -> sell(λ*1e6)
        end
    elseif material in metal
        nmetal = let spl = lookup_metal(material)
            function nmetal(λ)
                if λ > 1
                    throw(DomainError(λ, "Wavelength must be given in metres"))
                end
                return spl(λ*1e6)
            end
        end
        return nmetal
    else
        throw(DomainError(material, "Unknown material $material"))
    end
end

"""
    ref_index_fun(gases, P, T=roomtemp; lookup=nothing)

Get function which returns refractive index for gas mixture. `gases` is a `Tuple` of gas
identifiers (`Symbol`s) and `P` is a `Tuple` of equal length containing pressures.
"""
function ref_index_fun(gases::NTuple{N, Symbol}, P::NTuple{N, Number}, T=roomtemp; lookup=nothing) where N
    ngas = let funs=[χ1_fun(gi, Pi, T) for (gi, Pi) in zip(gases, P)]
        function ngas(λ)
            res = funs[1](λ)
            for ii in 2:length(gases)
                res += funs[ii](λ)
            end
            return sqrt(1 + res)
        end
    end
    return ngas
end

"""
    ref_index_fun(gases, T=roomtemp)

Get function which returns ref index for mixture as function of wavelength and densities.
"""
function ref_index_fun(gases::NTuple{N, Symbol}, T=roomtemp) where N
    let γs=[sellmeier_gas(gi) for gi in gases]
        function ngas(λ, densities::Vector{<:Number})
            χ1 = 0.0
            for (γi, di) in zip(γs, densities)
                χ1 += di * γi(λ*1e6)
            end
            sqrt(1 + χ1)
        end
    end
end

"""
    dispersion_func(order, n)

Get a function that calculates dispersion of order `order` for a refractive index given by
`n(λ)`.
"""
function dispersion_func(order, n)
    β(ω) = @. ω/c * real(n(wlfreq(ω)))
    βn(λ) = Maths.derivative(β, wlfreq(λ), order)
    return βn
end

"""
    dispersion_func(order, material, P=1, T=roomtemp; lookup=nothing)

Get a function to calculate dispersion. Arguments are the same as for [`dispersion`](@ref).
"""
function dispersion_func(order, material::Symbol, P=1.0, T=roomtemp; lookup=nothing)
    n = ref_index_fun(material, P, T, lookup=lookup)
    dispersion_func(order, n)
end

"""
    dispersion(order, material, λ, P=1.0, T=roomtemp; lookup=nothing)

Calculate the dispersion of order `order` of a given `material` at a wavelength `λ`.

For gases the pressure `P` (default:atmosphere) and the temperature `T` (default: room temp)
can also be specified. `lookup::Bool` determines whether a lookup table or a Sellmeier
expansion is used for the refractive index (default is material dependent).

# Examples
```jldoctest
julia> dispersion(2, :BK7, 400e-9) * 1e30 * 1e-3 # convert to fs^2/mm
122.03632107303108
```
"""
function dispersion(order, material::Symbol, λ, P=1.0, T=roomtemp; lookup=nothing)
    return dispersion_func(order, material, P, T; lookup=lookup).(λ)
end

"""
    fresnel(n2, θi; n1=1.0)

Calcualte reflection coefficients from Fresnel's equations.
"""
function fresnel(n2, θi; n1=1.0)
    θt = asin(n1*sin(θi)/n2)
    rs = (n1*cos(θi) - n2*cos(θt))/(n1*cos(θi) + n2*cos(θt))
    rp = (n2*cos(θi) - n1*cos(θt))/(n2*cos(θi) + n1*cos(θt))
    abs2(rs), angle(rs), abs2(rp), angle(rp)
end

"""
    γ3_gas(material::Symbol; source=nothing)

Calculate single-molecule third-order hyperpolarisability of a gas at given wavelength(s)
and at room temperature.
If `source` == `:Bishop`: Uses reference values to calculate γ
If `source` == `:Lehmeier` (default): Uses scaling factors to calculate χ3 at 1 bar and
scales by density to get to a single molecule i.e. the hyperpolarisability

References:
[1] Journal of Chemical Physics, AIP, 91, 3549-3551 (1989)
[2] Chemical Reviews, 94, 3-29 (1994)
[3] Optics Communications, 56(1), 67–72 (1985)
[4] Phys. Rev. A, vol. 42, 2578 (1990)
[5] Optics Letters Vol. 40, No. 24 (2015))
[6] Phys. Rev. A 2012, 85 (4), 043820. https://doi.org/10.1103/PhysRevA.85.043820.
"""
function γ3_gas(material::Symbol; source=nothing)
    # TODO: More Bishop/Shelton; Wahlstrand updated values.
    if source === nothing
        if material in (:He, :HeJ, :Ne, :Ar, :Kr, :Xe, :N2)
            source = :Lehmeier
        elseif material in (:H2, :CH4, :SF6)
            source = :Shelton
        elseif material in (:O2,)
            source = :Zahedpour
        elseif material in (:N2O,)
            source = :Wahlstrand
        else
            error("no default γ3 source for material: $material")
        end
    end
    dens = density(material, atm/bar, 273.15)
    if source == :Lehmeier
        # Table 1 in [3]
        if material in (:He, :HeJ)
            fac = 1
        elseif material == :Ne
            fac = 1.8
        elseif material == :Ar
            fac = 23.5
        elseif material == :Kr
            fac = 64.0
        elseif material == :Xe
            fac = 188.2
        elseif material == :N2
            fac = 21.1
        else
            throw(DomainError(material, "Lehmeier model does not include $material"))
        end
        return 4*fac*3.43e-28 / dens
    elseif source == :Shelton
        # ref [4], we use Table 1 to simply scale from
        # the paired gas (for which we use Lehmeier)
        # e.g. He for H2, N2 for CH4 or SF6
        if material == :H2
            return 15.77*γ3_gas(:He)
        elseif material == :CH4
            return 2.931*γ3_gas(:N2)
        elseif material == :SF6
            return 1.53*γ3_gas(:N2)
        else
            throw(DomainError(material, "Shelton model does not include $material"))
        end
    elseif source == :Zahedpour
        if material == :O2
            n0 = ref_index(:O2, 800e-9, atm/bar, roomtemp)
            ρ = density(:O2, atm/bar, roomtemp)
            n2 = 8.1e-24 # Table 1 in [5]
            return 4/3*ε_0*c*n0^2/ρ * n2
        else
            throw(DomainError(material, "Zahedpour model does not include $material"))
        end
    elseif source == :Wahlstrand
        if material == :N2O
            n0 = ref_index(:N2O, 800e-9, atm/bar, roomtemp)
            ρ = density(:N2O, atm/bar, roomtemp)
            n2 = 17.2e-24 # Table 1 in [6]
            return 4/3*ε_0*c*n0^2/ρ * n2
        else
            throw(DomainError(material, "Wahlstrand model does not include $material"))
        end
    else
        throw(DomainError(source, "Unkown γ3 model $source"))
    end
end

function χ3_gas(material::Symbol, P, T=roomtemp; source=nothing)
    return γ3_gas(material, source=source) .* density.(material, P, T)
end

function n2_gas(material::Symbol, P, T=roomtemp, λ=800e-9; source=nothing)
    n0 = ref_index(material, λ, P, T)
    return @. 3/4 * χ3_gas(material, P, T, source=source) / (ε_0*c*n0^2)
end

"""
    density(gas::Symbol, P, T=roomtemp)

Number density of `gas` [m^-3] at pressure `P` [bar] and temperature `T` [K].
"""
function density(gas::Symbol, P, T=roomtemp)
    P == 0 ? zero(P) : CoolProp.PropsSI("DMOLAR", "T", T, "P", bar*P, gas_str[gas])*N_A
end

"""
    pressure(gas, density, T=roomtemp)

Calculate the pressure in bar of the `gas` at number density `density` and temperature `T`.
"""
function pressure(gas, density, T=roomtemp)
    density == 0 ? zero(density) :
                   CoolProp.PropsSI("P", "T", T, "DMOLAR", density/N_A, gas_str[gas])/bar
end

"""
    densityspline(gas; Pmax, Pmin=0, N=2^10, T=roomtemp)

Create a `CSpline` interpolant for the density of the `gas` between pressures `Pmin` and
`Pmax` at temperature `T`. The spline is created using `N` samples.
"""
function densityspline(gas::Symbol; Pmax, Pmin=0, N=2^10, T=roomtemp)
    P = collect(range(Pmin, Pmax, length=N))
    ρ = density.(gas, P, T)
    Maths.CSpline(P, ρ)
end

"""
    ionisation_potential(material; unit=:SI)

Return the first ionisation potential of the `material` in a specific unit (default: SI).
Possible units are `:SI`, `:atomic` and `:eV`.
"""
function ionisation_potential(material; unit=:SI)
    if material in (:He, :HeJ)
        Ip = 0.9036
    elseif material == :Ne
        Ip = 0.7925
    elseif material == :Ar
        Ip = 0.5792
    elseif material == :Kr
        Ip = 0.5142
    elseif material == :Xe
        Ip = 0.4458
    elseif material == :H
        Ip = 0.5
    elseif material == :N2
        Ip = 0.5726
    elseif material == :H2
        Ip = 0.5669
    elseif material == :O2
        Ip = 0.443553
    elseif material == :CH4
        Ip = 0.4636
    elseif material == :N2O
        Ip = 0.474
    elseif material == :SF6
        Ip = 0.5
    else
        throw(DomainError(material, "Unknown material $material"))
    end

    if unit == :atomic
        return Ip
    elseif unit == :eV
        return @. 27.21138602*Ip
    elseif unit == :SI
        return @. 27.21138602*electron*Ip
    else
        throw(DomainError(unit, "Unknown unit $unit"))
    end
end

"""
    quantum_numbers(material)

Return the quantum numbers of the `material` for use in the PPT ionisation rate.
"""
function quantum_numbers(material)
    # Returns n, l, ion Z
    if material == :Ar
        return 3, 1, 1
    elseif material == :Ne
        return 2, 1, 1;
    elseif material == :Kr
        return 4, 1, 1
    elseif material == :Xe
        return 5, 1, 1
    elseif material in (:He, :HeJ)
        return 1, 0, 1
    elseif material == :O2
        return 2, 0, 0.53 # https://doi.org/10.1016/S0030-4018(99)00113-3
    elseif material == :N2
        return 2, 0, 0.9 # https://doi.org/10.1016/S0030-4018(99)00113-3
    else
        throw(DomainError(material, "Unknown material $material"))
    end
end

"""
    lookup_glass(material::Symbol)

Create a `CSpline` interpolant for look-up-table values of the refractive index.
"""
function lookup_glass(material::Symbol)
    if material == :SiO2
        data = data_glass(:SiO2)
        spl = Maths.BSpline(1e6*eV_to_m.(data[:, 1]), data[:, 2] + 1im * data[:, 3])
    else
        throw(DomainError(material, "Unknown metal $material"))
    end
    return spl
end

"""
    lookup_metal(material::Symbol)

Create a `CSpline` interpolant for look-up-table values of the refractive index.
"""
function lookup_metal(material::Symbol)
    data = data_metal(material)::Array{Float64,2}
    nspl = Maths.BSpline(data[:,1], data[:,2] .+ im.*data[:,3])
    return nspl
end

"""
    raman_parameters(material)

Get the Raman parameters for `material`.

# Fields
Fields in the returned named tuple must include:
- `kind::Symbol`: one of `:molecular` or ...

If `kind == :molecular` then the following must also be specified:
- `rotation::Symbol`: only `:nonrigid` or `:none` supported at present.
- `vibration::Symbol`: only `:sdo` or `:none` supported at present.

If `rotation == :nonrigid` then the following must also be specified:
- `B::Real`: the rotational constant [1/m]
- `Δα::Real`: molecular polarizability anisotropy [m^3]
- `qJodd::Integer`: nuclear spin parameter for odd `J`
- `qJeven::Integer`: nuclear spin parameter for even `J`
- `D::Real=0.0`: centrifugal constant [1/m]
Along with one of:
- `τ2r::Real`: coherence time [s]
- `Bρr::Real` : density dependent broadening coefficient [Hz/amagat]
If both `τ2r` and `Bρr` are specified, then `Bρr` takes precedence.
If `Bρr` is specified then we also need:
- `Aρr::Real` : self diffusion coefficient [Hz amagat]

If `vibration == :sdo` then the following must also be specified:
- `Ωv::Real`: vibrational frequency [rad/s]
- `dαdQ::Real`: isotropic averaged polarizability derivative [m^2]
- `μ::Real`: reduced molecular mass [kg]
Along with one of:
- `τ2v::Real`: coherence time [s]
- `Bρv::Real` : density dependent broadening coefficient [Hz/amagat]
If both `τ2v` and `Bρv` are specified, then `Bρv` takes precedence.
If `Bρv` is specified then we also need:
- `Aρv::Real` : self diffusion coefficient [Hz amagat]

# References
[1] Phys. Rev. A, 94, 023816 (2016)
[2] Phys. Rev. A, 85, 043820 (2012)
[3] Phys. Rev. A, 92, 063828 (2015)
[4] Journal of Raman Spectroscopy 2, 133 (1974)
[5] J. Phys. Chem., 91, 41 (1987)
[6] Applied Spectroscopy 23, 211 (1969)
[7] Phys. Rev. A, 34, 3, 1944 (1986)
[8] Can. J. Phys., 44, 4, 797 (1966)
[9] G. V. MIKHAtLOV, SOVIET PHYSICS JETP, vol. 36, no. 9, (1959).
[10] Phys. Rev. A, 33, 5, 3113 (1986)
[11] IEEE Journal of Quantum Electronics 1986, 22 (2), 332–336. https://doi.org/10.1109/JQE.1986.1072945.
[12] Phys. Rev. Lett. 1998, 81 (6), 1215–1218. https://doi.org/10.1103/PhysRevLett.81.1215.
[13] Optics Communications 1987, 64 (4), 393–397. https://doi.org/10.1016/0030-4018(87)90258-6.
[14] Science Advances 2020, 6 (34), eabb5375. https://doi.org/10.1126/sciadv.abb5375.
[15] Long, The Raman Effect; John Wiley & Sons, Ltd, 2002;
"""
function raman_parameters(material)
    if material == :N2
        rp = (kind = :molecular,
              rotation = :nonrigid,
              vibration = :sdo,
              B = 199.0, # [4]
              D = 5.74e-4, # [4]
              qJodd = 1, # [14] uses 3 as does [15]
              qJeven = 2, # [14] uses 6 as does [15]
              Δα = 6.7e-31, # [2] # note [14] use 1.86e-30
              # Bρr has a moderate dependence on J, which we ignore for now, taking J=8
              # [8] measured Bρr from 7 to 43 atm to be ~80e-3 cm^-1/atm,
              # which is translated to Hz/amg via
              # 80e-3*29979245800.0/(density(:N2, atm/bar)/amg)
              # giving ~2.6e9 Hz/amagat.
              # [7] gives ~3.3e9 Hz/amagat (measured at lower pressures), but they claim
              # their results are more accurate (of course!)
              Bρr = 3.3e9, # [7]
              Aρr = 0.0, # [7]
              dαdQ = 1.75e-20, # [6]
              Ωv = 2*π*2330.0*100.0*c, # [4]
              μ = (m_u*14.0067)^2/(2*m_u*14.0067),
              # For τ2v, [9] suggests pressure dependence is extremely weak up to 120 bar
              # [9] gives ~ 1.8 cm^-1, whereas Fig. 1 in [5] suggests something similar.
              # 1.8 cm^-1 = 0.054 THz
              # This gives a τ2v = 1/πΔν ~ 6 ps
              τ2v = 6e-12, # [5,9]
              )
    elseif material == :H2
        rp = (kind = :molecular,
              rotation = :nonrigid,
              vibration = :sdo,
              B = 5890.0, # [3]
              D = 5.0, # [3]
              qJodd = 3,
              qJeven = 1,
              Δα = 3e-31, # [3]
              Bρr = 114e6, # [7]
              Aρr = 6.15e6, # [7]
              dαdQ = 1.3e-20, # [3]
              Ωv = 2*π*124.5669e12,
              μ = (m_u*1.00784)^2/(2*m_u*1.00784),
              Bρv = 52.2e6, # [10]
              Aρv = 309e6, # [10]
              )
    elseif material == :D2
        rp = (kind = :molecular,
              rotation = :nonrigid,
              vibration = :sdo,
              B = 2930.0, # [3]
              D = 2.1, # [3]
              qJodd = 1,
              qJeven = 2,
              Δα = 3e-31, # [3]
              # TODO τ2r =
              dαdQ = 1.4e-20, # [3]
              # TODO Ωv =
              # TODO μ =
              # TODO τ2v =
              )
    elseif material == :O2
        rp = (kind = :molecular,
              rotation = :nonrigid,
              vibration = :sdo,
              B = 144.0, # [2]
              D = 0.0, # TODO
              qJodd = 1,
              qJeven = 0,
              Δα = 10.2e-31, # [2]
              # TODO τ2r =
              dαdQ = 1.46e-20, # [1]
              Ωv = 3e14, # [1]
              μ = 1.3e-26, # [1]
              # TODO τ2v =
              )
    elseif material == :N2O
        rp = (kind = :molecular,
              rotation = :nonrigid,
              vibration = :none, # TODO work out correct parameters here
              B = 41.0, # [2]
              D = 0.0, # TODO
              qJodd = 1, # [14]
              qJeven = 1, # [14]
              Δα = 28.1e-31, # [2] note that [14] uses twice this
              τ2r = 23.8e-12, # [14]
              # TODO dαdQ =
              Ωv = 2*π*1285*100.0*c,
              # TODO μ =
              # TODO τ2v =
             )
    elseif material == :CH4
        rp = (kind = :molecular,
              rotation = :none,
              vibration = :sdo,
              dαdQ = 1.04e-20, # [6]
              Ωv = 2*π*2914*100.0*c, # [6]
              μ = (1.00784*m_u)/4, #(m_u*12.0107*m_u*1.00784)/(m_u*12.0107 + m_u*1.00784),
              # TODO USE: (1) Taira, Y.; Ide, K.; Takuma, H. Accurate Measurement of the Pressure Broadening of the Ν1 Raman Line of CH4 in the 1–50 Atm Region by Inverse Raman Spectroscopy. Chemical Physics Letters 1982, 91 (4), 299–302. https://doi.org/10.1016/0009-2614(82)80160-7.
              #Bρv = X, # []
              #Aρv = X, # []
              τ2v = 28e-12, # [13]
             )
    elseif material == :SF6
        rp = (kind = :molecular,
                rotation = :none,
                vibration = :sdo,
                dαdQ = 1.23e-20, # [6]
                Ωv = 2*π*775*100.0*c, # [6]
                μ = (18.998403*m_u)/6,
                τ2v = 6.6e-12, # [13]
                )
    else
        throw(DomainError(material, "Unknown material $material"))
    end
    rp
end

"""
    lookup_mirror(type)

Create a `CSpline` interpolant for the complex-valued reflectivity of a mirror of `type`.
"""
function lookup_mirror(type)
    if type == :PC70
        # λ (nm), R(5deg) (%), R(19deg) (%)
        Rdat = CSV.File(joinpath(Utils.datadir(), "PC70_R.csv"))
        # λ (nm), GDD(5deg) (fs^2), GDD(19deg) (fs^2)
        GDDdat = CSV.File(joinpath(Utils.datadir(), "PC70_GDD.csv"))
        # Double sqrt creates average reflectivity per _reflection_ rather than per pair
        rspl = Maths.BSpline(Rdat.Wlgth*1e-9, sqrt.(sqrt.(Rdat.Rp5deg/100 .* Rdat.Rp19deg/100)))
        λGDD = GDDdat.Wlgth
        ω = wlfreq.(λGDD*1e-9)
        # average phase per _reflection_ rather than per pair
        ϕ = 1e-30/2 * Maths.cumtrapz(Maths.cumtrapz(GDDdat.GDDrp5deg.+GDDdat.GDDrp19deg, ω), ω)
        # ϕ has a large linear component - remove that
        ωfs = ω*1e-15
        ωfs0 = wlfreq(800e-9)*1e-15
        idcs =  2 .< ωfs .< 4 # large kinks at edge of frequency window confuse the fit
        p = Polynomials.fit(ωfs[idcs] .- ωfs0, ϕ[idcs], 5)
        p[2:end] = 0 # polynomials use 0-based indexing - only use constant and linear term
        ϕ .-= p.(ωfs .- ωfs0) # subtract linear part
        ϕspl = Maths.BSpline(λGDD*1e-9, ϕ)
        return λ -> rspl(λ) * exp(-1im*ϕspl(λ)) * Maths.planck_taper(
            λ, 400e-9, 450e-9, 1200e-9, 1300e-9)
    elseif type == :HD59
        dat = readdlm(joinpath(Utils.datadir(), "HD59_GDD.dat"); skipstart=2)
        λ = dat[:, 1] * 1e-9
        ω = wlfreq.(λ)
        GDD = dat[:, 2] .* 1e-30
        ϕ = Maths.cumtrapz(Maths.cumtrapz(GDD, ω), ω)
        ωfs = ω*1e-15
        ωfs0 = wlfreq(1030e-9)*1e-15
        p = Polynomials.fit(ωfs .- ωfs0, ϕ, 5)
        p[2:end] = 0 # polynomials use 0-based indexing - only use constant and linear term
        ϕ .-= p.(ωfs .- ωfs0) # subtract linear part
        ϕspl = Maths.BSpline(λ, ϕ)
        return λ -> exp(-1im*ϕspl(λ)) * Maths.planck_taper(
            λ, 993e-9, 1000e-9, 1060e-9, 1075e-9)
    elseif type == :PC147
        dat = readdlm(joinpath(Utils.datadir(), "PC147.txt"); skipstart=1)
        λR = dat[:, 1] * 1e-9
        R = dat[:, 2] # average reflectivity per mirror (complementary pair)
        rspl = Maths.BSpline(λR, sqrt.(R/100))
        λGDD = dat[:, 3] * 1e-9
        ω = wlfreq.(λGDD)
        GDD = dat[:, 4] .* 1e-30 # average GDD per mirror (complementary pair)
        ϕ = Maths.cumtrapz(Maths.cumtrapz(GDD, ω), ω)
        ωfs = ω*1e-15
        ωfs0 = wlfreq(1030e-9)*1e-15
        p = Polynomials.fit(ωfs .- ωfs0, ϕ, 5)
        p[2:end] = 0 # polynomials use 0-based indexing - only use constant and linear term
        ϕ .-= p.(ωfs .- ωfs0) # subtract linear part
        ϕspl = Maths.BSpline(λGDD, ϕ)
        return λ -> rspl(λ) * exp(-1im*ϕspl(λ)) * Maths.planck_taper(
            λ, 640e-9, 650e-9, 1350e-9, 1360e-9)
    elseif type == :HD120
        dat = readdlm(joinpath(Utils.datadir(), "HD120.csv"), ','; skipstart=1)
        λR = dat[:, 1] * 1e-9
        R = dat[:, 2] # reflectivity per mirror
        rspl = Maths.BSpline(λR, sqrt.(R/100))
        λGDD = dat[:, 3] * 1e-9
        ω = wlfreq.(λGDD)
        GDD = dat[:, 4] .* 1e-30 # GDD per mirror
        ϕ = Maths.cumtrapz(Maths.cumtrapz(GDD, ω), ω)
        ωfs = ω*1e-15
        ωfs0 = wlfreq(1030e-9)*1e-15
        p = Polynomials.fit(ωfs .- ωfs0, ϕ, 5)
        p[2:end] = 0 # polynomials use 0-based indexing - only use constant and linear term
        ϕ .-= p.(ωfs .- ωfs0) # subtract linear part
        ϕspl = Maths.BSpline(λGDD, ϕ)
        return λ -> rspl(λ) * exp(-1im*ϕspl(λ)) * Maths.planck_taper(
            λ, 880e-9, 900e-9, 1200e-9, 1220e-9)
    elseif type == :ThorlabsUMC
        # λ (nm), R(p) (%), R(s) (%)
        Rdat = CSV.File(joinpath(Utils.datadir(), "UCxx-15FS_R.csv"))
        # λ (nm), GD(p) (fs^2), GD(s) (fs^2)
        GDdat = CSV.File(joinpath(Utils.datadir(), "UCxx-15FS_GD.csv"))
        # Default to s-pol
        rspl = Maths.BSpline(Rdat.wl*1e-9, sqrt.(Rdat.Rs/100))
        λGD = GDdat.wl
        ω = wlfreq.(λGD*1e-9)
        # average phase per reflection, default to s-pol
        ϕ = Maths.cumtrapz(1e-15*GDdat.GDs, ω)
        # ϕ has a large linear component - remove that
        ωfs = ω*1e-15
        ωfs0 = wlfreq(800e-9)*1e-15
        idcs =  2 .< ωfs .< 4 # large kinks at edge of frequency window confuse the fit
        p = Polynomials.fit(ωfs[idcs] .- ωfs0, ϕ[idcs], 3)
        p[2:end] = 0 # polynomials use 0-based indexing - only use constant and linear term
        ϕ .-= p.(ωfs .- ωfs0) # subtract linear part
        ϕspl = Maths.BSpline(λGD*1e-9, ϕ)
        return λ -> rspl(λ) * exp(-1im*ϕspl(λ)) * Maths.planck_taper(
            λ, 640e-9, 650e-9, 1050e-9, 1100e-9)
    else
        throw(DomainError("Unknown mirror type $type"))
    end
end

end