text:
"""Module for the modelling of noise pollution."""
module Noises

using ArchGDAL
using GeoArrays
using Plots
using Shapefile

include(".\\Utils\\Functions.jl")

export run_noise

const agd = ArchGDAL
const ga = GeoArrays
const sf = Shapefile

repeat!(A::AbstractVector, count::Int64 ) = append!( A, repeat(A, count-1) )
# Functions to find the coordinates of the point resulting from the rotation of "(xp, yp)" by a angle "θ" around "(xc, yc)"
rotate_x( xp, yp, xc, yc, θ ) = round( Int64, (xp - xc)cos(deg2rad(θ)) - (yp - yc)sin(deg2rad(θ)) + xc )
rotate_y( xp, yp, xc, yc, θ ) = round( Int64, (xp - xc)sin(deg2rad(θ)) + (yp - yc)cos(deg2rad(θ)) + yc )
rotate_point( xp, yp, xc, yc, θ ) = θ == 0 ? (xp, yp) : ( rotate_x( xp, yp, xc, yc, θ ), rotate_y( xp, yp, xc, yc, θ ) )
distance( p0::Tuple{Float64, Float64}, pn::Tuple{Float64, Float64} ) =  √( sum( v -> v^2, pn .- p0 ) )

"""
    transmission_loss( r::Float64 )

Compute the transmission loss of a noise over `r` distance.
"""
function transmission_loss( radius::Float64 )
    return 20log10(radius)
end

"""
    atmospheric_absorpion_loss( radius::Float64, height_m::Float32, relative_humidity::Float64, temperature_k::Float64, frequency::Float64 )

Compute the attenuation of a sound at a distance `radius` due to the atmosphere.
"""
function atmospheric_absorption_loss( radius::Float64, height_m::Float32, relative_humidity::Float64, temperature_k::Float64, frequency::Float64 )
    # Calculate atmospheric absorption coefficient using ANSI S1.26-1995 standard
    # Convert elevation to atmospheric pressure
    atmospheric_pressure = 101.325( 1 - ( 2.25577 * 10^(-5) * height_m ) )^5.25588
    # Calculate derived values for subsequent equations
    p_atm_pressure = atmospheric_pressure / 101.325
    t_tr = temperature_k / 293.15
    # Convert relative humidity to molar concentration of water vapor
    C = ( -6.8346( 273.16 / temperature_k )^1.261 ) + 4.6151
    p_saturation_pressure = 10^C
    humidity = relative_humidity * p_saturation_pressure * p_atm_pressure^(-1)
    # Calculate relaxation frequency of O (equation 3)
    #   frO₂ = ( p_atm_pressure * ( (24 + 4.04e04) * humidity ) * (0.02 + humidity) ) / (0.391 + humidity)
    frO₂ = p_atm_pressure * ( 24 + ( 4.04e04humidity * ( (0.02 + humidity) / (0.391 + humidity) ) ) )
    # Calculate relaxation frequency of N (equation 4)
    frN₂ = p_atm_pressure * √t_tr * ( 9 + ( 280humidity * ℯ^( -4.170(t_tr^(-1/3) - 1) ) ) )
    # Calculate alpha (equation 5)
    term1 = 1.84 * 10^(-11) * p_atm_pressure^(-1) * √t_tr
    #   term2 = t_tr^(-2.5) * ( 0.01275 * ℯ^(-2239.1 / temp_k) * ( frO₂ / (frO₂^2 + freq^2) ) )
    term2 = t_tr^(-2.5)( 0.01275ℯ^(-2239.1 / temperature_k) / ( frO₂ + (frequency^2 / frO₂) ) )
    #   term3 = 0.1068 * ℯ^(-3352 / temp_k) * ( frN₂ / (frN₂^2 + freq^2) )
    term3 = t_tr^(-2.5) * ( 0.1068ℯ^(-3352 / temperature_k) / ( frN₂ + (frequency^2 / frN₂) ) )
    #   α = 8.686 * (frequency^2) * ( term1 + term2 + term3 )
    α = frequency^2 * (term1 + term2 + term3)
    return (α * radius) / 100
end

"""
    minmax( profile::AbstractArray{T, 1}, rel_h_src::Float32=0.0f0, rel_h_rec::Float32=0.0f0 )::Vector{T} where {T <: Tuple{Float64, Float32}}

Return the minimum and maximum peaks in `profile` coupled with their respective position in the vector.
"""
function minmax( profile::AbstractArray{T, 1}, rel_h_src::Float32=0.0f0, rel_h_rec::Float32=0.0f0 ) where {T <: Tuple{Float64, Float32}}
    if length(profile) <= 1
        return Tuple{Int64, Float32}[ (1, -rel_h_src), (1, -rel_h_rec) ]
    end
    # Test for the profile having a finite length; will be zero if directly overhead
    if abs( profile[1][1] - profile[end][1] ) < 10.0
        return Tuple{Int64, Float32}[ (1, rel_h_src), (length(profile), rel_h_rec) ]
    end
    # Define parameters A and B such that height of the source-receiver line
    # is given by z = Ax + B.
    x1 = profile[1][1]
    z1 = profile[1][2] + rel_h_src
    xn = profile[end][1]
    zn = profile[end][2] + rel_h_rec
    a = (zn - z1) / (xn - x1)
    b = z1 - a*x1

 #=
    # Look for the max and min
    prf = map( row -> ( row[1], row[2] - ( a*row[1] + b ) ), profile )
    sort!( prf, by=(row) -> row[2] )

    return [ prf[1], prf[end] ]
 =#
    min = ( 1, profile[1][2] )
    max = ( 1, profile[1][2] )
    for i in 1:length(profile)
        hgt = profile[i][2] - a*profile[i][1] + b
        if hgt > max[2]
            max = (i, hgt)
        end
        if hgt < min[2]
            min = (i, hgt)
        end
    end
    return Tuple{Int64, Float32}[min, max]
end

function delbaz( freq::Float64, flow_res::Float32 )
    dumr = 1.0 + 9.08 * (flow_res / freq)^0.75
    dumi = -11.9 * (flow_res / freq)^0.73
    return complex(dumr, dumi)
end

function subw( a::AbstractFloat, b::AbstractFloat, c::AbstractFloat, d::AbstractFloat, w::Complex )
    an = (a^2 - b^2  - d)^2 + (2.0 * a * b)^2
    r_i = c * (a^2 + b^2 + d)
    wr = b * r_i / an + w.re
    wi = a * r_i / an + w.im
    return complex(wr, wi)
end

function ww(t::Complex)
    a = abs(t.re)
    b = abs(t.im)

    if a <= 3.9 && b <= 3.0
        a1 = cos(2.0 * a * b)
        b1 = sin(2.0 * a * b)
        c1 = ℯ^( -2.0 * b * π / 0.8 ) - cos( 2.0 * a * π / 0.8 )
        d1 = sin( 2.0 * a * π / 0.8 )

        pq = 2.0 * ℯ^( -( a^2 + 2.0 * b * π / 0.8 - b^2 ) ) / ( c1^2 + d1^2 )
        p2 = pq * ( a1*c1 - b1*d1 )
        q2 = pq * ( a1*d1 - b1*c1 )

        aa = a^2
        bb = b^2
        ah = 0.8 * b / ( π * (aa + bb) )
        ak = 0.8 * a / ( π * (aa + bb) )

        for i in 1:5
            an = ( bb - aa + i^2 * 0.64 )^2 + 4.0 * aa * bb
            hk = 2.0 * 0.8 / π * ℯ^( -i^2 * 0.64 ) / an
            ah += b * hk * ( bb + aa + 0.64 * i^2 )
            ak += a * hk * ( bb + aa - 0.64 * i^2 )
        end

        if b < π/0.8
            ah += p2
            ak -= q2
        end

        if b == π/0.8
            ah += p2 / 2.0
            ak -= q2 / 2.0
        end

        w = complex(ah, ak)
    else
        cd = [
                 (0.4613135, 0.1901635),
                 (0.09999216, 1.7844927),
                 (0.002883894, 5.5253437)
             ]
        w = 0.0 + 0.0im
        for (c, d) in cd
            w = subw(a, b, c, d, w)
        end
    end

    if t.re < 0.0
        w = conj(w)
        a = -a
    end

    if t.im < 0.0
        wr = 2.0 * ℯ^(b^2 - a^2) * cos(2*a*b) - w.re
        wi = 2.0 * ℯ^(b^2 - a^2) * sin(2*a*b) - w.im
        w = complex(wr, wi)
    end

    return w
end

function qq( r::AbstractFloat, h::AbstractFloat, waveno::AbstractFloat, z::Complex )
    c = abs(h) / r
    n = (z.re * c + 1.0)^2 + (z.im * c)^2

    rr = ( ( c * abs2(z) ) - 1.0 ) / n
    ri = 2.0 * z.im * c / n

    nyr = z.re / abs2(z)
    nyi = -z.im / abs2(z)

    dumr = √(waveno * r / 4.0) * (nyr + c - nyi)
    dumi = √(waveno * r / 4.0) * (nyr + c + nyi)
    t = complex(dumr, dumi)

    w = ww(t)

    fr = 1.0 + √π * -imag(t*w)
    fi = √π * real(t*w)

    dumr = rr + (1.0 - rr) * fr + fi * ri
    dumi = ri + fi * (1.0 - rr) - ri * fr
    return z.re >= 1000.0 ? 1.0+0im : complex(dumr, dumi)
end

function qq2( d::AbstractFloat, src_h::AbstractFloat, rec_h::AbstractFloat, freq::AbstractFloat, z::Complex )
    waveno = 2.0 * π * freq / 340.0
    #   r1 = √( d^2 + (src_h - rec_h)^2 )
    r = √( d^2 + (src_h + rec_h)^2 )
    c = (src_h + rec_h) / r

    n = (z.re * c + 1.0)^2 + (z.im * c)^2
    rr = ( ( c * abs(z) )^2 - 1.0 ) / n
    ri = 2.0 * z.im * c / n

    nyr = z.re / abs2(z)
    nyi = -z.im / abs2(z)

    dumr = √(waveno * r / 4.0) * (nyr + c - nyi)
    dumi = √(waveno * r / 4.0) * (-nyr - c + nyi)
    t = complex(dumr, dumi)

    w = ww(-t)

    fr = 1.0 + √π * imag(t*w)
    fi = -√π * real(t*w)

    dumr = rr + (1.0 - rr) * fr + fi * ri
    dumi = ri + fi * (1.0 - rr) - ri * fr
    return complex(dumr, dumi)
end

"""
    egal( d1::Float64, d2::Float64, src_h::Float32, rec_h::Float32, src_flow_res::Float32, rec_flow_res::Float32, e_wind_vel::AbstractFloat, transition_height::AbstractFloat, turbulence::AbstractFloat, freq::Float64, dum::AbstractFloat )

Run level topography propagation model for spectral sound levels.
"""
function egal( d1::Float64, d2::Float64, src_h::Float32, rec_h::Float32, src_flow_res::Float32, rec_flow_res::Float32, e_wind_vel::Float64, transition_height::Float64,
               turbulence::Float64, freq::Float64 )
 #=
    arg1  = d1
    arg2  = d2
    arg3  = src_h
    arg4  = rec_h
    arg5  = src_flow_res
    arg6  = rec_flow_res
    arg7  = e_wind_vel
    arg8  = transition_height
    arg9  = turbulence
    arg10 = freq
    arg11 = dum
 =#
    src_rx = delbaz( freq, src_flow_res )
    rec_rx = delbaz( freq, rec_flow_res )

    rd = √( (d1 + d2)^2 + (src_h - rec_h)^2 )
    rr = √( (d1 + d2)^2 + (src_h + rec_h)^2 )

    la = 340.0 / freq
    fixed_speed = 2.0 * π * la^(-1)
    e_turbulence_scale = 0.0
    ar = 0.0

    kv = Float64[ 0.0, 0.0, 0.0 ]
    db = Complex{Float64}[ 0.0+0.0im, 0.0+0.0im, 0.0+0.0im ]

    if e_wind_vel != 0
        m = round( transition_height / ( la / 6.0 ) )
        fixed_speed = 2.0 * π * freq / (
                                 340.0 + (
                                             ( m-1 * la/6.0 + la/10.0 ) / 2.0
                                         ) * e_wind_vel/10.0
                             )
        e_turbulence_scale = 10e3 / freq

        eq(x, rrn, rrm ) = ℯ^complex( 0.0, -x*(rrn + rrm) )
        eq_cr( rrn, rrm ) = rrm * √( rrm * rrn * (rrn + rrm) )
        eq_jarr( kv, k, rrn, rrm; q1=nothing, q2=nothing ) =  ( eq(kv, rrn, rrm) - eq(k, rrn, rrm) ) * ( isnothing(q1) ? 1 : q1 ) * ( isnothing(q2) ? 1 : q2 ) / complex( eq_cr(rrn, rrm), 0.0 )

        for j in 1:m
            ha = (j-1) * (la/6.0) + (la/10.0)
            v = Float64[ e_wind_vel + turbulence / 10.0 * cos( ha * 2π/e_turbulence_scale ) ]
            push!(
                v,
                2e_wind_vel - v[1],
                e_wind_vel
            )

            for i in 1:3
                kv[i] = 2.0 * π * freq / ( 340.0 + (ha/2.0) * v[i]/10.0 )

                rr = Float64[
                    √( d1^2 + (src_h - ha)^2 ),
                    √( d1^2 + (src_h + ha)^2 ),
                    √( d2^2 + (ha - rec_h)^2 ),
                    √( d2^2 + (ha + rec_h)^2 )
                ]

                #   cr = eq_cr( rr[1], rr[3] )
                #   jarray[i] = ( eq( kv[i], rr[1], rr[3] ) - eq( fixed_speed, rr[1], rr[3] ) ) / (cr + 0.0im)
                #   db[i] += jarray[i]

                #   q2 = qq2( d2, ha, rec_h, freq, rec_rx )
                #   cr = eq_cr( rr[1], rr[4] )
                #   jarray[i] = ( eq( kv[i], rr[1], rr[4] ) - eq( fixed_speed, rr[1], rr[4] ) ) * q2 / (cr + 0.0im)
                #   db[i] += jarray[i]
                #
                #   q1 = qq2( d1, src_h, ha, freq, src_rx )
                #   cr = eq_cr( rr[2], rr[3] )
                #   jarray[i] = ( eq( kv[i], rr[2], rr[3] ) - eq( fixed_speed, rr[2], rr[3] ) ) * q1 / (cr + 0im)
                #   db[i] += jarray[i]

                #   cr = eq_cr( rr[2], rr[4] )
                #   jarray[i] = ( eq( kv[i], rr[2], rr[4] ) - eq( fixed_speed, rr[2], rr[4] ) ) * q1 * q2 / (cr + 0im)
                #   db[i] += jarray[i]

             #= LA SCRITTURA SOTTO E' EQUIVALENTE
                jarray = Complex{Float64}[
                    eq_jarr( kv[i], fixed_speed, rr[1], rr[3] ),
                    eq_jarr( kv[i], fixed_speed, rr[1], rr[4], q2=q2 ),
                    eq_jarr( kv[i], fixed_speed, rr[2], rr[3], q1=q1 ),
                    eq_jarr( kv[i], fixed_speed, rr[2], rr[4], q1=q1, q2=q2 )
                ]
                for val in jarray
                    db[i] += val
                end
             =#
                db[i] += eq_jarr( kv[i], fixed_speed, rr[1], rr[3] ) + eq_jarr( kv[i], fixed_speed, rr[1], rr[4], q2=q2 ) +
                    eq_jarr( kv[i], fixed_speed, rr[2], rr[3], q1=q1 ) + eq_jarr( kv[i], fixed_speed, rr[2], rr[4], q1=q1, q2=q2 )
            end
        end

        c = ℯ^complex( 0.0, π/4.0 )
        crh = (la / 6.0) * d2 * √(8.0 * π * k) / (16.0 * π^2)
        c *= crh + 0.0im

        db .*= c

    end
    q2 = qq2( d1+d2, src_h, rec_h, freq, rec_rx )
    c = ℯ^complex( 0.0, -fixed_speed*rr ) / complex( rr, 0.0 ) * q2
    ch = ℯ^complex( 0.0, -fixed_speed*rd ) / complex( rd, 0.0 )
    c += ch

    c = c / complex( 4.0*π, 0.0 )

    nff = 16.0 * π * 2.0 * rd^2
    if e_wind_vel != 0
        db .+= c
        ar = sum( x -> abs(x)^2, db )
        levturb = 4.3429 * log(ar * nff / 3.0)
        lnot = 4.3429 * log( abs(db[3])^2 * nff )
    else
        levturb = 0.0
        lnot = 4.3429 * log( abs(c)^2 * nff )
    end

    return (levturb, lnot)
