text:
#=

This file contains functions that set the quadrature streams for the RT model.

So far, there are the following quadrature types:

    - GaussQuadHemisphere
    - GaussQuadFullSphere
    - RadauQuad

=#

"""
$(FUNCTIONNAME)(::GaussQuadHemisphere,
                        Ltrunc::Int,
                        obs_geom::ObsGeometry,
                        pol_type,
                        arr_type)

Computes hemispheric quadrature points with Gauss method within [0-1] interval
SZA and all VZAs included with 0 weight
Returns computed quadrature points as [`QuadPoints`](@ref)
"""
function rt_set_streams(::GaussQuadHemisphere,
                        Ltrunc::Int,
                        obs_geom::ObsGeometry{FT},
                        pol_type,
                        arr_type) where {FT}

    @unpack sza, vza = obs_geom

    Nquad = (Ltrunc + 1) ÷ 2

    qp_μ, wt_μ = Scattering.gauleg(Nquad, 0.0, 1.0) # quadrature limits are 0.0-1.0

    # Solar angle
    μ₀ = cosd.(sza)

    # Assign zero-weights to remaining camera zenith angles
    qp_μ = unique(FT[qp_μ; cosd.(vza); μ₀]);
    n_eff = length(qp_μ) - length(wt_μ);
    wt_μ = FT[wt_μ; zeros(FT,n_eff)];
    Nquad = length(qp_μ);
    iμ₀ = nearest_point(qp_μ, μ₀);
    qp_μN = arr_type(reduce(vcat, (fill.(qp_μ, [pol_type.n]))))
    wt_μN = arr_type(reduce(vcat, (fill.(wt_μ, [pol_type.n]))))
    i_start = pol_type.n*(iμ₀-1) + 1
    return QuadPoints(μ₀, iμ₀, i_start, arr_type(qp_μ), arr_type(wt_μ), qp_μN, wt_μN, Nquad)
end

"""
$(FUNCTIONNAME)(::GaussQuadFullSphere,
                        Ltrunc::Int,
                        obs_geom::ObsGeometry,
                        pol_type,
                        arr_type)

Computes hemispheric quadrature points with Gauss method within [-1,1] to generate 2Nquad weigths and take the one from [0,1]
SZA and all VZAs included with 0 weight
Returns computed quadrature points as [`QuadPoints`](@ref)
"""
function rt_set_streams(::GaussQuadFullSphere,
                        Ltrunc::Int,
                        obs_geom::ObsGeometry{FT},
                        pol_type,
                        arr_type) where {FT}

    @unpack sza, vza = obs_geom

    Nquad = (Ltrunc + 1) ÷ 2
    qp_μ, wt_μ = gausslegendre(2Nquad) # quadrature limits are [-1,1]
    μ₀ = cosd.(sza)
    #qp_μ = unique([qp_μ; cosd.(vza)])
    # Assign zero-weights to remaining camera zenith angles
    qp_μ = unique(FT[qp_μ[Nquad + 1:end]; cosd.(vza); μ₀]);
    n_eff = length(qp_μ) - length(wt_μ[Nquad + 1:end]);
    wt_μ = FT[wt_μ[Nquad + 1:end]; zeros(FT,n_eff)];
    Nquad = length(qp_μ);

    iμ₀ = nearest_point(qp_μ, μ₀);
    qp_μN = arr_type(reduce(vcat, (fill.(qp_μ, [pol_type.n]))))
    wt_μN = arr_type(reduce(vcat, (fill.(wt_μ, [pol_type.n]))))
    i_start = pol_type.n*(iμ₀-1) + 1

    return QuadPoints(μ₀, iμ₀, i_start, arr_type(qp_μ), arr_type(wt_μ), qp_μN, wt_μN, Nquad)
end

