text:
__precompile__()
module DWARF
    using ObjFileBase
    using StructIO
    using AbstractTrees

    include("constants.jl")

    import Base: read, write, zero, bswap, isequal, show, print, hash, ==
    import Base: start, next, done

    import AbstractTrees: children, printnode

    abstract DWARFHeader
    abstract DWARFCUHeader <: DWARFHeader # Compilation Unit Header
    abstract DWARFTUHeader <: DWARFHeader # Type Unit Header
    abstract DWARFARHeader <: DWARFHeader # Address Range Header
    abstract DWARFPUBHeader <: DWARFHeader
    abstract PUBTableEntry
    abstract PUBTableSet

    export AbbrevTableEntry, AbbrevTableSet, ULEB128, SLEB128,
        attributes, dies

    module DWARF32
        using DWARF
        using StructIO

        @struct immutable CUHeader <: DWARF.DWARFCUHeader
            unit_length::UInt32
            version::UInt16
            debug_abbrev_offset::UInt32
            address_size::UInt8
        end align_packed

        @struct immutable TUHeader <: DWARF.DWARFTUHeader
            unit_length::UInt32
            version::UInt16
            debug_abbrev_offset::UInt32
            address_size::UInt8
            type_signature::UInt64
            type_offset::UInt32
        end align_packed

        @struct immutable ARHeader <: DWARF.DWARFARHeader
            unit_length::UInt32
            version::UInt16
            debug_info_offset::UInt32
            address_size::UInt8
            segment_size::UInt8
        end align_packed

        @struct immutable PUBHeader <: DWARF.DWARFPUBHeader
            unit_length::UInt32
            version::UInt16
            debug_info_offset::UInt32
            debug_info_length::UInt32
        end align_packed

        immutable PUBTableEntry <: DWARF.PUBTableEntry
            offset::UInt32
            name::String
        end

        immutable PUBTableSet <: DWARF.PUBTableSet
            header::PUBHeader
            entries::Array{PUBTableEntry,1}
        end

    end

    module DWARF64
        using DWARF
        using StructIO

        @struct immutable CUHeader <: DWARF.DWARFCUHeader
            unit_length::UInt64
            version::UInt16
            debug_abbrev_offset::UInt64
            address_size::UInt8
            type_offset::UInt64
        end align_packed

        @struct immutable TUHeader <: DWARF.DWARFTUHeader
            unit_length::UInt32
            version::UInt16
            debug_abbrev_offset::UInt32
            address_size::UInt8
            type_signature::UInt64
            type_offset::UInt64
        end align_packed

        @struct immutable ARHeader <: DWARF.DWARFARHeader
            unit_length::UInt32
            version::UInt16
            debug_info_offset::UInt64
            address_size::UInt8
            segment_size::UInt8
        end align_packed

        @struct immutable PUBHeader <: DWARF.DWARFPUBHeader
            unit_length::UInt32
            version::UInt16
            debug_info_offset::UInt64
            debug_info_length::UInt64
        end align_packed

        immutable PUBTableEntry <: DWARF.PUBTableEntry
            offset::UInt64
            name::String
        end

        immutable PUBTableSet <: DWARF.PUBTableSet
            header::PUBHeader
            entries::Array{PUBTableEntry,1}
        end
    end

    ### Display
    function show(io::IO, header::DWARFCUHeader)
        is64 = isa(header,DWARF64.CUHeader)
        printentry(io,"Length","0x",hex(header.unit_length))
        printentry(io,"Version",dec(header.version))
        printentry(io,"Abbrev Offset","0x",hex(header.debug_abbrev_offset))
        printentry(io,"Address Size","0x",hex(header.address_size))
        is64 && printentry(io,"Type Offset","0x",hex(header.type_offset))
    end

    ### LEB 128 types
    abstract LEB128{T}

    Base.convert{T<:LEB128}(::Type{T},x::Int64) = T(big(x))
    Base.convert{T<:Integer, S<:LEB128}(::Type{T}, x::S) = convert(T,x.val)

    immutable ULEB128{T} <: LEB128{T}
        val::T
    end

    immutable SLEB128{T} <: LEB128{T}
        val::T
    end

    show(io::IO,s::SLEB128) = show(io,s.val)
    hash(x::LEB128, h::UInt64) = hash(x.val, h)
    isequal(a::LEB128,b::LEB128) = (a.val==b.val)
    isequal(a::LEB128,b::Integer) = (a.val==b)
    isequal(a::Integer,b::LEB128) = (b==a.val)
    ==(a::LEB128,b::LEB128) = isequal(a,b)
    ==(a::LEB128,b::Integer) = isequal(a.val, b)
    ==(a::Integer,b::LEB128) = isequal(b, a.val)
    Base.zero{T<:LEB128}(::Type{T}) = convert(T,0)
    Base.hex(x::LEB128) = hex(x.val)

    function read{T}(io::IO, ::Type{ULEB128{T}})
        v = T(0)
        shift = 0
        while true
            c = read(io,UInt8)
            ((8*sizeof(T))-2 < shift) && throw(InexactError())
            v |= T(c&0x7f)<<shift
            if (c&0x80)==0 #is last bit
                break
            end
            shift+=7
        end
        ULEB128{T}(v)
    end
    read(io::IO, ::Type{ULEB128}) = read(io,ULEB128{UInt})

    function write(io::IO, x::ULEB128)
        x = x.val
        while true
            v = (x % UInt8) & 0x7f
            x >>= 7
            if x != 0
                v |= 0x80 # = ~0x7f
            end
            write(io,v)
            x == 0 && break
        end
    end

    function decode{T}(data::Array{UInt8,1}, offset, ::Type{ULEB128{T}})
        v = T(0)
        shift = 0
        i=0
        while true
            c = data[offset+i]
            i+=1
            ((8*sizeof(T))-2 < shift) && throw(InexactError())
            v |= T(c&0x7f)<<shift
            if (c&0x80)==0 #is last bit
                break
            end
            shift+=7
        end
        (offset+i,T(v))
    end

    function decode{T}(data::Array{UInt8,1}, offset, ::Type{SLEB128{T}})
        v = T(0)
        shift = 0
        c=0
        i=0
        while true
            c = data[offset+i]
            i+=1
            ((8*sizeof(T))-2 < shift) && throw(InexactError())
            v |= T(c&0x7f)<<shift
            shift+=7
            if (c&0x80)==0 #is last bit
                break
            end
        end
        if (c & 0x40) != 0
            v |= -(T(1)<<shift)
        end

        (offset+i,T(v))
    end

    function read{T}(io::IO, ::Type{SLEB128{T}})
        v = T(0)
        shift = 0
        c=0
        while(true)
            c = read(io,UInt8)
            ((8*sizeof(T))-2 < shift) && throw(InexactError())
            v |= T(c&0x7f)<<shift
            shift+=7
            if (c&0x80)==0 #is last bit
                break
            end
        end
        if (c & 0x40) != 0
            v |= -(T(1)<<shift)
        end

        SLEB128{T}(v)
    end
    read(io::IO, ::Type{SLEB128}) = read(io, SLEB128{Int})

    function write(io::IO, x::SLEB128)
        x = x.val
        more = true
        while more
            v = (x % UInt8) & 0x7f
            x >>= 7
            if ((x == 0 && ((v & 0x40) == 0)) || (x == -1 && ((v & 0x40) > 0)))
                more = false
            else
                v |= 0x80
            end
            write(io,v)
        end
    end

    const attr_color = :cyan

    immutable AttributeSpecification
        # These are both ULEB128s on disk, but we keep them as
        # integers in memory to avoid having to deal with BigInt
        name::UInt64
        form::UInt8
    end
    isequal(a::AttributeSpecification,b::AttributeSpecification) = (a.name == b.name)&&(a.form == b.form)
    ==(a::AttributeSpecification,b::AttributeSpecification) = isequal(a,b)

    function read(io::IO,::Type{AttributeSpecification},endianness::Symbol)
        name = read(io,ULEB128)
        form = read(io,ULEB128)
        AttributeSpecification(name,form)
    end

    function read(io::IO,header::DWARF.DWARFCUHeader,a::AttributeSpecification,endianness::Symbol)
        generic = read(io,Attributes.form2gattrT(a.form),header,a.name,a.form,endianness)
        # TODO: Return Actual Attributes
        generic
    end

    immutable AbbrevTableEntry
        # These are technically ULEB128 in the on disk format, but we keep them
        # as UInt64 in memory which should be large enough for anything ever
        # encoutered in reality
        code::UInt64
        tag::UInt64
        has_children::UInt8
        attributes::Array{AttributeSpecification,1}
    end

    immutable AbbrevTableSet
        entries::Array{AbbrevTableEntry,1}
    end
    zero(::Type{AbbrevTableEntry}) = AbbrevTableEntry(0,0,UInt8(0),Array(AttributeSpecification,0))
    const zero_entry = zero(AbbrevTableEntry)

    function readorskip
    end
    include("dies.jl")
    # Attributes
    module Attributes
        include("forms.jl")
    end
    using .Attributes

    module Expressions
        # TODO
        using DWARF
        using StructIO
        import StructIO: fix_endian

        type StateMachine{T}
            stack::Array{T,1}
            StateMachine(stack::Array{T,1}) = new(stack)
            StateMachine() = new(T[])
        end

        function operands(addr_type,opcode,opcodes,i,endianness)
            if opcode == DWARF.DW_OP_addr
                operand = fix_endian(reinterpret(addr_type,opcodes[i:(i+sizeof(addr_type)-1)])[1],endianness)
                i += sizeof(addr_type)
            elseif opcode == DWARF.DW_OP_const1u || opcode == DWARF.DW_OP_pick
                operand = convert(addr_type,opcodes[i])
                i += 1
            elseif opcode == DWARF.DW_OP_const1s
                operand = convert(addr_type,opcodes[i])
                # Yes, this is actually different from the above, since we need to sign extend properly
                push!(s.stack,convert(addr_type,fix_endian(reinterpret(Int8,opcodes[i]),endianness)))
                i += 1
            elseif opcode == DWARF.DW_OP_const2u # 1 2-byte constant
                operand = convert(addr_type,fix_endian(reinterpret(UInt16,opcodes[i:i+1])[1],endianness))
                i+=2
            elseif opcode == DWARF.DW_OP_const2s  || opcode == DWARF.DW_OP_bra || opcode == DWARF.DW_OP_skip # 1 2-byte constant
                operand = convert(addr_type,fix_endian(reinterpret(Int16,opcodes[i:i+1])[1],endianness))
                i+=2
            elseif opcode == DWARF.DW_OP_const4u # 1 4-byte constant
                operand = convert(addr_type,fix_endian(reinterpret(UInt32,opcodes[i:i+3])[1],endianness))
                i+=4
            elseif opcode == DWARF.DW_OP_const4u # 1 4-byte constant
                operand = convert(addr_type,fix_endian(reinterpret(Int32,opcodes[i:i+3])[1],endianness))
                i+=4
            elseif opcode == DWARF.DW_OP_const8u # 1 8-byte constant
                operand = convert(addr_type,fix_endian(reinterpret(UInt64,opcodes[i:i+7])[1],endianness))
                i+=8
            elseif opcode == DWARF.DW_OP_const8s # 1 8-byte constant
                operand = convert(addr_type,fix_endian(reinterpret(Int64,opcodes[i:i+7])[1],endianness))
                i+=8
            elseif opcode == DWARF.DW_OP_constu || opcode == DWARF.DW_OP_plus_uconst ||
                    opcode == DWARF.DW_OP_regx || opcode == DWARF.DW_OP_piece
                (i,operand) = DWARF.decode(opcodes,i,ULEB128{addr_type})
            elseif opcode == DWARF.DW_OP_consts || opcode == DWARF.DW_OP_fbreg ||
                opcode >= DWARF.DW_OP_breg0 && opcode <= DWARF.DW_OP_breg31
                (i,operand) = DWARF.decode(opcodes,i,SLEB128{typeof(signed(addr_type(0)))})
            elseif opcode == DWARF.DW_OP_bregx
                (i,reg) = DWARF.decode(opcodes,i,ULEB128{UInt})
                (i,offset) = DWARF.decode(opcodes,i,SLEB128{addr_type})
                operand = (reg,offset)
            elseif opcode == DWARF.DW_OP_bit_piece
                (i,reg) = DWARF.decode(opcodes,i,ULEB128{UInt})
                (i,offset) = DWARF.decode(opcodes,i,ULEB128{addr_type})
                operand = (reg,offset)
            else
                return (i,)
            end
            return (i,operand)
        end

        function evaluate_generic_instruction{T}(s::StateMachine{T},opcodes,i,
                getreg_func::Function,getword_func,addr_func,endianness::Symbol)
            opcode = opcodes[i]
            i+=1
            if opcode == DWARF.DW_OP_deref
                addr = pop!(s.stack)
                push!(s.stack,getword_func(addr))
            elseif in(opcode,(DWARF.DW_OP_const1u,DWARF.DW_OP_const1s,DWARF.DW_OP_const2u,
                              DWARF.DW_OP_const2s,DWARF.DW_OP_const4u,DWARF.DW_OP_const4s,DWARF.DW_OP_const8u,
                              DWARF.DW_OP_const8s,DWARF.DW_OP_constu,DWARF.DW_OP_consts))
                (i,val) = operands(T,opcode,opcodes,i,endianness)
                push!(s.stack,convert(T,val))
            elseif opcode == DWARF.DW_OP_addr
                (i,val) = operands(T,opcode,opcodes,i,endianness)
                push!(s.stack,addr_func(convert(T,val)))
            elseif opcode == DWARF.DW_OP_dup
                push!(s.stack,s.stack[length(s.stack)-1])
            elseif opcode == DWARF.DW_OP_pick
                (i,val) = operands(T,opcode,opcodes,i,endianness)
                push!(s.stack,s.stack[val])
            elseif opcode == DWARF.DW_OP_swap
                top = length(s.stack)
                val = s.stack[top]
                s.stack[top] = s.stack[top-1]
                s.stack[top-1] = val
            elseif opcode == DWARF.DW_OP_rot
                top = length(s.stack)
                val = s.stack[top]
                s.stack[top] = s.stack[top-1]
                s.stack[top-1] = s.stack[top-2]
                s.stack[top-2] = val
            elseif opcode == DWARF.DW_OP_xderef
                error("Unimplemented")
            elseif opcode == DWARF.DW_OP_abs
                push!(s.stack,abs(signed(pop!(s.stack))))
            elseif opcode == DWARF.DW_OP_and
                push!(s.stack,pop!(s.stack)&pop!(s.stack))
            elseif opcode == DWARF.DW_OP_div
                top = pop!(s.stack)
                push!(s.stack,pop!(s.stack)/top)
            elseif opcode == DWARF.DW_OP_minus
                top = pop!(s.stack)
                push!(s.stack,pop!(s.stack)-top)
            elseif opcode == DWARF.DW_OP_minus
                top = pop!(s.stack)
                push!(s.stack,pop!(s.stack)%top)
            elseif opcode == DWARF.DW_OP_plus
                push!(s.stack,pop!(s.stack)+pop!(s.stack))
            elseif opcode == DWARF.DW_OP_mul
                push!(s.stack,pop!(s.stack)*pop!(s.stack))
            elseif opcode == DWARF.DW_OP_neg
                push!(s.stack,-(signed(pop!(s.stack))))
            elseif opcode == DWARF.DW_OP_not
                push!(s.stack,~(pop!(s.stack)))
            elseif opcode == DWARF.DW_OP_plus_uconst
                (i,val) = DWARF.decode(opcodes,i,ULEB128{UInt})
                push!(s.stack,pop!(s.stack)+UInt(val))
            elseif opcode == DWARF.DW_OP_shl
                top = pop!(s.stack)
                push!(s.stack,pop!(s.stack)<<top)
            elseif opcode == DWARF.DW_OP_shr
                top = pop!(s.stack)
                push!(s.stack,pop!(s.stack)>>top)
            elseif opcode == DWARF.DW_OP_shra
                top = pop!(s.stack)
                push!(s.stack,convert(T,signed(pop!(s.stack))>>top))
            elseif opcode == DWARF.DW_OP_xor
                top = pop!(s.stack)
                push!(s.stack,pop!(s.stack)$top)
            elseif opcode == DWARF.DW_OP_le
                top = pop!(s.stack)
                push!(s.stack,convert(T,pop!(s.stack)<=top))
            elseif opcode == DWARF.DW_OP_ge
                top = pop!(s.stack)
                push!(s.stack,convert(T,pop!(s.stack)>=top))
            elseif opcode == DWARF.DW_OP_eq
                top = pop!(s.stack)
                push!(s.stack,convert(T,pop!(s.stack)==top))
            elseif opcode == DWARF.DW_OP_lt
                top = pop!(s.stack)
                push!(s.stack,convert(T,pop!(s.stack)<top))
            elseif opcode == DWARF.DW_OP_gt
                top = pop!(s.stack)
                push!(s.stack,convert(T,pop!(s.stack)>top))
            elseif opcode == DWARF.DW_OP_ne
                top = pop!(s.stack)
                push!(s.stack,convert(T,pop!(s.stack)!=top))
            elseif opcode == DWARF.DW_OP_skip
                i += fix_endian(reinterrept(Int16,opcodes[i:i+1])[1],endianness)
            elseif opcode == DWARF.DW_OP_bra
                (i,skip) = operands(T,opcode,opcodes,i,endianness)
                if(pop!(s.stack) != 0)
                    i += skip
                end
            elseif opcode == DWARF.DW_OP_call2 || opcode == DWARF.DW_OP_call4 || opcode == DWARF.DW_OP_call_ref
                error("Unimplemented")
            elseif opcode >= DWARF.DW_OP_lit1 && opcode <= DWARF.DW_OP_lit31
                push!(s.stack,opcode-DWARF.DW_OP_lit1+1)
            elseif opcode >= DWARF.DW_OP_breg0 && opcode <= DWARF.DW_OP_breg31
                (i,offset) = operands(T,opcode,opcodes,i,endianness)
                push!(s.stack,Int(getreg_func(opcode-DWARF.DW_OP_breg0)) + offset)
            elseif opcode == DWARF.DW_OP_fbreg
                (i,offset) = operands(T,opcode,opcodes,i,endianness)
                val = getreg_func(opcode)
                push!(s.stack,val + offset%UInt64)
            elseif opcode == DWARF.DW_OP_bregx
                (i,(val,offset)) = operands(T,opcode,opcodes,i,endianness)
                push!(s.stack,getreg_func(val) + offset)
            elseif opcode == DWARF.DW_OP_call_frame_cfa
                push!(s.stack,getreg_func(DWARF.DW_OP_call_frame_cfa))
                i -= 1
            elseif opcode == DWARF.DW_OP_nop
                #NOP
            elseif opcode == DWARF.DW_OP_stack_value
                # Handled explicitly elsewhere
            else
                return (i-1,false)
            end
            (i,true)
        end

        function evaluate_generic{T}(s::StateMachine{T},opcodes::Array{UInt8,1},getreg_func::Function,getword_func,addr_func,endianness::Symbol)
            i=1
            while i <= length(opcodes)
                (opcodes[i] == DWARF.DW_OP_stack_value) &&
                    return true
                i,res = evaluate_generic_instruction(s,opcodes,i,getreg_func,getword_func,addr_func,endianness)
                if !res
                    error("Unrecognized Opcode ",opcodes[i])
                end
            end
            return false
        end

        function op_name(opcode)
            if DWARF.DW_OP_lit1 <= opcode <= DWARF.DW_OP_lit31
                return string("DW_OP_lit",1+opcode-DWARF.DW_OP_lit1)
            elseif DWARF.DW_OP_reg0 <= opcode <= DWARF.DW_OP_reg31
                return string("DW_OP_reg",opcode-DWARF.DW_OP_reg0)
            elseif DWARF.DW_OP_breg0 <= opcode <= DWARF.DW_OP_breg31
                return string("DW_OP_breg",opcode-DWARF.DW_OP_breg0)
            else
                return DWARF.DW_OP[opcode]
            end
        end

        function print_expression(io::IO, addr_type, opcodes::Array{UInt8,1},endianness::Symbol)
            i = 1
            while i <= length(opcodes)
                opcode = opcodes[i]
                i += 1
                ops = operands(addr_type, opcode, opcodes, i, endianness)
                i = ops[1]
                print_with_color(:blue, io, op_name(opcode))
                print(io," ")
                if length(ops) > 1
                    operand = ops[2]
                    if !isa(operand,Tuple)
                        operand = (operand,)
                    end
                    print(io,join(map(repr,operand)," ")," ")
                end
            end
        end

        immutable RegisterLocation
            i::Int32
        end

        immutable MemoryLocation{T}
            i::T
        end

        function evaluate_simple_location{T}(s::StateMachine{T},opcodes::Array{UInt8,1},
                getreg_func::Function,getword_func,addr_func,endianness::Symbol)
            i=1
            opcode = opcodes[i]
            if opcode >= DWARF.DW_OP_reg0 && opcode <= DWARF.DW_OP_reg31
                return RegisterLocation(opcode-DWARF.DW_OP_reg0)
            elseif opcode == DWARF.DW_OP_regx
                (i,val) = DWARF.decode(opcodes,i+1,ULEB128{UInt})
                return RegisterLocation(val)
            else
                if evaluate_generic(s,opcodes,getreg_func,getword_func,addr_func,endianness)
                    return last(s.stack)
                else
                    return MemoryLocation{T}(last(s.stack))
                end
            end
        end
    end

    # The line table is encoded as a state machine program
    # operating on a register machine, whose register represent the
    # values the debugger needs to know about the current source location
    module LineTableSupport
        using StructIO

        import ..ULEB128, ..SLEB128, ..DWARF
        import Base: ==
        using ObjFileBase: printfield

        immutable HeaderStub{T}
            length::T
            version::UInt16
            header_length::T
            minimum_instruction_length::UInt8
            maximum_operations_per_instruction::UInt8
            default_is_stmt::UInt8
            line_base::Int8
            line_range::Int8
            opcode_base::UInt8
        end

        immutable FileEntry
            name::String
            dir_idx::UInt
            timestamp::UInt
            filelength::UInt
        end

        Base.isequal(x::FileEntry,y::FileEntry) =
            (x.name == y.name && x.dir_idx == y.dir_idx && x.timestamp == y.timestamp && x.filelength == y.filelength)
        ==(x::FileEntry, y::FileEntry) = isequal(x,y)

        function Base.read(io::IO,::Type{FileEntry})
            s = readstring(io)
            if endof(s) == 0
                return FileEntry(s,0,0,0)
            end
            return FileEntry(s,read(io,ULEB128{UInt}),read(io,ULEB128{UInt}),read(io,ULEB128{UInt}))
        end

        function unpack(io,::Type{HeaderStub})
            T = UInt32
            length = read(io,T)
            if length == 0xffffffff
                T = UInt64
                length = read(io,T)
            end
            version = read(io,UInt16)
            header_length = read(io,T)
            minimum_instruction_length = read(io,UInt8)
            maximum_operations_per_instruction = version >= 4 ? read(io,UInt8) : 1
            HeaderStub{T}(length,version,header_length,minimum_instruction_length,maximum_operations_per_instruction,
                read(io,UInt8),read(io,Int8),read(io,Int8),read(io,UInt8))
        end

        immutable Header{T}
            stub::HeaderStub{T}
            standard_opcode_lengths::Vector{UInt8}
            include_directories::Vector{String}
            file_names::Vector{FileEntry}
        end

        function readstring(io)
            ret = Array(UInt8,0)
            while true
                c = read(io,UInt8)
                c == 0 && break
                push!(ret,c)
            end
            String(ret)
        end

        function read_header(io)
            stub = unpack(io,HeaderStub)
            standard_opcode_lengths = Array(UInt8,max(0,stub.opcode_base-1))
            read!(io,standard_opcode_lengths)
            include_directories = String[]
            while true
                s = readstring(io)
                endof(s) == 0 && break
                push!(include_directories,s)
            end
            file_names = FileEntry[]
            while true
                f = read(io,FileEntry)
                endof(f.name) == 0 && break
                push!(file_names,f)
            end
            Header{header_type(stub)}(stub,standard_opcode_lengths,include_directories,file_names)
        end

        header_type{T}(h::Header{T}) = T
        header_type{T}(h::HeaderStub{T}) = T

        immutable RegisterState
            address::Int
            op_index::Int
            file::Int
            line::Int
            column::UInt
            is_stmt::Bool
            basic_block::Bool
            end_sequence::Bool
            prologue_end::Bool
            epilogue_begin::Bool
            isa::Int
            discriminator::Int

            # Initial Register state as defined by DWARF standard
            function RegisterState(default_is_stmt::Bool)
                new(Int(0),Int(0),Int(1),Int(1),Int(0),default_is_stmt,
                    false,false,false,false,Int(0),Int(0))
            end

            function RegisterState(x::RegisterState;
                    address = x.address, op_index = x.op_index, file = x.file, line = x.line,
                    column = x.column, is_stmt = x.is_stmt, basic_block = x.basic_block,
                    end_sequence = x.end_sequence, prologue_end = x.prologue_end,
                    epilogue_begin = x.epilogue_begin, isa = x.isa, discriminator = x.discriminator)
                new(address,op_index,file,line,column,is_stmt,basic_block,end_sequence,prologue_end,
                    epilogue_begin,isa,discriminator)
            end

            RegisterState(address,op_index,file,line,column,is_stmt,basic_block,end_sequence,
                    prologue_end,epilogue_begin,isa,discriminator) =
                new(address,op_index,file,line,column,is_stmt,basic_block,end_sequence,
                    prologue_end,epilogue_begin,isa,discriminator)
        end

        Base.isequal(x::RegisterState,y::RegisterState) = (
            x.address == y.address &&
            x.op_index == y.op_index &&
            x.file == y.file &&
            x.line == y.line &&
            x.column == y.column &&
            x.is_stmt == y.is_stmt &&
            x.basic_block == y.basic_block &&
            x.end_sequence == y.end_sequence &&
            x.prologue_end == y.prologue_end &&
            x.epilogue_begin == y.epilogue_begin &&
            x.isa == y.isa &&
            x.discriminator == y.discriminator)
        ==(x::RegisterState, y::RegisterState) = isequal(x,y)

        type StateMachine
            header::Header
            state::RegisterState
        end

        function pc_adv!(m::StateMachine,op_adv)
            m.state = RegisterState(m.state,
                address = m.state.address + m.header.stub.minimum_instruction_length *
                    div(m.state.op_index + op_adv,m.header.stub.maximum_operations_per_instruction),
                op_index = m.state.op_index +
                    mod(m.state.op_index + op_adv,m.header.stub.maximum_operations_per_instruction))
        end

        function pcl_adv!(m::StateMachine,op; advance_line = true)
            adj_opc = op - m.header.stub.opcode_base
            pc_adv!(m,div(adj_opc,m.header.stub.line_range))
            if advance_line
                m.state = RegisterState(m.state,
                    line = m.state.line + m.header.stub.line_base + mod(adj_opc,m.header.stub.line_range))
            end
        end

        const DW_LNE_special_adv = 0x80
        const DW_LNS_OFF = 0x20

        function decode_next(io, header)
            op = read(io,UInt8)
            if op == 0
                # Extended opcode
                pos = position(io)
                size::UInt = read(io,ULEB128{UInt})
                ex_op = read(io,UInt8)
                if ex_op == DWARF.DW_LNE_end_sequence
                    position(io) > pos+size+1 && error("Malformed extended instruction")
                    return (DW_LNS_OFF+ex_op,())
                elseif ex_op == DWARF.DW_LNE_set_address
                    addrsize = pos+size - position(io) + 1
                    if addrsize == 4
                        return (DW_LNS_OFF+ex_op,(read(io,UInt32),))
                    elseif addrsize == 8
                        return (DW_LNS_OFF+ex_op,(read(io,UInt64),))
                    else
                        error("Unsupported target address size $addrsize")
                    end
                elseif ex_op == DWARF.DW_LNE_define_file
                    return (DW_LNS_OFF+ex_op,(read(io,FileEntry),))
                elseif ex_op == DWARF.DW_LNE_set_discriminator
                    return (DW_LNS_OFF+ex_op,(read(io,ULEB128),))
                else
                    error("Unrecognized extended opcode $ex_op")
                end
                position(io) > pos+size+1 && error("Malformed extended instruction (op=$ex_op, pos=$pos, size=$size, iopos=$(position(io)))")
            elseif op < header.stub.opcode_base
                # standard opcode
                if op == DWARF.DW_LNS_copy || op == DWARF.DW_LNS_negate_stmt || op == DWARF.DW_LNS_set_basic_block ||
                        op == DWARF.DW_LNS_const_add_pc || op == DWARF.DW_LNS_set_prologue_end ||
                        op == DWARF.DW_LNS_set_epilogue_begin
                    if header.standard_opcode_lengths[op] != 0
                        error("Malformed Instruction")
                    end
                    return (op,())
                elseif op == DWARF.DW_LNS_advance_pc || op == DWARF.DW_LNS_set_column || op == DWARF.DW_LNS_set_file ||
                        op == DWARF.DW_LNS_set_isa
                    if header.standard_opcode_lengths[op] != 1
                        error("Malformed Instruction")
                    end
                    return (op,(read(io,ULEB128),))
                elseif op == DWARF.DW_LNS_advance_line
                    if header.standard_opcode_lengths[op] != 1
                        error("Malformed Instruction")
                    end
                    return (op,(read(io,SLEB128),))
                elseif op == DWARF.DW_LNS_fixed_advance_pc
                    if header.standard_opcode_lengths[op] != 1
                        error("Malformed Instruction")
                    end
                    return (op,(read(io,UInt16),))
                else
                    error("Unknown opcode")
                end
            else
                # special opcode
                # decode the opcode
                return (DW_LNE_special_adv,(op,))
            end
        end

        function step(io,m::StateMachine)
            opcode, operands = decode_next(io, m.header)
            if opcode == DW_LNS_OFF+DWARF.DW_LNE_end_sequence
                m.state = RegisterState(m.state,end_sequence = true)
                ret = (true,m.state)
                m.state = RegisterState(m.header.stub.default_is_stmt > 0)
                return ret
            elseif opcode == DW_LNS_OFF+DWARF.DW_LNE_set_address
                m.state = RegisterState(m.state,address = operands[1])
            elseif opcode == DW_LNS_OFF+DWARF.DW_LNE_define_file
                push!(m.file_names,operands[1])
            elseif opcode == DW_LNS_OFF+DWARF.DW_LNE_set_discriminator
                m.state = RegisterState(m.state,discriminator = operands[1])
            elseif opcode == DWARF.DW_LNS_copy
                ret = (true,m.state)
                m.state = RegisterState(m.state,
                    discriminator = 0,
                    basic_block = false,
                    prologue_end = false,
                    epilogue_begin = false)
                return ret
            elseif opcode == DWARF.DW_LNS_advance_pc
                pc_adv!(m,UInt(operands[1]))
            elseif opcode == DWARF.DW_LNS_advance_line
                m.state = RegisterState(m.state,line = m.state.line + Int(operands[1]))
            elseif opcode == DWARF.DW_LNS_set_file
                m.state = RegisterState(m.state,file = operands[1])
            elseif opcode == DWARF.DW_LNS_set_column
                m.state = RegisterState(m.state,column = operands[1])
            elseif opcode == DWARF.DW_LNS_negate_stmt
                m.state = RegisterState(m.state,is_stmt = !m.state.is_stmt)
            elseif opcode == DWARF.DW_LNS_set_basic_block
                m.state = RegisterState(m.state,basic_block = true)
            elseif opcode == DWARF.DW_LNS_const_add_pc
                pcl_adv!(m,255; advance_line = false)
            elseif opcode == DWARF.DW_LNS_fixed_advance_pc
                m.state = RegisterState(m.state,
                    address = m.state.address + operands[1],
                    op_index = 0)
            elseif opcode == DWARF.DW_LNS_set_prologue_end
                if m.header.standard_opcode_lengths[DWARF.DW_LNS_set_prologue_end] != 0
                    error("Malformed Instruction")
                end
                m.state = RegisterState(m.state, prologue_end = false)
            elseif opcode == DWARF.DW_LNS_set_epilogue_begin
                m.state = RegisterState(m.state, epilogue_begin = false)
            elseif opcode == DWARF.DW_LNS_set_isa
                m.state = RegisterState(m.state, isa = operands[1])
            elseif opcode == DW_LNE_special_adv
                pcl_adv!(m,operands[1])
                ret = (true, m.state)
                m.state = RegisterState(m.state, basic_block = false, prologue_end = false, epilogue_begin = false,
                    discriminator = 0)
                return ret
            else
                error("Unhandled opcode")
            end
            return (false,m.state)
        end

        function state_step(io,m::StateMachine)
            stop = false
            state = m.state
            while !stop
                (stop, state) = step(io,m)
            end
            return state
        end

        # Iterate over the line table
        immutable LineTable
            io::IO
            header::Header
            start::Int
            LineTable(io,header,start) = new(io,header,start)
            LineTable(io) = (pos=position(io);new(io,read_header(io),pos))
        end

        function dump_program(out::IO, x::LineTable)
            seek(x.io,x.start + x.header.stub.header_length + 2*sizeof(header_type(x.header)) + sizeof(UInt16) +
                (header_type(x.header) == UInt64 ? sizeof(UInt32) : 0))
            while position(x.io) <= (x.start + x.header.stub.length)
                opcode, operands = decode_next(x.io, x.header)
                print_with_color(:blue, out, opcode == DW_LNE_special_adv ? "DW_LNE_special_adv" :
                    opcode > DW_LNS_OFF ? DWARF.DW_LNE[opcode-DW_LNS_OFF] : DWARF.DW_LNS[opcode])
                for operand in operands
                    print(out, ' ')
                    if opcode == DWARF.DW_LNE_set_address+DW_LNS_OFF
                        print(out, "0x", hex(operand))
                    elseif isa(operand, ULEB128) || isa(operand, SLEB128)
                        print(out, typeof(operand).parameters[1](operand))
                    else
                        print(out, operand)
                    end
                end
                println(out)
            end
        end

        function Base.showcompact(io::IO, state::RegisterState)
            printfield(io, string("0x",hex(state.address)), 18); print(io,' ')
            printfield(io, state.line, 6); print(io,' ')
            printfield(io, state.column, 6); print(io,' ')
            printfield(io, state.file, 6); print(io,' ')
            printfield(io, state.isa, 3); print(io,' ')
            printfield(io, state.discriminator, 13, align = :right); print(io,' ')
            state.is_stmt && print(io, "is_stmt ")
            state.end_sequence && print(io, "end_sequence ")
            state.prologue_end && print(io, "prologue_end ")
            state.epilogue_begin && print(io, "epilogue_begin ")
            println(io)
        end

        # Follow llvm-dwarfdump format
        function dump_table(io::IO, x::LineTable)
            println(io, "Address            Line   Column File   ISA Discriminator Flags")
            println(io, "------------------ ------ ------ ------ --- ------------- -------------")
            for state in x
                showcompact(io, state)
            end
        end

        import Base: start, next, done

        function start(x::LineTable)
            # The header length is the number of bytes after the header length
            # The fields from the start are length (header_type), version (UInt16) and header_length(header_type)
            # plus an additional UInt32 if the header_type is UInt64
            seek(x.io,x.start + x.header.stub.header_length + 2*sizeof(header_type(x.header)) + sizeof(UInt16) +
                (header_type(x.header) == UInt64 ? sizeof(UInt32) : 0))
            StateMachine(x.header,RegisterState(x.header.stub.default_is_stmt > 0))
        end
        next(x::LineTable,m::StateMachine) = (state_step(x.io,m),m)
        done(x::LineTable,m::StateMachine) = position(x.io) > (x.start + x.header.stub.length)
    end

    import .LineTableSupport.LineTable
    export LineTable

    zero(::Type{AttributeSpecification}) = AttributeSpecification(0,0)

    const tag_color = :blue

    function tag_name(tag)
        get(DW_TAG, tag, "Unknown tag $tag")
    end

    immutable ARTableEntry{S,T}
        segment::S
        address::T
        length::T
    end

    immutable LocationListEntry{T}
        first::T
        last::T
        data::Array{UInt8,1}
    end

    immutable LocationList{T}
        entries::Vector{LocationListEntry{T}}
    end

    function show{T}(io::IO, x::LocationListEntry{T})
        print(io,repr(x.first)," - ", repr(x.last), ": ")
        Expressions.print_expression(io, T, x.data,:NativeEndian)
        println(io)
    end

    function show(io::IO, x::LocationList)
        println(io,"Location List:")
        for e in x.entries
            show(io, e)
        end
    end

    isequal{S,T}(a::ARTableEntry{S,T},b::ARTableEntry{S,T}) = (a.segment==b.segment)&&(a.address==b.address)&&(a.length==b.length)
    =={S,T}(a::ARTableEntry{S,T},b::ARTableEntry{S,T}) = isequal(a,b)

    zero{S,T}(::Type{ARTableEntry{S,T}}) = ARTableEntry{S,T}(zero(S),zero(T),zero(T))

    const DEBUG_SECTIONS = ObjFileBase.DEBUG_SECTIONS

    @struct immutable InitialLength
        val::UInt32
    end

    for (t,t32,t64) in ((:(DWARF.DWARFCUHeader),:(DWARF32.CUHeader),:(DWARF32.CUHeader)),
                      (:(DWARF.DWARFTUHeader),:(DWARF32.TUHeader),:(DWARF32.TUHeader)),
                      (:(DWARF.DWARFARHeader),:(DWARF32.ARHeader),:(DWARF32.ARHeader)),
                      (:(DWARF.DWARFPUBHeader),:(DWARF32.PUBHeader),:(DWARF32.PUBHeader)))
        @eval function read(io::IO,::Type{($t)},endianness::Symbol)
                # XXX: Is there a better way to do this?
                l = unpack(io,InitialLength,endianness)
                if l.val<0xfffffff0 #Is a 32bit DWARF record
                    skip(io,-sizeof(InitialLength)) #Try to rewind
                    return unpack(io,($t32),endianness)
                elseif l.val == 0xffffffff
                    return unpack(io,($t64),endianness)
                else
                    error("Unkown Compilation Unit Header Type")
                end
            end
        @eval read(io::IO,::Type{$t}) = read(io,$t,:NativeEndian)
    end

    immutable Zero
    end

    bswap(::Zero) = Zero()
    read(io::IO,::Type{Zero}) = Zero()
    zero(::Type{Zero}) = Zero()

    function size_to_inttype(size)
        if size == 0
            return Zero
        elseif size == 1
            return UInt8
        elseif size == 2
            return UInt16
        elseif size == 4
            return UInt32
        elseif size == 8
            return UInt64
        elseif size == 16
            return UInt128
        else
            error("Unsupported size unit $size")
        end
    end

    immutable ARTableSet{S,T}
        header::DWARFARHeader
        entries::Array{ARTableEntry{S,T},1}
    end

    type ARTable
        sets::Array{ARTableSet,1}
    end

    # aranges tables
    function read(io::IO,::Type{ARTableSet},endianness::Symbol)
        header = read(io,DWARFARHeader,endianness)
        t = ARTableEntry{size_to_inttype(header.segment_size),size_to_inttype(header.address_size)}
        table = ARTableSet(header,Array(t,0))
        entry_size = sizeof(t)
        while true
            skip(io,position(io)%entry_size) # Align to a boundary that is a multiple of the entry size
            r = unpack(io,t,dummy_dict,align_packed,endianness)
            if r==zero(t)
                break
            end
            push!(table.entries,r)
        end
        table
    end
    read(io::IO,::Type{ARTableSet}) = read(io,ARTableSet,:NativeEndian)

    # pubnames/pubtypes
    type PUBTable
        sets::Array{PUBTableSet,1}
    end

    function read{T<:PUBTableEntry}(io::IO,::Type{T},endianness::Symbol)
        offset = read(io,T.types[1])
        StructIO.needs_bswap(endianness) && (offset = bswap(offset))
        if offset != 0
            name = strip(readuntil(io,'\0'),'\0')
        else
            name = ""
        end
        T(offset,name)
    end

    function read(io::IO,::Type{PUBTableSet},endianness::Symbol)
        header = read(io,DWARFPUBHeader,endianness)
        if typeof(header) == DWARF32.PUBHeader
            t=DWARF32.PUBTableSet(header,Array(DWARF32.PUBTableEntry,0))
        elseif typeof(header) == DWARF64.PUBHeader
            t=DWARF64.PUBTableSet(header,Array(DWARF64.PUBTableEntry,0))
        end
        while true
            entry = read(io,eltype(t.entries),endianness)
            if entry.offset == 0
                break
            end
            push!(t.entries,entry)
        end
        t
    end
    read(io::IO,::Type{PUBTableSet}) = read(io,PUBTableSet,:NativeEndian)

    function read(io::IO,::Type{AbbrevTableEntry},endianness::Symbol)
        code = read(io,ULEB128)
        if code != 0
            tag = read(io,ULEB128)
            has_children = read(io,UInt8)
            ret = AbbrevTableEntry(code,tag,has_children,Array(AttributeSpecification,0))
            while true
                e = read(io,AttributeSpecification,endianness)
                if(e == zero(AttributeSpecification))
                    break
                end
                push!(ret.attributes,e)
            end
            return ret
        else
            return zero(AbbrevTableEntry)
        end
    end

    function read(io::IO,::Type{AbbrevTableSet},endianness::Symbol)
        ret = AbbrevTableSet(Array(AbbrevTableEntry,0))
        while true
            e = read(io,AbbrevTableEntry,endianness)
            if e.code == 0
                break
            end
            push!(ret.entries,e)
        end
        ret
    end

    function read{T}(io::IO,::Type{LocationList{T}})
        ret = LocationListEntry{T}[]
        while true
            pos = position(io)
            start = read(io, T)
            last = read(io, T)
            if start == typemax(T)
                error("Base address selection not implemented")
            elseif start == 0 && last == 0
                break
            else
                length = read(io, UInt16)
                push!(ret,LocationListEntry{T}(start,last,read(io,length)))
            end
        end
        LocationList{T}(ret)
    end

    include("navigate.jl")
    include("cfi.jl")
    include("precompile.jl")