end

"""
    varysurf( dists::AbstractArray{T1, 1}, ground_type::AbstractArray{T2, 1}, src_h::Float32, rec_h::Float32, soft_atten::AbstractFloat, hard_atten::AbstractFloat ) where {T1 <: Float64, T2 <:Float32}

Run mixed terrain propagation model.
"""
function varysurf( dists::AbstractArray{T1, 1}, ground_type::AbstractArray{T2, 1}, src_h::Float32, rec_h::Float32, soft_atten::AbstractFloat, hard_atten::AbstractFloat ) where {T1 <: Float64, T2 <:Float32}
    drefl = src_h * dists[end] / (src_h + rec_h)

    srcInf = ground_type[1] == Soft ?
             src_h > 3.0 ?
                 (20.0 * src_h / 3.0 - 10.0) * src_h :
                 10.0 * src_h :
             30.0 * src_h

    recInf = ground_type[end] == Soft ?
             rec_h > 3.0 ?
                 (20.0 * rec_h / 3.0 - 10.0) * rec_h :
                 10.0 * rec_h :
             30.0 * rec_h

    srcInf = min( srcInf, 0.7*drefl )
    recInf = min( recInf, 0.7*(dists[end]-drefl) )
    recInf = dists[end] - recInf

    Δl = sum = denom = 0
    for i in 1:length(dists)
        w1 = w2 = 0
        if dists[i] > srcInf && dists[i] < recInf
            w1 = dists[i] < drefl ? rec_h : src_h
            w2, Δl = ground_type[i] == Hard ? (0.5, hard_atten) : (1.0, soft_atten)
        end
        sum += w1 * w2 * Δl
        denom += w1 * w2
    end

    return sum / denom
end

function fres( y::AbstractFloat )
    c = 0.797885
    x = c * y
    f = (1.0 + 0.962x) / (2.0 + 1.792x + 3.104x^2 )
    g = 1.0 / (2.0 + 4.142x + 3.492x^2 + 6.67x^3 )
    si = sin( (x / c)^2 )
    co = cos( (x / c)^2 )
    return complex( (-f * si + g * co)/c, (f * co + g * si)/c )
end

function diffraction!( aalast::Ref{Float32}, r1::AbstractFloat, a::AbstractFloat, al2::AbstractFloat, pm::AbstractFloat, any::AbstractFloat, k::AbstractFloat )
    df = -ℯ^complex(0.0, k * r1 + π / 4.0) / complex(r1, 0.0)
    tangent = tan( (π + pm * al2) / (2.0 * any) )
    aa = tangent != 0 ? ( 1.0 / tangent / (2.0 * any) / √(2.0 * π * k * a) ) : aalast[]
    aalast[] = aa
    df *= complex(aa, 0.0)
    n = pm < -0.9 ?
            al2 > π + any * π ? 1 :
                al2 > π - any * π ? 0 : -1 :
            al2 > any * π - π ? 1 : 0
    xv = 2.0 * k * a * cos( (2.0 * n * any * π - al2) / 2.0 )^2
    return df * ℯ^complex(0.0, -xv) * fres(√xv) * complex(0.0, -2.0√xv)
end

function calc_mirror( locs::Float64, points::AbstractArray{T}; source::Bool ) where {T <: Tuple{Float64, Float32}}
    i, j, l = source ? (3,2,1) : (3,4,5)

    # Distance from image to top of hill
    Δx, Δy = source ? points[i]-locs : ( locs[1]-points[i][1], points[i][2]-locs[2] )
    rr = √( Δx^2 + Δy^2 )
    # Angle from image to top of hill
    θi = atan( Δy, Δx )
    # Angle of the hillside
    Δxh, Δyh = source ? points[i]-points[j] : ( points[j][1]-points[i][1], points[i][2]-points[j][2] )
    θh = atan( Δyh, Δxh )

    if θi < θh
        # Angle of the flat
        Δxf, Δyf = source ? points[j]-points[l] : ( points[l][1]-points[j][1], points[j][2]-points[l][2] )
        θf = atan( Δyf, Δxf )
        # Angle of the image path relative to the flat
        θ = θi - θf
        # Net image source-receiver height, as needed by QQ
        Δz = rr * sin(θ)
        rr = abs( rr * cos(θ) )
        return (rr, Δz)
    else
        return nothing
    end
end

"""
    bakkernn( hills::AbstractArray{T, 1}, src_loc::AbstractArray{T,}, rec_loc::AbstractArray{T, 1}, src_flow_res::Float32, ber_flow_res::Float32, rec_flow_res::Float32, freq::Float64 ) where {T <: Tuple{Float64, Float32}}

Run hill topogrphy propagation model for spectral sound levels.
"""
function bakkernn( hills::AbstractArray{T, 1}, src_loc::AbstractArray{T, 1}, rec_loc::AbstractArray{T, 1}, src_flow_res::Float32, ber_flow_res::Float32, rec_flow_res::Float32, freq::Float64 ) where {T <: Tuple{Float64, Float64}}
    # Delany-Bazley under source, berm and reciver
    dbs = delbaz.( freq, [ src_flow_res, ber_flow_res, rec_flow_res ] )
    waveno = 2.0 * π * freq / 340.0

 # Mirror Source
    # Distance from image source to top of hill
    Δx, Δy = hills[3] .- src_loc[2]
    rr = √( Δx^2 + Δy^2 )
    # Angle from image to top of hill
    θi = atan( Δy, Δx )
    # Angle of the hillside
    Δxh, Δyh = hills[3] .- hills[2]
    θh = atan( Δyh, Δxh )

    qs = 0.0 + 0.0im
    if θi < θh
        # Angle of the flat
        Δxf, Δyf = hills[2] .- hills[1]
        θf = atan( Δyf, Δxf )
        # Angle of the image path relative to the flat
        θ = θi - θf
        # Net image source to receiver height, as needed by QQ
        Δz = rr * sin(θ)

        rr = abs( rr * cos(θ) )
        qs = qq( rr, Δz, waveno, dbs[1] )
    end
    #   rr_Δz = calc_mirror( src_loc[2], hills, source=true )
    #   qs = isnothing(rr_Δz) ? 0.0 + 0.0im : qq( rr_Δz..., waveno, dbs[1]  )

 # Mirror Receiver
    # Distance from image receiver to top of hill
    Δx = rec_loc[2][1] - hills[3][1]
    Δy = hills[3][2] - rec_loc[2][2]
    rr = √( Δx^2 + Δy^2 )
    # Angle from image to top of hill
    θi = atan( Δy, Δx )
    # Angle of the hillside
    Δxh = hills[4][1] - hills[3][1]
    Δyh = hills[3][2] - hills[4][2]
    θh = atan( Δyh, Δxh )

    qr = 0.0 + 0.0im
    if θi < θh
        # Angle of the flat
        Δxf = hills[5][1] - hills[4][1]
        Δyf = hills[4][2] - hills[5][2]
        θf = atan( Δyf, Δxf )
        # Angle of the image path relative to the flat
        θ = θi - θf
        # Net image source to receiver height, as needed by QQ
        Δz = rr * sin(θ)
        rr = abs( rr * cos(θ) )
        qr = qq( rr, Δz, waveno, dbs[1] )
    end
    #   rr_Δz = calc_mirror( rec_loc[2], hills, source=false )
    #   qr = isnothing(rr_Δz) ? 0.0 + 0.0im : qq( rr_Δz..., waveno, dbs[1]  )

 # Wedge Angle
    Δx, Δz = hills[3] .- hills[2]
    θ1 = atan( Δx, Δz )
    Δx = hills[4][1] - hills[3][1]
    Δz = hills[3][2] - hills[4][2]
    θ2 = atan( Δx, Δz )
    θ = θ1 + θ2

    pt = 0.0 + 0.0im
    f0dir = 0
    f1dir = 0
    f0refl = 0
    f1refl = 0
    for i in 1:4
        src_i = ( i == 1 || i == 3 ) ? 1 : 2
        rec_i = i <= 2 ? 1 : 2

        relev_refl_factor = i == 2 ? qs :
                                i == 3 ? qr :
                                    i == 4 ? qs * qr : 1.0 + 0.0im

        # Get length and angle of path from source to top of wedge
        Δx, Δz = hills[3] .- src_loc[src_i]
        # Distance
        rh0 = √( Δx^2 + Δz^2 )
        # Angle, clockwise from straight down
        θh = atan( Δx, Δz )
        f0 = θh - θ1

        Δx = rec_loc[rec_i][1] - hills[3][1]
        Δz = hills[3][2] - rec_loc[rec_i][2]
        rh1 = √( Δx^2 + Δz^2 )
        # Angle, counterclockwise from straight down
        θ = atan( Δx, Δz )
        f1 = 2π - θ1 - θh

        if i == 1
            f0dir = f0
            f1dir = f1
        end
        if i == 4
            f0refl = f0
            f1refl = f1
        end

        tot_propag_path = rh0 + rh1

        if f0 > 0 && ( (f1 + θ) < 2.0π )
            h_over_wedgeleg = sin(f0) * tot_propag_path
            wedge_impedence1 = qq( tot_propag_path, h_over_wedgeleg, waveno, dbs[2] )
            h_over_wedgeleg = sin( 2 * π - f1 - θ ) * tot_propag_path
            wedge_impedence2 = qq( tot_propag_path, h_over_wedgeleg, waveno, dbs[2] )

            a = rh0 * rh1 / tot_propag_path
            any = 2.0 - θ / π
            aalast = Ref{Float32}(0.0f0)
            dif1 = diffraction!( aalast, tot_propag_path, a, f1-f0, -1.0, any, waveno )
            dif2 = diffraction!( aalast, tot_propag_path, a, f1+f0, -1.0, any, waveno ) * wedge_impedence1
            dif3 = diffraction!( aalast, tot_propag_path, a, f1+f0, 1.0, any, waveno ) * wedge_impedence2
            dif4 = diffraction!( aalast, tot_propag_path, a, f1-f0, 1.0, any, waveno ) * wedge_impedence1 * wedge_impedence2
            pl = dif1 + dif2 + dif3 + dif4

            # NON SONO CERTO I DUE MODI SIANO EQUIVALENTI
            #   pl = diffraction( tot_propag_path, a, f1-f0, -1.0, any, waveno ) +
            #        diffraction( tot_propag_path, a, f1+f0, -1.0, any, waveno ) * wedge_impedence1 +
            #        diffraction( tot_propag_path, a, f1+f0, 1.0, any, waveno ) * wedge_impedence2 +
            #        diffraction( tot_propag_path, a, f1-f0, 1.0, any, waveno ) * wedge_impedence1 * wedge_impedence2

            #   pl = sum( diffraction.(
            #               tot_propag_path,
            #               a,
            #               [ f1-f0, f1+f0, f1+f0, f1-f0 ],
            #               [  -1.0,  -1.0,   1.0,   1.0 ],
            #               any,
            #               waveno
            #             ) .* [ 1.0, wedge_impedence1, wedge_impedence2, wedge_impedence1*wedge_impedence2 ] )
            pl *= relev_refl_factor
            pt += pl
        end
    end

 # Direct path source to receiver
    if π + f0dir - f1dir > 0
        Δx, Δz = rec_loc[1] .- src_loc[1]
        rd = √( Δx^2 + Δz^2 )
        po = ℯ^complex(0.0, waveno*rd) / complex(rd, 0.0)
        pt += po
    end
 # Path mirrored source to receiver
    if π + f0refl - f1dir > 0
        Δx, Δz = rec_loc[1] .- src_loc[2]
        rd = √( Δx^2 + Δz^2 )
        po = ℯ^complex(0.0, waveno*rd) * qq( rd, Δz, waveno, dbs[1] ) / complex(rd, 0.0)
        pt += po
    end
 # Path source to mirrored receiver
    if π + f0dir - f1refl > 0
        Δx, Δz = rec_loc[2] .- src_loc[1]
        rd = √( Δx^2 + Δz^2 )
        po = ℯ^complex(0.0, waveno*rd) * qq( rd, Δz, waveno, dbs[1] ) / complex(rd, 0.0)
        pt += po
    end

    Δx, Δz = rec_loc[1] .- src_loc[1]
    rd = √( Δx^2 + Δz^2 )
    level = 4.34log( (rd * abs(pt))^2 )

    return level
end

"""
    dal( first_second_dist::Float64, second_third_dist::Float64, src_h::Float32, rec_h::Float32, src_slope_α::Float64, flow_res1::Float32, flow_res2::Float32, freq::Float64 )

Run valley topography propagation model for spectral sound levels.
"""
function dal( first_second_dist::Float64, second_third_dist::Float64, src_h::Float32, rec_h::Float32, src_slope_α::Float64, flow_res1::Float32, flow_res2::Float32, freq::Float64 )
    # Delany-Bazley for source and receiver leg
    dbs = delbaz.( freq, Float32[flow_res1, flow_res2] )
    waveno = 2.0 * π * freq / 340.0

    calc_r( ra, rb, α ) = √( ra^2 + rb^2 - 2.0 * ra * rb * cos(α) )

    rh0 = √( src_h^2 + first_second_dist^2 )
    rh1 = √( rec_h^2 + second_third_dist^2 )
    r1 = rh0 + rh1
    f0 = atan(src_h, first_second_dist)
    f1 = 2.0 * π - src_slope_α - atan(rec_h, second_third_dist)

    rd = calc_r( rh0, rh1, f0-f1 )
    direct_field = ℯ^complex(0.0, waveno*rd) / complex(rd, 0.0)

    h_over_wedgeleg = sin(f0) * r1
    wedge_impedence1 = qq( r1, h_over_wedgeleg, waveno, dbs[1] )
    h_over_wedgeleg = sin(2.0 * π - f1 - src_slope_α) * r1
    wedge_impedence2 = qq( r1, h_over_wedgeleg, waveno, dbs[2] )

    a = rh0 * rh1 / r1
    any = 2.0 - src_slope_α / π
    aalast = Ref{Float32}(0.0)
    pl = diffraction!( aalast, r1, a, f1-f0, -1.0, any, waveno ) +
        ( diffraction!( aalast, r1, a, f1+f0, -1.0, any, waveno ) * wedge_impedence1 ) +
        ( diffraction!( aalast, r1, a, f1+f0, 1.0, any, waveno ) * wedge_impedence2 ) +
        ( diffraction!( aalast, r1, a, f1-f0, 1.0, any, waveno ) * wedge_impedence1 * wedge_impedence2 )
 #=
    pl = sum( diffraction!.(
                Ref(aalast),
                r1,
                a,
                Float64[ f1-f0, f1+f0, f1+f0, f1-f0 ],
                Float64[  -1.0,  -1.0,   1.0,   1.0 ],
                any,
                waveno
              ) .* Float64[ 1.0, wedge_impedence1, wedge_impedence2, wedge_impedence1*wedge_impedence2 ] )
 =#
    if f1-f0 < π
        pl += direct_field
    end

    if f1+f0 < π
        rs = calc_r( rh0, rh1, f1+f0 )
        q = qq( rs, src_h+rh1*sin(π-f1), waveno, dbs[1] )
        ps = ℯ^complex(0.0, waveno*rs) / complex(rs, 0.0) * q
        pl += ps
    end

    θ = atan(rec_h, first_second_dist)

    if ( f1 - f0 + 2.0 * θ ) < π
        rr = calc_r( rh0, rh1, f1-f0+2.0*θ )
        q = qq(rr, rec_h+rh0*sin(f0 + src_slope_α - π), waveno, dbs[2] )
        pr = ℯ^complex(0.0, waveno*rr) / complex(rr, 0.0) * q
        pl += pr
    end

    if ( f1 + f0 + 2.0 * θ ) < π
        rb = calc_r( rh0, rh1, f1+f0+2.0*θ )
        q1 = qq(rb, src_h+rh1*sin(2.0 * src_slope_α - 3.0 * π + f1), waveno, dbs[1] )
        q2 = qq(rb, src_h+rh0*sin(-f0 + src_slope_α - π), waveno, dbs[2] )
        pb = ℯ^complex(0.0, waveno*rb) / complex(rb, 0.0) * q1 * q2
        pl += pb
    end

    return 4.34log( (rd*abs(pl))^2 )