"""
$(FUNCTIONNAME)(::RadauQuad,
                        Ltrunc::Int,
                        obs_geom::ObsGeometry,
                        pol_type,
                        arr_type)

Computes hemispheric quadrature points with Gauss-Radau method in two intervals between [0,SZA] and [SZA,1] to include SZA as full quadrature point for DNI
SZA included as full weighted node, VZAs included with 0 weight
Returns computed quadrature points as [`QuadPoints`](@ref)
"""
function rt_set_streams(::RadauQuad,
                        Ltrunc::Int,
                        obs_geom::ObsGeometry{FT},
                        pol_type,
                        arr_type) where {FT}

    @unpack obs_alt, sza, vza, vaz = obs_geom

    # Ltrunc + 1 = number of spherical coefficients considered (+1 for l=0)
    # quadtype = 'r' for Radau (with DNI), 'g' for Gauss (no DNI)
    # sza = single solar zenith angle
    # lza = array of or single line-of-sight zenith angle
    Nquad = (Ltrunc + 1) ÷ 2

    tqp_μ₀, twt_μ₀ = gaussradau(Nquad)
    qp_μ₀ = -reverse(tqp_μ₀)
    wt_μ₀ =  reverse(twt_μ₀)
    μ₀ = cosd(sza) # check for degree to radian conversion

    if μ₀ ∈ qp_μ₀

        qp_μ = zeros(FT, Nquad)
        wt_μ = zeros(FT, Nquad)

        for i = 1:Nquad
            qp_μ[i] = (1 + qp_μ₀[i]) / 2
            wt_μ[i] = wt_μ₀[i]
        end

        Ncam = length(vza)
        μ = cosd.(vza) # check for degree to radian conversion

        # Screen out duplicate camera zenith angles
        qp_μ = unique([qp_μ; cosd.(vza)])

        # Assign zero-weights to remaining camera zenith angles
        wt_μ = FT[wt_μ; zeros(length(qp_μ) - length(wt_μ))]

        Nquad = length(qp_μ)

    else

        qp_μ = zeros(FT, 2Nquad)
        wt_μ = zeros(FT, 2Nquad)

        for i = 1:Nquad
            qp_μ[i] = (μ₀ + μ₀ * qp_μ₀[i]) / 2
            wt_μ[i] = μ₀ * wt_μ₀[i] / 2
            qp_μ[Nquad + i] = ((1 + μ₀) + (1 - μ₀) * qp_μ₀[i]) / 2
            wt_μ[Nquad + i] = (1 - μ₀) * wt_μ₀[i] / 2
        end

        Ncam = length(vza)
        μ = cosd.(vza) # check for degree to radian conversion

        # Screen out duplicate camera zenith angles
        qp_μ = unique([qp_μ; cosd.(vza)])

        # Assign zero-weights to remaining camera zenith angles
        wt_μ = FT[wt_μ; zeros(length(qp_μ) - length(wt_μ))]
        Nquad = length(qp_μ)

    end

    iμ₀ = nearest_point(qp_μ, μ₀);
    qp_μN = arr_type(reduce(vcat, (fill.(qp_μ, [pol_type.n]))))
    wt_μN = arr_type(reduce(vcat, (fill.(wt_μ, [pol_type.n]))))
    i_start = pol_type.n*(iμ₀-1) + 1
    return QuadPoints(μ₀, iμ₀, i_start, arr_type(qp_μ), arr_type(wt_μ), qp_μN, wt_μN, Nquad)
end


text:
using ApproxFun, Test, BlockArrays, BlockBandedMatrices, SpecialFunctions
    import ApproxFun: testspace, testtransforms, testmultiplication,
                      testbandedoperator, testcalculus, Block
    import SpecialFunctions: factorial

@testset "Fourier" begin
@testset "Fourier transforms" begin
    for d in (PeriodicSegment(0.1,0.5),Circle(1.0+im,2.0))
        testtransforms(CosSpace(d);minpoints=2)
        testtransforms(SinSpace(d))

        testtransforms(Taylor(d))
        testtransforms(Hardy{false}(d))

        testspace(Laurent(d);hasintegral=false)
        testspace(Fourier(d);hasintegral=false)
    end

    P = ApproxFun.plan_transform(Fourier(),4)
    v = randn(4) .+ im.*randn(4)
    @test P*v == P*real(v) + im*(P*imag(v))
end

d = PeriodicSegment()

@test sum(Fun(1,CosSpace())) ≈ 2π
@test sum(Fun(SinSpace(),[1])) == 0

f=Fun(t->cos(t)+cos(3t),CosSpace)
@test f(0.1) ≈ cos(0.1)+cos(3*0.1)
@test (f*f-Fun(t->(cos(t)+cos(3t))^2,CosSpace)).coefficients|>norm <100eps()

f=Fun(exp,Taylor(Circle()))
@test f(exp(0.1im)) ≈ exp(exp(0.1im))
@test f(1.0) ≈ exp(1.0)
g=Fun(z->1/(z-0.1),Hardy{false}(Circle()))
@test (f(1.)+g(1.)) ≈ (exp(1.) + 1/(1-.1))

## Periodic
f=Fun(x->exp(-10sin((x-.1)/2)^2),Laurent)
@test f(.5) ≈ (Conversion(space(f),Fourier(domain(f)))*f)(.5)
@test f(.5) ≈ Fun(f,Fourier)(.5)

Γ=Circle(1.1,2.2)
z=Fun(Fourier(Γ))
@test space(z)==Fourier(Γ)
@test z(1.1+2.2exp(0.1im)) ≈ 1.1+2.2exp(0.1im)

@test Fun(Laurent(0..2π),[1,1.,1.])(0.1) ≈ 1+2cos(0.1)
@test Fun(Laurent(-1..1),[1,1.,1.])(0.1) ≈ 1+2cos(π*(0.1+1))
@test Fun(Laurent(0..1),[1,1.,1.])(0.1) ≈ 1+2cos(2π*0.1)