end #module


text:
using DiffTests
using Test

n = rand()
x, y = rand(5, 5), rand(26)
A, B = rand(5, 5), rand(5, 5)

# f returns Number

for f in DiffTests.NUMBER_TO_NUMBER_FUNCS
    @test isa(f(n), Number)
end

for f in DiffTests.VECTOR_TO_NUMBER_FUNCS
    @test isa(f(y), Number)
end

for f in DiffTests.MATRIX_TO_NUMBER_FUNCS
    @test isa(f(x), Number)
end

for f in DiffTests.TERNARY_MATRIX_TO_NUMBER_FUNCS
    @test isa(f(A, B, x), Number)
end

# f returns Array

for f in DiffTests.NUMBER_TO_ARRAY_FUNCS
    @test isa(f(n), Array)
end

for f in DiffTests.ARRAY_TO_ARRAY_FUNCS
    @test isa(f(A), Array)
    @test isa(f(y), Array)
end

for f in DiffTests.MATRIX_TO_MATRIX_FUNCS
    @test isa(f(A), Array)
end

for f in DiffTests.BINARY_MATRIX_TO_MATRIX_FUNCS
    @test isa(f(A, B), Array)
end

# f! returns Nothing

for f! in DiffTests.INPLACE_ARRAY_TO_ARRAY_FUNCS
    @test isa(f!(y, x), Nothing)