end

"""
    onCut( distances::AbstractArray{T1, 1}, heights::AbstractArray{T2, 1}, impdcs::AbstractArray{T2}, src_h::Float32, rec_h::Float32, nfreq::Int64, freqs::AbstractArray{T1, 1} ) where {T1 <: Float64, T2 <: Float32}

Return the attenuation of a sound along a terrain cut reppresented by vectors `distances`, `heights` and `impedences`.
"""
function oneCut( distances::AbstractArray{T1, 1}, heights::AbstractArray{T2, 1}, impdcs::AbstractArray{T2}, src_h::Float32, rec_h::Float32, nfreq::Int64, freqs::AbstractArray{T1, 1} ) where {T1 <: Float64, T2 <: Float32}

    ihard = 0
    isoft = 0
    flow_ress = Float32[0.0f0, 0.0f0]
    atten = 0.0
    attenuations = Float64[]
    profile = Tuple{Float64, Float32}[]
    flowpr = Tuple{Float64, Float32}[]
    ignd = Int64[]

 # ====================================================== Section to Process Profile =====================================================================

    for i in 1:length(distances)
        push!( profile, (distances[i], heights[i]) )
        push!( flowpr, (distances[i], impdcs[i]) )
        if flowpr[i][2] <= 1000.0
            push!( ignd, 0 )
            isoft += 1
            flow_ress[1] += flowpr[i][2]
        else
            push!( ignd, 1 )
            ihard += 1
            flow_ress[2] += flowpr[i][2]
        end
    end
    flow_ress[1] /= max(isoft, 1)
    flow_ress[2] /= max(ihard, 1)

    if flow_ress[1] == 0.0
        flow_ress[1] = 200.0
    end
    if flow_ress[2] == 0.0
        flow_ress[2] = 10.0^6
    end

 # =======================================================================================================================================================

    # Points of interest:
    #                          min near src ||   max  || min near rec
    #                                x   z  || x   z  || x   z
    points = Tuple{Int64, Float64}[ (0, 0.0), (0, 0.0), (0, 0.0) ]
    # Max point (second point)
    points[2] = minmax( profile, src_h, rec_h )[2]
    # Min point near source (first point)
    points[1] = minmax( profile[ 1:points[2][1] ], src_h, 0.0f0 )[1]
    # Min point near receiver (third point)
    ntemp = length(profile) - points[2][1]
    #   xs = [ point[1] for point in profile[ points[2][1]+1:ntemp ] ]
    #   res = minmax( xs, 0.0, rec_h )[1]
    res = minmax( profile[ points[2][1]+1:ntemp ], 0.0f0, rec_h )[1]
    points[3] = ( points[2][1] - 1 + res[1], res[2] )

    hillxz = Tuple{Float64, Float64}[ profile[1] ]
    klocs = Int64[1]
    nmm = 1
    for i in 1:3
        push!( hillxz, profile[ points[i][1] ] )
        if points[i][1] != klocs[nmm]
            nmm += 1
            if length(klocs) >= nmm
                klocs[nmm]= points[i][1]
            else
                push!( klocs, points[i][1] )
            end
        end
    end

    push!( hillxz, profile[end] )
    if klocs[nmm] < length(profile)
        nmm += 1
        if length(klocs) >= nmm
            klocs[nmm] = length(profile)
        else
            push!( klocs, length(profile) )
        end
    end

    if points[1][1] == points[2][1]
        hillxz[2] = hillxz[1]
    end

    if points[3][1] == points[2][1]
        hillxz[4] = hillxz[5]
    end

    if hillxz[1][1] == hillxz[2][1]
        dx, dy = hillxz[3] .- hillxz[1]
        hillxz[2] = hillxz[1] .+ 0.1 .* (dx, dy)
    end

    if hillxz[4][1] == hillxz[5][1]
        dx, dy = hillxz[5] .- hillxz[3]
        hillxz[4] = hillxz[3] .+ 0.1 .* (dx, dy)
    end

 # ================================================= Profile and Supporting Stuff Established ============================================================

  # ----------------------------------------------- Level Model -----------------------------------------------------

    if nmm <= 2
        ax = profile[1][1]
        ay = 0.0
        az = profile[1][2]
        ox = profile[end][1]
        oy = 0.0
        oz = profile[end][2]

        d = √( (ax-ox)^2 + (ay-oy)^2 )

        for j in 1:nfreq
            if ihard > 0
    # NON MI E' CHIARO IL SENSO DI "duml"
                atten = egal( d/2, d/2, src_h, rec_h, flow_ress[2], flow_ress[2], 0.0, 0.0, 0.0, freqs[j])[2]
                atten = attenh
            end
            if isoft > 0
    # NON MI E' CHIARO IL SENSO DI "duml"
                attens = egal( d/2, d/2, src_h, rec_h, flow_ress[1], flow_ress[1], 0.0, 0.0, 0.0, freqs[j])[2]
                atten = attens
            end
            if ihard > 0 && isoft > 0
                atten = varysurf( distances, ignd, src_h, rec_h, attens, attenh )
            end
            push!( attenuations, atten )
        end
    end

  # ------------------------------------------------- Hill Model ----------------------------------------------------

    zz2 = 0.0
    zcrit = 0.0
    if nmm == 3
        # Total distance dist and distance to second point dsl
        dist = profile[ klocs[3] ][1] - profile[ klocs[1] ][1]
        dsl = profile[ klocs[2] ][1] - profile[ klocs[1] ][1]
        # Altitudes at the three points
        zz1 = profile[ klocs[1] ][2]
        zz2 = profile[ klocs[2] ][2]
        zz3 = profile[ klocs[3] ][2]
        zcrit = zz1 + (zz3-zz1) * dsl / dist
    end

    if nmm > 3 || ( nmm == 3 && zz2 >= zcrit )
        # Set up the source and receiver locations, normal to the corresponding plateaus
        cosθ = 1.0
        sinθ = 0.0
        Δx, Δz = hillxz[2] .- hillxz[1]
        if Δx > 0
            θ = atan(Δz, Δx)
            cosθ = cos(θ)
            sinθ = sin(θ)
        end

        # Source location is hs above the start of the terrain cut
        srcloc = [ hillxz[1] .+ (0, src_h) ]
        # Reflect the original source image
        push!( srcloc, srcloc[1] .+ (2*src_h*cosθ).*(sinθ, -cosθ) )

        cosθ = 1.0
        sinθ = 0.0
        Δx, Δz = hillxz[5] .- hillxz[4]
        if Δx > 0
            θ = atan(Δz, Δx)
            cosθ = cos(θ)
            sinθ = sin(θ)
        end

        # Right over the end of the cut receiver
        recloc = [ hillxz[5] .+ (0, rec_h) ]
        # reflect the original receiver image
        push!( recloc, recloc[1] .+ (2*rec_h*cosθ).*(sinθ, -cosθ) )

        for j in 1:nfreq
            if ihard > 0
                attenh = bakkernn( hillxz, srcloc, recloc, flow_ress[2], flow_ress[2], flow_ress[2], freqs[j] )
                atten = attenh
            end
            if isoft > 0
                attens = bakkernn( hillxz, srcloc, recloc, flow_ress[1], flow_ress[1], flow_ress[1], freqs[j] )
                atten = attens
            end
            if ihard > 0 && isoft > 0
                atten = varysurf( distances, ignd, src_h, rec_h, attens, attenh )
            end
            push!( attenuations, atten )
        end
    end

  # ------------------------------------------------ Valley Model ---------------------------------------------------

    if nmm == 3 && zz2 < zcrit
        diff1 = profile[ klocs[2] ] .- profile[ klocs[1] ]
        diff2 = profile[ klocs[3] ] .- profile[ klocs[2] ]

        α1 = atan( diff1[2], diff1[1] )
        α2 = atan( diff2[2], diff2[1] )
        α = α2 - α1 + π
        # d0, d1 = @. √sum( [diff1^2, diff2^2] )
        d0 = @. √( diff1[1]^2 + diff1[2]^2 )
        d1 = @. √( diff2[1]^2 + diff2[2]^2 )

        for j in 1:nfreq
            attenh = 0.0
            attens = 0.0
            if ihard > 0
                attenh = dal( d0, d1, src_h, rec_h, α, flow_ress[2], flow_ress[2], freqs[j] )
                atten = attenh
            end
            if isoft > 0
                attens = dal( d0, d1, src_h, rec_h, α, flow_ress[1], flow_ress[1], freqs[j] )
                atten = attens
            end
            if ihard > 0 && isoft > 0
                atten = varysurf( distances, ignd, src_h, rec_h, attens, attenh )
            end
            push!( attenuations, atten )
        end
    end
    return attenuations
end

# Based on "https://www.geeksforgeeks.org/dda-line-generation-algorithm-computer-graphics/"
"""
    DDA( dtm::GeoArrays.GeoArray{Float32}, impedences::GeoArrays.GeoArray{Float32}, r0::Int64, c0::Int64, rn::Int64, cn::Int64 )

Digital Differential Analyzer, rasterizes a line from indexes (`r0`, `c0`) to (`rn`,`cn`) returning all the cells in `dtm` and `impedences` (The two rasters must
contain the heights and the resistivity concerning the same terrain) crossed by the line.
"""
function DDA( dtm::AbstractArray{Float32}, impedences::AbstractArray{Float32}, r0::Int64, c0::Int64, rn::Int64, cn::Int64 )
    Δr = rn - r0
    Δc = cn - c0
    steps = max( abs(Δr), abs(Δc) )
    r_inc = Δr / steps
    c_inc = Δc / steps
    r = Float64(r0)
    c = Float64(c0)
    heigths_profile = Float32[]
    impedences_profile = Float32[]
    coords_profile = Tuple{Float64, Float64}[]
    for i in 1:steps
        rint, cint = round.(Int64, [r, c])
        push!( heigths_profile, dtm[rint, cint][1] )
        push!( impedences_profile, impedences[rint, cint] )
        push!( coords_profile, Tuple{Float64, Float64}(ga.coords(dtm, Int64[rint, cint])) )
        r += r_inc
        c += c_inc
    end
    return heigths_profile, impedences_profile, coords_profile
end

#=
function ground_loss( dtm::GeoArrays.GeoArray{Float32}, x0::Float64, y0::Float64, frequency::Int64, noData::Float32 )
    # Cell dimensions
    Δx::Float64, Δy::Float64 = ( ga.coords( dtm, size(dtm)[1:2] ) .- ga.coords( dtm, [1,1] ) ) ./ size(dtm)[1:2]
    # Source cell
    r0, c0 = ga.indices( dtm, [x0, y0] )
    # Source height
    h0::Float32 = dtm[r0, c0][1]
    # Maximum radius according to transmission loss
    max_radius = ceil(10^(dB / 20))
    # Number of cell fitting the radius of the area of effect of the sound
    cell_num = ceil( Int64, max_radius / max(Δx, Δy) )
    # Limits of the area
    row_begin = r0 - cell_num
    row_end = r0 + cell_num
    col_begin = c0 - cell_num
    col_end = c0 + cell_num
    # Trivial profiles
    # Profile containing the heights
    heights_profiles = push!(
        Vector{Vector{Float32}}(),
        dtm[ r0, c0:-1:col_begin ], # x axis first half
        dtm[ r0, c0:col_end ], # x axis second half
        dtm[ r0:-1:row_begin, c0 ], # y axis first half
        dtm[ r0:row_end, c0 ], # y axis second half
        [ dtm[r0+i, c0-i][1] for i in 0:-1:-cell_num ], # first diagonal first half
        [ dtm[r0+i, c0-i][1] for i in 0:cell_num ], # first diagonal second half
        [ dtm[r0+i, c0+i][1] for i in 0:-1:-cell_num ], # second diagonal first half
        [ dtm[r0+i, c0+i][1] for i in 0:cell_num ] # second diagonal second half
    )
    # Profile containing the positions
    coords_profiles = push!(
        Vector{Vector{Tuple{Float64, Float64}}}(),
        # x and y axes
        [ Tuple{Float64, Float64}(ga.coords(dtm, [r0, col])) for col in c0:-1:col_begin ],
        [ Tuple{Float64, Float64}(ga.coords(dtm, [r0, col])) for col in c0:col_end ],
        [ Tuple{Float64, Float64}(ga.coords(dtm, [row, c0])) for row in r0:-1:row_begin ],
        [ Tuple{Float64, Float64}(ga.coords(dtm, [row, c0])) for row in r0:row_end ],
        # diagonals
        [ Tuple{Float64, Float64}(ga.coords(dtm, [r0+i, c0-i])) for i in 0:-1:-cell_num ],
        [ Tuple{Float64, Float64}(ga.coords(dtm, [r0+i, c0-i])) for i in 0:cell_num ],
        [ Tuple{Float64, Float64}(ga.coords(dtm, [r0+i, c0+i])) for i in 0:-1:-cell_num ],
        [ Tuple{Float64, Float64}(ga.coords(dtm, [r0+i, c0+i])) for i in 0:cell_num ]
    )
    # Vector containing the indexes of the points on first quadrant of the border of the area of interest
    endpoints = vcat(
        [ (row_begin, col) for col in c0+1:col_end-1 ],
        [ (row, col_end) for row in row_begin+1:r0 ]
    )
    r, c = (0, 0)
    dists = nothing
    # Matrix with the resulting intenisty levels on the area of interest
    intenisty_matrix = fill( noData, row_end - row_begin + 1, col_end - col_begin + 1 )

 # PLACEHOLDER PER I PROFILI OTTENUTI DAL RASTER DELLE IMPEDENZE
    impdcs = zeros(Float32, length(heights_profiles[1]))

    # Compute and insert the values for the trivial profiles( x and y axes and diagonal profiles )
    @inbounds for (heights, coords) in zip(heights_profiles, coords_profiles)
        # Vector holding the distances from the source to the current point
        dists = [ distance((x0, y0), coords[1]) ]
        for j in 2:length(heights)
            # Add distance of the current point from source to the distances vector
            push!( dists, distance((x0, y0), coords[j]) )
            # Current cell
            r, c = ga.indices( dtm, [ coords[j]... ] )
            # If the cell should have a value
            if dtm[r, c] != noData
                intenisty_matrix[ ( (r, c) .- (row_begin, col_begin) .+ 1 )... ] = onCut( view(dists, 1:j), view(heights, 1:j), view(impdcs, 1:j), h0, heights[j], 1, frequency )[end]
            end
        end
    end
    # Compute the values for the remainder of the area
    @inbounds for point in endpoints
        for α in [0, 90, 180, 270]
            # Arrays of the heigths and the respective coordnates
            heights, impedences, coords = DDA( dtm, r0, c0, rotate_point( point..., r0, c0, α )... )
            # Compute array of the distances of each point of the profile from the source
            dists = [ distance((x0, y0), coords[1]) ]
            # Array of the resulting attenuations for each point of a single profile
            for j in 2:length(heights)
                push!( dists, distance((x0, y0), coords[j]) )
                r, c = ga.indices( dtm, [ coords[j]... ] )
                if dtm[r, c] != noData
                    intenisty_matrix[ ( (r, c) .- (row_begin, col_begin) .+ 1 )... ] = onCut( view(dists, 1:j), view(heights, 1:j), view(impdcs, 1:j), h0, heights[j], 1, frequency )[end]
                end
            end
        end
    end
    return intenisty_matrix