@test abs(Fun(cos,Circle())(exp(0.1im))-cos(exp(0.1im)))<100eps()
@test abs(Fun(cos,Circle())'(exp(0.1im))+sin(exp(0.1im)))<100eps()
@test abs(Fun(cos,Circle())'(exp(0.1im))+Fun(sin,Circle())(exp(0.1im)))<100eps()

@test norm(Fun(x->Fun(cos,Fourier(-π .. π),20)(x),20)-Fun(cos,20)) <100eps()
@test norm(Fun(x->Fun(cos,Fourier(-π .. π))(x))-Fun(cos)) <100eps()
@test norm(Fun(cos,Fourier)'+Fun(sin,Fourier)) < 100eps()
@test norm(Fun(x->Fun(cos,Laurent)(x))-Fun(cos)) <100eps()
@test norm(Fun(cos,Laurent)'+Fun(sin,Laurent)) < 100eps()
@test norm(Fun(cos,Circle())'+Fun(sin,Circle()))<100eps()

for f in (Fun(θ->sin(sin(θ)),SinSpace()),Fun(θ->cos(θ)+cos(3θ),CosSpace()),
            Fun(θ->sin(sin(θ)),Fourier()),Fun(θ->cos(θ)+cos(3θ),CosSpace()))
    @test norm(integrate(f)'-f)<10eps()
end

let f=Fun(exp,Circle())
    @test component(f,1)(exp(0.1im)) ≈ exp(exp(0.1im))
    @test f(exp(0.1im)) ≈ exp(exp(0.1im))
    @test norm(f'-f)<100eps()
    @test norm(integrate(f)+1-f)<100eps()
end

let f=Fun(x->exp(-10sin((x-.1)/2)^2),Fourier)
    @test real(f)(.1) ≈ f(.1)
end

@test (Fun(z->sin(z)*cos(1/z),Circle())*Fun(z->exp(z)*airyai(1/z),Circle()))(exp(.1im)) ≈
                (z->sin(z)*cos(1/z)*exp(z)*airyai(1/z))(exp(.1im))

## Calculus

let f=Fun(t->cos(t),CosSpace)
    D=Derivative(space(f))
    @test (D*f)(.1) ≈ -sin(.1)
    @test f'(.1) ≈ -sin(.1)
end

let f=Fun(t->sin(t),SinSpace)
    D=Derivative(space(f))
    @test (D*f)(.1) ≈ cos(.1)
    @test f'(.1) ≈ cos(.1)
end

let f=Fun(cos,Fourier)
    @test norm((Derivative(space(f))^2)*f+f)<100eps()
end

## Taylor

@test Fun(Taylor())  == Fun(Taylor(),[0.,1.])

@test Fun(Taylor())(1.0) ≈ 1.0
@test Fun(Taylor(Circle(0.1,2.2)))(1.0) ≈ 1.0
@test Fun(Taylor(Circle(0.1+0.1im,2.2)))(1.0) ≈ 1.0

@test Multiplication(Fun(Taylor()),Taylor())[1:3,1:3] == [0. 0. 0.; 1. 0. 0.; 0. 1. 0.]

for d in (Circle(),Circle(0.5),Circle(-0.1,2.))
    S=Taylor(d)
    D=Derivative(S)
    ef=Fun(exp,S)
    @test norm((D*ef-ef).coefficients)<4000eps()
    @test norm((D^2*ef-ef).coefficients)<200000eps()
    u=[Evaluation(S,0.),D-I]\[1.;0.]
    @test norm((u-ef).coefficients)<200eps()
    @test norm((Integral(S)*Fun(exp,S)+ef.coefficients[1]-ef).coefficients)<100eps()

    f=Fun(z->exp(1/z)-1,Hardy{false}(d))
    df=Fun(z->-1/z^2*exp(1/z),Hardy{false}(d))
    @test norm((Derivative()*f-df).coefficients)<1000eps()
    @test norm((Derivative()^2*f-df').coefficients)<100000eps()
    @test norm((f'-df).coefficients)<1000eps()
end

d=Circle()
S=Taylor(d)
D=Derivative(S)
D-I
ef=Fun(exp,S)
@test norm((D*ef-ef).coefficients)<1000eps()
@test norm((D^2*ef-ef).coefficients)<100000eps()
u=[Evaluation(S,0.),D-I]\[1.;0.]

# check's Derivative constructor works
D=Derivative(Taylor(PeriodicSegment()))

## Multiplication

s=Fun(t->(sin(t)+sin(2t))*cos(sin(t)),SinSpace)
b=Fun(t->(sin(t)+sin(3t)),SinSpace)

@test (s*s)(.1) ≈ s(.1)^2
@test (s*b)(.1) ≈ s(.1)*b(.1)

s=Fun(t->(cos(t)+cos(2t))*cos(cos(t)),CosSpace)
b=Fun(t->(1+cos(t)+cos(3t)),CosSpace)

@test (s*s)(.1) ≈ s(.1)^2
@test (s*b)(.1) ≈ s(.1)*b(.1)

s=Fun(t->(cos(t)+cos(2t))*cos(cos(t)),CosSpace)
b=Fun(t->(sin(t)+sin(3t)),SinSpace)

@test (s*b)(.1) ≈ s(.1)*b(.1)

s=Fun(t->(sin(t)+sin(2t))*cos(sin(t)),SinSpace)
b=Fun(t->(1+cos(t)+cos(3t)),CosSpace)

@test (s*b)(.1) ≈ s(.1)*b(.1)

a=Fun(t->exp(cos(t)+sin(t)),Fourier)
b=Fun(t->sin(t)+cos(3t)+1,Fourier)

@test (a*b)(.1) ≈ a(.1)*b(.1)

a=Fun(t->exp(cos(t)),CosSpace)
b=Fun(t->sin(t)+cos(3t)+1,Fourier)

@test (a*b)(.1) ≈ a(.1)*b(.1)

a=Fun(t->sin(sin(t)),SinSpace)
b=Fun(t->sin(t)+cos(3t)+1,Fourier)

@test (a*b)(.1) ≈ a(.1)*b(.1)

# Check bug in off centre Circle
c2=-0.1+.2im;r2=0.3;
d2=Circle(c2,r2)
z=Fun(identity,d2)

@test z(-0.1+.2im+0.3*exp(0.1im)) ≈ (-0.1+.2im+0.3*exp(0.1im))

# false Circle
@test Fun(exp,Fourier(Circle(0.,1.,false)))(exp(0.1im)) ≈ exp(exp(.1im))
@test Fun(exp,Laurent(Circle(0.,1.,false)))(exp(0.1im)) ≈ exp(exp(.1im))

## Reverse orientation

f=Fun(z->1/z,Taylor(1/Circle()))
@test f(exp(0.1im)) ≈ exp(-0.1im)

## exp(z)

z=Fun(identity,Circle())
cfs=exp(z).coefficients[1:2:end]
for k=1:length(cfs)
    @test abs(cfs[k]-1/factorial(1.0(k-1))) ≤ 1E-10
end

##  Norms

@test sum(Fun(CosSpace(),[1.]))/(2π) ≈ 1.
@test sum(Fun(CosSpace(),[0.,1.])^2)/(2π) ≈ 0.5
@test sum(Fun(CosSpace(),[0.,0.,1.])^2)/(2π) ≈ 0.5
@test sum(Fun(CosSpace(),[0.,0.,0.,1.])^2)/(2π) ≈ 0.5

@test sum(Fun(SinSpace(),[0.,1.])^2)/(2π) ≈ 0.5
@test sum(Fun(SinSpace(),[0.,0.,1.])^2)/(2π) ≈ 0.5
@test sum(Fun(SinSpace(),[0.,0.,0.,1.])^2)/(2π) ≈ 0.5

## Bug in multiplicaiton

@test Fun(SinSpace(),Float64[])^2 == Fun(SinSpace(),Float64[])
@test Fun(Fourier(),[1.])^2 ≈ Fun(Fourier(),[1.])

B=Evaluation(Laurent(0..2π),0,1)
@test B*Fun(sin,domainspace(B)) ≈ 1.0

## Conversion between reverse

C = Conversion(SinSpace()⊕CosSpace(),Fourier())
@test C[Block(1), Block(1)] ≈ [0 1; 1 0]
@test ApproxFun.defaultgetindex(C, Block.(1:2), Block.(1:2)) isa AbstractMatrix

testbandedoperator(C)

## Diagonal Derivative
D = Derivative(Laurent())
@test isdiag(D)

## Test bug in multiplication

y = Fun(Circle())
@test (y^2) ≈ Fun(z->z^2,domain(y))

## Test Multiplication
mySin = Fun(Fourier(),[0,1.0])
A = Multiplication(mySin,Fourier())
@test A.op[1,1] == 0

mySin = Fun(Fourier(),[0,1])
A = Multiplication(mySin,Fourier())
@test A.op[1,1] == 0

@testset "Negatively oriented circles" begin
    f1 = Fun(z -> exp(1/z), Circle(0.0,0.2))
    f̃1 = Fun(z -> exp(1/z), Circle(0.0,0.2,false))
    f̃2 = Fun(z -> exp(1/z), Circle(0.0,0.3,false))

    @test f1(0.2exp(0.1im)) ≈ exp(1/(0.2exp(0.1im)))
    @test f̃1(0.2exp(0.1im)) ≈ exp(1/(0.2exp(0.1im)))
    @test f̃2(0.3exp(0.1im)) ≈ exp(1/(0.3exp(0.1im)))

    @test sum(f1) ≈ -sum(f̃1)
    @test sum(f̃1) ≈ sum(f̃2)
end

@testset "Fourier inplace" begin
    S = Fourier()

    x = [1.,2,3,4,5]
    y = similar(x)
    z = similar(x)
    P = ApproxFun.plan_transform(S, x)
    P! = ApproxFun.plan_transform!(S, x)
    mul!(y, P, x)
    @test x ≈ [1.,2,3,4,5]
    mul!(z, P!, x)
    @test x ≈ [1.,2,3,4,5]
    @test y ≈ z ≈ P*x ≈ P!*copy(x)

    P = ApproxFun.plan_itransform(S, x)
    P! = ApproxFun.plan_itransform!(S, x)
    mul!(y, P, x)
    @test x ≈ [1.,2,3,4,5]
    mul!(z, P!, x)
    @test x ≈ [1.,2,3,4,5]
    @test y ≈ z ≈ P*x ≈ P!*copy(x)
end
end

ext: jl


text:
using CSV
using DataFrames

DD226_DIR="/home/graham_s/OU/DD226/docs/"
MD_DIR = "$DD226_DIR/sections/"
BIB_FILE = "$DD226_DIR/data/DD226.bib"
META_FILE = "$DD226_DIR/data/DD226-metadata.yaml"
OUT_DIR = "/var/www/ou/stb/"
INCLUDE_DIR = "/home/graham_s/VirtualWorlds/projects/ou/stb.jl/web/includes/"
PANDOC_DIR = "/home/graham_s/pandoc_data/"
NullableString = Union{Missing,AbstractString}

const DEFAULT_OPTS = Dict(
    "data-dir"          => PANDOC_DIR,
    "css"               =>
            [
                "/css/ou-clone.css",
                "/css/oustb.css"
            ],
    "include-in-header" => [ "$INCLUDE_DIR/ou-js-headers.html" ],
    "from"              => "markdown+yaml_metadata_block",
    "section-divs"      => true,
    "standalone"        => true,
    # "number-sections"   => true,
    "bibliography"      => BIB_FILE,
    "metadata-file"     => META_FILE,
    "default-image-extension"=>"svg",
    "csl"               => "$PANDOC_DIR/chicago-note-bibliography.csl",
    "to"                => "html5",
    # FIXME this is broken "filter"            => "pandoc-citeproc",
    "template"          => "$INCLUDE_DIR/ou-template.html"

);

function makeoptarray( opts :: Dict, format=missing ) :: AbstractArray
    out = []
    oldformat = opts["to"]
    if ! ismissing(format)
         opts["to"] = format
    end
    s = ""
    for (k,v) in opts
        if v === true
            v = ""
        end
        if length(k) == 1
            k = "-$k"
        else
            k = "--$k"
        end
        if ! (typeof( v ) <: Array)
            v = [v]
        end
        for ve in v
            push!( out, k )
            if ve !== ""
                push!( out, ve )
            end
        end
    end
    # opts["to"] = oldformat
    out
end

function addone(
    pos,
    count,
    title::String,
    output :: NullableString,
    form::NullableString,
    content::AbstractString,
    model::NullableString,
    prev_content::NullableString,
    next_content::NullableString )

    opts = deepcopy( DEFAULT_OPTS )
    includes = []
    if ! ismissing(form)
        push!( includes, "$(INCLUDE_DIR)/$(form).html" )
    end
    if ! ismissing(output)
        push!( includes, "$(INCLUDE_DIR)/$(output).html" )
    end
    if ! ismissing(model)
        opts["include-after-body"] = [ "$(INCLUDE_DIR)/run-$model-js.html" ]
    end
    # if pos > 1
    #     opts["number-offset"]=(pos-1)
    # end
    opts["include-before-body"] = includes
    opts["o"] = "$(OUT_DIR)$(content).html"
    # opts["metadata"] = "title:$title"
    links = []

    if ! ismissing( prev_content )
        push!( links, "prev:$(prev_content).html")
    end
    if ! ismissing( next_content )
        push!( links, "next:$(next_content).html")
    end
    opts["variable"] = links
    optsarr = makeoptarray( opts )

    push!(optsarr, "$(MD_DIR)$(content).md" )
    push!(optsarr, "$(MD_DIR)footnotes.md" )

    cmd=`/usr/bin/pandoc $optsarr`
    println( cmd )
    # println( cmd )
    run( `$cmd` )
end

df = CSV.File( "$INCLUDE_DIR/ou-files.csv") |> DataFrame

# addone( 1, 2, "Introduction",missing,missing,"intro",missing,missing)

npages = size( df )[1]
outfile = "/home/graham_s/tmp/alltext.md"
outf = open( outfile,"w")
wholething = ""
for i in 1:npages
    prev = missing
    if i > 1
        prev = df[i-1,:content]
    end
    next = missing
    if i < npages
        next = df[i+1,:content]
    end
    page = df[i,:]
    addone( i, npages, page.title, page.output, page.form, page.content, page.model, prev, next )

    filecont = read("$(MD_DIR)$(page.content).md")
    write( outf, "##$(page.title)\n\n")
    write( outf, filecont )
end

close( outf )


text:
function PotentialSimulationSetupRB(ssd::SolidStateDetector{T, :cartesian}, grid::Grid{T, 3, :cartesian} = Grid(ssd),
                potential_array::Union{Missing, Array{T, 3}} = missing; weighting_potential_contact_id::Union{Missing, Int} = missing) where {T}#::PotentialSimulationSetupRB{T} where {T}
    is_weighting_potential::Bool = !ismissing(weighting_potential_contact_id)

    @inbounds begin
        begin # Geometrical weights of the Axes
            nx, ny, nz = size(grid)

            # X-axis
            ax_x::Vector{T} = collect(grid.axes[1]) # real grid points/ticks -> the potential at these ticks are going to be calculated
            Δx::Vector{T} = diff(ax_x) # difference between real ticks
            x_ext::Vector{T} = get_extended_ticks(grid.axes[1])
            Δx_ext::Vector{T} = diff(x_ext)
            Δx_ext_inv::Vector{T} = inv.(Δx_ext)
            mpx::Vector{T} = midpoints(x_ext)
            mpx_inv::Vector{T} = inv.(mpx)
            Δmpx::Vector{T} = diff(mpx)
            Δmpx_inv::Vector{T} = inv.(Δmpx)
            Δhmpxr::Vector{T} = mpx[2:end] - ax_x # distances between midpoints and real grid points (half distances -> h), needed for weights wrr, wrl, ... & dV (volume element)
            Δhmpxl::Vector{T} = ax_x - mpx[1:end - 1]
            wxr::Vector{T} = Δmpx_inv .* Δhmpxr # weights for epislon_r adding
            wxl::Vector{T} = Δmpx_inv .* Δhmpxl
            wx::Array{T, 2} = zeros(T, 4, length(wxr) + 1)
            wx[1, 1:length(wxr)] = wxr
            wx[2, 1:length(wxr)] = wxl
            wx[3, 1:length(wxr)] = Δmpx
            wx[4, :] = Δx_ext_inv
            gw_x::GeometricalCartesianAxisWeights{T} = GeometricalCartesianAxisWeights{T}( wx ) # Weights needed for Field Simulation loop
            x_ext_mid::T = (x_ext[end - 1] - x_ext[2]) / 2
            grid_boundary_factor_x_right::T = abs((x_ext[end - 1] - x_ext_mid) / (x_ext[end] - x_ext_mid))
            grid_boundary_factor_x_left::T = abs((x_ext[2] - x_ext_mid) / (x_ext[1] - x_ext_mid))

            # Y-axis
            ax_y::Vector{T} = collect(grid.axes[2]) # real grid points/ticks -> the potential at these ticks are going to be calculated
            Δy::Vector{T} = diff(ax_y) # difference between real ticks
            y_ext::Vector{T} = get_extended_ticks(grid.axes[2])
            Δy_ext::Vector{T} = diff(y_ext)
            Δy_ext_inv::Vector{T} = inv.(Δy_ext)
            mpy::Vector{T} = midpoints(y_ext)
            mpy_inv::Vector{T} = inv.(mpy)
            Δmpy::Vector{T} = diff(mpy)
            Δmpy_inv::Vector{T} = inv.(Δmpy)
            Δhmpyr::Vector{T} = mpy[2:end] - ax_y # distances between midpoints and real grid points (half distances -> h), needed for weights wrr, wrl, ... & dV (volume element)
            Δhmpyl::Vector{T} = ax_y - mpy[1:end - 1]
            wyr::Vector{T} = Δmpy_inv .* Δhmpyr # weights for epislon_r adding
            wyl::Vector{T} = Δmpy_inv .* Δhmpyl
            wy::Array{T, 2} = zeros(T, 4, length(wyr) + 1)
            wy[1, 1:length(wyr)] = wyr
            wy[2, 1:length(wyr)] = wyl
            wy[3, 1:length(wyr)] = Δmpy
            wy[4, :] = Δy_ext_inv
            gw_y::GeometricalCartesianAxisWeights{T} = GeometricalCartesianAxisWeights{T}( wy ) # Weights needed for Field Simulation loop
            y_ext_mid::T = (y_ext[end - 1] - y_ext[2]) / 2
            grid_boundary_factor_y_right::T = abs((y_ext[end - 1] - y_ext_mid) / (y_ext[end] - y_ext_mid))
            grid_boundary_factor_y_left::T = abs((y_ext[2] - y_ext_mid) / (y_ext[1] - y_ext_mid))

            # Z-axis
            ax_z::Vector{T} = collect(grid.axes[3]) # real grid points/ticks -> the potential at these ticks are going to be calculated
            Δz::Vector{T} = diff(ax_z) # difference between real ticks
            z_ext::Vector{T} = get_extended_ticks(grid.axes[3])
            Δz_ext::Vector{T} = diff(z_ext)
            Δz_ext_inv::Vector{T} = inv.(Δz_ext)
            mpz::Vector{T} = midpoints(z_ext)
            mpz_inv::Vector{T} = inv.(mpz)
            Δmpz::Vector{T} = diff(mpz)
            Δmpz_inv::Vector{T} = inv.(Δmpz)
            Δhmpzr::Vector{T} = mpz[2:end] - ax_z # distances between midpoints and real grid points (half distances -> h), needed for weights wrr, wrl, ... & dV (volume element)
            Δhmpzl::Vector{T} = ax_z - mpz[1:end - 1]
            wzr::Vector{T} = Δmpz_inv .* Δhmpzr # weights for epislon_r adding
            wzl::Vector{T} = Δmpz_inv .* Δhmpzl
            wz::Array{T, 2} = zeros(T, 4, length(wzr) + 1)
            wz[1, 1:length(wzr)] = wzr
            wz[2, 1:length(wzr)] = wzl
            wz[3, 1:length(wzr)] = Δmpz
            wz[4, :] = Δz_ext_inv
            gw_z::GeometricalCartesianAxisWeights{T} = GeometricalCartesianAxisWeights{T}( wz ) # Weights needed for Field Simulation loop
            z_ext_mid::T = (z_ext[end - 1] - z_ext[2]) / 2
            grid_boundary_factor_z_right::T = abs((z_ext[end - 1] - z_ext_mid) / (z_ext[end] - z_ext_mid))
            grid_boundary_factor_z_left::T = abs((z_ext[2] - z_ext_mid) / (z_ext[1] - z_ext_mid))

            geom_weights::NTuple{3, AbstractGeometricalAxisWeights{T}} = (gw_x, gw_y, gw_z) # Weights needed for Field Simulation loop
            grid_boundary_factors::NTuple{3, NTuple{2, T}} = (  (grid_boundary_factor_x_left, grid_boundary_factor_x_right),
                                                                (grid_boundary_factor_y_left, grid_boundary_factor_y_right),
                                                                (grid_boundary_factor_z_left, grid_boundary_factor_z_right))
        end

        # detector_material_ϵ_r::T = ssd.material_detector.ϵ_r
        # environment_material_ϵ_r::T = ssd.material_environment.ϵ_r

        bulk_is_ptype::Bool = ssd.semiconductors[1].bulk_type == :ptype ? true : false
        contact_bias_voltages::Vector{T} = T[contact.potential for contact in ssd.contacts]
        minimum_applied_potential::T = minimum(contact_bias_voltages)
        maximum_applied_potential::T = maximum(contact_bias_voltages)
        bias_voltage::T = maximum_applied_potential - minimum_applied_potential
        depletion_handling_potential_limit::T = -bias_voltage
        sor_consts::Vector{T} = [1]

        ϵ = Array{T, 3}(undef, length(mpx), length(mpy), length(mpz))
        ρ_tmp = Array{T, 3}(undef, length(mpx), length(mpy), length(mpz))
        for iz in 1:size(ϵ, 3)
            pos_z::T = mpz[iz]
            for iy in 1:size(ϵ, 2)
                pos_y::T = mpy[iy]
                for ix in 1:size(ϵ, 1)
                    pos_x = mpx[ix]
                    pt::CartesianPoint{T} = CartesianPoint{T}(pos_x, pos_y, pos_z)
                    ρ_tmp[ix, iy, iz]::T, ϵ[ix, iy, iz]::T = get_ρ_and_ϵ(pt, ssd)
                end
            end
        end
        ϵ0_inv::T = inv(ϵ0)
        ρ_tmp *= ϵ0_inv

        volume_weights::Array{T, 4} = RBExtBy2Array(T, grid)
        ρ::Array{T, 4} = RBExtBy2Array(T, grid)
        for iz in range(2, stop = length(z_ext) - 1)
            inz::Int = iz - 1
            pos_z::T = z_ext[iz]
            for iy in range(2, stop = length(y_ext) - 1)
                iny::Int = iy - 1
                pos_y::T = y_ext[iy]
                for ix in range(2, stop = length(x_ext) - 1)
                    inx::Int = ix - 1;
                    irbx::Int = rbidx(inx)
                    pos_x::T = x_ext[ix]

                    rbi::Int = iseven(inx + iny + inz) ? rb_even::Int : rb_odd::Int

                    ρ_cell::T = 0
                    if !is_weighting_potential
                        ρ_cell += ρ_tmp[ ix,  iy,  iz] * wxr[inx] * wyr[iny] * wzr[inz]
                        ρ_cell += ρ_tmp[ ix,  iy, inz] * wxr[inx] * wyr[iny] * wzl[inz]
                        ρ_cell += ρ_tmp[ ix, iny,  iz] * wxr[inx] * wyl[iny] * wzr[inz]
                        ρ_cell += ρ_tmp[ ix, iny, inz] * wxr[inx] * wyl[iny] * wzl[inz]

                        ρ_cell += ρ_tmp[inx,  iy,  iz] * wxl[inx] * wyr[iny] * wzr[inz]
                        ρ_cell += ρ_tmp[inx,  iy, inz] * wxl[inx] * wyr[iny] * wzl[inz]
                        ρ_cell += ρ_tmp[inx, iny,  iz] * wxl[inx] * wyl[iny] * wzr[inz]
                        ρ_cell += ρ_tmp[inx, iny, inz] * wxl[inx] * wyl[iny] * wzl[inz]
                    end

                    # right weight in x: wxr
                    wxr_eps::T = ϵ[  ix,  iy, inz + 1] * wyr[iny] * wzr[inz]
                    wxr_eps   += ϵ[  ix, iny, inz + 1] * wyl[iny] * wzr[inz]
                    wxr_eps   += ϵ[  ix,  iy, inz ]    * wyr[iny] * wzl[inz]
                    wxr_eps   += ϵ[  ix, iny, inz ]    * wyl[iny] * wzl[inz]
                    # left weight in x: wxr
                    wxl_eps::T = ϵ[ inx,  iy, inz + 1] * wyr[iny] * wzr[inz]
                    wxl_eps   += ϵ[ inx, iny, inz + 1] * wyl[iny] * wzr[inz]
                    wxl_eps   += ϵ[ inx,  iy, inz ]    * wyr[iny] * wzl[inz]
                    wxl_eps   += ϵ[ inx, iny, inz ]    * wyl[iny] * wzl[inz]
                    # right weight in y: wyr
                    wyr_eps::T = ϵ[ inx,  iy, inz + 1] * wxl[inx] * wzr[inz]
                    wyr_eps   += ϵ[  ix,  iy, inz + 1] * wxr[inx] * wzr[inz]
                    wyr_eps   += ϵ[ inx,  iy, inz ]    * wxl[inx] * wzl[inz]
                    wyr_eps   += ϵ[  ix,  iy, inz ]    * wxr[inx] * wzl[inz]
                    # left weight in y: wyl
                    wyl_eps::T = ϵ[ inx, iny, inz + 1] * wxl[inx] * wzr[inz]
                    wyl_eps   += ϵ[  ix, iny, inz + 1] * wxr[inx] * wzr[inz]
                    wyl_eps   += ϵ[ inx, iny, inz ]    * wxl[inx] * wzl[inz]
                    wyl_eps   += ϵ[  ix, iny, inz ]    * wxr[inx] * wzl[inz]
                    # right weight in z: wzr
                    wzr_eps::T = ϵ[  ix,  iy, inz + 1] * wxr[inx] * wyr[iny]
                    wzr_eps   += ϵ[  ix, iny, inz + 1] * wxr[inx] * wyl[iny]
                    wzr_eps   += ϵ[ inx,  iy, inz + 1] * wxl[inx] * wyr[iny]
                    wzr_eps   += ϵ[ inx, iny, inz + 1] * wxl[inx] * wyl[iny]
                    # left weight in z: wzr
                    wzl_eps::T = ϵ[  ix,  iy,    inz ] * wxr[inx] * wyr[iny]
                    wzl_eps   += ϵ[  ix, iny,    inz ] * wxr[inx] * wyl[iny]
                    wzl_eps   += ϵ[ inx,  iy,    inz ] * wxl[inx] * wyr[iny]
                    wzl_eps   += ϵ[ inx, iny,    inz ] * wxl[inx] * wyl[iny]

                    volume_weight::T = wxr_eps * Δx_ext_inv[ ix] * Δmpy[iny] * Δmpz[inz]
                    volume_weight   += wxl_eps * Δx_ext_inv[inx] * Δmpy[iny] * Δmpz[inz]
                    volume_weight   += wyr_eps * Δy_ext_inv[ iy] * Δmpx[inx] * Δmpz[inz]
                    volume_weight   += wyl_eps * Δy_ext_inv[iny] * Δmpx[inx] * Δmpz[inz]
                    volume_weight   += wzr_eps * Δz_ext_inv[ iz] * Δmpx[inx] * Δmpy[iny]
                    volume_weight   += wzl_eps * Δz_ext_inv[inz] * Δmpx[inx] * Δmpy[iny]

                    volume_weights[ irbx, iy, iz, rbi ] = inv(volume_weight)

                    dV::T = Δmpx[inx] * Δmpy[iny] * Δmpz[inz]
                    ρ[ irbx, iy, iz, rbi ] = dV * ρ_cell
                end
            end
        end

        potential::Array{T, 3} = ismissing(potential_array) ? zeros(T, size(grid)...) : potential_array
        pointtypes::Array{PointType, 3} = ones(PointType, size(grid)...)
        set_pointtypes_and_fixed_potentials!( pointtypes, potential, grid, ssd, weighting_potential_contact_id = weighting_potential_contact_id  )
        rbpotential::Array{T, 4}  = RBExtBy2Array( potential, grid )
        rbpointtypes::Array{T, 4} = RBExtBy2Array( pointtypes, grid )
        potential = clear(potential)
        pointtypes = clear(pointtypes)
    end # @inbounds

    fssrb::PotentialSimulationSetupRB{T, 3, 4, :cartesian} = PotentialSimulationSetupRB{T, 3, 4, :cartesian}(
        grid,
        rbpotential,
        rbpointtypes,
        volume_weights,
        ρ,
        ϵ,
        geom_weights,
        sor_consts,
        bias_voltage,
        maximum_applied_potential,
        minimum_applied_potential,
        depletion_handling_potential_limit,
        bulk_is_ptype,
        grid_boundary_factors
     )
    return fssrb
end

function ElectricPotentialArray(fssrb::PotentialSimulationSetupRB{T, 3, 4, :cartesian})::Array{T, 3} where {T}
    pot::Array{T, 3} = Array{T, 3}(undef, size(fssrb.grid))
    for iz in axes(pot, 3)
        irbz::Int = iz + 1
        for iy in axes(pot, 2)
            irby::Int = iy + 1
            idxsum::Int = iz + iy
            for ix in axes(pot, 1)
                irbx::Int = rbidx(ix)
                rbi::Int = iseven(idxsum + ix) ? rb_even::Int : rb_odd::Int
                pot[ix, iy, iz] = fssrb.potential[ irbx, irby, irbz, rbi ]
            end
        end
    end
    return pot
end

function PointTypeArray(fssrb::PotentialSimulationSetupRB{T, 3, 4, :cartesian})::Array{PointType, 3} where {T}
    pointtypes::Array{PointType, 3} = zeros(PointType, size(fssrb.grid))
    for iz in axes(pointtypes, 3)
        irbz::Int = iz + 1
        for iy in axes(pointtypes, 2)
            irby::Int = iy + 1
            idxsum::Int = iz + iy
            for ix in axes(pointtypes, 1)
                irbx::Int = rbidx(ix)
                rbi::Int = iseven(idxsum + ix) ? rb_even::Int : rb_odd::Int
                pointtypes[ix, iy, iz] = fssrb.pointtypes[irbx, irby, irbz, rbi ]
            end
        end
    end
    return pointtypes
end

function ChargeDensityArray(fssrb::PotentialSimulationSetupRB{T, 3, 4, :cartesian})::Array{T} where {T}
    ρ::Array{T, 3} = zeros(T, size(fssrb.grid))
    for iz in axes(ρ, 3)
        irbz::Int = iz + 1
        Δmpz::T = fssrb.geom_weights[3].weights[3, iz]
        for iy in axes(ρ, 2)
            irby::Int = iy + 1
            idxsum::Int = iz + iy
            Δmpy::T = fssrb.geom_weights[2].weights[3, iy]
            Δmpzy::T = Δmpz * Δmpy
            for ix in axes(ρ, 1)
                irbx::Int = rbidx(ix)
                rbi::Int = iseven(idxsum + ix) ? rb_even::Int : rb_odd::Int
                dV::T =  Δmpzy * fssrb.geom_weights[1].weights[3, ix]

                ρ[ix, iy, iz] = fssrb.ρ[irbx, irby, irbz, rbi ] / dV
            end
        end
    end
    return ρ
end