end

for f! in DiffTests.INPLACE_NUMBER_TO_ARRAY_FUNCS
    @test isa(f!(y, n), Nothing)
end

ext: jl


text:
#https://en.wikipedia.org/wiki/Rotation_matrix#In_three_dimensions

function euler_rotation_matrix(θ)
	a = θ[1]
	b = θ[2]
	c = θ[3]

	[cos(a) * cos(b) (cos(a) * sin(b) * sin(c) - sin(a) * cos(c)) (cos(a) * sin(b) * cos(c) + sin(a) * sin(c));
	 sin(a) * cos(b) (sin(a) * sin(b) * sin(c) + cos(a) * cos(c)) (sin(a) * sin(b) * cos(c) - cos(a) * sin(c));
	 -sin(b) cos(b) * sin(c) cos(b) * cos(c)]
end


text:
module OpticalProperties

using MyPhysicalConstants, Interpolations, Parameters, QuadGK

# Functions
export convert_prop, skin_depth,
       permittivity, refractive_index,eps_si,eps_sin,
       resistivity

# Types
export OptProp, ElectricalProperties,
       Model, Polariton,Cst,
       Bruggeman, MaxwellGarnett,
       SimpleMixingPar,SimpleMixingPerp,
       MaterialFile, Sellmeier,
       ResistivityFile, MobilityModel,
       Cbn, Sic,
       Si_n_doped,Si_p_doped