end

# function noise_level( dtm::GeoArrays.GeoArray{Float32}, terrain_impedence::GeoArrays.GeoArray{Float32}, x0::Float64, y0::Float64, relative_humidity::Float64, temperature_k::FLoat64, frequencies::Union{Int64, Vector{Int64}}, noData::Float32 )
function noise_level( dtm::GeoArrays.GeoArray{Float32}, x0::Float64, y0::Float64, relative_humidity::Float64, temperature_k::FLoat64, frequency::Int64, noData::Float32 )
    # Cell dimensions
    Δx::Float64, Δy::Float64 = ( ga.coords( dtm, size(dtm)[1:2] ) .- ga.coords( dtm, [1,1] ) ) ./ size(dtm)[1:2]
    # Source cell
    r0, c0 = ga.indices( dtm, [x0, y0] )
    # Source height
    h0::Float32 = dtm[r0, c0][1]
    # Maximum radius according to transmission loss
    max_radius = ceil(10^(dB / 20))
    # Number of cell fitting the radius of the area of effect of the sound
    cell_num = ceil( Int64, max_radius / max(Δx, Δy) )
    # Limits of the area
    row_begin = r0 - cell_num
    row_end = r0 + cell_num
    col_begin = c0 - cell_num
    col_end = c0 + cell_num
    # Trivial profiles
    # Profile containing the heights
    heights_profiles = push!(
        Vector{Vector{Float32}}(),
        dtm[ r0, c0:-1:col_begin ], # x axis first half
        dtm[ r0, c0:col_end ], # x axis second half
        dtm[ r0:-1:row_begin, c0 ], # y axis first half
        dtm[ r0:row_end, c0 ], # y axis second half
        [ dtm[r0+i, c0-i][1] for i in 0:-1:-cell_num ], # first diagonal first half
        [ dtm[r0+i, c0-i][1] for i in 0:cell_num ], # first diagonal second half
        [ dtm[r0+i, c0+i][1] for i in 0:-1:-cell_num ], # second diagonal first half
        [ dtm[r0+i, c0+i][1] for i in 0:cell_num ] # second diagonal second half
    )
    # Profile containing the positions
    coords_profiles = push!(
        Vector{Vector{Tuple{Float64, Float64}}}(),
        # x and y axes
        [ Tuple{Float64, Float64}(ga.coords(dtm, [r0, col])) for col in c0:-1:col_begin ],
        [ Tuple{Float64, Float64}(ga.coords(dtm, [r0, col])) for col in c0:col_end ],
        [ Tuple{Float64, Float64}(ga.coords(dtm, [row, c0])) for row in r0:-1:row_begin ],
        [ Tuple{Float64, Float64}(ga.coords(dtm, [row, c0])) for row in r0:row_end ],
        # diagonals
        [ Tuple{Float64, Float64}(ga.coords(dtm, [r0+i, c0-i])) for i in 0:-1:-cell_num ],
        [ Tuple{Float64, Float64}(ga.coords(dtm, [r0+i, c0-i])) for i in 0:cell_num ],
        [ Tuple{Float64, Float64}(ga.coords(dtm, [r0+i, c0+i])) for i in 0:-1:-cell_num ],
        [ Tuple{Float64, Float64}(ga.coords(dtm, [r0+i, c0+i])) for i in 0:cell_num ]
    )
    # Vector containing the indexes of the points on first quadrant of the border of the area of interest
    endpoints = vcat(
        [ (row_begin, col) for col in c0+1:col_end-1 ],
        [ (row, col_end) for row in row_begin+1:r0 ]
    )
    r, c = (0, 0)
    dists = nothing
    # Matrix with the resulting intenisty levels on the area of interest
    intenisty_matrix = fill( noData, row_end - row_begin + 1, col_end - col_begin + 1 )

 # PLACEHOLDER PER I PROFILI OTTENUTI DAL RASTER DELLE IMPEDENZE
    impdcs = zeros(Float32, length(heights_profiles[1]))

    # Compute and insert the values for the trivial profiles( x and y axes and diagonal profiles )
    @inbounds for (heights, coords) in zip(heights_profiles, coords_profiles)
        # Vector holding the distances from the source to the current point
        dists = [ distance((x0, y0), coords[1]) ]
        for j in 2:length(heights)
            # Add distance of the current point from source to the distances vector
            push!( dists, distance((x0, y0), coords[j]) )
            # Current cell
            r, c = ga.indices( dtm, [ coords[j]... ] )
            # If the cell should have a value
            if dtm[r, c] != noData
                intenisty_matrix[ ( (r, c) .- (row_begin, col_begin) .+ 1 )... ] =
                    # Loss of intensity due to propagation
                    transmission_loss( dists[j] ) -
                    # Loss of intensity due to the atmosphere
                    atmospheric_loss( dists[j], heights[j], relative_humidity, temperature_k, frequency ) -
                    # Loss of intenisty due to the terrain
                    onCut( view(dists, 1:j), view(heights, 1:j), view(impdcs, 1:j), h0, heights[j], 1, [dB] )[end]
            end
        end
    end
    # Compute the values for the remainder of the area
    @inbounds for point in endpoints
        for α in [0, 90, 180, 270]
            # Arrays of the heigths and the respective coordnates
            heights, impedences, coords = DDA( dtm, r0, c0, rotate_point( point..., r0, c0, α )... )
            # Compute array of the distances of each point of the profile from the source
            dists = [ distance((x0, y0), coords[1]) ]
            # Array of the resulting attenuations for each point of a single profile
            for j in 2:length(heights)
                push!( dists, distance((x0, y0), coords[j]) )
                r, c = ga.indices( dtm, [ coords[j]... ] )
                if dtm[r, c] != noData
                    intenisty_matrix[ ( (r, c) .- (row_begin, col_begin) .+ 1 )... ] =
                        # Loss of intensity due to propagation
                        transmission_loss( dists[j] ) -
                        # Loss of intensity due to the atmosphere
                        atmospheric_loss( dists[j], heights[j], relative_humidity, temperature_k, frequency ) -
                        # Loss of intenisty due to the terrain
                        onCut( view(dists, 1:j), view(heights, 1:j), view(impedences, 1:j), h0, heights[j], 1, [dB] )[end]
                end
            end
        end
    end
    return intenisty_matrix
end
=#

"""
    run_noise( dem_file::AbstractString, terrain_impedences::AbstractString, source_file::AbstractString, temperature_K::Float64, relative_humidity::Float64, dB::Float64, frequency::Float64 )

Create and save as `output_path` a raster containing the results of model of dispersion of airborne pollutants.Run a simulation of plumes of turbidity induced by dredging.

#Arguments
- `dem_file::AbstractString`: path to the raster containing the height of the terrain in each cell.
- `terrain_impedence_file::AbstractString`: path to the raster containing the terrain impedences.
- `source_file::AbstractString`: path to the shapefile containing the source point of the noise.
- `temperature_K::Float64`: outside mean temperature in kelvin.
- `relative_humidity::Float64`: relative humidity.
- `dB::Float64`: sound intenisty in decibel.
- `frequency::Float64`: frequency of the sound in hertz.
"""
# function noise_level( dtm::GeoArrays.GeoArray{Float32}, x0::Float64, y0::Float64, relative_humidity::Float64, temperature_k::FLoat64, frequency::Int64, noData::Float32 )
function run_noise(; dem_file::AbstractString, terrain_impedences_file::String="", source_file::String, temperature_K::Float64, relative_humidity::Float64,
                     intensity_dB::Float64, frequency::Float64, output_path::String=".\\noise_otput_model.tiff" )
 # Input rasters and source point
    noData = -9999.0f0
    dtm = replace( ga.read(dem_file), missing => noData )
    impedences = isempty(terrain_impedences_file) ? fill( 0.0f0, size(dtm) ) : replace( ga.read(terrain_impedences_file), missing => noData )
    src = sf.Table(source_file)
 # Coordinates of the source
    x0 = src.geometry[1].x
    y0 = src.geometry[1].y
 # Source cell
    r0, c0 = ga.indices(dtm, [x0, y0])
 # Source height
    h0 = dtm[r0, c0][1]
 # Dimensions of a cell
    Δx, Δy = ( ga.coords(dtm, size(dtm)[1:2]) .- ga.coords(dtm, [1,1]) ) ./ size(dtm)[1:2]
 # Maximum radius according to transmission loss
    max_radius = ceil(10.0^((intensity_dB - 32.0) / 20.0))
 # Number of cell fitting the radius of the area of effect of the sound
    cell_num = ceil( Int64, max_radius / max(Δx, Δy) )
 # Limits of the area
    row_begin = r0 - cell_num
    row_end = r0 + cell_num
    col_begin = c0 - cell_num
    col_end = c0 + cell_num
 # Trivial profiles
  # Heights and impedences profiles
    heights_profiles = Vector{Float32}[ Float32[] for i in 1:8 ]
    impedences_profiles = deepcopy(heights_profiles)
    coords_profiles = Vector{Tuple{Float64, Float64}}[ Tuple{Float64, Float64}[] for i in 1:8 ]
    # The index will be used to insert the profiles in the correct places
    i = 1
   # x axis
    @inbounds for int in StepRange{Int64, Int64}[ c0:-1:col_begin, c0:1:col_end ]
        for c in int
            push!(heights_profiles[i], dtm[r0, c][1])
            push!(impedences_profiles[i], impedences[r0, c][1])
            push!( coords_profiles[i], Tuple{Float64, Float64}(ga.coords(dtm, Int64[r0, c])) )
        end
        i += 1
    end
   # y axis
    @inbounds for int in StepRange{Int64, Int64}[ r0:-1:row_begin, r0:1:row_end ]
        for r in int
            push!(heights_profiles[i], dtm[r, c0][1])
            push!(impedences_profiles[i], impedences[r, c0][1])
            push!( coords_profiles[i], Tuple{Float64, Float64}(ga.coords(dtm, Int64[r, c0])) )
        end
        i += 1
    end
   # diagonals
    @inbounds for int in StepRange{Int64, Int64}[ 0:-1:-cell_num, 0:1:cell_num ]
        for j in int
            push!(heights_profiles[i], dtm[r0 + j, c0 - j][1])
            push!(heights_profiles[i+2], dtm[r0 + j, c0 + j][1])
            push!(impedences_profiles[i], impedences[r0 + j, c0 - j][1])
            push!(impedences_profiles[i+2], impedences[r0 + j, c0 + j][1])
            push!( coords_profiles[i], Tuple{Float64, Float64}(ga.coords(dtm, Int64[r0 + j, c0 - j])) )
            push!( coords_profiles[i+2], Tuple{Float64, Float64}(ga.coords(dtm, Int64[r0 + j, c0 + j])) )
        end
        i += 1
    end
 # Vector containing the indexes of the points on first quadrant of the border of the area of interest
    endpoints = vcat(
        Tuple{Int64, Int64}[ (row_begin, col) for col in c0+1:col_end-1 ],
        Tuple{Int64, Int64}[ (row, col_end) for row in row_begin+1:r0 ]
    )
 # Matrix with the resulting intenisty levels on the area of interest
    intensity_matrix = Float32[ noData for i in 1:(row_end - row_begin + 1), j in 1:(col_end - col_begin + 1) ]
 # Compute and insert the values for the trivial profiles( x and y axes and diagonal profiles )
    @inbounds for (heights, impdcs, coords) in zip(heights_profiles, impedences_profiles, coords_profiles)
     # Vector holding the distances from the source to the current point
        dists = Float64[ distance((x0, y0), coords[1]) ]
     # The following cycle allows to compute the attenuations on each cell of the profiles
        for j in 2:length(heights)
         # Add distance of the current point from source to the distances vector
            push!( dists, distance((x0, y0), coords[j]) )
         # Current cell
            r, c = ga.indices(dtm, [coords[j]...])
         # If the cell should have a value
            if dtm[r, c] != noData
             # Attenuation due to terrain
                ground_loss = oneCut( view(dists, 1:j), view(heights, 1:j), view(impdcs, 1:j), h0, heights[j], 1, Float64[frequency] )[end]
             # `ground_loss` returns a value smaller than 0 for invisible/unreachable terrain
                if ground_loss >= 0.0
                 # Resulting sound intensity on the cell
                    intensity = intensity_dB - transmission_loss(dists[j]) - ground_loss
                 # Intensity values below 32 dB are ignored
                    if intensity >= 32.0
                        intensity_matrix[ ( (r, c) .- (row_begin, col_begin) .+ 1 )... ] = Float32(intensity)
                    end
                end
            end
        end
    end
 # Compute the values for the remainder of the area
    @inbounds for point in endpoints
        for α in Int64[0, 90, 180, 270]
         # Arrays of the heigths and the respective coordnates
            heights, impdcs, coords = DDA( dtm, impedences, r0, c0, rotate_point( point..., r0, c0, α )... )
         # Array of the distances of each point of the profile from the source
            dists = Float64[ distance((x0, y0), coords[1]) ]
         # Compute the array of the resulting attenuations for each point of a single profile
            for j in 2:length(heights)
             # Add to the distance vector as required at the j-th iteration
                push!( dists, distance((x0, y0), coords[j]) )
             # Current cell
                r, c = ga.indices(dtm, [coords[j]...])
             # If the cell should have a value
                if dtm[r, c] != noData
                 # Attenuation due to terrain
                    ground_loss = oneCut( view(dists, 1:j), view(heights, 1:j), view(impdcs, 1:j), h0, heights[j], 1, Float64[frequency] )[end]
                 # `ground_loss` returns a value smaller than 0 for invisible/unreachable terrain
                    if ground_loss >= 0.0
                     # Resulting sound intensity on the cell
                        intensity = intensity_dB - transmission_loss(dists[j]) - ground_loss
                     # Intensity values below 32 dB are ignored
                        if intensity >= 32.0
                            intensity_matrix[ ( (r, c) .- (row_begin, col_begin) .+ 1 )... ] = Float32(intensity)
                        end
                    end
                end
            end
        end
    end
    geotransform = agd.getgeotransform(dtm)
    geotransform[[1, 4]] = Functions.toCoords(geotransform, row_begin, col_begin)
    Functions.writeRaster(intensity_matrix, agd.getdriver("GTiff"), geotransform, dtm.crs.val, noData, output_path)
end

end # module

max_stars_repo_name: DavideFavaro/Tirocinio

max_issues_repo_name: DavideFavaro/Tirocinio

max_forks_repo_name: DavideFavaro/Tirocinio


text:
# Generated functions
using Dictionaries
function retrieve_code_info(sigtypes, world = Base.get_world_counter())
  S = Tuple{map(s -> Core.Compiler.has_free_typevars(s) ? typeof(s.parameters[1]) : s, sigtypes)...}
  _methods = Base._methods_by_ftype(S, -1, world)
  isempty(_methods) && @error("method $(sigtypes) does not exist, may-be run it once")
  type_signature, raw_static_params, method = _methods[1] # method is the same as we would get by invoking methods(+, (Int, Int)).ms[1]

  # this provides us with the CodeInfo
  method_instance = Core.Compiler.specialize_method(method, type_signature, raw_static_params, false)
  code_info = Core.Compiler.retrieve_code_info(method_instance)
end

struct Calls
    stamps::Vector{Float64} # contains the time stamps
    event::Vector{Symbol}  # name of the function that is being recorded
    startstop::Vector{Symbol} # if the time stamp corresponds to start or to stop
    i::Ref{Int}
end

function Calls(n::Int)
    Calls(Vector{Float64}(undef, n+1), Vector{Symbol}(undef, n+1), Vector{Symbol}(undef, n+1), Ref{Int}(0))
end

global const to = Calls(100)

function Base.show(io::IO, calls::Calls)
    for i in 1:calls.i[]
        println(io, calls.stamps[i] - calls.stamps[1],"  ", calls.startstop[i],"  ",calls.event[i])
    end
end

function Base.push!(calls::Calls, s::Symbol, ev::Symbol)
    n = calls.i[] = calls.i[] + 1
    n > length(calls.stamps) && return
    calls.event[n] = ev
    calls.startstop[n] = s
    calls.stamps[n] = time()
end

reset!(calls::Calls) = calls.i[] = 0

struct Context{T<:Union{Nothing, Vector{Symbol}}}
    functions::T
end
Context() = Context(nothing)

overdubbable(ex::Expr) = false
overdubbable(ex) = false
timable(ex::Expr) = ex.head == :call
timable(ex) = false

rename_args(ex, slot_vars, ssa_vars) = ex
rename_args(ex::Expr, slot_vars, ssa_vars) = Expr(ex.head, rename_args(ex.args, slot_vars, ssa_vars)...)
rename_args(args::AbstractArray, slot_vars, ssa_vars) = map(a -> rename_args(a, slot_vars, ssa_vars), args)
rename_args(r::Core.ReturnNode, slot_vars, ssa_vars) = Core.ReturnNode(rename_args(r.val, slot_vars, ssa_vars))
rename_args(a::Core.SlotNumber, slot_vars, ssa_vars) = slot_vars[a.id]
rename_args(a::Core.SSAValue, slot_vars, ssa_vars) = ssa_vars[a.id]

assigned_vars(ex) = []
assigned_vars(ex::Expr) = assigned_vars(ex.args)
assigned_vars(args::AbstractArray) = mapreduce(assigned_vars, vcat, args)
assigned_vars(r::Core.ReturnNode) = assigned_vars(r.val)
assigned_vars(a::Core.SlotNumber) = []
assigned_vars(a::Core.SSAValue) = [a.id]

exportname(ex::GlobalRef) = ex.name
exportname(ex::Expr) = ex.args[1]

overdub(ctx::Context, f::Core.IntrinsicFunction, args...) = f(args...)

@generated function overdub(ctx::Context, f::F, args...) where {F}
    ci = retrieve_code_info((F, args...))
    slot_vars = Dict(enumerate(ci.slotnames))
    # ssa_vars = Dict(i => gensym(:left) for i in 1:length(ci.code))
    ssa_vars = Dict(i => Symbol(:L, i) for i in 1:length(ci.code))
    used = assigned_vars(ci.code) |> distinct
    exprs = []
    for i in 1:length(args)
        push!(exprs, Expr(:(=), ci.slotnames[i+1], :(args[$(i)])))
    end
    for (i, ex) in enumerate(ci.code)
        ex = rename_args(ex, slot_vars, ssa_vars)
        if ex isa Core.ReturnNode
            push!(exprs, Expr(:return, ex.val))
            continue
        end
        if timable(ex)
            fname = exportname(ex)
            fname = :(Symbol($(fname)))
            push!(exprs, Expr(:call, :push!, :to, :(:start), fname))
            ex = overdubbable(ex) ? Expr(:call, :overdub, :ctx, ex.args...) : ex
            ex = i ∈ used ? Expr(:(=) , ssa_vars[i], ex) : ex
            push!(exprs, ex)
            push!(exprs, Expr(:call, :push!, :to, :(:stop), fname))
        else
            ex = i ∈ used ? Expr(:(=) , ssa_vars[i], ex) : ex
            push!(exprs, ex)
        end
    end
    r = Expr(:block, exprs...)
    @show r
    # println("  ")
    r
end

function foo(x, y)
   z = x * y
   z + sin(y)
end
reset!(to)
overdub(Context(), foo, 1.0, 1.0)
to

reset!(to)
overdub(ctx, Base.Math.sin_kernel, 1.0)

ext: jl


text:
using Match

# Возвращает вынесенный за скобки константный множитель
# Подразумевается, что выражение уже упрощено
consts(x) = (1, x)
function consts(ex::Expr)
    _ex = copy(ex)
    c = consts!(_ex)
    @match length(_ex.args[2:end]) begin
        0 => (c, 1)
        1 => (c, _ex.args[2])
        _ => (c, simplify(_ex))
    end
end

consts!(x) = 1
function consts!(ex::Expr)
    op = ex.args[1]
    args = ex.args[2:end]

    if op == :*
        constants, new_args = partition(args, x -> x isa Number)
        ex.args = [op, new_args...]

        inner_consts = [consts!(a) for a ∈ new_args]
        return reduce(*, vcat(constants, inner_consts); init=1)
    elseif op == :/
        #TODO: improve (?)
        if args[1] isa Number
            ex.args[2] = 1
            return args[1] / consts!(args[2])
        elseif args[2] isa Number
            ex.args[3] = 1
            return consts!(args[1]) / args[2]
        end
    end

    1
end


text:
using Distributions, StatsFuns, StatsBase

include("portfolio_initialization_template.jl")

const numVariables = 1000 + 1
const numAssets = numVariables-1

#*============ MODEL PARAMETERS ================
const mean = Array{Float64}(numAssets)
const stdev = Array{Float64}(numAssets)
for i = 1:numAssets
	mean[i] = 1.05 + 0.3*(numAssets-i)/(numAssets-1)
	stdev[i] = (0.05 + 0.6*(numAssets-i)/(numAssets-1))/3
end

const sigma = zeros(Float64, numAssets, numAssets)
for i = 1:numAssets
	sigma[i,i] = stdev[i]
end
#*===============================================

# compute probability that constraints aren't satisfied
function computeRiskLevel(y_opt::Array{Float64})

	riskLevel::Float64 = 1.0 - normcdf((sum(mean[i]*y_opt[i] for i = 1:numAssets) - y_opt[numAssets+1])/norm(sigma*(y_opt[1:numAssets])))

	return riskLevel
end

# evaluate the constraint values for given decision vector x
# and realization of the random variables xi
function evaluateConstraint(x::Array{Float64},xi::Array{Float64})

	return evaluateConstraint(x,xi,1.0)
end