# Constants
export Cu,SiO2,Si,SiN,Vacuum,
       Au,Au_latella,Al,
       W,Ti,TiW,TiW_v2,
       TiN,TiN_cethil,
       pSi_masetti,nSi_masetti,
       pSi_sze,nSi_sze,
       SiC,cBN,Si_cst

abstract type AbstractMaterial end
abstract type OptProp <: AbstractMaterial end
abstract type ElectricalProperties <: AbstractMaterial end

# Generic models
struct Model <: OptProp
    eps0   :: Float64
    wp     :: Float64
    w0     :: Float64
    gamma0 :: Float64
    gamma1 :: Float64
end
Model(eps0,wp,w0,gamma0) = Model(eps0,wp,w0,gamma0,0.0)

gamma1(mfp=1.0,a=0.0; vf=0.0) = a*vf/mfp
# Gold
Au(mfp=1.0,a=0.0) = Model(9.4,13584.25e12,0.0,109.96e12,gamma1(mfp,a;vf=vf_au))
Au_latella(mfp,a) = Model(1.0,1.37e16,0.0,5.32e13,gamma1(mfp,a;vf=vf_au))

#Aluminium
Al(mfp=1.0,a=0.0) = Model(1.0,2.24e16,0.0,1.22e14,gamma1(mfp,a ; vf=vf_ag))

"Drude model of Tungsten: fitted from data obtained from Ordal et al."
const W = Model(6.0, 6.2e15, 0.0, 2.6e14,0.0)

"Drude model of Titanium: fitted from data obtained from Ordal et al."
const Ti = Model(1.0, 3.82e15, 0.0, 719.6e11, 0.0)

"Drude model of Titanium Tungsten: obtained from Varpula et al."
const TiW    = Model(1.0,sqrt(7.75e5*1e15/8.85e-12), 0.0, 1.0e15, 0.0)
const TiW_v2 = Model(1.0,sqrt(5.0e5*1e15/8.85e-12), 0.0, 1.0e15, 0.0)