# evaluate the constraint values for given decision vector x,
# realization of the random variables xi, and constraint scalings
function evaluateConstraint(x::Array{Float64},xi::Array{Float64},scalings::Float64)

	con::Float64 = (x[numVariables] - xi'*x[1:numVariables-1])/scalings

	return con
end

# evaluate constraint gradients
function evaluateConstraintGradients(x::Array{Float64},xi::Array{Float64})

	return evaluateConstraintGradients(x,xi,1.0)
end

# evaluate constraint gradients
function evaluateConstraintGradients(x::Array{Float64},xi::Array{Float64},scalings::Float64)

	grad_con = zeros(Float64,numVariables)
	for j = 1:numVariables-1
		grad_con[j] = -xi[j]
	end
	grad_con[numVariables] = 1
	grad_con /= scalings

	return grad_con
end

# generate random samples of the stock returns
function generateRandomSamples(numSamples::Int64)

	xi = Array{Float64}(numAssets,numSamples)
	for i = 1:numAssets
		xi[i,:] = rand(Normal(mean[i],stdev[i]), numSamples)
	end

	return xi
end

# project on to the set {x[1:n] >= 0, sum(x[1:n]) = 1, x[n+1] >= target}
function projectOntoFeasibleSet(z::Array{Float64},target::Float64,proj_params::Array{Float64})

	const n::Int64 = numVariables-1
	y = deepcopy(z)

	v = Array{Float64}(1)
	v[1] = y[1]
	w = Array{Float64}(0)
	rho::Float64 = y[1] - 1.0

	for j = 2:n
		if(y[j] > rho)
			rho = rho + (y[j] - rho)/(size(v,1) + 1)
			if(rho > y[j] - 1.0)
				push!(v,y[j])
			else
				append!(w,v)
				v = [y[j]]
				rho = y[j] - 1.0
			end
		end
	end

	if(size(w,1) > 0)
		for j = 1:size(w,1)
			if(w[j] > rho)
				push!(v,w[j])
				rho = rho + (w[j] - rho)/size(v,1)
			end
		end
	end

	numDeleted::Int64 = 1

	while(numDeleted > 0)
		initialCard::Int64 = size(v,1)
		numDeleted = 0
		for j = 1:initialCard
			if(v[j-numDeleted] <= rho)
				rho = rho + (rho - v[j-numDeleted])/(size(v,1)-1)
				deleteat!(v,j-numDeleted)
				numDeleted += 1
			end
		end
	end

	x_proj = zeros(numVariables)
	for j = 1:n
		x_proj[j] = max(0,y[j] - rho)
	end
	x_proj[numVariables] = max(target,z[numVariables])

	return x_proj, proj_params
end

# get the stochastic gradient of the approximation
function getStochasticGradient(x::Array{Float64},numSamples::Int64,mu::Float64,tau::Float64,scalings::Float64)

	xi = generateRandomSamples(numSamples)

	stochgrad = zeros(numVariables)

	for samp = 1:numSamples
		con::Float64 = evaluateConstraint(x,xi[:,samp],scalings)

		grad_con = evaluateConstraintGradients(x,xi[:,samp],scalings)

		baseFactor1::Float64 = tau/((sqrt(mu)*exp(0.5*tau*con) + sqrt(1.0/mu)*exp(-0.5*tau*con))^2)

		stochgrad += baseFactor1*grad_con/numSamples
	end

	return stochgrad
end

# get the stochastic gradient of the approximation
function getStochasticGradientForFixedSample(x::Array{Float64},xi::SubArray{Float64,2,Array{Float64,2},Tuple{Base.Slice{Base.OneTo{Int64}},UnitRange{Int64}},true},mu::Float64,tau::Float64,scalings::Float64)

	numSamples::Int64 = size(xi,2)

	stochgrad = zeros(numVariables)

	for samp = 1:numSamples
		con::Float64 = evaluateConstraint(x,xi[:,samp],scalings)

		grad_con = evaluateConstraintGradients(x,xi[:,samp],scalings)

		baseFactor1::Float64 = tau/((sqrt(mu)*exp(0.5*tau*con) + sqrt(1.0/mu)*exp(-0.5*tau*con))^2)

		stochgrad += baseFactor1*grad_con/numSamples
	end

	return stochgrad

end

# estimate constraint scalings
function estimateScalings(x::Array{Float64},numSamples::Int64)

	scaling_tol::Float64 = 1E-06

	xi = generateRandomSamples(numSamples)

	constraint_values = zeros(numSamples)
	for samp = 1:numSamples
		constraint_values[samp] = evaluateConstraint(x,xi[:,samp])
	end

	scalings::Float64 = median(abs.(constraint_values))

	scalings = max(scalings,scaling_tol)

	return scalings
end

# estimate Lipschitz constant of gradient of approximation
function estimateCompositeLipschitzConstant(x_ref::Array{Float64},numSamplesForLipVarEst::Int64,numGradientSamples::Int64,batchSize::Int64,mu::Float64,tau::Float64,scalings::Float64,objectiveBound::Float64)

	proj_params = zeros(Float64,3)

	sample_radius::Float64 = norm(x_ref)/10.0

	xi = generateRandomSamples(batchSize)

	L_avg::Float64 = 0.0
	for batch = 1:batchSize
		xi_start::Int64 = batch
		xi_end::Int64 = batch

		L_max::Float64 = 0.0
		for samp = 1:numSamplesForLipVarEst
			x_1 = pickPointOnSphere(x_ref,sample_radius)
			x_1,~ = projectOntoFeasibleSet(x_1,objectiveBound,proj_params)

			x_2 = pickPointOnSphere(x_ref,sample_radius)
			x_2,~ = projectOntoFeasibleSet(x_2,objectiveBound,proj_params)

			@views grad_1 = getStochasticGradientForFixedSample(x_1,xi[:,xi_start:xi_end],mu,tau,scalings)
			@views grad_2 = getStochasticGradientForFixedSample(x_2,xi[:,xi_start:xi_end],mu,tau,scalings)

			L_est::Float64 = norm(grad_2 - grad_1)/norm(x_2 - x_1)

			if(L_est > L_max)
				L_max = L_est
			end
		end

		L_avg += L_max/batchSize
	end

	return L_avg
end

# estimate step length factor
function estimateStepLength(x_ref::Array{Float64},numGradientSamples::Int64,
							numSamplesForLipVarEst::Int64,batchSize::Int64,maxNumIterations::Int64,
							minNumReplicates::Int64,mu::Float64,tau::Float64,
							scalings::Float64,objectiveBound::Float64)

	proj_params = zeros(Float64,3)

	sample_radius::Float64 = norm(x_ref)/10.0

	rho_est::Float64 = estimateCompositeLipschitzConstant(x_ref,numSamplesForLipVarEst,numGradientSamples,batchSize,mu,tau,scalings,objectiveBound)

	variance_est::Float64 = 0.0
	for samp = 1:numSamplesForLipVarEst
		x_1 = pickPointOnSphere(x_ref,sample_radius)
		x_1,~ = projectOntoFeasibleSet(x_1,objectiveBound,proj_params)

		xi = generateRandomSamples(numGradientSamples*batchSize)

		variance_1::Float64 = 0.0
		for batch = 1:batchSize
			xi_start::Int64 = (batch-1)*numGradientSamples + 1
			xi_end::Int64 = batch*numGradientSamples

			@views grad_1 = getStochasticGradientForFixedSample(x_1,xi[:,xi_start:xi_end],mu,tau,scalings)

			variance_1 += norm(grad_1)^2/batchSize
		end

		if(variance_1 > variance_est)
			variance_est = variance_1
		end
	end

	gamma::Float64 = 1/sqrt(rho_est*variance_est)

	stepLength::Float64 = gamma/sqrt((maxNumIterations+1.0)*minNumReplicates)

	return stepLength
end

# given a reference point, pick a random point within a given radius
function pickPointOnSphere(x_orig::Array{Float64},radius::Float64)

	if(radius < 1E-09)
		error("Sampling radius in pickPointOnSphere is almost zero! Set manually.")
	end

	n::Int64 = size(x_orig,1)
	unifRand = rand(1)
	factor1::Float64 = radius*(unifRand[1])^(1.0/n)
	normRand = rand(Normal(0,1), n)
	factor2::Float64 = norm(normRand)

	x_samp = x_orig + factor1*normRand/factor2

	return x_samp
end

max_stars_repo_name: rohitkannan/SA-for-CCP

max_issues_repo_name: rohitkannan/SA-for-CCP

max_forks_repo_name: rohitkannan/SA-for-CCP


text:
#!/usr/bin/env julia

# Copyright 2022 John T. Foster
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#    http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
module project2

using MPI
using RecipesBase
using Plots

struct Grid{T}
    dx::T
    dy::T
    xrange::NTuple{2, T}
    yrange::NTuple{2, T}
    u::Matrix{T}
end

function Grid(nx::Integer=10, ny::Integer=10,
              xrange::NTuple{2, Real}=(0.0, 1.0),
              yrange::NTuple{2, Real}=(0.0, 1.0))
    xmin, xmax = xrange
    ymin, ymax = yrange
    dx = (xmax - xmin) / (nx - 1)
    dy = (ymax - ymin) / (ny - 1)
    t = typeof(dx)
    u = zeros(t, (nx, ny))
    Grid{t}(dx, dy, xrange, yrange, u)
end

function set_boundary_condition!(g::Grid, f::Function=(x,y) -> 0.0;
        sides::Union{String, Tuple{Vararg{String}}})

    xmin, xmax = g.xrange
    ymin, ymax = g.yrange
    nx = size(g.u)[2]
    ny = size(g.u)[1]

    if sides == "all"
        sides = ("bottom", "top", "left", "right")
    elseif typeof(sides) == String
        sides = (sides,)
    end

    for side in sides

        if side == "left"
            g.u[:, 1] = f.(xmin, LinRange(ymin, ymax, ny))
        elseif side == "right"
            g.u[:, end] = f.(xmax, LinRange(ymin, ymax, ny))
        elseif side == "bottom"
            g.u[1, :] = f.(LinRange(xmin, xmax, nx), ymin)
        elseif side == "top"
            g.u[end, :] = f.(LinRange(xmin, xmax, nx), ymax)
        end
    end

end

function reset!(g::Grid)
    g.u .= zero(typeof(g.dx))
end

function iterate!(g::Grid)

    u = g.u
    nx, ny = size(u)

    dx2, dy2 = g.dx^2, g.dy^2

    err = 0.0

    for i = 2:(nx-1)
        for j = 2:(ny-1)

            tmp = u[i, j]

            u[i,j] = ((u[i-1, j] + u[i+1, j]) * dy2 +
                      (u[i, j-1] + u[i, j+1]) * dx2) / (dx2 + dy2) / 2

            diff = u[i,j] - tmp

            err += diff * diff
        end
    end

    sqrt(err)
end

function solve!(g::Grid; max_iterations::Integer=10000,
                         tolerance::AbstractFloat=1.0e-12,
                         quiet::Bool=false)

    for i=1:max_iterations
        error = iterate!(g)
        if error < tolerance
            if !quiet
                println("Solution converged in $(i) iterations")
            end
            break
        end
    end
end

function partition_grid(g::Union{Grid, Nothing}, comm::MPI.Comm)::Grid
    # Return Grid{t}(dx, dy, xrange, yrange, my_u) on each rank
end

function solve!(my_g::Grid, comm::MPI.Comm; max_iterations::Integer=10000,
                                            tolerance::AbstractFloat=1.0e-12,
                                            quiet::Bool=false)
   #Returns nothing, should update my_g in place
end

function get_solution!(g::Union{Grid, Nothing}, my_g::Grid, comm)
   #Returns nothing, should update g in place
end

@recipe function f(g::Grid)
    xmin, xmax = g.xrange
    ymin, ymax = g.yrange
    nx = size(g.u)[2]
    ny = size(g.u)[1]
    LinRange(xmin, xmax, nx), LinRange(ymin, ymax, ny), g.u
end

function run_serial()
    g = Grid(101, 101)
    set_boundary_condition!(g, (x, y) -> 10, sides=("top", "bottom"))
    solve!(g)

    ENV["GKSwstype"]="nul"
    p = contourf(g, aspect_ratio=:equal, xlims=g.xrange);
    savefig(p, "laplace_serial.png");
end

function run_parallel()
    MPI.Init()
    comm = MPI.COMM_WORLD
    rank = MPI.Comm_rank(comm)

    ENV["GKSwstype"]="nul"

    if rank == 0
        g = Grid(101, 101)
        set_boundary_condition!(g, (x, y) -> 10, sides=("left", "right"))
    else
        g = nothing
    end
    my_g = partition_grid(g, comm)
    solve!(my_g, comm)
    get_solution!(g, my_g, comm)
    if rank == 0
        p = contourf(g, aspect_ratio=:equal, xlims=g.xrange);
        savefig(p, "laplace_parallel.png");
    end
end

export run_serial, run_parallel, set_boundary_condition!, partition_grid, Grid, solve!, get_solution!

end


text:
__precompile__()
module SimpleIntegrals
export integral
using ArgCheck

integral(xs, ys; window=nothing) = _integral(xs, ys, window)
_integral(xs, ys, ::Nothing) = trapezoid(xs, ys)

function _integral(xs, ys,window)
    @argcheck length(window) == 2
    a,b = window
    trapezoid(xs,ys,a,b)
end

function compute_trapezoid_return_type(::Type{T}, ::Type{S}) where {T,S}
    p = (zero(T) * zero(S)) / 2
    typeof(p + p)
end

function trapezoid(xs, ys)
    @argcheck issorted(xs)
    index = eachindex(xs, ys)[2:end]
    T = compute_trapezoid_return_type(eltype(xs), eltype(ys))
    trapezoid_kernel(T,xs,ys,index)
end

function trapezoid_kernel(::Type{T},xs, ys, index) where {T}
    ret = zero(T)::T
    @simd for i in index
        Δ::T = (xs[i] - xs[i-1]) * (ys[i] + ys[i-1])
        @inbounds ret += Δ
    end
    ret::T / 2
end
function trapezoid_kernel(::Type{T},xs::AbstractRange, ys, index) where {T}
    ret = zero(T)
    ret += ys[first(index)-1]/2
    ret -= ys[last(index)]/2
    ret += sum(view(ys, index))
    ret *= step(xs)
    ret
end

function trapezoid(xs, ys, a,b)
    @argcheck !isempty(xs)
    @argcheck first(xs) <= a <= b <= last(xs)
    @argcheck eachindex(xs) == eachindex(ys) DimensionMismatch
    @argcheck issorted(xs)
    T = compute_trapezoid_return_type(eltype(xs), eltype(ys))
    length(xs) == 1 && return zero(T)
    i1 = searchsortedfirst(xs, a)
    i2 = searchsortedlast(xs, b)
    @check get(xs, i1-1, -Inf*oneunit(a)) <= a <= xs[i1]
    @check xs[i2] <= b <= get(xs, i2+1, Inf*oneunit(b))
    if i2 < i1
         ya = linterpol(a, xs[i1-1], xs[i1], ys[i1-1], ys[i1])
         yb = linterpol(b, xs[i2], xs[i2+1], ys[i2], ys[i2+1])
         return (b-a) * (yb + ya) /2
    else
        inner_index = (i1+1):i2
        ret = trapezoid_kernel(T, xs, ys, inner_index)
        if i1 != first(eachindex(xs))
            x0 = a
            @assert xs[i1-1] <= x0 <= xs[i1]
            y0 = linterpol(x0, xs[i1-1], xs[i1], ys[i1-1], ys[i1])
            ret += (xs[i1] - x0) * (ys[i1] + y0)/2
        end
        if i2 != last(eachindex(xs))
            x0 = b
            @assert xs[i2] <= x0 <= xs[i2+1]
            y0 = linterpol(x0, xs[i2], xs[i2+1], ys[i2], ys[i2+1])
            ret += (x0 - xs[i2]) * (y0 + ys[i2])/2
        end
    end
    ret
end

function linterpol(x,x1,x2,y1,y2)
    @argcheck x1 <= x <= x2
    w1 = (x2-x)/(x2-x1)
    w2 = 1 - w1
    w1*y1 + w2*y2
end

end # module


text:
struct SanchezLacombeParam <: EoSParam
    Mw::SingleParam{Float64}
    r::SingleParam{Float64}
    v_ast::SingleParam{Float64}
    omega::SingleParam{Float64}
    epsilon_ast::PairParam{Float64}
end

abstract type SanchezLacombeModel <: SAFTModel end
@newmodel SanchezLacombe SanchezLacombeModel SanchezLacombeParam
const SL = SanchezLacombe

const k_B = 1.380649e-23 # m2 kg s-2 K-1
#im supposing that this is the expression for total helmholtz energy, including ideal terms
function eos(model::SanchezLacombeModel,V,T,z=SA[1.0])
    rᵢ = model.params.r.values
    ϵᵢⱼ∗ = model.params.epsilon_ast.values
    ϵᵢ∗ = model.params.epsilon_ast.diagvalues
    ω = model.params.omega.values
    v∗ = v_mix(model,V,T,z)
    Nᵢ = z
    N = ∑(Nᵢ)
    r = dot(Nᵢ,rᵢ)/N
    ϵ∗ = eps_mix(model,V,T,z)
    T∗ = ϵ∗/k_B
    T̃ = T/T∗
    V∗ = r*N*v∗
    ṽ = V/V∗
    ṽinv = 1/ṽ
    res = (ṽ-1)*log(one(ṽ)-ṽ) + log(ṽinv)/r
    #this step could be optimized
    #also, this fails in the case of Ni[i] = 0
    ϕᵢ = rᵢ.*Nᵢ ./ (r*N)
    res += sum(log.(ϕᵢ ./ω) .* ϕᵢ ./ rᵢ)
    res *= T̃
    res -= ṽinv
    return res
end

function eps_mix(model::SanchezLacombeModel,V,T,z=SA[1.0])
    rᵢ = model.params.r.values
    ϵᵢⱼ∗ = model.params.epsilon_ast.values
    ϵᵢ∗ = model.params.epsilon_ast.diagvalues
    Nᵢ = z
    N = ∑(Nᵢ)
    r = dot(Nᵢ,rᵢ)/N
    Σei = ∑(rᵢ[i]*Nᵢ[i]*ϵᵢ∗[i] for i in 1:length(z))/(r*N)
    return Σei
end

function v_mix(model::SanchezLacombeModel,V,T,z=SA[1.0])
    return 1
end


text:
#
# @description Implements the Blowfish block cipher algorithm.
#
# @author Alejandro Claro (alejandro.claro@gmail.com)
#
# Copyright 2017 All rights reserved.
# Use of this source code is governed by a MIT-style license that can be found in the LICENSE file.
#
export BlowfishCipher
export BLOWFISH_BLOCK_SIZE, BLOWFISH_LEGAL_KEY_SIZES
export block_size, key_size, encrypt, decrypt

# @description Defines the Blowfish cipher data struture.
struct BlowfishCipher <: BlockCipher
  parray::Vector{Vector{UInt8}}
  sboxes::Vector{Vector{UInt32}}
  key_size::UInt64

  # @description Initialize the Blowfish cipher.
  #
  # @param {Vector{UInt8}} key The symmetric key.
  function BlowfishCipher(key::Vector{UInt8})
    p, sb = compute_blowfish_boxes(key)
    return new(p, sb, length(key))
  end
end

# @description The DES cipher block size in bytes.
const BLOWFISH_BLOCK_SIZE = UInt64(8)

# @description The DES cipher legal key sizes in bytes.
const BLOWFISH_LEGAL_KEY_SIZES = UInt64[ x for x in 4 : 56 ]

# @description Gets the cipher block size.
#
# @param {BlowfishCipher} self The cipher data struture.
#
# @return {Integer} The block size in bytes.
function block_size(self::BlowfishCipher)
  return BLOWFISH_BLOCK_SIZE
end

# @description Gets the cipher key size.
#
# @param {BlowfishCipher} self The cipher data struture.
#
# @return {Integer} The key size in bytes.
function key_size(self::BlowfishCipher)
  return self.key_size
end

# @description Encrypts the given plaintext block with the key set at initialization.
#
# @param {BlowfishCipher} self      The cipher data struture.
# @param {Vector{UInt8}}  plaintext The message to encrypt.
#
# @return {Vector{UInt8}} The resulting ciphertext block.
function encrypt(self::BlowfishCipher, plaintext::Vector{UInt8})
  if length(plaintext) != block_size(self)
    throw(ArgumentError("Invalid input plaintext length."))
  end

  return transform_blowfish_block(plaintext, self.sboxes, self.parray)
end

# @description Decrypts the given ciphertext block with the key set at initialization.
#
# @param {BlowfishCipher} self       The cipher data struture.
# @param {Vector{UInt8}}  ciphertext The message to decrypt.
#
# @return {Vector{UInt8}} The resulting plaintext block.
function decrypt(self::BlowfishCipher, ciphertext::Vector{UInt8})
  if length(ciphertext) != block_size(self)
    throw(ArgumentError("Invalid input ciphertext length."))
  end

  return transform_blowfish_block(ciphertext, self.sboxes, reverse(self.parray))
end

# PRIVATE IMPLEMENTATION #######################################################

# @description Blowfish P-array values derived from the hexadecimal digits of pi.
const BLOWFISH_PARRAY = [
  UInt8[ 0x24, 0x3F, 0x6A, 0x88 ],
  UInt8[ 0x85, 0xA3, 0x08, 0xD3 ],
  UInt8[ 0x13, 0x19, 0x8A, 0x2E ],
  UInt8[ 0x03, 0x70, 0x73, 0x44 ],
  UInt8[ 0xA4, 0x09, 0x38, 0x22 ],
  UInt8[ 0x29, 0x9F, 0x31, 0xD0 ],
  UInt8[ 0x08, 0x2E, 0xFA, 0x98 ],
  UInt8[ 0xEC, 0x4E, 0x6C, 0x89 ],
  UInt8[ 0x45, 0x28, 0x21, 0xE6 ],
  UInt8[ 0x38, 0xD0, 0x13, 0x77 ],
  UInt8[ 0xBE, 0x54, 0x66, 0xCF ],
  UInt8[ 0x34, 0xE9, 0x0C, 0x6C ],
  UInt8[ 0xC0, 0xAC, 0x29, 0xB7 ],
  UInt8[ 0xC9, 0x7C, 0x50, 0xDD ],
  UInt8[ 0x3F, 0x84, 0xD5, 0xB5 ],
  UInt8[ 0xB5, 0x47, 0x09, 0x17 ],
  UInt8[ 0x92, 0x16, 0xD5, 0xD9 ],
  UInt8[ 0x89, 0x79, 0xFB, 0x1B ]
]

# @Description BlowFish S-Box 1 values derived from the hexadecimal digits of pi.
const BLOWFISH_SBOX1 = UInt32[
  0xD1310BA6, 0x98DFB5AC, 0x2FFD72DB, 0xD01ADFB7,
  0xB8E1AFED, 0x6A267E96, 0xBA7C9045, 0xF12C7F99,
  0x24A19947, 0xB3916CF7, 0x0801F2E2, 0x858EFC16,
  0x636920D8, 0x71574E69, 0xA458FEA3, 0xF4933D7E,
  0x0D95748F, 0x728EB658, 0x718BCD58, 0x82154AEE,
  0x7B54A41D, 0xC25A59B5, 0x9C30D539, 0x2AF26013,
  0xC5D1B023, 0x286085F0, 0xCA417918, 0xB8DB38EF,
  0x8E79DCB0, 0x603A180E, 0x6C9E0E8B, 0xB01E8A3E,
  0xD71577C1, 0xBD314B27, 0x78AF2FDA, 0x55605C60,
  0xE65525F3, 0xAA55AB94, 0x57489862, 0x63E81440,
  0x55CA396A, 0x2AAB10B6, 0xB4CC5C34, 0x1141E8CE,
  0xA15486AF, 0x7C72E993, 0xB3EE1411, 0x636FBC2A,
  0x2BA9C55D, 0x741831F6, 0xCE5C3E16, 0x9B87931E,
  0xAFD6BA33, 0x6C24CF5C, 0x7A325381, 0x28958677,
  0x3B8F4898, 0x6B4BB9AF, 0xC4BFE81B, 0x66282193,
  0x61D809CC, 0xFB21A991, 0x487CAC60, 0x5DEC8032,
  0xEF845D5D, 0xE98575B1, 0xDC262302, 0xEB651B88,
  0x23893E81, 0xD396ACC5, 0x0F6D6FF3, 0x83F44239,
  0x2E0B4482, 0xA4842004, 0x69C8F04A, 0x9E1F9B5E,
  0x21C66842, 0xF6E96C9A, 0x670C9C61, 0xABD388F0,
  0x6A51A0D2, 0xD8542F68, 0x960FA728, 0xAB5133A3,
  0x6EEF0B6C, 0x137A3BE4, 0xBA3BF050, 0x7EFB2A98,
  0xA1F1651D, 0x39AF0176, 0x66CA593E, 0x82430E88,
  0x8CEE8619, 0x456F9FB4, 0x7D84A5C3, 0x3B8B5EBE,
  0xE06F75D8, 0x85C12073, 0x401A449F, 0x56C16AA6,
  0x4ED3AA62, 0x363F7706, 0x1BFEDF72, 0x429B023D,
  0x37D0D724, 0xD00A1248, 0xDB0FEAD3, 0x49F1C09B,
  0x075372C9, 0x80991B7B, 0x25D479D8, 0xF6E8DEF7,
  0xE3FE501A, 0xB6794C3B, 0x976CE0BD, 0x04C006BA,
  0xC1A94FB6, 0x409F60C4, 0x5E5C9EC2, 0x196A2463,
  0x68FB6FAF, 0x3E6C53B5, 0x1339B2EB, 0x3B52EC6F,
  0x6DFC511F, 0x9B30952C, 0xCC814544, 0xAF5EBD09,
  0xBEE3D004, 0xDE334AFD, 0x660F2807, 0x192E4BB3,
  0xC0CBA857, 0x45C8740F, 0xD20B5F39, 0xB9D3FBDB,
  0x5579C0BD, 0x1A60320A, 0xD6A100C6, 0x402C7279,
  0x679F25FE, 0xFB1FA3CC, 0x8EA5E9F8, 0xDB3222F8,
  0x3C7516DF, 0xFD616B15, 0x2F501EC8, 0xAD0552AB,
  0x323DB5FA, 0xFD238760, 0x53317B48, 0x3E00DF82,
  0x9E5C57BB, 0xCA6F8CA0, 0x1A87562E, 0xDF1769DB,
  0xD542A8F6, 0x287EFFC3, 0xAC6732C6, 0x8C4F5573,
  0x695B27B0, 0xBBCA58C8, 0xE1FFA35D, 0xB8F011A0,
  0x10FA3D98, 0xFD2183B8, 0x4AFCB56C, 0x2DD1D35B,
  0x9A53E479, 0xB6F84565, 0xD28E49BC, 0x4BFB9790,
  0xE1DDF2DA, 0xA4CB7E33, 0x62FB1341, 0xCEE4C6E8,
  0xEF20CADA, 0x36774C01, 0xD07E9EFE, 0x2BF11FB4,
  0x95DBDA4D, 0xAE909198, 0xEAAD8E71, 0x6B93D5A0,
  0xD08ED1D0, 0xAFC725E0, 0x8E3C5B2F, 0x8E7594B7,
  0x8FF6E2FB, 0xF2122B64, 0x8888B812, 0x900DF01C,
  0x4FAD5EA0, 0x688FC31C, 0xD1CFF191, 0xB3A8C1AD,
  0x2F2F2218, 0xBE0E1777, 0xEA752DFE, 0x8B021FA1,
  0xE5A0CC0F, 0xB56F74E8, 0x18ACF3D6, 0xCE89E299,
  0xB4A84FE0, 0xFD13E0B7, 0x7CC43B81, 0xD2ADA8D9,
  0x165FA266, 0x80957705, 0x93CC7314, 0x211A1477,
  0xE6AD2065, 0x77B5FA86, 0xC75442F5, 0xFB9D35CF,
  0xEBCDAF0C, 0x7B3E89A0, 0xD6411BD3, 0xAE1E7E49,
  0x00250E2D, 0x2071B35E, 0x226800BB, 0x57B8E0AF,
  0x2464369B, 0xF009B91E, 0x5563911D, 0x59DFA6AA,
  0x78C14389, 0xD95A537F, 0x207D5BA2, 0x02E5B9C5,
  0x83260376, 0x6295CFA9, 0x11C81968, 0x4E734A41,
  0xB3472DCA, 0x7B14A94A, 0x1B510052, 0x9A532915,
  0xD60F573F, 0xBC9BC6E4, 0x2B60A476, 0x81E67400,
  0x08BA6FB5, 0x571BE91F, 0xF296EC6B, 0x2A0DD915,
  0xB6636521, 0xE7B9F9B6, 0xFF34052E, 0xC5855664,
  0x53B02D5D, 0xA99F8FA1, 0x08BA4799, 0x6E85076A
]

# @Description BlowFish S-Box 2 values derived from the hexadecimal digits of pi.
const BLOWFISH_SBOX2 = UInt32[
  0x4B7A70E9, 0xB5B32944, 0xDB75092E, 0xC4192623,
  0xAD6EA6B0, 0x49A7DF7D, 0x9CEE60B8, 0x8FEDB266,
  0xECAA8C71, 0x699A17FF, 0x5664526C, 0xC2B19EE1,
  0x193602A5, 0x75094C29, 0xA0591340, 0xE4183A3E,
  0x3F54989A, 0x5B429D65, 0x6B8FE4D6, 0x99F73FD6,
  0xA1D29C07, 0xEFE830F5, 0x4D2D38E6, 0xF0255DC1,
  0x4CDD2086, 0x8470EB26, 0x6382E9C6, 0x021ECC5E,
  0x09686B3F, 0x3EBAEFC9, 0x3C971814, 0x6B6A70A1,
  0x687F3584, 0x52A0E286, 0xB79C5305, 0xAA500737,
  0x3E07841C, 0x7FDEAE5C, 0x8E7D44EC, 0x5716F2B8,
  0xB03ADA37, 0xF0500C0D, 0xF01C1F04, 0x0200B3FF,
  0xAE0CF51A, 0x3CB574B2, 0x25837A58, 0xDC0921BD,
  0xD19113F9, 0x7CA92FF6, 0x94324773, 0x22F54701,
  0x3AE5E581, 0x37C2DADC, 0xC8B57634, 0x9AF3DDA7,
  0xA9446146, 0x0FD0030E, 0xECC8C73E, 0xA4751E41,
  0xE238CD99, 0x3BEA0E2F, 0x3280BBA1, 0x183EB331,
  0x4E548B38, 0x4F6DB908, 0x6F420D03, 0xF60A04BF,
  0x2CB81290, 0x24977C79, 0x5679B072, 0xBCAF89AF,
  0xDE9A771F, 0xD9930810, 0xB38BAE12, 0xDCCF3F2E,
  0x5512721F, 0x2E6B7124, 0x501ADDE6, 0x9F84CD87,
  0x7A584718, 0x7408DA17, 0xBC9F9ABC, 0xE94B7D8C,
  0xEC7AEC3A, 0xDB851DFA, 0x63094366, 0xC464C3D2,
  0xEF1C1847, 0x3215D908, 0xDD433B37, 0x24C2BA16,
  0x12A14D43, 0x2A65C451, 0x50940002, 0x133AE4DD,
  0x71DFF89E, 0x10314E55, 0x81AC77D6, 0x5F11199B,
  0x043556F1, 0xD7A3C76B, 0x3C11183B, 0x5924A509,
  0xF28FE6ED, 0x97F1FBFA, 0x9EBABF2C, 0x1E153C6E,
  0x86E34570, 0xEAE96FB1, 0x860E5E0A, 0x5A3E2AB3,
  0x771FE71C, 0x4E3D06FA, 0x2965DCB9, 0x99E71D0F,
  0x803E89D6, 0x5266C825, 0x2E4CC978, 0x9C10B36A,
  0xC6150EBA, 0x94E2EA78, 0xA5FC3C53, 0x1E0A2DF4,
  0xF2F74EA7, 0x361D2B3D, 0x1939260F, 0x19C27960,
  0x5223A708, 0xF71312B6, 0xEBADFE6E, 0xEAC31F66,
  0xE3BC4595, 0xA67BC883, 0xB17F37D1, 0x018CFF28,
  0xC332DDEF, 0xBE6C5AA5, 0x65582185, 0x68AB9802,
  0xEECEA50F, 0xDB2F953B, 0x2AEF7DAD, 0x5B6E2F84,
  0x1521B628, 0x29076170, 0xECDD4775, 0x619F1510,
  0x13CCA830, 0xEB61BD96, 0x0334FE1E, 0xAA0363CF,
  0xB5735C90, 0x4C70A239, 0xD59E9E0B, 0xCBAADE14,
  0xEECC86BC, 0x60622CA7, 0x9CAB5CAB, 0xB2F3846E,
  0x648B1EAF, 0x19BDF0CA, 0xA02369B9, 0x655ABB50,
  0x40685A32, 0x3C2AB4B3, 0x319EE9D5, 0xC021B8F7,
  0x9B540B19, 0x875FA099, 0x95F7997E, 0x623D7DA8,
  0xF837889A, 0x97E32D77, 0x11ED935F, 0x16681281,
  0x0E358829, 0xC7E61FD6, 0x96DEDFA1, 0x7858BA99,
  0x57F584A5, 0x1B227263, 0x9B83C3FF, 0x1AC24696,
  0xCDB30AEB, 0x532E3054, 0x8FD948E4, 0x6DBC3128,
  0x58EBF2EF, 0x34C6FFEA, 0xFE28ED61, 0xEE7C3C73,
  0x5D4A14D9, 0xE864B7E3, 0x42105D14, 0x203E13E0,
  0x45EEE2B6, 0xA3AAABEA, 0xDB6C4F15, 0xFACB4FD0,
  0xC742F442, 0xEF6ABBB5, 0x654F3B1D, 0x41CD2105,
  0xD81E799E, 0x86854DC7, 0xE44B476A, 0x3D816250,
  0xCF62A1F2, 0x5B8D2646, 0xFC8883A0, 0xC1C7B6A3,
  0x7F1524C3, 0x69CB7492, 0x47848A0B, 0x5692B285,
  0x095BBF00, 0xAD19489D, 0x1462B174, 0x23820E00,
  0x58428D2A, 0x0C55F5EA, 0x1DADF43E, 0x233F7061,
  0x3372F092, 0x8D937E41, 0xD65FECF1, 0x6C223BDB,
  0x7CDE3759, 0xCBEE7460, 0x4085F2A7, 0xCE77326E,
  0xA6078084, 0x19F8509E, 0xE8EFD855, 0x61D99735,
  0xA969A7AA, 0xC50C06C2, 0x5A04ABFC, 0x800BCADC,
  0x9E447A2E, 0xC3453484, 0xFDD56705, 0x0E1E9EC9,
  0xDB73DBD3, 0x105588CD, 0x675FDA79, 0xE3674340,
  0xC5C43465, 0x713E38D8, 0x3D28F89E, 0xF16DFF20,
  0x153E21E7, 0x8FB03D4A, 0xE6E39F2B, 0xDB83ADF7
]

# @Description BlowFish S-Box 3 values derived from the hexadecimal digits of pi.
const BLOWFISH_SBOX3 = UInt32[
  0xE93D5A68, 0x948140F7, 0xF64C261C, 0x94692934,
  0x411520F7, 0x7602D4F7, 0xBCF46B2E, 0xD4A20068,
  0xD4082471, 0x3320F46A, 0x43B7D4B7, 0x500061AF,
  0x1E39F62E, 0x97244546, 0x14214F74, 0xBF8B8840,
  0x4D95FC1D, 0x96B591AF, 0x70F4DDD3, 0x66A02F45,
  0xBFBC09EC, 0x03BD9785, 0x7FAC6DD0, 0x31CB8504,
  0x96EB27B3, 0x55FD3941, 0xDA2547E6, 0xABCA0A9A,
  0x28507825, 0x530429F4, 0x0A2C86DA, 0xE9B66DFB,
  0x68DC1462, 0xD7486900, 0x680EC0A4, 0x27A18DEE,
  0x4F3FFEA2, 0xE887AD8C, 0xB58CE006, 0x7AF4D6B6,
  0xAACE1E7C, 0xD3375FEC, 0xCE78A399, 0x406B2A42,
  0x20FE9E35, 0xD9F385B9, 0xEE39D7AB, 0x3B124E8B,
  0x1DC9FAF7, 0x4B6D1856, 0x26A36631, 0xEAE397B2,
  0x3A6EFA74, 0xDD5B4332, 0x6841E7F7, 0xCA7820FB,
  0xFB0AF54E, 0xD8FEB397, 0x454056AC, 0xBA489527,
  0x55533A3A, 0x20838D87, 0xFE6BA9B7, 0xD096954B,
  0x55A867BC, 0xA1159A58, 0xCCA92963, 0x99E1DB33,
  0xA62A4A56, 0x3F3125F9, 0x5EF47E1C, 0x9029317C,
  0xFDF8E802, 0x04272F70, 0x80BB155C, 0x05282CE3,
  0x95C11548, 0xE4C66D22, 0x48C1133F, 0xC70F86DC,
  0x07F9C9EE, 0x41041F0F, 0x404779A4, 0x5D886E17,
  0x325F51EB, 0xD59BC0D1, 0xF2BCC18F, 0x41113564,
  0x257B7834, 0x602A9C60, 0xDFF8E8A3, 0x1F636C1B,
  0x0E12B4C2, 0x02E1329E, 0xAF664FD1, 0xCAD18115,
  0x6B2395E0, 0x333E92E1, 0x3B240B62, 0xEEBEB922,
  0x85B2A20E, 0xE6BA0D99, 0xDE720C8C, 0x2DA2F728,
  0xD0127845, 0x95B794FD, 0x647D0862, 0xE7CCF5F0,
  0x5449A36F, 0x877D48FA, 0xC39DFD27, 0xF33E8D1E,
  0x0A476341, 0x992EFF74, 0x3A6F6EAB, 0xF4F8FD37,
  0xA812DC60, 0xA1EBDDF8, 0x991BE14C, 0xDB6E6B0D,
  0xC67B5510, 0x6D672C37, 0x2765D43B, 0xDCD0E804,
  0xF1290DC7, 0xCC00FFA3, 0xB5390F92, 0x690FED0B,
  0x667B9FFB, 0xCEDB7D9C, 0xA091CF0B, 0xD9155EA3,
  0xBB132F88, 0x515BAD24, 0x7B9479BF, 0x763BD6EB,
  0x37392EB3, 0xCC115979, 0x8026E297, 0xF42E312D,
  0x6842ADA7, 0xC66A2B3B, 0x12754CCC, 0x782EF11C,
  0x6A124237, 0xB79251E7, 0x06A1BBE6, 0x4BFB6350,
  0x1A6B1018, 0x11CAEDFA, 0x3D25BDD8, 0xE2E1C3C9,
  0x44421659, 0x0A121386, 0xD90CEC6E, 0xD5ABEA2A,
  0x64AF674E, 0xDA86A85F, 0xBEBFE988, 0x64E4C3FE,
  0x9DBC8057, 0xF0F7C086, 0x60787BF8, 0x6003604D,
  0xD1FD8346, 0xF6381FB0, 0x7745AE04, 0xD736FCCC,
  0x83426B33, 0xF01EAB71, 0xB0804187, 0x3C005E5F,
  0x77A057BE, 0xBDE8AE24, 0x55464299, 0xBF582E61,
  0x4E58F48F, 0xF2DDFDA2, 0xF474EF38, 0x8789BDC2,
  0x5366F9C3, 0xC8B38E74, 0xB475F255, 0x46FCD9B9,
  0x7AEB2661, 0x8B1DDF84, 0x846A0E79, 0x915F95E2,
  0x466E598E, 0x20B45770, 0x8CD55591, 0xC902DE4C,
  0xB90BACE1, 0xBB8205D0, 0x11A86248, 0x7574A99E,
  0xB77F19B6, 0xE0A9DC09, 0x662D09A1, 0xC4324633,
  0xE85A1F02, 0x09F0BE8C, 0x4A99A025, 0x1D6EFE10,
  0x1AB93D1D, 0x0BA5A4DF, 0xA186F20F, 0x2868F169,
  0xDCB7DA83, 0x573906FE, 0xA1E2CE9B, 0x4FCD7F52,
  0x50115E01, 0xA70683FA, 0xA002B5C4, 0x0DE6D027,
  0x9AF88C27, 0x773F8641, 0xC3604C06, 0x61A806B5,
  0xF0177A28, 0xC0F586E0, 0x006058AA, 0x30DC7D62,
  0x11E69ED7, 0x2338EA63, 0x53C2DD94, 0xC2C21634,
  0xBBCBEE56, 0x90BCB6DE, 0xEBFC7DA1, 0xCE591D76,
  0x6F05E409, 0x4B7C0188, 0x39720A3D, 0x7C927C24,
  0x86E3725F, 0x724D9DB9, 0x1AC15BB4, 0xD39EB8FC,
  0xED545578, 0x08FCA5B5, 0xD83D7CD3, 0x4DAD0FC4,
  0x1E50EF5E, 0xB161E6F8, 0xA28514D9, 0x6C51133C,
  0x6FD5C7E7, 0x56E14EC4, 0x362ABFCE, 0xDDC6C837,
  0xD79A3234, 0x92638212, 0x670EFA8E, 0x406000E0
]

# @Description BlowFish S-Box 4 values derived from the hexadecimal digits of pi.
const BLOWFISH_SBOX4 = UInt32[
  0x3A39CE37, 0xD3FAF5CF, 0xABC27737, 0x5AC52D1B,
  0x5CB0679E, 0x4FA33742, 0xD3822740, 0x99BC9BBE,
  0xD5118E9D, 0xBF0F7315, 0xD62D1C7E, 0xC700C47B,
  0xB78C1B6B, 0x21A19045, 0xB26EB1BE, 0x6A366EB4,
  0x5748AB2F, 0xBC946E79, 0xC6A376D2, 0x6549C2C8,
  0x530FF8EE, 0x468DDE7D, 0xD5730A1D, 0x4CD04DC6,
  0x2939BBDB, 0xA9BA4650, 0xAC9526E8, 0xBE5EE304,
  0xA1FAD5F0, 0x6A2D519A, 0x63EF8CE2, 0x9A86EE22,
  0xC089C2B8, 0x43242EF6, 0xA51E03AA, 0x9CF2D0A4,
  0x83C061BA, 0x9BE96A4D, 0x8FE51550, 0xBA645BD6,
  0x2826A2F9, 0xA73A3AE1, 0x4BA99586, 0xEF5562E9,
  0xC72FEFD3, 0xF752F7DA, 0x3F046F69, 0x77FA0A59,
  0x80E4A915, 0x87B08601, 0x9B09E6AD, 0x3B3EE593,
  0xE990FD5A, 0x9E34D797, 0x2CF0B7D9, 0x022B8B51,
  0x96D5AC3A, 0x017DA67D, 0xD1CF3ED6, 0x7C7D2D28,
  0x1F9F25CF, 0xADF2B89B, 0x5AD6B472, 0x5A88F54C,
  0xE029AC71, 0xE019A5E6, 0x47B0ACFD, 0xED93FA9B,
  0xE8D3C48D, 0x283B57CC, 0xF8D56629, 0x79132E28,
  0x785F0191, 0xED756055, 0xF7960E44, 0xE3D35E8C,
  0x15056DD4, 0x88F46DBA, 0x03A16125, 0x0564F0BD,
  0xC3EB9E15, 0x3C9057A2, 0x97271AEC, 0xA93A072A,
  0x1B3F6D9B, 0x1E6321F5, 0xF59C66FB, 0x26DCF319,
  0x7533D928, 0xB155FDF5, 0x03563482, 0x8ABA3CBB,
  0x28517711, 0xC20AD9F8, 0xABCC5167, 0xCCAD925F,
  0x4DE81751, 0x3830DC8E, 0x379D5862, 0x9320F991,
  0xEA7A90C2, 0xFB3E7BCE, 0x5121CE64, 0x774FBE32,
  0xA8B6E37E, 0xC3293D46, 0x48DE5369, 0x6413E680,
  0xA2AE0810, 0xDD6DB224, 0x69852DFD, 0x09072166,
  0xB39A460A, 0x6445C0DD, 0x586CDECF, 0x1C20C8AE,
  0x5BBEF7DD, 0x1B588D40, 0xCCD2017F, 0x6BB4E3BB,
  0xDDA26A7E, 0x3A59FF45, 0x3E350A44, 0xBCB4CDD5,
  0x72EACEA8, 0xFA6484BB, 0x8D6612AE, 0xBF3C6F47,
  0xD29BE463, 0x542F5D9E, 0xAEC2771B, 0xF64E6370,
  0x740E0D8D, 0xE75B1357, 0xF8721671, 0xAF537D5D,
  0x4040CB08, 0x4EB4E2CC, 0x34D2466A, 0x0115AF84,
  0xE1B00428, 0x95983A1D, 0x06B89FB4, 0xCE6EA048,
  0x6F3F3B82, 0x3520AB82, 0x011A1D4B, 0x277227F8,
  0x611560B1, 0xE7933FDC, 0xBB3A792B, 0x344525BD,
  0xA08839E1, 0x51CE794B, 0x2F32C9B7, 0xA01FBAC9,
  0xE01CC87E, 0xBCC7D1F6, 0xCF0111C3, 0xA1E8AAC7,
  0x1A908749, 0xD44FBD9A, 0xD0DADECB, 0xD50ADA38,
  0x0339C32A, 0xC6913667, 0x8DF9317C, 0xE0B12B4F,
  0xF79E59B7, 0x43F5BB3A, 0xF2D519FF, 0x27D9459C,
  0xBF97222C, 0x15E6FC2A, 0x0F91FC71, 0x9B941525,
  0xFAE59361, 0xCEB69CEB, 0xC2A86459, 0x12BAA8D1,
  0xB6C1075E, 0xE3056A0C, 0x10D25065, 0xCB03A442,
  0xE0EC6E0E, 0x1698DB3B, 0x4C98A0BE, 0x3278E964,
  0x9F1F9532, 0xE0D392DF, 0xD3A0342B, 0x8971F21E,
  0x1B0A7441, 0x4BA3348C, 0xC5BE7120, 0xC37632D8,
  0xDF359F8D, 0x9B992F2E, 0xE60B6F47, 0x0FE3F11D,
  0xE54CDA54, 0x1EDAD891, 0xCE6279CF, 0xCD3E7E6F,
  0x1618B166, 0xFD2C1D05, 0x848FD2C5, 0xF6FB2299,
  0xF523F357, 0xA6327623, 0x93A83531, 0x56CCCD02,
  0xACF08162, 0x5A75EBB5, 0x6E163697, 0x88D273CC,
  0xDE966292, 0x81B949D0, 0x4C50901B, 0x71C65614,
  0xE6C6C7BD, 0x327A140A, 0x45E1D006, 0xC3F27B9A,
  0xC9AA53FD, 0x62A80F00, 0xBB25BFE2, 0x35BDD2F6,
  0x71126905, 0xB2040222, 0xB6CBCF7C, 0xCD769C2B,
  0x53113EC0, 0x1640E3D3, 0x38ABBD60, 0x2547ADF0,
  0xBA38209C, 0xF746CE76, 0x77AFA1C5, 0x20756060,
  0x85CBFE4E, 0x8AE88DD8, 0x7AAAF9B0, 0x4CF9AA7E,
  0x1948C25C, 0x02FB8A8C, 0x01C36AE4, 0xD6EBE1F9,
  0x90D4F869, 0xA65CDEA0, 0x3F09252D, 0xC208E69F,
  0xB74E6132, 0xCE77E25B, 0x578FDFE3, 0x3AC372E6
]

# @description The full set of Blowfish initial S-Boxes.
const BLOWFISH_SBOXES = [ BLOWFISH_SBOX1, BLOWFISH_SBOX2, BLOWFISH_SBOX3, BLOWFISH_SBOX4 ]

# @description Computes the p-array and s-boxes for Blowfish algorithm.
#
# @param key The symmetric key.
#
# @return {Vector{Vector{UInt8}}, Vector{Vector{UInt32}}} The P-array and the S-boxes.
function compute_blowfish_boxes(key::Vector{UInt8})
  if length(key) ∉ BLOWFISH_LEGAL_KEY_SIZES
    throw(ArgumentError("Invalid key length. $(length(key) * 8) bits provided."))
  end

  parray = deepcopy(BLOWFISH_PARRAY)
  sboxes = deepcopy(BLOWFISH_SBOXES)
  block  = zeros(UInt8, BLOWFISH_BLOCK_SIZE)
  index  = 0

  for j in 1 : length(parray)
    for k in 1 : length(parray[j])
      parray[j][k] = parray[j][k] ⊻ key[(index % length(key)) + 1]
      index += 1
    end
  end

  for j in 1 : 2 : length(parray)
    block = transform_blowfish_block(block, sboxes, parray)
    parray[j]     = block[1 : 4]
    parray[j + 1] = block[5 : 8]
  end

  for j in 1 : length(sboxes)
    for k in 1 : 2 : length(sboxes[j])
      block  = transform_blowfish_block(block, sboxes, parray)
      values = reinterpret(UInt32, reverse(block))
      sboxes[j][k]     = values[2]
      sboxes[j][k + 1] = values[1]
    end
  end

  return parray, sboxes
end

# @description Performs the Blowfish's block transformation.
#
# @param {Vector{UInt8}}  block  The block bytes (8 bytes).
# @param {Vector{UInt32}} sboxes The derived S-boxes.
# @param {Vector{UInt8}}  parray The derived P-array.
#
# @return {Array{UInt8}} The transformed block.
function transform_blowfish_block(block::Vector{UInt8}, sboxes::Vector{Vector{UInt32}}, parray::Vector{Vector{UInt8}})
  @assert length(block) == 8

  ln = block[1 : 4]
  rn = block[5 : 8]

  for i in 1 : 16
    ln = map(xor, ln, parray[i])
    rn = map(xor, rn, compute_blowfish_round(ln, sboxes))
    ln, rn = rn, ln
  end

  ln = map(xor, ln, parray[17])
  rn = map(xor, rn, parray[18])

  return vcat(rn, ln)
end

# @description Computes the Blowfish's f(*) function for the given half-block and S-boxes.
#
# @param {Vector{UInt8}} x      The half-block bytes (4 bytes).
# @param {Vector{UInt8}} sboxes The derived S-boxes.
#
# @return {Array{UInt8}} The round function value.
function compute_blowfish_round(x::Vector{UInt8}, sboxes::Vector{Vector{UInt32}})
  @assert length(x) == 4

  h = sboxes[1][x[1] + 1] + sboxes[2][x[2] + 1]
  y = h ⊻ sboxes[3][x[3] + 1]
  z = y + sboxes[4][x[4] + 1]

  return reverse(reinterpret(UInt8, [ z ]))
end


text:
using Laszip
using Test

# write your own tests here
#@test 1 == 1

out = [0.0 0 0; 1 0 1; 0 1 1; 1 1 2; 0 0 2];
dat2las("lixo.laz", out);

in = las2dat("lixo.laz");

@test in == out

# Remove garbage
rm("lixo.laz")

max_stars_repo_name: joa-quim/Laszip

max_issues_repo_name: joa-quim/Laszip

max_forks_repo_name: joa-quim/Laszip


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/FEMBase.jl/blob/master/LICENSE

abstract type AbstractAnalysis end
abstract type AbstractResultsWriter end

mutable struct Analysis{A<:AbstractAnalysis}
    name :: String
    problems :: Vector{Problem}
    fields :: Dict{String, AbstractField}
    results_writers :: Vector{AbstractResultsWriter}
    properties :: A
end

"""
    Analysis(A, analysis_name)

Create a new analysis of type `A`, where `A` is a subtype of `AbstractAnalysis`.
Analysis can be e.g. `Linear` for linear quasistatic analysis, `Nonlinear` for
nonlinear quasistatic analysis or `Modal` for natural frequency analysis.

# Examples

```julia
analysis = Analysis(Linear, "linear quasistatic analysis of beam structure")
```
"""
function Analysis(::Type{A}, name::String="$A Analysis") where A<:AbstractAnalysis
    analysis = Analysis{A}(name, [], Dict(), [], A())
    @info("Creating a new analysis of type $A with name `$name`.")
    return analysis
end

function add_problem!(analysis::Analysis, problems::Problem...)
    for problem in problems
        @info("Adding problem `$(problem.name)` to analysis `$(analysis.name)`.")
        push!(analysis.problems, problem)
    end
    return nothing
end

add_problems!(analysis::Analysis, problems::Union{Vector, Tuple}) = add_problem!(analysis, problems...)

"""
    add_problems!(analysis, problem...)

Add problem(s) to analysis.

# Examples

Add two problems, `beam` and `bc`, to linear quasistatic analysis:
```julia
beam = Problem(Beam, "beam structure", 6)
bc = Problem(Dirichlet, "fix", 6, "displacement")
analysis = Analysis(Linear, "linear quasistatic analysis")
add_problems!(analysis, beam, bc)
```
"""
add_problems!(analysis::Analysis, problems...) = add_problem!(analysis, problems...)

function get_problems(analysis::Analysis)
    return analysis.problems
end

function get_problem(analysis::Analysis, problem_name)
    problems = filter(p -> p.name == problem_name, get_problems(analysis))
    if length(problems) != 1
        error("Several problem with name $problem_name found from analysis $(analysis.name).")
    end
    return first(problems)
end

function add_results_writer!(analysis::Analysis, writer::W) where W<:AbstractResultsWriter
    push!(analysis.results_writers, writer)
    return nothing
end

function get_results_writers(analysis::Analysis)
    return analysis.results_writers
end

function run!(::Analysis{A}) where A<:AbstractAnalysis
    @info("This is a placeholder function for running an analysis $A for a set of problems.")
end

function write_results!(::Analysis{A}, ::W) where {A<:AbstractAnalysis, W<:AbstractResultsWriter}
    @info("Writing the results of analysis $A is not supported by a results writer $W")
    return nothing
end

function write_results!(analysis)
    results_writers = get_results_writers(analysis)
    if isempty(results_writers)
        @info("No result writers attached to the analysis $(analysis.name). " *
              "In order to get results of the analysis stored to the disk, one " *
              "must attach some results writer to the analysis using " *
              "add_results_writer!, e.g. xdmf_writer = Xdmf(\"results\"); " *
              "add_results_writer!(analysis, xdmf_writer)")
        return nothing
    end
    for results_writer in results_writers
        write_results!(analysis, results_writer)
    end
    return nothing
end


text:
getcolumn(x::AbstractVector{T}, colidx::Int) where {T} =
   colidx === 1 ? x : throw(DomainError("colidx ($colidx) > 1"))

getcolumn(x::AbstractArray{T,2}, colidx::Int) where {T} =
   1 <= colidx <= 2 ? x[:,colidx] : throw(DomainError("!(1 <= colidx ($colidx) <= 2)"))

function getcolumn(x::AbstractArray{T,3}, colidx::Tuple{Int,Int,Int}) where {T}
   if colidx[1] == 0
      x[:, colidx[2], colidx[3]]
   elseif colidx[2] == 0
      x[colidx[1], :, colidx[3]]
   elseif colidx[3] == 0
      x[colidx[1], colidx[2], :]
   else
      throw(DomainError("bad colidx ($colidx) <zero not found>"))
   end
end

# get DataFrame column from the column id as string or symbol
getcolumn(df::DataFrames.DataFrame, str::String) = getindex(df, Symbol(str))
getcolumn(df::DataFrames.DataFrame, sym::Symbol) = df[sym]

# get TimeArray column from the column id as string or symbol
getcolumn(ta::TimeSeries.TimeArray, sym::Symbol) = getindex(ta, string(sym)).values
getcolumn(ta::TimeSeries.TimeArray, str::String) = getindex(ta, str).values


text:
using PkgPage
using Pkg
using Test
import Franklin
using FranklinUtils
using NodeJS
const P = PkgPage
const F = Franklin

F.newmodule("Utils")
Base.eval(Main.Utils, quote using PkgPage end)

include("latex/basic.jl")
include("latex/environments.jl")

include("instantiate.jl")