"Drude model of Titanium Nitride obtained from Varpula et al."
const TiN = Model(1.0,sqrt(7.57e5/1.5e-15/8.85e-12), 0.0, 1.0/(1.5e-15), 0.0)

"""
    gamma0_TiN_cethil(thickness)

    Drude model for TiN based on reflectance measurements performed at CETHIL 04/2021.
# Arguments
* `thickness`: layer thickness in meter
"""
function gamma0_TiN_cethil(thickness)
    th = thickness*1e9
    if th<=20
        gamma = 0.0087e15*th + 0.8939e15
    else
        gamma = 1.13e15
    end
end

"""
    TiN_cethil(thickness)

    Returns the Model object for a Drude model with the parameters obtained from reflectance measurements performed at CETHIL 04/2021.
# Arguments
* `thickness`: layer thickness in meter
"""
TiN_cethil(thickness) = Model(1.0, 6.93e15,0.0,gamma0_TiN_cethil(thickness), 0.0)

struct Polariton{T} <: OptProp
    eps_fin :: T
    w_lo    :: T # rad/s
    w_to    :: T # rad/s
    gamma   :: T # rad/s
end
SiC = Polariton(6.7,1.827e14,1.495e14,8.971e11)
cBN = Polariton(4.46,2.451e14,1.985e14,9.934e11)

struct LorentzTerm{T} <: OptProp
    f      :: T
    wp     :: T
    w0     :: T
    gamma0 :: T
end

function susceptibility(model :: LorentzTerm, w)
    @unpack f, wp, w0, gamma0 = model
    return f*wp^2/(w^2 - w0^2 + im*w*gamma0)
end

struct DrudeTerm{T} <: OptProp
    wp     :: T
    gamma0 :: T
    gamma1 :: T
end

DrudeTerm(wp, gamma0) = DrudeTerm(wp, gamma0, 0.0)

function susceptibility(model :: DrudeTerm, w)
    @unpack wp,gamma0, gamma1 = model
    return wp^2/(w*(w + im*(gamma0+gamma1)))
end

struct DrudeLorentzSerie{T} <: OptProp
    eps0    :: T
    drude   :: DrudeTerm{T}
    lorentz :: Vector{LorentzTerm{T}}
end

function permittivity(model :: DrudeLorentzSerie, w)
    @unpack eps0, drude, lorentz = model
    ϵ = eps0 - susceptibility(drude, w)
    for term in lorentz
        ϵ -= susceptibility(term, w)
    end
    return ϵ
end

struct Sellmeier{T} <: OptProp
    values :: Array{Tuple{T,T},1}
end

const Si  = Sellmeier([(10.6684293 , 0.301516485), (0.0030434748 , 1.13475115), (1.54133408 , 1104.0)])
const SiN = Sellmeier([(3.0249 , 0.1353406), (40314.0 , 1239.842)])

# Constant permittivity
struct Cst <: OptProp
    val :: Complex{Float64}
end

Cst()  = Cst(1.0+im*0.0)
const Vacuum = Cst()
const Si_cst = Cst(11.7 + im*0.0)

# Types for effective medium models
struct Bruggeman{T, U , V} <: OptProp
   phase1  :: T
   phase2  :: U
   volfrac :: V
end

struct MaxwellGarnett{T, U , V} <: OptProp
   particle :: T
   matrix   :: U
   volfrac  :: V
end

struct SimpleMixingPar{T,U,V} <: OptProp
   phase1  :: T
   phase2  :: U
   volfrac :: V
end

struct SimpleMixingPerp{T,U,V} <: OptProp
   phase1  :: T
   phase2  :: U
   volfrac :: V
end

# From file containing measurements
struct MaterialFile{T, U, V} <: OptProp
    re :: V
    im :: V
end

# Outer constructor to dispatch on options
MaterialFile(T, U, re :: V, im :: V) where V = MaterialFile{T,U,V}(re,im)

function Base.show(io:: IO, mf :: MaterialFile{T,U}) where {T,U}
    println(io, "property: $T")
    println(io, "parameter: $U")
end

struct ResistivityFile{T} <: ElectricalProperties
    rho :: T
end

struct MobilityModel{T} <: ElectricalProperties
    μ1   :: T
    μ2   :: T
    μmax :: T
    cr   :: T
    cs   :: T
    α    :: T
    β    :: T
    pc   :: T
    charge :: T
end

const nSi_masetti = MobilityModel(68.5,56.1,1414.0,9.2e16,3.41e20,0.711,1.98,0.0,0.0)
const pSi_masetti = MobilityModel(44.9,29.0,470.5,2.23e17,6.10e20,0.719,2.0,9.23e16,1.0)

# Create interpolations
include("interpolation_objects.jl")
include("SiN.jl")

#To be deprecated
include("singleton_types.jl")

"""
    mobility(material,N)

Compute the electric mobility for a material in cm^2 V^{-1} s^{-1}.

# Arguments
* `material :: MobilityModel` : the material for which you need the resisitivity
* `N :: Real` : electron or hole concentration in cm^{-3}

# Example
```julia
julia> mobility(nSi,1e15)
4.593579841190333
```
"""

function mobility(material :: ElectricalProperties, N)
    @unpack μ1, μ2, μmax, cr, cs, α, β, pc, charge = material
    return μ1*exp(-charge*pc/N) + (μmax - (1.0 - charge)*μ1)/(1 + (N/cr)^α) - μ2/(1 + (cs/N)^β)
end

"""
    resistivity(material,N)

Compute the electric resistivity for a material in Ohm cm.

# Arguments
* `material :: ResistivityFile` : the material for which you need the resisitivity
* `N :: Real` : impurity concentration in cm^{-3}

# Example
```julia
julia> resistivity(nSi,1e15)
4.593579841190333
```
"""
resistivity(material :: ResistivityFile, N) = material.rho(N)
resistivity(material :: ElectricalProperties, N) = 1/(N*electron*mobility(material,N))

"""
Instances of doped Silicon. N varies between 3e19 and 5e20 cm^{-3}
rho is the dc resistivity
"""

Si_n_doped(material::ElectricalProperties, N) = Model(11.7,sqrt(N*1e6*electron^2/(0.27*m0)/epsilon0),0.0,N*1e4*electron^2*resistivity(material,N)/(0.27*m0),0.0)
Si_p_doped(material::ElectricalProperties,N) = Model(11.7,sqrt(N*1e6*electron^2/(0.34*m0)/epsilon0),0.0,N*1e4*electron^2*resistivity(material,N)/(0.34*m0),0.0)

"""
    permittivity(material,w)

Compute the dielectric permittivity for a material.

# Arguments
* `material :: OptProp` : is a type for the material you want to pass
* `w :: Real` : pulsation in rad s^(-1)

# Example
```julia
julia> permittivity(SiC,1e13)
6.805820438080644 + 0.002847538251166107im
```
"""

function permittivity(material::Model,w) :: Complex{Float64}
     eps0   = material.eps0
     wp     = material.wp
     w0     = material.w0
     gamma0 = material.gamma0
     gamma1 = material.gamma1
    return  eps0 + wp^2/(w0*w0 - w*w - im*w*(gamma0+gamma1))
end

function permittivity(material::Polariton,w) :: Complex{Float64}
    eps_fin = material.eps_fin
    w_lo    = material.w_lo
    w_to    = material.w_to
    gamma   = material.gamma
    return  eps_fin*(w^2-w_lo^2 + im*gamma*w)/(w^2-w_to^2 + im*gamma*w)
end

function permittivity(material::Cst,w) :: Complex{Float64}
    return material.val
end

function permittivity(material::Bruggeman,w) :: Complex{Float64}
    eps1 = permittivity(material.phase1,w)
    eps2 = permittivity(material.phase2,w)
    p    = sqrt(eps1)/sqrt(eps2)
    b    = ((3.0*material.volfrac-1.0)*(1.0/p - p)+p)/4.0
    z    = b+sqrt(b*b + 0.5)
    return transpose(z*sqrt(eps2)*sqrt(eps1))
end

function permittivity(material::MaxwellGarnett,w) :: Complex{Float64}
    epsm = permittivity(material.matrix,w)
    epsp = permittivity(material.particle,w)
    f    = material.volfrac
    num  = 3.0*epsm + (1.0 + 2.0*f)*(epsp - epsm)
    den  = 3.0*epsm + (1.0 - f)*(epsp - epsm)
    return epsm*num/den
end

function permittivity(material::SimpleMixingPar,w)
    eps1 = permittivity(material.phase1,w)
    eps2 = permittivity(material.phase2,w)
    f    = material.volfrac
    return f*eps1 + (1-f)*eps2
end

function permittivity(material::SimpleMixingPerp,w)
    eps1 = permittivity(material.phase1,w)
    eps2 = permittivity(material.phase2,w)
    f    = material.volfrac
    return 1/(f/eps1 + (1-f)/eps2)
end

function permittivity(material :: MaterialFile{:refr_ind, :wavelength}, w)
    lamb = 2.0*pi*c0/w*1e6 # wavelength in microns
    return (material.re(lamb)^2 - material.im(lamb)^2) + im*2*material.im(lamb)*material.re(lamb)
end

function permittivity(material :: MaterialFile{:permittivity, :frequency}, w)
    return material.re(w) + im*material.im(w)
end

function permittivity(material :: Sellmeier, w)
    lamb   = 2.0*pi*c0/w*1e6
    eps    = 1.0 + im*0.0
    for val in material.values
        eps += val[1]/(1-(val[2]/lamb)^2)
    end
    return eps
end

# Refractive index
refractive_index(material :: OptProp, w) = sqrt(permittivity(material,w))

# Skin depth
skin_depth(material :: OptProp , w ) = c0/imag(refractive_index(material,w))/w

end # module


text:
#
# types.jl -
#
# Defines camera models and methods specific to given camera models.  In
# particular, the addresses for registers required to be able to stream images
# from a CoaXPress camera.
#
#------------------------------------------------------------------------------
#
# This file is part of the `Phoenix.jl` package which is licensed under the MIT
# "Expat" License.
#
# Copyright (C) 2017-2019, Éric Thiébaut (https://github.com/emmt/Phoenix.jl).
# Copyright (C) 2016, Éric Thiébaut & Jonathan Léger.
#

# Generic camera model.
struct GenericCameraModel <: CameraModel end

function _openhook(cam::Camera{GenericCameraModel})
    if false
        # Set default configuration (as in `default.pcf`).
        cam[PHX_BOARD_VARIANT]      = PHX_DIGITAL
        cam[PHX_CAM_TYPE]           = PHX_CAM_AREASCAN_ROI
        cam[PHX_CAM_SRC_DEPTH]      = 8
        cam[PHX_CAM_SRC_COL]        = PHX_CAM_SRC_MONO
        cam[PHX_CAM_ACTIVE_XOFFSET] = 0
        cam[PHX_CAM_ACTIVE_YOFFSET] = 0
        cam[PHX_CAM_ACTIVE_XLENGTH] = 640
        cam[PHX_CAM_ACTIVE_YLENGTH] = 480
        cam[PHX_CAM_HTAP_DIR]       = PHX_CAM_HTAP_LEFT
        cam[PHX_CAM_HTAP_TYPE]      = PHX_CAM_HTAP_LINEAR
        cam[PHX_CAM_HTAP_NUM]       = 1
        cam[PHX_CAM_VTAP_DIR]       = PHX_CAM_VTAP_TOP
        cam[PHX_CAM_VTAP_TYPE]      = PHX_CAM_VTAP_LINEAR
        cam[PHX_CAM_VTAP_NUM]       = 1
        cam[PHX_COMMS_DATA]         = PHX_COMMS_DATA_8
        cam[PHX_COMMS_STOP]         = PHX_COMMS_STOP_1
        cam[PHX_COMMS_PARITY]       = PHX_COMMS_PARITY_NONE
        cam[PHX_COMMS_SPEED]        = 9600
        cam[PHX_COMMS_FLOW]         = PHX_COMMS_FLOW_NONE
    end
    return PHX_OK
end

struct ActiveSiliconTestModel <: CameraModel end

_starthook(cam::Camera{ActiveSiliconTestModel}) =
    exec(cam, RegisterCommand{UInt32}(0x6000, 0))

_stophook(cam::Camera{ActiveSiliconTestModel}) =
    exec(cam, RegisterCommand{UInt32}(0x6000, 1))

struct AdimecOpal2000cModel <: CameraModel end

_starthook(cam::Camera{AdimecOpal2000cModel}) =
    exec(cam, RegisterCommand{UInt32}(0x8204, 1))

_stophook(cam::Camera{AdimecOpal2000cModel}) =
    exec(cam, RegisterCommand{UInt32}(0x8208, 1))

struct AdimecOpal2000mModel <: CameraModel end

_starthook(cam::Camera{AdimecOpal2000mModel}) =
    exec(cam, RegisterCommand{UInt32}(0x8204, 1))

_stophook(cam::Camera{AdimecOpal2000mModel}) =
    exec(cam, RegisterCommand{UInt32}(0x8208, 1))

struct AdimecQuartz4A180Model <: CameraModel end

_starthook(cam::Camera{AdimecQuartz4A180Model}) =
    exec(cam, RegisterCommand{UInt32}(0x8204, 1))

_stophook(cam::Camera{AdimecQuartz4A180Model}) =
    exec(cam, RegisterCommand{UInt32}(0x8208, 1))

struct E2vEliixa16kColorModel <: CameraModel end

_starthook(cam::Camera{E2vEliixa16kColorModel}) =
    exec(cam, RegisterCommand{UInt32}(0x700C, 0))

_stophook(cam::Camera{E2vEliixa16kColorModel}) =
    exec(cam, RegisterCommand{UInt32}(0x7010, 0))

struct E2vEliixa16kMonoModel <: CameraModel end

_starthook(cam::Camera{E2vEliixa16kMonoModel}) =
    exec(cam, RegisterCommand{UInt32}(0x700C, 0))

_stophook(cam::Camera{E2vEliixa16kMonoModel}) =
    exec(cam, RegisterCommand{UInt32}(0x7010, 0))

struct ImperxModel <: CameraModel end

_starthook(cam::Camera{ImperxModel}) =
    exec(cam, RegisterCommand{UInt32}(0x10000020, 1))

_stophook(cam::Camera{ImperxModel}) =
    exec(cam, RegisterCommand{UInt32}(0x10000024, 1))

struct IoIndustries4M140Model <: CameraModel end

_starthook(cam::Camera{IoIndustries4M140Model}) =
    exec(cam, RegisterCommand{UInt32}(0x10000020, 1))

_stophook(cam::Camera{IoIndustries4M140Model}) =
    exec(cam, RegisterCommand{UInt32}(0x10000024, 1))

struct ISVI_C25M_CXP_Model <: CameraModel end

_starthook(cam::Camera{ISVI_C25M_CXP_Model}) =
    exec(cam, RegisterCommand{UInt32}(0x6050, 1))

_stophook(cam::Camera{ISVI_C25M_CXP_Model}) =
    exec(cam, RegisterCommand{UInt32}(0x6050, 0))

struct ISVI_M12M_CXP_Model <: CameraModel end

_starthook(cam::Camera{ISVI_M12M_CXP_Model}) =
    exec(cam, RegisterCommand{UInt32}(0x6050, 1))

_stophook(cam::Camera{ISVI_M12M_CXP_Model}) =
    exec(cam, RegisterCommand{UInt32}(0x6050, 0))

struct ISVI_M25M_CXP_Model <: CameraModel end

_starthook(cam::Camera{ISVI_M25M_CXP_Model}) =
    exec(cam, RegisterCommand{UInt32}(0x6050, 1))

_stophook(cam::Camera{ISVI_M25M_CXP_Model}) =
    exec(cam, RegisterCommand{UInt32}(0x6050, 0))

struct JAI_SP_5000M_Model <: CameraModel end

_starthook(cam::Camera{JAI_SP_5000M_Model}) =
    exec(cam, RegisterCommand{UInt32}(0x200B0, 1))

_stophook(cam::Camera{JAI_SP_5000M_Model}) =
    exec(cam, RegisterCommand{UInt32}(0x200B4, 1))

# struct MikrotronMC408xModel <: CameraModel end
#
# _starthook(cam::Camera{MikrotronMC408xModel}) =
#     exec(cam, RegisterCommand{UInt32}(0x8204, 1))
#
# _stophook(cam::Camera{MikrotronMC408xModel}) =
#     exec(cam, RegisterCommand{UInt32}(0x8208, 1))
include("MikrotronMC408x.jl")
import .MikrotronMC408x: MikrotronMC408xModel
using .MikrotronMC408x
MikrotronMC408x.@exportmethods

struct OptronisCL4000Model <: CameraModel end

_starthook(cam::Camera{OptronisCL4000Model}) =
    exec(cam, RegisterCommand{UInt32}(0x601C, 1))

_stophook(cam::Camera{OptronisCL4000Model}) =
    exec(cam, RegisterCommand{UInt32}(0x601C, 0))

struct OptronisCL25000Model <: CameraModel end

_starthook(cam::Camera{OptronisCL25000Model}) =
    exec(cam, RegisterCommand{UInt32}(0x601C, 1))

_stophook(cam::Camera{OptronisCL25000Model}) =
    exec(cam, RegisterCommand{UInt32}(0x601C, 0))


text:
#module optical_flow

#export get_flows, compute_warping_operator

using PyCall
using Suppressor
using LinearOperators

is_initiated = false

function __init__()
    is_initiated = true
    py"""
    # reference: https://github.com/pathak22/pyflow.git
    import numpy as np
    import pyflow

    def py_flow(img1,img2,alpha=0.012, ratio=0.75, minWidth=20, nOuterFPIterations=7, nInnerFPIterations=1, nSORIterations=30, colType=1):
        img1 = np.array(img1)
        img2 = np.array(img2)
        height, width = img1.shape

        im1 = np.array(img1).reshape(height,width,1).copy(order='C')
        im2 = np.array(img2).reshape(height,width,1).copy(order='C')

        u, v, im2W = pyflow.coarse2fine_flow(
            im1, im2, alpha, ratio, minWidth, nOuterFPIterations, nInnerFPIterations,
            nSORIterations, colType)

        return u, v, im2W
    """
end

function get_flows(u::Array{T,3})::Array{T,4} where {T<:AbstractFloat}
    @suppress begin
        #HACK!
        if !is_initiated
            __init__()
        end
        height, width, frames = size(u)
        result = zeros(height, width, 2, frames)
        for i in 1:frames-1
            flow_x, flow_y, im2Warped = py"py_flow"(u[:,:,i],u[:,:,i+1])
            result[:,:,1,i] = flow_x
            result[:,:,2,i] = flow_y
        end
        return result
    end
end

using SparseArrays
using LinearOperators

function sub2ind(array_shape, rows, cols)
    out = (cols.-1).*array_shape[1] .+ rows
    return convert.(Int64, out)
end

function compute_warping_operator(flow)
    # ref: https://github.com/HendrikMuenster/flexBox/blob/79abc7285703911cca2653434cca5bcefc79c722/operators/generateMatrix/warpingOperator.m
    """
    u: gray image
    v: flow [H, W, 2]
    """
    H, W, dims = size(flow)
    shape = (H,W)
    XX, YY = repeat(collect(1:W)', H, 1), repeat(collect(1:H), 1, W)

    targetPoint2Mat = XX + flow[:,:,1]
    targetPoint1Mat = YY + flow[:,:,2]

    x1 = floor.(targetPoint1Mat).-1
    x2 = x1 .+ 1
    x3 = x1 .+ 2
    x4 = x1 .+ 3

    y1 = floor.(targetPoint2Mat).-1
    y2 = y1 .+ 1
    y3 = y1 .+ 2
    y4 = y1 .+ 3

    v2 = targetPoint1Mat .- x2
    v1 = targetPoint2Mat .- y2

    indicator = ones(H,W)
    indicator = (x1 .> 0).*indicator
    indicator = (x4 .<= W).*indicator
    indicator = (y1 .> 0).*indicator
    indicator = (y4 .<= H).*indicator

    indicator = indicator .> 0
    v1 = v1[indicator]
    v2 = v2[indicator]

    idxList = sub2ind(shape, YY, XX)
    idxList_ = collect(Iterators.flatten(idxList))

    #find the indices
    list_J  = vcat(
        idxList_[sub2ind(shape, x2[indicator], y2[indicator])],
        idxList_[sub2ind(shape, x3[indicator], y2[indicator])],
        idxList_[sub2ind(shape, x2[indicator], y3[indicator])],
        idxList_[sub2ind(shape, x3[indicator], y3[indicator])])

    # list_J  = vcat(
    #     idxList_[sub2ind(shape, x1[indicator], y1[indicator])],
    #     idxList_[sub2ind(shape, x2[indicator], y1[indicator])],
    #     idxList_[sub2ind(shape, x3[indicator], y1[indicator])],
    #     idxList_[sub2ind(shape, x4[indicator], y1[indicator])],
    #     idxList_[sub2ind(shape, x1[indicator], y2[indicator])],
    #     idxList_[sub2ind(shape, x2[indicator], y2[indicator])],
    #     idxList_[sub2ind(shape, x3[indicator], y2[indicator])],
    #     idxList_[sub2ind(shape, x4[indicator], y2[indicator])],
    #     idxList_[sub2ind(shape, x1[indicator], y3[indicator])],
    #     idxList_[sub2ind(shape, x2[indicator], y3[indicator])],
    #     idxList_[sub2ind(shape, x3[indicator], y3[indicator])],
    #     idxList_[sub2ind(shape, x4[indicator], y3[indicator])],
    #     idxList_[sub2ind(shape, x1[indicator], y4[indicator])],
    #     idxList_[sub2ind(shape, x2[indicator], y4[indicator])],
    #     idxList_[sub2ind(shape, x3[indicator], y4[indicator])],
    #     idxList_[sub2ind(shape, x4[indicator], y4[indicator])])

    list_J = collect(Iterators.flatten(list_J))

    list_val  = ( (1 .-v1) .* (1 .-v2))
    list_val = vcat(list_val,(v2 .* (1 .-v1)))
    list_val = vcat(list_val,(v1 .* (1 .-v2)))
    list_val = vcat(list_val,(v1 .* v2))
    # list_val =vcat(
    #     (v1.*v2.*(v1 .- 1).^2 .*(v2 .- 1).^2)/4,
    #     -1 .*(v1.*(v1 .- 1).^2 .*(3 .*v2.^3 .- 5 .*v2.^2 .+ 2))./4,
    #     -1 .*(v1.*v2.*(v1 .- 1).^2 .*(- 3 .*v2.^2 + 4 .*v2 .+ 1))./4,
    #     -1 .*(v1.*v2.^2 .*(v1 .- 1).^2 .*(v2 .- 1))./4,
    #     -1 .*(v2.*(v2 .- 1).^2 .*(3 .*v1.^3 .- 5 .*v1.^2 .+ 2))./4,
    #     ((3 .*v1.^3 .- 5 .*v1.^2 .+ 2).*(3 .*v2.^3 .- 5 .*v2.^2 .+ 2))./4,
    #     (v2.*(- 3 .*v2.^2 .+ 4 .*v2 .+ 1).*(3 .*v1.^3 .- 5 .*v1.^2 .+ 2))./4,
    #     (v2.^2 .*(v2 .- 1).*(3 .*v1.^3 .- 5 .*v1.^2 .+ 2))./4,
    #     -1 .*(v1.*v2.*(v2 .- 1).^2 .*(- 3 .*v1.^2 .+ 4 .*v1 .+ 1))./4,
    #     (v1.*(- 3 .*v1.^2 .+ 4 .*v1 .+ 1).*(3 .*v2.^3 .- 5 .*v2.^2 .+ 2))./4,
    #     (v1.*v2.*(- 3 .*v1.^2 .+ 4 .*v1 .+ 1).*(- 3 .*v2.^2 .+ 4 .*v2 .+ 1))./4,
    #     (v1.*v2.^2 .*(v2 .- 1).*(- 3 .*v1.^2 .+ 4 .*v1 .+ 1))./4,
    #     -1 .*(v1.^2 .*v2.*(v1 .- 1).*(v2 .- 1).^2)./4,
    #     (v1.^2 .*(v1 .- 1).*(3 .*v2.^3 .- 5 .*v2.^2 .+ 2))./4,
    #     (v1.^2 .*v2.*(v1 .- 1).*(- 3 .*v2.^2 .+ 4 .*v2 .+ 1))./4,
    #     (v1.^2 .*v2.^2 .*(v1 .- 1).*(v2 .- 1))./4
    # )

    list_val = collect(Iterators.flatten(list_val))

    list_I = idxList[indicator]
    #list_I =collect(Iterators.flatten(repeat(list_I,1,16)))
    list_I =collect(Iterators.flatten(repeat(list_I,1,4)))

    # list_I_py, list_J_py, list_val_py = py"py_warping_operator"(flow)
    # list_I_py = list_I_py.+1
    # list_J_py = list_J_py.+1
    # W_mat_py = sparse(list_I_py, list_J_py, list_val_py, H*W, H*W)

    W_mat = dropzeros!(sparse(list_I, list_J, list_val, H*W, H*W))

    return LinearOperator(W_mat)
end

#end

max_stars_repo_name: QPImaging/xfromProjections

max_issues_repo_name: QPImaging/xfromProjections

max_forks_repo_name: QPImaging/xfromProjections


text:
export get_E_all
#ENERGIES
function get_E_all(x)
    #Units are STILL not clear
    V = zero(Float64)
    E_au_au = get_V_au_au(x, V)
    E_coup = get_E_coup(x)
    E_neutral = get_E_neutral(x) + E_au_au
    E_ion = get_E_ion(x) + E_au_au
    hp = [E_neutral, E_ion, E_coup]
    return hp
end

#---------------------------Coupling----------------------------------

function get_E_coup(x::Array{Float64,2})
    E_au_n_coup = get_E_au_n_coup(x)
    E_au_o_coup = get_E_au_o_coup(x)
    return E_au_n_coup + E_au_o_coup
end

#---------------------------Au-N Coupling----------------------------------
function get_E_au_n_coup_loop!(x, V, i)
    delta_xn = @. x[i+2, :] - x[1, :]
    delta_xn[1] = @.delta_xn[1] - cell[1]*round(delta_xn[1]/cell[1])
    delta_xn[2] = @.delta_xn[2] - cell[2]*round(delta_xn[2]/cell[2])

    cutoff_test = @. abs(delta_xn) - cutoff
    bool_cutoff_test = @. cutoff_test <= 0.0
    if sum(bool_cutoff_test) != 0.0
        norm_delta_xn = sqrt(delta_xn' * delta_xn)
        if norm_delta_xn <= cutoff
            V_iterate = coup_a_N * (
            1.0/ (1.0 + coup_b_N * exp(coup_β_N*norm_delta_xn)) - coup_cutoff_N
            )
            V = V + V_iterate
        end
    end
    return V
end

function get_E_au_n_coup(x)
    V = 0.0
    for i in 1:N
        V = get_E_au_n_coup_loop!(x, V, i)
    end
    return V
end
#---------------------------Au-O Coupling----------------------------------
function get_E_au_o_coup_loop!(x, V, i)
    delta_xo = @. x[i+2, :] - x[2, :]
    delta_xo[1] = @.delta_xo[1] - cell[1]*round(delta_xo[1]/cell[1])
    delta_xo[2] = @.delta_xo[2] - cell[2]*round(delta_xo[2]/cell[2])

    cutoff_test = @. abs(delta_xo) - cutoff
    bool_cutoff_test = @. cutoff_test <= 0.0
    if sum(bool_cutoff_test) != 0.0
        norm_delta_xo = sqrt(delta_xo' * delta_xo)
        if norm_delta_xo <= cutoff
            V_coup = coup_a_O /(1.0 + coup_b_O * exp(coup_β_O*norm_delta_xo))
            V_cut = coup_cutoff_O
            V = V + V_coup - V_cut
        end
    end
    return V
end

function get_E_au_o_coup(x)
    V = 0.0
    for i in 1:N
        V = get_E_au_o_coup_loop!(x, V, i)
    end
    return V
end
#---------------------------Ionic----------------------------------

function get_E_ion(x)
    E_n_o_ion= get_E_n_o_ion(x[1:2, :])
    E_au_n_ion = get_E_au_n_ion(x)
    E_au_o_ion = get_E_au_o_ion(x)
    zcom = zcom_NO(x[1:2, :])
    E_image = get_E_image_ionic(zcom)
    return E_n_o_ion + E_au_n_ion+ E_au_o_ion + E_image + KI
end

#---------------------------Ionic N O ----------------------------------

function get_E_n_o_ion(x_no::Array{Float64,2})
    r_no = x_no[1, :] - x_no[2,:]
    norm_r_no = sqrt(r_no' * r_no)
    return F_i * (1 - exp(-δ_i*(norm_r_no - rNO_e_i)))^2
end

#---------------------------Ionic Au O ----------------------------------
function get_E_au_o_ion_loop!(x::Array{Float64,2}, V::Float64, i::Int)
    delta_xo = @. x[i+2, :] - x[2, :]
    delta_xo[1] = @.delta_xo[1] - cell[1]*round(delta_xo[1]/cell[1])
    delta_xo[2] = @.delta_xo[2] - cell[2]*round(delta_xo[2]/cell[2])

    cutoff_test = @. abs(delta_xo) - cutoff
    bool_cutoff_test = @. cutoff_test <= 0.0
    if sum(bool_cutoff_test) != 0.0
        norm_delta_xo = sqrt(delta_xo' * delta_xo)
        if norm_delta_xo <= cutoff
            V_iterate = A_n * (
            exp(-α_i * norm_delta_xo)
            - exp_alpha_i_cutoff
            )
            V = V + V_iterate
        end
    end
    return V
end

function get_E_au_o_ion(x::Array{Float64,2})
    V = 0.0
    for i in 1:N
        V = get_E_au_o_ion_loop!(x, V, i)
    end
    return V
end

#---------------------------Ionic Au N ----------------------------------

function get_E_au_n_ion_loop!(x::Array{Float64,2}, V::Float64, i::Int)
    delta_xn = @. x[i+2, :] - x[1, :]
    delta_xn[1] = @.delta_xn[1] - cell[1]*round(delta_xn[1]/cell[1])
    delta_xn[2] = @.delta_xn[2] - cell[2]*round(delta_xn[2]/cell[2])

    cutoff_test = @. abs(delta_xn) - cutoff
    bool_cutoff_test = @. cutoff_test <= 0.0
    if sum(bool_cutoff_test) != 0.0
        norm_delta_xn = sqrt(delta_xn' * delta_xn)
        if norm_delta_xn <= cutoff
            r_no = x[1, :] - x[2,:]
            norm_r_no = sqrt(r_no' * r_no)
            cos_theta = (x[1, 3] - x[2, 3])/norm_r_no
            V_iterate = @. B_i * (
            exp(-2*β_i *(norm_delta_xn - rN_sur_e_i)) -
            exp_beta_i_cutoff1 -
            2*cos_theta^2 * exp(-β_i *(norm_delta_xn - rN_sur_e_i))
            - exp_beta_i_cutoff2
            )
            V = V + V_iterate
        end
    end
    return V
end

function get_E_au_n_ion(x::Array{Float64,2})
    V = 0.0
    for i in 1:N
        V = get_E_au_n_ion_loop!(x, V, i)
    end
    return V
end

#---------------------------Ionic Image Potential----------------------------------
function get_E_image_ionic(zcom::Float64)
    return -D_i/sqrt(C_i^2 + ((zcom - z0_i))^2)
end
#---------------------------Ionic N-O----------------------------------
function get_E_n_o_ion(x_no::Array{Float64,2})
    r_no = x_no[1, :] - x_no[2,:]
    norm_r_no = sqrt(r_no' * r_no)
    return F_i * (1 - exp(-δ_i*(norm_r_no - rNO_e_i)))^2
end

#---------------------------Neutral----------------------------------

function get_E_neutral(x::Array{Float64,2})
    E_n_o_neutral = get_E_n_o_neutral(x[1:2, :])
    E_au_n_neutral = get_E_au_n_neutral(x)
    E_au_o_neutral = get_E_au_o_neutral(x)
    return E_n_o_neutral + E_au_n_neutral + E_au_o_neutral
end

#---------------------------Neutral N-O----------------------------------

function get_E_n_o_neutral(x_no::Array{Float64,2})
    r_no = @. x_no[1, :] - x_no[2,:]
    norm_r_no = sqrt(r_no' * r_no)
    return F_n * (1 - exp(-δ_n*(norm_r_no - r_0_NO)))^2
end

#---------------------------Neutral AU_O----------------------------------
function get_E_au_o_neutral_loop!(x, i, V)
    delta_xn = @. x[i+2, :] - x[2, :]
    delta_xn[1] = @.delta_xn[1] - cell[1]*round(delta_xn[1]/cell[1])
    delta_xn[2] = @.delta_xn[2] - cell[2]*round(delta_xn[2]/cell[2])

    cutoff_test = @. abs(delta_xn) - cutoff
    bool_cutoff_test = @. cutoff_test <= 0.0
    if sum(bool_cutoff_test) != 0.0
        norm_delta_xn = sqrt(delta_xn' * delta_xn)
        if norm_delta_xn <= cutoff
            V = V + A_n * (exp(-α_n*norm_delta_xn) - exp_alpha_n_cutoff)
        end
    end
    return V
end

function get_E_au_o_neutral(x)
    V = 0.0
    for i in 1:N
        V = get_E_au_o_neutral_loop!(x, i, V)
    end
    return V
end

#---------------------------Neutral AU_N----------------------------------
# function get_E_au_n_neutral_loop!(x::Array{Float64,2}, i::Int, V::Float64)
#     delta_xn = @. x[i+2, :] - x[1, :]
#     delta_xn[1] = @.delta_xn[1] - cell[1]*round(delta_xn[1]/cell[1])
#     delta_xn[2] = @.delta_xn[2] - cell[2]*round(delta_xn[2]/cell[2])
#
#     cutoff_test = @. abs(delta_xn) - cutoff
#     bool_cutoff_test = @. cutoff_test <= 0.0
#     if sum(bool_cutoff_test) != 0.0
#         norm_delta_xn = sqrt(delta_xn' * delta_xn)
#         if norm_delta_xn <= cutoff
#             V = V + B_n * (exp(-β_n*norm_delta_xn) - exp_beta_n_cutoff)
#         end
#     end
#     return V
# end
#
# function get_E_au_n_neutral(x::Array{Float64,2})
#     V = 0.0
#     for i in 1:N
#         V = get_E_au_n_neutral_loop!(x, i, V)
#     end
#     return V
# end
function get_E_au_n_neutral_loop!(x::Array{Float64,2}, i::Int, V::Float64)
    delta_xn = @. x[i+2, :] - x[1, :]
    delta_xn[1] = @.delta_xn[1] - cell[1]*round(delta_xn[1]/cell[1])
    delta_xn[2] = @.delta_xn[2] - cell[2]*round(delta_xn[2]/cell[2])

    cutoff_test = @. abs(delta_xn) - cutoff
    bool_cutoff_test = @. cutoff_test <= 0.0
    if sum(bool_cutoff_test) != 0.0
        norm_delta_xn = sqrt(delta_xn' * delta_xn)
        if norm_delta_xn <= cutoff
            V = V + B_n * (exp(-β_n*norm_delta_xn) - exp_beta_n_cutoff)
        end
    end
    return V
end

function get_E_au_n_neutral(x::Array{Float64,2})
    V = 0.0
    for i in 1:N
        V = get_E_au_n_neutral_loop!(x, i, V)
    end
    return V
end

#---------------------------GOLD AU AU LATTICE----------------------------------
function V_au_au_loop_if!(x::Array{Float64,2}, V::Float64, i::Int, j::Int)
    if nn[i, j] != 0
        m::Int = nn[i, j]+2
        xm = x[m, :]
        xi = x[i+2, :]
        temp = xm .- xi
        temp .= temp./cell
#        @views temp =  (x[m, :] .- x[i+2, :])./cell
        temp_floor = floor.(temp .+ 0.5)
        temp =  temp .- temp_floor
        r =  temp .* cell .- r0[:, j]

        if j == 1 || j== 7
            V = V + dot(r, d1_new * r)
        elseif j == 2 || j== 8
            V = V + dot(r, d2_new * r)
        elseif j == 3 || j== 9
            V = V + dot(r, d3_new * r)
        elseif j == 4 || j== 10
            V = V + dot(r, d4_new * r)
        elseif j == 5|| j== 11
            V = V + dot(r, d5_new * r)
        elseif j == 6 || j== 12
            V = V + dot(r, d6_new * r)
        end
    end
    return V
end
function V_au_au_loop!(x::Array{Float64,2}, V::Float64)
    for i in 1:N
        for j in 1:12
            V = V_au_au_loop_if!(x, V, i, j)
        end
    end
    return V
end

function get_V_au_au(x::Array{Float64,2}, V)
    V = V_au_au_loop!(x, V)
    return V/4.0
end


text:
# some functors (#367)
struct MyModel
    a::Int
end
@model function (f::MyModel)(x)
    m ~ Normal(f.a, 1)
    return x ~ Normal(m, 1)
end
struct MyZeroModel end
@model function (::MyZeroModel)(x)
    m ~ Normal(0, 1)
    return x ~ Normal(m, 1)
end

@testset "model.jl" begin
    @testset "convenience functions" begin
        model = gdemo_default

        # sample from model and extract variables
        vi = VarInfo(model)
        s = vi[@varname(s)]
        m = vi[@varname(m)]

        # extract log pdf of variable object
        lp = getlogp(vi)

        # log prior probability
        lprior = logprior(model, vi)
        @test lprior ≈ logpdf(InverseGamma(2, 3), s) + logpdf(Normal(0, sqrt(s)), m)

        # log likelihood
        llikelihood = loglikelihood(model, vi)
        @test llikelihood ≈ loglikelihood(Normal(m, sqrt(s)), [1.5, 2.0])

        # log joint probability
        ljoint = logjoint(model, vi)
        @test ljoint ≈ lprior + llikelihood
        @test ljoint ≈ lp
    end

    @testset "rng" begin
        model = gdemo_default

        for sampler in (SampleFromPrior(), SampleFromUniform())
            for i in 1:10
                Random.seed!(100 + i)
                vi = VarInfo()
                model(Random.GLOBAL_RNG, vi, sampler)
                vals = DynamicPPL.getall(vi)

                Random.seed!(100 + i)
                vi = VarInfo()
                model(Random.GLOBAL_RNG, vi, sampler)
                @test DynamicPPL.getall(vi) == vals
            end
        end
    end

    @testset "defaults without VarInfo, Sampler, and Context" begin
        model = gdemo_default

        Random.seed!(100)
        s, m = model()

        Random.seed!(100)
        @test model(Random.GLOBAL_RNG) == (s, m)
    end

    @testset "nameof" begin
        @model function test1(x)
            m ~ Normal(0, 1)
            return x ~ Normal(m, 1)
        end
        @model test2(x) = begin
            m ~ Normal(0, 1)
            x ~ Normal(m, 1)
        end
        function test3 end
        @model function (::typeof(test3))(x)
            m ~ Normal(0, 1)
            return x ~ Normal(m, 1)
        end
        function test4 end
        @model function (a::typeof(test4))(x)
            m ~ Normal(0, 1)
            return x ~ Normal(m, 1)
        end

        @test nameof(test1(rand())) == :test1
        @test nameof(test2(rand())) == :test2
        @test nameof(test3(rand())) == :test3
        @test nameof(test4(rand())) == :test4

        # callables
        @test nameof(MyModel(3)(rand())) == Symbol("MyModel(3)")
        @test nameof(MyZeroModel()(rand())) == Symbol("MyZeroModel()")
    end

    @testset "Internal methods" begin
        model = gdemo_default

        # sample from model and extract variables
        vi = VarInfo(model)

        # Second component of return-value of `evaluate!!` should
        # be a `DynamicPPL.AbstractVarInfo`.
        evaluate_retval = DynamicPPL.evaluate!!(model, vi, DefaultContext())
        @test evaluate_retval[2] isa DynamicPPL.AbstractVarInfo

        # Should not return `AbstractVarInfo` when we call the model.
        call_retval = model()
        @test !any(map(x -> x isa DynamicPPL.AbstractVarInfo, call_retval))
    end

    @testset "rand" begin
        model = gdemo_default

        Random.seed!(1776)
        s, m = model()
        sample_namedtuple = (; s=s, m=m)
        sample_dict = Dict(:s => s, :m => m)

        # With explicit RNG
        @test rand(Random.seed!(1776), model) == sample_namedtuple
        @test rand(Random.seed!(1776), NamedTuple, model) == sample_namedtuple
        @test rand(Random.seed!(1776), Dict, model) == sample_dict

        # Without explicit RNG
        Random.seed!(1776)
        @test rand(model) == sample_namedtuple
        Random.seed!(1776)
        @test rand(NamedTuple, model) == sample_namedtuple
        Random.seed!(1776)
        @test rand(Dict, model) == sample_dict
    end
end

max_stars_repo_name: TuringLang/MicroPPL

max_issues_repo_name: TuringLang/MicroPPL

max_forks_repo_name: TuringLang/MicroPPL


text:
using ActorInterfaces.Classic
using QuickActors
import ActorInterfaces.Implementation.Tick
using Test

#abstract type MyMessage end
#Classic.SendStyle(::Type{<:MyMessage}) = Sendable()

mutable struct Counter
    counter::Int
end

struct Increment end #<: MyMessage end

@ctx function (me::Counter)(::Increment)
    me.counter += 1
end

struct Spawner end

struct SpawnTree #<: MyMessage
    childcount::UInt8
    depth::UInt8
end

@ctx function (me::Spawner)(msg::SpawnTree)
    if msg.depth > 0
        for i = 1:msg.childcount
            child = spawn(Spawner())
            send(child, SpawnTree(msg.childcount, msg.depth - 1))
        end
    end
    return nothing
end

# struct NoMessage end

# mutable struct RaceTester
#     arrived::Bool
#     RaceTester() = new(false)
# end
# struct RacingMessage end
# Classic.SendStyle(::Type{RacingMessage}) = Racing()

# function Classic.onmessage(me::Actor{RaceTester}, msg::RacingMessage, ::Racing)
#     behavior(me).arrived = true
# end
# function Classic.onmessage(me::Actor{RaceTester}, msg::RacingMessage)
#     error("Invalid delivery of Racing!")
# end

@testset "Empty Scheduler" begin
    s = QuickScheduler()
    @test isempty(s.msgs)
    @test isempty(s.actorcache)
    while Tick.tick!(s) end
    @test isempty(s.msgs)
    @test isempty(s.actorcache)
end

@testset "Spawning from the outside" begin
    s = QuickScheduler()
    counter = Counter(0)
    counteraddr = spawn!(s, counter)
    @test length(s.actorcache) == 1
    @test counteraddr isa Addr
    @test isempty(s.msgs)
end

@testset "Sending messages from the outside" begin
    s = QuickScheduler()
    counter = Counter(0)
    counteraddr = spawn!(s, counter)
    send!(s, counteraddr, Increment())
    @test counter.counter == 0
    @test !isempty(s.msgs)
    while Tick.tick!(s) end
    @test counter.counter == 1
    @test isempty(s.msgs)
    send!(s, counteraddr, Increment())
    send!(s, counteraddr, Increment())
    send!(s, counteraddr, Increment())
    while Tick.tick!(s) end
    @test counter.counter == 4
end

const TREE_HEIGHT = 19
const TREE_SIZE = 2^(TREE_HEIGHT + 1) - 1

@testset "Tree: Spawning children and sending messages to them" begin
    s = QuickScheduler()
    root = Spawner()
    rootaddr = spawn!(s, root)
    send!(s, rootaddr, SpawnTree(2, TREE_HEIGHT))
    println("Building a tree of $TREE_SIZE actors and delivering the same amount of messages")
    @time while Tick.tick!(s) end
    @test length(s.actorcache) == TREE_SIZE
end

# @testset "Sending NonSendable" begin
#     s = QuickScheduler()
#     root = Spawner()
#     rootaddr = spawn!(s, root)
#     @test_throws Classic.SendingNonSendable send(s.actorcache[rootaddr], rootaddr, NoMessage())
# end

# @testset "Sending Racing" begin
#     s = QuickScheduler()
#     tester = RaceTester()
#     testeraddr = spawn!(s, tester)
#     send!(s, testeraddr, RacingMessage())
#     while Tick.tick!(s) end
#     @test tester.arrived == true
# end

struct BecamePinger
    depth::Int
end

struct BecamePonger
    depth::Int
end

struct BecamePing end # <: MyMessage end
struct BecamePong end #<: MyMessage end

@ctx function (me::BecamePinger)(msg::BecamePong)
    depth = me.depth
    if depth > 0
        become(BecamePonger(depth))
        send(self(), BecamePing())
    end
    return nothing
end

@ctx function (me::BecamePonger)(msg::BecamePing)
    depth = me.depth
    if depth > 0
        become(BecamePinger(depth - 1))
        send(self(), BecamePong())
    end
    return nothing
end

const PING_COUNT = 5_00_000

@testset "Become" begin
    s = QuickScheduler()
    becamer = BecamePinger(PING_COUNT)
    becameraddr = spawn!(s, becamer)
    send!(s, becameraddr, BecamePong())
    println("Becoming $(PING_COUNT * 2) times to a different type and delivering the same amount of messages")
    @time while Tick.tick!(s) end
    @test s.actorcache[becameraddr].actor.behavior.depth == 0
end


text:
using Revise

includet("src/Spring.jl")
IArray = InternalArray

# 1. Solving the equation using the euler method ------------------------------------------

# We are asked to find solutions for the time interval [t0, t0+40].
# I want to execute the following with a step of 1e-3,
# so I'll set an appropiate number of iterations.
const step = 1e-3
const nit = Int(40/step)  # round to avoid small rounding errors

prob = SpringProblem(p0 = 0.8, v0 = 0, m = 1, k = 1.5, b = 0.3, A = 0.4, ω = 2.4)

x = eulermethod(prob, step, nit)

# 2. Plotting the result ------------------------------------------------------------------

timeline = [prob.t0 + i*step for i in 0:length(x)-1]
position = plot(timeline, IArray(x, 1), label = "movement")
speed = plot(timeline, IArray(x, 2), label = "speed")
plot(position, speed)

savefig("media/spring_std.png")

max_stars_repo_name: useredsa/numeric-differential-equations

max_issues_repo_name: useredsa/numeric-differential-equations

max_forks_repo_name: useredsa/numeric-differential-equations