text:
module Cards

export Suit, Card, Hand, ♣, ♢, ♡, ♠, .., deal, points

import Base: *, |, &

"""
Encode a suit as a 2-bit value (low bits of a `UInt8`):

- 0 = ♣ (clubs)
- 1 = ♢ (diamonds)
- 2 = ♡ (hearts)
- 3 = ♠ (spades)

The suits have global constant bindings: `♣`, `♢`, `♡`, `♠`.
"""
struct Suit
    i::UInt8
    Suit(s::Integer) = 0 ≤ s ≤ 3 ? new(s) :
        throw(ArgumentError("invalid suit number: $s"))
end

char(s::Suit) = Char(0x2663-s.i)
Base.string(s::Suit) = string(char(s))
Base.show(io::IO, s::Suit) = print(io, char(s))

const ♣ = Suit(0)
const ♢ = Suit(1)
const ♡ = Suit(2)
const ♠ = Suit(3)

const suits = [♣, ♢, ♡, ♠]

"""
Encode a playing card as a 6-bit integer (low bits of a `UInt8`):

- low bits represent rank from 0 to 15
- high bits represent suit (♣, ♢, ♡ or ♠)

Ranks are assigned as follows:

- numbered cards (2 to 10) have rank equal to their number
- jacks, queens and kings have ranks 11, 12 and 13
- there are low and high aces with ranks 1 and 14
- there are low and high jokers with ranks 0 and 15

This allows any of the standard orderings of cards ranks to be
achieved simply by choosing which aces or which jokers to use.
There are a total of 64 possible card values with this scheme,
represented by `UInt8` values `0x00` through `0x3f`.
"""
struct Card
    value::UInt8
end

function Card(r::Integer, s::Integer)
    0 ≤ r ≤ 15 || throw(ArgumentError("invalid card rank: $r"))
    return Card(((s << 4) % UInt8) | (r % UInt8))
end
Card(r::Integer, s::Suit) = Card(r, s.i)

suit(c::Card) = Suit((0x30 & c.value) >>> 4)
rank(c::Card) = (c.value & 0x0f) % Int8

function Base.show(io::IO, c::Card)
    r = rank(c)
    if 1 ≤ r ≤ 14
        r == 10 && print(io, '1')
        print(io, "1234567890JQKA"[r])
    else
        print(io, '\U1f0cf')
    end
    print(io, suit(c))
end

*(r::Integer, s::Suit) = Card(r, s)

for s in "♣♢♡♠", (r,f) in zip(11:14, "JQKA")
    ss, sc = Symbol(s), Symbol("$f$s")
    @eval (export $sc; const $sc = Card($r,$ss))
end

"""
Represent a hand (set) of cards using a `UInt64` bit set.
"""
struct Hand <: AbstractSet{Card}
    cards::UInt64
    Hand(cards::UInt64) = new(cards)
end

bit(c::Card) = one(UInt64) << c.value
bits(s::Suit) = UInt64(0xffff) << 16(s.i)

function Hand(cards)
    hand = Hand(zero(UInt64))
    for card in cards
        card isa Card || throw(ArgumentError("not a card: $repr(card)"))
        i = bit(card)
        hand.cards & i == 0 || throw(ArgumentError("duplicate cards are not supported"))
        hand = Hand(hand.cards | i)
    end
    return hand
end

Base.in(c::Card, h::Hand) = (bit(c) & h.cards) != 0
Base.length(h::Hand) = count_ones(h.cards)
Base.isempty(h::Hand) = h.cards == 0
Base.lastindex(h::Hand) = length(h)

function Base.iterate(h::Hand, s::UInt8 = trailing_zeros(h.cards) % UInt8)
    (h.cards >>> s) == 0 && return nothing
    c = Card(s); s += true
    c, s + trailing_zeros(h.cards >>> s) % UInt8
end

function Base.unsafe_getindex(h::Hand, i::UInt8)
    card, s = 0x0, 0x5
    while true
        mask = 0xffff_ffff_ffff_ffff >> (0x40 - (0x1<<s) - card)
        card += UInt8(i > count_ones(h.cards & mask) % UInt8) << s
        s > 0 || break
        s -= 0x1
    end
    return Card(card)
end
Base.unsafe_getindex(h::Hand, i::Integer) = Base.unsafe_getindex(h, i % UInt8)

function Base.getindex(h::Hand, i::Integer)
    @boundscheck 1 ≤ i ≤ length(h) || throw(BoundsError(h,i))
    return Base.unsafe_getindex(h, i)
end

function Base.show(io::IO, hand::Hand)
    if isempty(hand) || !get(io, :compact, false)
        print(io, "Hand([")
        for card in hand
            print(io, card)
            (bit(card) << 1) ≤ hand.cards && print(io, ", ")
        end
        print(io, "])")
    else
        for suit in suits
            s = hand & suit
            isempty(s) && continue
            show(io, suit)
            for card in s
                r = rank(card)
                if r == 10
                    print(io, '\u2491')
                elseif 1 ≤ r ≤ 14
                    print(io, "1234567890JQKA"[r])
                else
                    print(io, '\U1f0cf')
                end
            end
        end
    end
end

a::Hand | b::Hand = Hand(a.cards | b.cards)
a::Hand | c::Card = Hand(a.cards | bit(c))
c::Card | h::Hand = h | c

a::Hand & b::Hand = Hand(a.cards & b.cards)
h::Hand & s::Suit = Hand(h.cards & bits(s))
s::Suit & h::Hand = h & s

Base.intersect(s::Suit, h::Hand) = h & s
Base.intersect(h::Hand, s::Suit) = intersect(s::Suit, h::Hand)

*(rr::OrdinalRange{<:Integer}, s::Suit) = Hand(Card(r,s) for r in rr)
..(r::Integer, c::Card) = (r:rank(c))*suit(c)
..(a::Card, b::Card) = suit(a) == suit(b) ? rank(a)..b :
    throw(ArgumentError("card ranges need matching suits: $a vs $b"))

const deck = Hand(Card(r,s) for s in suits for r = 2:14)

Base.empty(::Type{Hand}) = Hand(zero(UInt64))

@eval Base.rand(::Type{Hand}) = Hand($(deck.cards) & rand(UInt64))

function deal!(counts::Vector{<:Integer}, hands::AbstractArray{Hand}, offset::Int=0)
    for rank = 2:14, suit = 0:3
        while true
            hand = rand(1:4)
            if counts[hand] > 0
                counts[hand] -= 1
                hands[offset + hand] |= Card(rank, suit)
                break
            end
        end
    end
    return hands
end

deal() = deal!(fill(13, 4), fill(empty(Hand), 4))

function deal(n::Int)
    counts = fill(0x0, 4)
    hands = fill(empty(Hand), 4, n)
    for i = 1:n
        deal!(fill!(counts, 13), hands, 4(i-1))
    end
    return permutedims(hands)
end

function points(hand::Hand)
    p = 0
    for rank = 11:14, suit = 0:3
        card = Card(rank, suit)
        p += (rank-10)*(card in hand)
    end
    return p
end

end # Cards


text:
using GLVisualize, GeometryTypes, Colors, GLAbstraction
using CompScienceMeshes

const cm = vec(mapslices(
    r->Colors.RGBA{Float32}(r[1],r[2],r[3], 1.0),
    readcsv(Pkg.dir("CompScienceMeshes","examples","cm.csv")),
    [2]))

function patch(mesh::Mesh, amp = nothing)

    if amp == nothing
        amp = Float32[v[end] for v in mesh.vertices]
    else
        vtc, nc = vertextocellmap(mesh)
        amp = Float32[mean(amp[vtc[i,1:nc[i]]]) for i in 1:numvertices(mesh)]
    end

    @show length(amp)

    m, M = extrema(amp)

    #cm = map(x -> RGBA{Float32}(x, 1.0), colormap("RdBu"))
    nc = length(cm)
    at = Float32[(nc-1)*(a-m)/(M-m) for a in amp]

    @assert length(at) == numvertices(mesh)

    vertices = Point3f0[v for v in mesh.vertices]
    faces = GLTriangle[f .- 1 for f in mesh.faces]

    #colors = RGBA{U8}[RGBA{U8}(rand(), rand(), rand(), 1.) for i=1:5]
    #attribute_id = rand(0f0:4f0, length(faces))

    sphere_mesh = GLNormalAttributeMesh(
        vertices=vertices, faces=faces,
        attributes=cm, attribute_id=at)

end

## Example
# m = meshsphere(1.0, 0.2)
# sphere_mesh = patch(m)
#
error("stop")

##
fcr, geo = facecurrents(u, RT)
mp = patch(geo,real.(norm.(fcr)))

##
window = glscreen()
moveright = translationmatrix(Vec3f0(0,0,0))
_view(visualize(mp, model=moveright), window)
renderloop(window)

ext: jl


text:
"""
    ShowMakie([; kwargs...]) <: ShowBackend

A backend for showing block data that uses
[Makie.jl](https://github.com/JuliaPlots/Makie.jl) figures for
visualization.

Keyword arguments are passed through to the constructed `Figure`s.
"""
struct ShowMakie <: ShowBackend
    size::Tuple{Int,Int}
    kwargs::Any
end
ShowMakie(sz = (500, 500); kwargs...) = ShowMakie(sz, kwargs)


text:
__precompile__()

module Insane

include("parser.jl")
include("util.jl")
include("codegen.jl")
include("repl.jl")

export @λ_str, insane, @insane_load

macro λ_str(code)
    insane(">($code)") |> esc
end

function insane(code)
    parser = InsaneParser(code)
    parse_space!(parser)
    expr = parse_expr!(parser)
    codegen!(expr, [])
end

macro insane_load(path)
    code = open(readstring, joinpath(Base.source_dir(), path))
    insane(">($code)") |> esc
end

__init__() = if isdefined(Base, :active_repl)
    run_insane_REPL()
end

end


text:
# demonstrates the use of sixelplot to add a plot
# and overriding the axes limits

using SixelGraphics

x = linspace(0,4*pi)
y1 = sin(x)
y2 = cos(x)

s = sixelplot(x,y1, typ='b', xlim=[-2,14])

sixelplot(s,x,y2, typ='b')


text:
import Base.*

export SWCutoff, SplineCutoff, Shift, HS, C0Shift, C1Shift, C2Shift

# this file is include from Potentials.jl
# i.e. it is part of JuLIP.Potentials
# it implements three types of cutoff potentials.

######################## Stillinger-Weber type cutoff

"""
`cutsw(r, Rc, Lc)`

Implementation of the C^∞ Stillinger-Weber type cut-off potential
    1.0 / ( 1.0 + exp( lc / (Rc-r) ) )

`d_cutinf` implements the first derivative
"""
@inline cutsw(r::T, Rc, Lc) where {T} =
    T(1) / ( T(1) + exp( Lc / ( max(Rc-r, T(0)) + T(1e-2) ) ) )

"derivative of `cutsw`"
@inline function cutsw_d(r::T, Rc, Lc) where {T}
    t = T(1) / ( max(Rc-r, T(0)) + T(1e-2) )    # a numerically stable (Rc-r)^{-1}
    e = T(1) / (T(1) + exp(Lc * t))             # compute exponential only once
    return - Lc * (T(1) - e) * e * t^2
end

"""
`type SWCutoff`: SW type cut-off potential with C^∞ regularity
     1.0 / ( 1.0 + exp( Lc / (Rc-r) ) )

This is not very optimised: one could speed up `evaluate_d` significantly
by avoiding multiple evaluations.

### Parameters

* `Rc` : cut-off radius
* `Lc` : scale
"""
struct SWCutoff{T} <: PairPotential
    Rc::T
    Lc::T
    e0::T
end

@pot SWCutoff

evaluate(p::SWCutoff, r::Number) = p.e0 * cutsw(r, p.Rc, p.Lc)
evaluate_d(p::SWCutoff, r::Number) = p.e0 * cutsw_d(r, p.Rc, p.Lc)
cutoff(p::SWCutoff) = p.Rc

# simplified constructor to ensure compatibility
SWCutoff(Rc, Lc) = SWCutoff(Rc, Lc, one(typeof(Rc)))

# kw-constructor (original SW parameters, pair term)
SWCutoff(; Rc=1.8, Lc=1.0, e0=1.0) = SWCutoff(Rc, Lc, e0)

######################## Shift-Cutoff:

"""
`Shift{ORD}` : a shift-cutoff function

It is not technically a cut-off function but a cut-off operator. Let
`f(r)` with real arguments and range, then
```
g = f * Shift(o, rcut)
```
creates a new function `g(r) = (f(r) - p(r)) * (r < rcut)` where
`p(r)` is a polynomial of order `o` such that all derivatives of `g` up
to order `o` vanish at `r = rcut`. There order must be an integer
between -1 and 2. Choosing `o=-1` creates a Heaviside function.

### Equivalent constructors:
```
HS(r::Float64) = Shift(-1, r)
C0Shift(r::Float64) = Shift(0, r)
C1Shift(r::Float64) = Shift(1, r)
C2Shift(r::Float64) = Shift(2, r)
```

### Example:
```
lj = LennardJones()  # standard lennad-jones potential
V = lj * C2Shift(2.5)
```
Now `V` is a C2-continuous `PairPotential` with support (0, 2.5]."""
struct Shift{ORD, TV, T} <: PairPotential
   ord::Val{ORD}
   V::TV
   rcut::T
   Vcut::T
   dVcut::T
   ddVcut::T
end

@pot Shift

const HS{TV} = Shift{-1, TV}
const C0Shift{TV} = Shift{0, TV}
const C1Shift{TV} = Shift{1, TV}
const C2Shift{TV} = Shift{2, TV}

cutoff(p::Shift) = p.rcut

# the basic constructors
"see documentation for `Shift`"
HS(r::AbstractFloat) = Shift(-1, r)
"see documentation for `Shift`"
C0Shift(r::AbstractFloat) = Shift(0, r)
"see documentation for `Shift`"
C1Shift(r::AbstractFloat) = Shift(1, r)
"see documentation for `Shift`"
C2Shift(r::AbstractFloat) = Shift(2, r)
Shift(o::Int, r::AbstractFloat) = Shift(Val(o), nothing, r, 0.0, 0.0, 0.0)
Shift(V, p::Shift{-1}) = Shift(p.ord, V, p.rcut, 0.0, 0.0, 0.0)
Shift(V, p::Shift{0}) = Shift(p.ord, V, p.rcut, V(p.rcut), 0.0, 0.0)
Shift(V, p::Shift{1}) = Shift(p.ord, V, p.rcut, V(p.rcut), (@D V(p.rcut)), 0.0)
Shift(V, p::Shift{2}) = Shift(p.ord, V, p.rcut, V(p.rcut), (@D V(p.rcut)), (@DD V(p.rcut)))
*(V::PairPotential, p::Shift{ORD, Nothing}) where {ORD} = Shift(V, p)
*(p::Shift{ORD, Nothing}, V::PairPotential) where {ORD} = Shift(V, p)

@inline evaluate(p::Shift{-1}, r::Number) = r < p.rcut ? p.V(r) : 0.0
@inline evaluate_d(p::Shift{-1}, r::Number) = r < p.rcut ? (@D p.V(r)) : 0.0
@inline evaluate_dd(p::Shift{-1}, r::Number) = r < p.rcut ? (@DD p.V(r)) : 0.0

@inline evaluate(p::Shift{0}, r::Number) = r < p.rcut ? (p.V(r) - p.Vcut) : 0.0
@inline evaluate_d(p::Shift{0}, r::Number) = r < p.rcut ? (@D p.V(r)) : 0.0
@inline evaluate_dd(p::Shift{0}, r::Number) = r < p.rcut ? (@DD p.V(r)) : 0.0

@inline evaluate(p::Shift{1}, r::Number) = r >= p.rcut ? 0.0 :
      (p.V(r) - p.Vcut - p.dVcut * (r - p.rcut))
@inline evaluate_d(p::Shift{1}, r::Number) = r < p.rcut ? ((@D p.V(r)) - p.dVcut) : 0.0
@inline evaluate_dd(p::Shift{1}, r::Number) = r < p.rcut ? (@DD p.V(r)) : 0.0

@inline evaluate(p::Shift{2}, r::Number) = r >= p.rcut ? 0.0 :
   (p.V(r) - p.Vcut - p.dVcut * (r - p.rcut) - 0.5 * p.ddVcut * (r-p.rcut)^2)
@inline evaluate_d(p::Shift{2}, r::Number) = r >= p.rcut ? 0.0 :
   ((@D p.V(r)) - p.dVcut - p.ddVcut * (r - p.rcut))
@inline evaluate_dd(p::Shift{2}, r::Number) = r >= p.rcut ? 0.0 :
   ((@DD p.V(r)) - p.ddVcut)

######################## Quintic Spline cut-off:

# TODO: switch to general types ?

@inline function fcut(r, r0, r1)
    s = 1.0 - (r-r0) / (r1-r0)
    return (s >= 1.0) + (0.0 < s < 1.0) * (@fastmath 6.0 * s^5 - 15.0 * s^4 + 10.0 * s^3)
end

"Derivative of `fcut`; see documentation of `fcut`."
@inline function fcut_d(r, r0, r1)
    s = 1-(r-r0) / (r1-r0)
    return -(0 < s < 1) * (@fastmath (30*s^4 - 60 * s^3 + 30 * s^2) / (r1-r0))
end

"Second order derivative of `fcut`; see documentation of `fcut`."
function fcut_dd(r, r0, r1)
  s = 1-(r-r0) / (r1-r0)
  return (0 < s < 1) * (@fastmath (120*s^3 - 180 * s^2 + 60 * s) / (r1-r0)^2)
end

"""
`SplineCutoff` : Piecewise quintic C^{2,1} regular polynomial.

Parameters:

* r0 : inner cut-off radius
* r1 : outer cut-off radius.
"""
mutable struct SplineCutoff <: PairPotential
   r0::Float64
   r1::Float64
end

@pot SplineCutoff

@inline evaluate(p::SplineCutoff, r::Number) = fcut(r, p.r0, p.r1)
@inline evaluate_d(p::SplineCutoff, r::Number) = fcut_d(r, p.r0, p.r1)
evaluate_dd(p::SplineCutoff, r::Number) = fcut_dd(r, p.r0, p.r1)
cutoff(p::SplineCutoff) = p.r1
Base.string(p::SplineCutoff) = "SplineCutoff(r0=$(p.r0), r1=$(p.r1))"

# ============ DRAFT: cos-based cutoff ============

@inline function coscut(r::T, r0, r1) where {T <: AbstractFloat}
   if r <= r0
      return T(1.0)
   elseif r > r1
      return T(0.0)
   else
      return @fastmath T(0.5) * (cos( π * (r-r0) / (r1-r0) ) + T(1.0))
   end
end

@inline function coscut_d(r::T, r0, r1) where {T <: AbstractFloat}
   if r0 < r < r1
      return @fastmath - π/(2*(r1-r0)) * sin( π * (r-r0) / (r1-r0) )
   else
      return T(0.0)
   end
end


text:
function get_date(obj::TimeStamped)
    return jcall(obj, "getDate", AbsoluteDate, ())
end


text:
################################################################################
#
#   IMPLEMENTATIONS OF DIFFERENT LATTICE CONSTRUCTION RELATED STUFF
#
#   STRUCTURE OF THE FILE
#
#   1) CONSTRUCTION BASED ON UNITCELLS (2D & 3D)
#   - Periodic Boundaries
#   - Open Boundaries
#   - Semi-Periodic Boundaries
#   - General Boundaries
#
#   2) CONSTRUCTION BASED ON BOND DISTANCE (2D & 3D)
#
#   3) CONSTRUCTION BASED ON SHAPE (2D & 3D)
#   - General Shape
#   - Sphere
#   - Box
#
################################################################################

#-----------------------------------------------------------------------------------------------------------------------------
#
#   INDIVIDUAL FUNCTIONS FOR LATTICES COMPOSED FROM UNITCELLS
#   BUILDING FUNCTIONS FOR NORMAL PLACEMENT OF CELLS
#
#   Parameters that have to be passed:
#   - Unitcell object from which lattice is build
#   - Array with integers indicating the extent of the lattice along this lattice vector
#
#   Functions will start with most specific and will be generalized later on
#
#-----------------------------------------------------------------------------------------------------------------------------

#-----------------------------------------------------------------------------------------------------------------------------
#
#   Construction of PERIODIC Lattices (i.e. lattices, where all directions have periodic boundary conditions)
#
#-----------------------------------------------------------------------------------------------------------------------------

# FOR 2D AND 3D
function getLatticePeriodic2D(unitcell::Unitcell, repetition_array::Array{Int64}; save::Bool=false, load::Bool=false)

    # extract the cardinal directions of the lattice from the array
    N_a1 = abs(repetition_array[1])
    N_a2 = abs(repetition_array[2])

    # generate the filename of the output
    if contains(unitcell.filename, FOLDER_UNITCELLS)
        filename = replace(unitcell.filename, FOLDER_UNITCELLS, FOLDER_LATTICES)
        filename = replace(filename, ".jld", "_periodic_$(N_a1)_$(N_a2).jld")
        filename = replace(filename, "_unitcell_", "_lattice_")
    else
        filename = "$(FOLDER_LATTICES)$(split(unitcell.filename, FOLDER_UNITCELLS[end])[end])"
        filename = replace(filename, ".jld", "_periodic_$(N_a1)_$(N_a2).jld")
        filename = replace(filename, "_unitcell_", "_lattice_")
    end

    # check if can be loaded
    if isfile(filename) && load
        return Lattice(filename)
    end

    # build the lattice

    # load the data from the unit cell
    uc_basis            = unitcell.basis
    uc_connections      = unitcell.connections
    uc_lattice_vectors  = unitcell.lattice_vectors

    # GENERATE NEW POSITIONS
	positions = Array{Float64, 1}[]
    positions_indices = Int64[]

	# define the index function to get the index of an element in the array
	function index(i,j,alpha)
		return size(uc_basis,1)*((i-1)*N_a2+j-1) + alpha
	end

    # define positions array to be filled
	for element in 1:N_a1*N_a2*size(uc_basis, 1)
		push!(positions, uc_basis[1])
        push!(positions_indices, 1)
	end
	# set all positions to their correct values
	for i in 1:N_a1
	for j in 1:N_a2
    for a in 1:size(uc_basis, 1)
		positions[index(i,j,a)] = uc_basis[a] + i*uc_lattice_vectors[1] + j*uc_lattice_vectors[2]
        positions_indices[index(i,j,a)] = a
	end
	end
    end

    # GENERATE NEW CONNECTIONS
	connections = Array{Any,1}[]

	# iterate over all unit cells
	for i in 1:N_a1
	for j in 1:N_a2
		# add all connections for unitcell (i,j)
		for connection in uc_connections
			# calculate the index from where the connection goes
			index_from = index(i,j,connection[1])
			# calculate the aimed unitcell
			i_to = i + connection[4][1]
			j_to = j + connection[4][2]
			# check if the connection goes around in a1 direction
            offset_a1 = 0
			while i_to < 1
                i_to += N_a1
                offset_a1 += -1
            end
            while i_to > N_a1
                i_to -= N_a1
                offset_a1 += 1
			end
			# check if the connection goes around in a2 direction
            offset_a2 = 0
			while j_to < 1
                j_to += N_a2
                offset_a2 += -1
            end
            while j_to > N_a2
                j_to -= N_a2
                offset_a2 += 1
			end
			# get the index to where the connection goes
			index_to = index(i_to, j_to, connection[2])
            # generate a new connection
            # format [<from index>; <to index>; <strength>; (<lattice displaced by lattice vector j>)]
            connection_new = [index_from; index_to; connection[3]; (offset_a1, offset_a2)]
			# register as connection
			push!(connections, connection_new)
		end
	end
	end

    # generate new lattice vectors, now spanning the WHOLE lattice
    lattice_vectors = Array{Float64, 1}[]
    a1 = uc_lattice_vectors[1] .* N_a1
    a2 = uc_lattice_vectors[2] .* N_a2
    push!(lattice_vectors, a1)
    push!(lattice_vectors, a2)

    # save everything to a Lattice object
    lattice = Lattice(
        unitcell,
        [N_a1, N_a2],
        lattice_vectors,
        positions,
        positions_indices,
        connections,
        filename
        )
    # save the lattice object
    if save
        saveLattice(lattice)
    end

    # return the lattice
    return lattice

end
function getLatticePeriodic3D(unitcell::Unitcell, repetition_array::Array{Int64}; save::Bool=false, load::Bool=false)

    # extract the cardinal directions of the lattice from the array
    N_a1 = abs(repetition_array[1])
    N_a2 = abs(repetition_array[2])
    N_a3 = abs(repetition_array[3])

    # generate the filename of the output
    if contains(unitcell.filename, FOLDER_UNITCELLS)
        filename = replace(unitcell.filename, FOLDER_UNITCELLS, FOLDER_LATTICES)
        filename = replace(filename, ".jld", "_periodic_$(N_a1)_$(N_a2)_$(N_a3).jld")
        filename = replace(filename, "_unitcell_", "_lattice_")
    else
        filename = "$(FOLDER_LATTICES)$(split(unitcell.filename, FOLDER_UNITCELLS[end])[end])"
        filename = replace(filename, ".jld", "_periodic_$(N_a1)_$(N_a2)_$(N_a3).jld")
        filename = replace(filename, "_unitcell_", "_lattice_")
    end

    # check if can be loaded
    if isfile(filename) && load
        return Lattice(filename)
    end

    # build the lattice

    # load the data from the unit cell
    uc_basis            = unitcell.basis
    uc_connections      = unitcell.connections
    uc_lattice_vectors  = unitcell.lattice_vectors

    # GENERATE NEW POSITIONS
	positions = Array{Float64, 1}[]
    positions_indices = Int64[]

	# define the index function to get the index of an element in the array
	function index(i,j,k,alpha)
		return size(uc_basis,1)*((i-1)*N_a2*N_a3 + (j-1)*N_a3 + (k-1)) + alpha
	end

    # define positions array to be filled
	for element in 1:N_a1*N_a2*N_a3*size(uc_basis, 1)
		push!(positions, uc_basis[1])
        push!(positions_indices, 1)
	end
	# set all positions to their correct values
	for i in 1:N_a1
	for j in 1:N_a2
	for k in 1:N_a3
    for a in 1:size(uc_basis, 1)
		positions[index(i,j,k,a)] = uc_basis[a] .+ i.*uc_lattice_vectors[1] .+ j.*uc_lattice_vectors[2] .+ k.*uc_lattice_vectors[3]
        positions_indices[index(i,j,k,a)] = a
	end
	end
    end
    end

    # GENERATE NEW CONNECTIONS
	connections = Array{Any,1}[]

	# iterate over all unit cells
	for i in 1:N_a1
	for j in 1:N_a2
	for k in 1:N_a3
		# add all connections for unitcell (i,j)
		for connection in uc_connections
			# calculate the index from where the connection goes
			index_from = index(i,j,k,Int(connection[1]))
			# calculate the aimed unitcell
			i_to = Int(i + connection[4][1])
			j_to = Int(j + connection[4][2])
			k_to = Int(k + connection[4][3])
			# check if the connection goes around in a1 direction
            offset_a1 = 0
			while i_to < 1
                i_to += N_a1
                offset_a1 += -1
            end
            while i_to > N_a1
                i_to -= N_a1
                offset_a1 += 1
			end
			# check if the connection goes around in a2 direction
            offset_a2 = 0
			while j_to < 1
                j_to += N_a2
                offset_a2 += -1
            end
            while j_to > N_a2
                j_to -= N_a2
                offset_a2 += 1
			end
			# check if the connection goes around in a3 direction
            offset_a3 = 0
			while k_to < 1
                k_to += N_a3
                offset_a3 += -1
            end
            while k_to > N_a3
                k_to -= N_a3
                offset_a3 += 1
			end
			# get the index to where the connection goes
			index_to = index(i_to, j_to, k_to, Int(connection[2]))
            # generate a new connection
            # format [<from index>; <to index>; <strength>; (<lattice displaced by lattice vector j>)]
            connection_new = [index_from; index_to; connection[3]; (offset_a1, offset_a2, offset_a3)]
			# register as connection
			push!(connections, connection_new)
		end
	end
	end
    end

    # generate new lattice vectors, now spanning the WHOLE lattice
    lattice_vectors = Array{Float64, 1}[]
    a1 = uc_lattice_vectors[1] .* N_a1
    a2 = uc_lattice_vectors[2] .* N_a2
    a3 = uc_lattice_vectors[3] .* N_a3
    push!(lattice_vectors, a1)
    push!(lattice_vectors, a2)
    push!(lattice_vectors, a3)

    # save everything to a Lattice object
    lattice = Lattice(
        unitcell,
        [N_a1, N_a2, N_a3],
        lattice_vectors,
        positions,
        positions_indices,
        connections,
        filename
        )
    # save the lattice object
    if save
        saveLattice(lattice)
    end

    # return the lattice object
    return lattice

end

"""
    getLatticePeriodic(unitcell::Unitcell, repetition_array::Array{Int64} [; save::Bool, load::Bool])

Function to construct a finite lattice with periodic boundary conditions on all sides out of a `Unitcell` object.
The number of unitcells that are put together in each elementery direction is passed in the
`repetition_array`.

Additionally, the newly created `Lattice` object can directly be saved. If this has been done before, passing a `load=true`
will allow to load the object instead of creating it again.

Note that this function works for both 2D and 3D unitcells.

# Examples

```julia-repl
julia> lattice = getLatticePeriodic(unitcell, [10, 20])
LatticePhysics.Lattice(...)

julia> lattice = getLatticePeriodic(unitcell, [10, 20], load=true)
LatticePhysics.Lattice(...)
```
"""
function getLatticePeriodic(unitcell::Unitcell, repetition_array::Array{Int64}; save::Bool=false, load::Bool=false)

    # check how many lattice vectors the unitcell has
    number_lv = size(unitcell.lattice_vectors,1)

    # determine which method to call
    if number_lv == 1
        println("Periodic lattices in 1D not implemented yet")
        return
    elseif number_lv == 2
        return getLatticePeriodic2D(unitcell, repetition_array, save=save, load=load)
    elseif number_lv == 3
        return getLatticePeriodic3D(unitcell, repetition_array, save=save, load=load)
    else
        println("Periodic lattices in dimensions larger 3D not implemented yet")
        return
    end

end

export getLatticePeriodic

#-----------------------------------------------------------------------------------------------------------------------------
#
#   Construction of OPEN Lattices (i.e. lattices, where all directions have open boundary conditions)
#
#-----------------------------------------------------------------------------------------------------------------------------

# FOR 2D AND 3D
function getLatticeOpen2D(unitcell::Unitcell, repetition_array::Array{Int64}; save::Bool=false, load::Bool=false)

    # extract the cardinal directions of the lattice from the array
    N_a1 = abs(repetition_array[1])
    N_a2 = abs(repetition_array[2])

    # generate the filename of the output
    if contains(unitcell.filename, FOLDER_UNITCELLS)
        filename = replace(unitcell.filename, FOLDER_UNITCELLS, FOLDER_LATTICES)
        filename = replace(filename, ".jld", "_open_$(N_a1)_$(N_a2).jld")
        filename = replace(filename, "_unitcell_", "_lattice_")
    else
        filename = "$(FOLDER_LATTICES)$(split(unitcell.filename, FOLDER_UNITCELLS[end])[end])"
        filename = replace(filename, ".jld", "_open_$(N_a1)_$(N_a2).jld")
        filename = replace(filename, "_unitcell_", "_lattice_")
    end

    # check if can be loaded
    if isfile(filename) && load
        return Lattice(filename)
    end

    # build the lattice

    # load the data from the unit cell
    uc_basis            = unitcell.basis
    uc_connections      = unitcell.connections
    uc_lattice_vectors  = unitcell.lattice_vectors

    # GENERATE NEW POSITIONS
	positions = Array{Float64, 1}[]
    positions_indices = Int64[]

	# define the index function to get the index of an element in the array
	function index(i,j,alpha)
		return size(uc_basis,1)*((i-1)*N_a2+j-1) + alpha
	end

    # define positions array to be filled
	for element in 1:N_a1*N_a2*size(uc_basis, 1)
		push!(positions, uc_basis[1])
		push!(positions_indices, 1)
	end
	# set all positions to their correct values
	for i in 1:N_a1
	for j in 1:N_a2
    for a in 1:size(uc_basis, 1)
		positions[index(i,j,a)] = uc_basis[a] + i*uc_lattice_vectors[1] + j*uc_lattice_vectors[2]
        positions_indices[index(i,j,a)] = a
	end
	end
    end

    # GENERATE NEW CONNECTIONS
	connections = Array{Any,1}[]

	# iterate over all unit cells
	for i in 1:N_a1
	for j in 1:N_a2
		# add all connections for unitcell (i,j)
		for connection in uc_connections
			# calculate the index from where the connection goes
			index_from = index(i,j,connection[1])
			# calculate the aimed unitcell
			i_to = i + connection[4][1]
			j_to = j + connection[4][2]
			# check if the connection goes around in a1 direction
            offset_a1 = 0
			while i_to < 1
                i_to += N_a1
                offset_a1 += -1
            end
            while i_to > N_a1
                i_to -= N_a1
                offset_a1 += 1
			end
			# check if the connection goes around in a2 direction
            offset_a2 = 0
			while j_to < 1
                j_to += N_a2
                offset_a2 += -1
            end
            while j_to > N_a2
                j_to -= N_a2
                offset_a2 += 1
			end
            # check if it is periodic, if yes: ignore
            if (offset_a1, offset_a2) != (0,0)
                continue
            end
			# get the index to where the connection goes
			index_to = index(i_to, j_to, connection[2])
            # generate a new connection
            # format [<from index>; <to index>; <strength>; (<lattice displaced by lattice vector j>)]
            connection_new = [index_from; index_to; connection[3]; (offset_a1, offset_a2)]
			# register as connection
			push!(connections, connection_new)
		end
	end
	end

    # generate new lattice vectors, now spanning the WHOLE lattice (zero lattice vectors indicating no periodicity = finite)
    lattice_vectors = Array{Float64, 1}[]

    # save everything to a Lattice object
    lattice = Lattice(
        unitcell,
        [N_a1, N_a2],
        lattice_vectors,
        positions,
        positions_indices,
        connections,
        filename
        )
    # save the lattice object
    if save
        saveLattice(lattice)
    end

    # return the lattice
    return lattice

end
function getLatticeOpen3D(unitcell::Unitcell, repetition_array::Array{Int64}; save::Bool=false, load::Bool=false)

    # extract the cardinal directions of the lattice from the array
    N_a1 = abs(repetition_array[1])
    N_a2 = abs(repetition_array[2])
    N_a3 = abs(repetition_array[3])

    # generate the filename of the output
    if contains(unitcell.filename, FOLDER_UNITCELLS)
        filename = replace(unitcell.filename, FOLDER_UNITCELLS, FOLDER_LATTICES)
        filename = replace(filename, ".jld", "_open_$(N_a1)_$(N_a2)_$(N_a3).jld")
        filename = replace(filename, "_unitcell_", "_lattice_")
    else
        filename = "$(FOLDER_LATTICES)$(split(unitcell.filename, FOLDER_UNITCELLS[end])[end])"
        filename = replace(filename, ".jld", "_open_$(N_a1)_$(N_a2)_$(N_a3).jld")
        filename = replace(filename, "_unitcell_", "_lattice_")
    end

    # check if can be loaded
    if isfile(filename) && load
        return Lattice(filename)
    end

    # build the lattice

    # load the data from the unit cell
    uc_basis            = unitcell.basis
    uc_connections      = unitcell.connections
    uc_lattice_vectors  = unitcell.lattice_vectors

    # GENERATE NEW POSITIONS
	positions = Array{Float64, 1}[]
    positions_indices = Int64[]

	# define the index function to get the index of an element in the array
	function index(i,j,k,alpha)
		return size(uc_basis,1)*((i-1)*N_a2*N_a3 + (j-1)*N_a3 + (k-1)) + alpha
	end

    # define positions array to be filled
	for element in 1:N_a1*N_a2*N_a3*size(uc_basis, 1)
		push!(positions, uc_basis[1])
		push!(positions_indices, 1)
	end
	# set all positions to their correct values
	for i in 1:N_a1
	for j in 1:N_a2
	for k in 1:N_a3
    for a in 1:size(uc_basis, 1)
		positions[index(i,j,k,a)] = uc_basis[a] + i*uc_lattice_vectors[1] + j*uc_lattice_vectors[2] + k*uc_lattice_vectors[3]
        positions_indices[index(i,j,k,a)] = a
	end
	end
    end
    end

    # GENERATE NEW CONNECTIONS
	connections = Array{Any,1}[]

	# iterate over all unit cells
	for i in 1:N_a1
	for j in 1:N_a2
	for k in 1:N_a3
		# add all connections for unitcell (i,j)
		for connection in uc_connections
			# calculate the index from where the connection goes
			index_from = index(i,j,k,connection[1])
			# calculate the aimed unitcell
			i_to = i + connection[4][1]
			j_to = j + connection[4][2]
			k_to = k + connection[4][3]
			# check if the connection goes around in a1 direction
            offset_a1 = 0
			while i_to < 1
                i_to += N_a1
                offset_a1 += -1
            end
            while i_to > N_a1
                i_to -= N_a1
                offset_a1 += 1
			end
			# check if the connection goes around in a2 direction
            offset_a2 = 0
			while j_to < 1
                j_to += N_a2
                offset_a2 += -1
            end
            while j_to > N_a2
                j_to -= N_a2
                offset_a2 += 1
			end
			# check if the connection goes around in a3 direction
            offset_a3 = 0
			while k_to < 1
                k_to += N_a3
                offset_a3 += -1
            end
            while k_to > N_a3
                k_to -= N_a3
                offset_a3 += 1
			end
            # check if it is periodic, if yes: ignore
            if (offset_a1, offset_a2, offset_a3) != (0,0,0)
                continue
            end
			# get the index to where the connection goes
			index_to = index(i_to, j_to, k_to, connection[2])
            # generate a new connection
            # format [<from index>; <to index>; <strength>; (<lattice displaced by lattice vector j>)]
            connection_new = [index_from; index_to; connection[3]; (offset_a1, offset_a2, offset_a3)]
			# register as connection
			push!(connections, connection_new)
		end
	end
	end
    end

    # generate new lattice vectors, now spanning the WHOLE lattice (zero lattice vectors indicating no periodicity = finite)
    lattice_vectors = Array{Float64, 1}[]

    # save everything to a Lattice object
    lattice = Lattice(
        unitcell,
        [N_a1, N_a2, N_a3],
        lattice_vectors,
        positions,
        positions_indices,
        connections,
        filename
        )
    # save the lattice object
    if save
        saveLattice(lattice)
    end

    # return the lattice object
    return lattice

end

"""
    getLatticeOpen(unitcell::Unitcell, repetition_array::Array{Int64} [; save::Bool, load::Bool])

Function to construct a finite lattice with open boundary conditions
(bonds are cut when connecting to sites outside the wanted area) on all sides out of a `Unitcell` object.
The number of unitcells that are put together in each elementery direction is passed in the
`repetition_array`.

Additionally, the newly created `Lattice` object can directly be saved. If this has been done before, passing a `load=true`
will allow to load the object instead of creating it again.

Note that this function works for both 2D and 3D unitcells.

# Examples

```julia-repl
julia> lattice = getLatticeOpen(unitcell, [10, 20])
LatticePhysics.Lattice(...)

julia> lattice = getLatticeOpen(unitcell, [10, 20], load=true)
LatticePhysics.Lattice(...)
```
"""
function getLatticeOpen(unitcell::Unitcell, repetition_array::Array{Int64}; save::Bool=false, load::Bool=false)

    # check how many lattice vectors the unitcell has
    number_lv = size(unitcell.lattice_vectors,1)

    # determine which method to call
    if number_lv == 1
        println("Open lattices in 1D not implemented yet")
        return
    elseif number_lv == 2
        return getLatticeOpen2D(unitcell, repetition_array, save=save, load=load)
    elseif number_lv == 3
        return getLatticeOpen3D(unitcell, repetition_array, save=save, load=load)
    else
        println("Open lattices in dimensions larger 3D not implemented yet")
        return
    end

end

export getLatticeOpen

#-----------------------------------------------------------------------------------------------------------------------------
#
#   Construction of SEMIPERIODIC Lattices
#   (i.e. lattices, where some directions have open and some have periodic boundary conditions)
#
#   negative numbers denote periodic directions
#   positive numbers denote open directions
#
#-----------------------------------------------------------------------------------------------------------------------------

# FOR 2D AND 3D
function getLatticeSemiperiodic2D(unitcell::Unitcell, repetition_array::Array{Int64}; save::Bool=false, load::Bool=false)

    # extract the cardinal directions of the lattice from the array
    N_a1 = repetition_array[1]
    N_a2 = repetition_array[2]

    # check if correct method is called
    if N_a1*N_a2 > 0
        # something is wrong, either flake or totally periodic
        println("wrong method called, either no periodic direction or all directions periodic desired, but can only deliver semiperiodic!")
        return
    end

    # generate the filename of the output
    if contains(unitcell.filename, FOLDER_UNITCELLS)
        filename = replace(unitcell.filename, FOLDER_UNITCELLS, FOLDER_LATTICES)
        filename = replace(filename, ".jld", "_semiperiodic_$(N_a1)_$(N_a2).jld")
        filename = replace(filename, "_unitcell_", "_lattice_")
    else
        filename = "$(FOLDER_LATTICES)$(split(unitcell.filename, FOLDER_UNITCELLS[end])[end])"
        filename = replace(filename, ".jld", "_semiperiodic_$(N_a1)_$(N_a2).jld")
        filename = replace(filename, "_unitcell_", "_lattice_")
    end

    # check if can be loaded
    if isfile(filename) && load
        return Lattice(filename)
    end

    # build the lattice

    # load the data from the unit cell
    uc_basis            = unitcell.basis
    uc_connections      = unitcell.connections
    uc_lattice_vectors  = unitcell.lattice_vectors

    # GENERATE NEW POSITIONS
	positions = Array{Float64, 1}[]
    positions_indices = Int64[]

    # check which direction periodic
    if N_a1 < 0
        N_a1 *= -1
        periodic_direction = 1
    else
        N_a2 *= -1
        periodic_direction = 2
    end

	# define the index function to get the index of an element in the array
	function index(i,j,alpha)
		return size(uc_basis,1)*((i-1)*N_a2+j-1) + alpha
	end

    # define positions array to be filled
	for element in 1:N_a1*N_a2*size(uc_basis, 1)
		push!(positions, uc_basis[1])
        push!(positions_indices, 1)
	end
	# set all positions to their correct values
	for i in 1:N_a1
	for j in 1:N_a2
    for a in 1:size(uc_basis, 1)
		positions[index(i,j,a)] = uc_basis[a] + i*uc_lattice_vectors[1] + j*uc_lattice_vectors[2]
        positions_indices[index(i,j,a)] = a
	end
	end
    end

    # GENERATE NEW CONNECTIONS
	connections = Array{Any,1}[]

	# iterate over all unit cells
	for i in 1:N_a1
	for j in 1:N_a2
		# add all connections for unitcell (i,j)
		for connection in uc_connections
			# calculate the index from where the connection goes
			index_from = index(i,j,connection[1])
			# calculate the aimed unitcell
			i_to = i + connection[4][1]
			j_to = j + connection[4][2]
			# check if the connection goes around in a1 direction
            offset_a1 = 0
			while i_to < 1
                i_to += N_a1
                offset_a1 += -1
            end
            while i_to > N_a1
                i_to -= N_a1
                offset_a1 += 1
			end
			# check if the connection goes around in a2 direction
            offset_a2 = 0
			while j_to < 1
                j_to += N_a2
                offset_a2 += -1
            end
            while j_to > N_a2
                j_to -= N_a2
                offset_a2 += 1
			end
            # check if it outside the open boundary
            if periodic_direction == 1 && offset_a2 != 0
                continue
            elseif periodic_direction == 2 && offset_a1 != 0
                continue
            end
			# get the index to where the connection goes
			index_to = index(i_to, j_to, connection[2])
            # generate a new connection
            # format [<from index>; <to index>; <strength>; (<lattice displaced by lattice vector j>)]
            if periodic_direction == 1
                connection_new = [index_from; index_to; connection[3]; (offset_a1)]
            else
                connection_new = [index_from; index_to; connection[3]; (offset_a2)]
            end
			# register as connection
			push!(connections, connection_new)
		end
	end
	end

    # generate new lattice vectors, now spanning the WHOLE lattice (zero lattice vectors indicating no periodicity = finite)
    lattice_vectors = Array{Float64, 1}[]
    if periodic_direction == 1
        a1 = uc_lattice_vectors[1] .* N_a1
        push!(lattice_vectors, a1)
    else
        a2 = uc_lattice_vectors[2] .* N_a2
        push!(lattice_vectors, a2)
    end

    # save everything to a Lattice object
    lattice = Lattice(
        unitcell,
        [N_a1, N_a2],
        lattice_vectors,
        positions,
        positions_indices,
        connections,
        filename
        )
    # save the lattice object
    if save
        saveLattice(lattice)
    end

    # return the lattice
    return lattice
end
function getLatticeSemiperiodic3D(unitcell::Unitcell, repetition_array::Array{Int64}; save::Bool=false, load::Bool=false)

    # extract the cardinal directions of the lattice from the array
    N_a1 = repetition_array[1]
    N_a2 = repetition_array[2]
    N_a3 = repetition_array[3]

    # check which directions are periodic
    periodic_directions = []
    if N_a1 < 0
        push!(periodic_directions, 1)
    end
    if N_a2 < 0
        push!(periodic_directions, 2)
    end
    if N_a3 < 0
        push!(periodic_directions, 3)
    end

    # check if correct method called
    if length(periodic_directions) in [0,3]
        # something is wrong, either flake or totally periodic
        println("wrong method called, either no periodic direction or all directions periodic desired, but can only deliver semiperiodic!")
        return
    end

    # generate the filename of the output
    if contains(unitcell.filename, FOLDER_UNITCELLS)
        filename = replace(unitcell.filename, FOLDER_UNITCELLS, FOLDER_LATTICES)
        filename = replace(filename, ".jld", "_semiperiodic_$(N_a1)_$(N_a2)_$(N_a3).jld")
        filename = replace(filename, "_unitcell_", "_lattice_")
    else
        filename = "$(FOLDER_LATTICES)$(split(unitcell.filename, FOLDER_UNITCELLS[end])[end])"
        filename = replace(filename, ".jld", "_semiperiodic_$(N_a1)_$(N_a2)_$(N_a3).jld")
        filename = replace(filename, "_unitcell_", "_lattice_")
    end

    # check if can be loaded
    if isfile(filename) && load
        return Lattice(filename)
    end

    # build the lattice

    # load the data from the unit cell
    uc_basis            = unitcell.basis
    uc_connections      = unitcell.connections
    uc_lattice_vectors  = unitcell.lattice_vectors

    # GENERATE NEW POSITIONS
	positions = Array{Float64, 1}[]
	positions_indices = Int64[]

    # turn all directions to positive numbers
    N_a1 = abs(N_a1)
    N_a2 = abs(N_a2)
    N_a3 = abs(N_a3)

	# define the index function to get the index of an element in the array
	function index(i,j,k,alpha)
		return size(uc_basis,1)*((i-1)*N_a2*N_a3 + (j-1)*N_a3 + (k-1)) + alpha
	end

    # define positions array to be filled
	for element in 1:N_a1*N_a2*N_a3*size(uc_basis, 1)
		push!(positions, uc_basis[1])
        push!(positions_indices, 1)
	end
	# set all positions to their correct values
	for i in 1:N_a1
	for j in 1:N_a2
	for k in 1:N_a3
    for a in 1:size(uc_basis, 1)
		positions[index(i,j,k,a)] = uc_basis[a] + i*uc_lattice_vectors[1] + j*uc_lattice_vectors[2] + k*uc_lattice_vectors[3]
        positions_indices[index(i,j,k,a)] = a
	end
	end
    end
    end

    # GENERATE NEW CONNECTIONS
	connections = Array{Any,1}[]

	# iterate over all unit cells
	for i in 1:N_a1
	for j in 1:N_a2
	for k in 1:N_a3
		# add all connections for unitcell (i,j)
		for connection in uc_connections
			# calculate the index from where the connection goes
			index_from = index(i,j,k,connection[1])
			# calculate the aimed unitcell
			i_to = i + connection[4][1]
			j_to = j + connection[4][2]
			k_to = k + connection[4][3]
			# check if the connection goes around in a1 direction
            offset_a1 = 0
			while i_to < 1
                i_to += N_a1
                offset_a1 += -1
            end
            while i_to > N_a1
                i_to -= N_a1
                offset_a1 += 1
			end
			# check if the connection goes around in a2 direction
            offset_a2 = 0
			while j_to < 1
                j_to += N_a2
                offset_a2 += -1
            end
            while j_to > N_a2
                j_to -= N_a2
                offset_a2 += 1
			end
			# check if the connection goes around in a3 direction
            offset_a3 = 0
			while k_to < 1
                k_to += N_a3
                offset_a3 += -1
            end
            while k_to > N_a3
                k_to -= N_a3
                offset_a3 += 1
			end
            # offsets
            offsets = []
            # check if the offsets are correct
            if offset_a1 != 0 && !(1 in periodic_directions)
                continue
            elseif 1 in periodic_directions
                push!(offsets, offset_a1)
            end
            if offset_a2 != 0 && !(2 in periodic_directions)
                continue
            elseif 2 in periodic_directions
                push!(offsets, offset_a2)
            end
            if offset_a3 != 0 && !(3 in periodic_directions)
                continue
            elseif 3 in periodic_directions
                push!(offsets, offset_a3)
            end
			# get the index to where the connection goes
			index_to = index(i_to, j_to, k_to, connection[2])
            # generate a new connection
            # format [<from index>; <to index>; <strength>; (<lattice displaced by lattice vector j>)]
            if length(offsets) == 1
                connection_new = [index_from; index_to; connection[3]; offsets[1]]
            else # offset has to have 2 entries
                connection_new = [index_from; index_to; connection[3]; (offsets[1], offsets[2])]
            end
			# register as connection
			push!(connections, connection_new)
		end
	end
	end
    end

    # generate new lattice vectors, now spanning the WHOLE lattice
    lattice_vectors = Array{Float64, 1}[]
    a1 = uc_lattice_vectors[1] .* N_a1
    a2 = uc_lattice_vectors[2] .* N_a2
    a3 = uc_lattice_vectors[3] .* N_a3
    if 1 in periodic_directions
        push!(lattice_vectors, a1)
    end
    if 2 in periodic_directions
        push!(lattice_vectors, a2)
    end
    if 3 in periodic_directions
        push!(lattice_vectors, a3)
    end

    # save everything to a Lattice object
    lattice = Lattice(
        unitcell,
        [N_a1, N_a2, N_a3],
        lattice_vectors,
        positions,
        positions_indices,
        connections,
        filename
        )
    # save the lattice object
    if save
        saveLattice(lattice)
    end

    # return the lattice object
    return lattice
end

"""
    getLatticeSemiperiodic(unitcell::Unitcell, repetition_array::Array{Int64} [; save::Bool, load::Bool])

Function to construct a finite lattice with a mix of open and periodic boundary conditions out of a `Unitcell` object.
The number of unitcells that are put together in each elementery direction is passed in the
`repetition_array`, where negative entries denote periodic boundary conditions whereas positive entries denote open boundary conditions.
Note that this function can only handle a mix of periodic and open but not purely open or periodic boundary conditions.

Also, the newly created `Lattice` object can directly be saved. If this has been done before, passing a `load=true`
will allow to load the object instead of creating it again.

Note that this function works for both 2D and 3D unitcells.

# Examples

```julia-repl
julia> lattice = getLatticeSemiperiodic(unitcell, [10, -20])
LatticePhysics.Lattice(...)

julia> lattice = getLatticeSemiperiodic(unitcell, [10, -20], load=true)
LatticePhysics.Lattice(...)
```
"""
function getLatticeSemiperiodic(unitcell::Unitcell, repetition_array::Array{Int64}; save::Bool=false, load::Bool=false)

    # check how many lattice vectors the unitcell has
    number_lv = size(unitcell.lattice_vectors,1)

    # determine which method to call
    if number_lv == 1
        println("Semiperiodic lattices in 1D not implemented yet")
        return
    elseif number_lv == 2
        return getLatticeSemiperiodic2D(unitcell, repetition_array, save=save, load=load)
    elseif number_lv == 3
        return getLatticeSemiperiodic3D(unitcell, repetition_array, save=save, load=load)
    else
        println("Semiperiodic lattices in dimensions larger 3D not implemented yet")
        return
    end

end

export getLatticeSemiperiodic

"""
    getLattice(unitcell::Unitcell, repetition_array::Array{Int64} [; save::Bool, load::Bool])
    getLattice(unitcell::Unitcell, repetitions::Int64 [; save::Bool, load::Bool])

Function to construct a finite lattice with a mix of open and periodic boundary conditions out of a `Unitcell` object.
The number of unitcells that are put together in each elementery direction is passed in the
`repetition_array`, where negative entries denote periodic boundary conditions whereas positive entries denote open boundary conditions.
Note that this function CAN handle any mix of periodic and open including purely open or periodic boundary conditions.

Also, the newly created `Lattice` object can directly be saved. If this has been done before, passing a `load=true`
will allow to load the object instead of creating it again.

Note that this function works for both 2D and 3D unitcells.

# Examples

```julia-repl
julia> lattice = getLattice(unitcell, [10, -20])
LatticePhysics.Lattice(...)

julia> lattice = getLattice(unitcell, [10, -20], load=true)
LatticePhysics.Lattice(...)

julia> lattice = getLattice(unitcell, -20, load=true)
LatticePhysics.Lattice(...)
```
"""
function getLattice(unitcell::Unitcell, repetition_array::Array{Int64}; save=false, load=false)

    # check how many periodic directions
    number_pd = 0
    for N in repetition_array
        if N < 0
            number_pd = number_pd + 1
        end
    end

    # check how many open directions
    number_od = 0
    for N in repetition_array
        if N > 0
            number_od = number_od + 1
        end
    end

    # check if all directions are either periodic or open
    if number_od + number_pd != length(repetition_array)
        println("Some directions are not specified correctly! Abort...")
        return
    end

    # check which method to call
    if number_pd == 0
        # no periodc directions present, OPEN
        return getLatticeOpen(unitcell, repetition_array, save=save, load=load)
    elseif number_od == 0
        # no open directions present, PERIODIC (but multiply all ranges by -1)
        return getLatticePeriodic(unitcell, repetition_array.*-1, save=save, load=load)
    else
        # both periodic and open directions present, SEMIPERIODIC
        return getLatticeSemiperiodic(unitcell, repetition_array, save=save, load=load)
    end

end
function getLattice(unitcell::Unitcell, repetitions::Int64; save=false, load=false)
    # just parse through to more general method
    return getLattice(unitcell, repetitions.*ones(Int64, size(unitcell.lattice_vectors,1)), save=save, load=load)
end

export getLattice

#-----------------------------------------------------------------------------------------------------------------------------
#
#   BUILDING FUNCTIONS FOR LATTICES THAT ARE BUILD BY BOND DISTANCE FROM AN ORIGIN SITE
#   These lattices have of course open boundaries
#
#   Parameters that have to be passed:
#   - Unitcell object from which lattice is build
#   - Integer with extent of the lattice from the origin site outwards
#   - OPTIONAL: origin site index
#
#-----------------------------------------------------------------------------------------------------------------------------

# generate a lattice in 2D by bond distance (open boundary conditions)
function getLatticeByBondDistance2D(unitcell::Unitcell, bonddistance::Int64; origin::Int64=1, load::Bool=false, save::Bool=false)

    # generate the filename of the output
    if contains(unitcell.filename, FOLDER_UNITCELLS)
        filename = replace(unitcell.filename, FOLDER_UNITCELLS, FOLDER_LATTICES)
        filename = replace(filename, ".jld", "_by_bonddistance_$(bonddistance)_from_$(origin).jld")
        filename = replace(filename, "_unitcell_", "_lattice_")
    else
        filename = "$(FOLDER_LATTICES)$(split(unitcell.filename, FOLDER_UNITCELLS[end])[end])"
        filename = replace(filename, ".jld", "_by_bonddistance_$(bonddistance)_from_$(origin).jld")
        filename = replace(filename, "_unitcell_", "_lattice_")
    end

    # check if the filename already exists
    if isfile(filename) && load
        # return the lattice
        return Lattice(filename)
    end

    # load the data from the unit cell
    uc_basis            = unitcell.basis
    uc_connections      = unitcell.connections
    uc_lattice_vectors  = unitcell.lattice_vectors

    # arrays for new positions and connections
    positions   = []
    connections = Array{Any,1}[]

    # checklist for all sites that are checked if added etc
    checklist   = []
    # format for sites on checklist: ([pos_x, pos_y], [uc_i, uc_j], index_in_UC, bd_current)

    # push the origin site to the checklist
    push!(
        checklist, (uc_basis[origin], [0,0], origin, 0)
    )

    # iterate while the checklist is not empty
    while size(checklist,1) > 0

        # get the item that was on the checklist for longest
        item_to_handle = shift!(checklist)

        # check if the item is already in the positions list
        found = false
        for p in positions
            if item_to_handle[2] == p[2] && item_to_handle[3] == p[3]
                # if yes, continue
                found = true
                break
            end
        end
        if found
            continue
        end

        # if not, push it into
        push!(positions, item_to_handle)
        # for all connections
        index_from = size(positions, 1)

        # insert all connections to sites that are already inside the positions list and all other sites into the checklist
        for c in uc_connections
            # check if correct connections
            if c[1] != item_to_handle[3]
                continue
            end
            # search for the other element of the connection
            i_to = item_to_handle[2][1] + c[4][1]
            j_to = item_to_handle[2][2] + c[4][2]
            a_to = c[2]
            index_to = -1
            for (index,item_handled) in enumerate(positions)
                if item_handled[2] == [i_to, j_to] && item_handled[3] == a_to
                    # make a new connection
                    index_to = index
                    # break the loop
                    break
                end
            end
            # determine whether the element is already inside the list
            if index_to > 0
                # element is at index index_to an can be linked
                connection_new = [index_from; index_to; c[3]; (0, 0)]
                # check if the connetion is already added
                if connection_new in connections
                    continue
                end
                # register as connection
                push!(connections, connection_new)
            else
                # element not in list yet, maybe should be added
                if item_to_handle[4] < bonddistance
                    # format for sites on checklist: ([pos_x, pos_y], [uc_i, uc_j], index_in_UC, bd_current)
                    push!(checklist, (
                        (i_to * uc_lattice_vectors[1]) .+ (j_to * uc_lattice_vectors[2]) .+ uc_basis[c[2]],
                        [i_to, j_to],
                        c[2],
                        item_to_handle[4]+1
                    ))
                end
            end
        end

    end

    # change the format of positions
    positions_TMP = positions

    # erase positions
    positions = Array{Float64, 1}[]
    positions_indices = Int64[]

    # insert the real positions
    for p in positions_TMP
        push!(positions, p[1])
        push!(positions_indices, p[3])
    end

    # insert missing connections (if (i to j) is present, insert (j to i))
    for c in connections
        c_proposed = [c[2]; c[1]; c[3]; (0, 0)]
        if !(c_proposed in connections)
            push!(connections, c_proposed)
        else
            #println("connection already there")
        end
    end

    # save everything to a Lattice object
    lattice = Lattice(
            unitcell,
            Int64[],
            Array{Float64, 1}[],
            positions,
            positions_indices,
            connections,
            filename
        )
    # save the lattice object
    if save
        saveLattice(lattice)
    end

    # return the lattice object
    return lattice

end
# generate a lattice in 3D by bond distance (open boundary conditions)
function getLatticeByBondDistance3D(unitcell::Unitcell, bonddistance::Int64; origin::Int64=1, load::Bool=false, save::Bool=false)

    # generate the filename of the output
    if contains(unitcell.filename, FOLDER_UNITCELLS)
        filename = replace(unitcell.filename, FOLDER_UNITCELLS, FOLDER_LATTICES)
        filename = replace(filename, ".jld", "_by_bonddistance_$(bonddistance)_from_$(origin).jld")
        filename = replace(filename, "_unitcell_", "_lattice_")
    else
        filename = "$(FOLDER_LATTICES)$(split(unitcell.filename, FOLDER_UNITCELLS[end])[end])"
        filename = replace(filename, ".jld", "_by_bonddistance_$(bonddistance)_from_$(origin).jld")
        filename = replace(filename, "_unitcell_", "_lattice_")
    end

    # check if the filename already exists
    if isfile(filename) && load
        # return the lattice
        return Lattice(filename)
    end

    # load the data from the unit cell
    uc_basis            = unitcell.basis
    uc_connections      = unitcell.connections
    uc_lattice_vectors  = unitcell.lattice_vectors

    # arrays for new positions and connections
    positions   = []
    connections = Array{Any, 1}[]

    # checklist for all sites that are checked if added etc
    checklist   = []
    # format for sites on checklist: ([pos_x, pos_y, pos_z], [uc_i, uc_j, uc_k], index_in_UC, bd_current)

    # push the origin site to the checklist
    push!(
        checklist, (uc_basis[origin], [0,0,0], origin, 0)
    )

    # iterate while the checklist is not empty
    while size(checklist,1) > 0

        # get the item that was on the checklist for longest
        item_to_handle = shift!(checklist)

        # check if the item is already in the positions list
        found = false
        for p in positions
            if item_to_handle[2] == p[2] && item_to_handle[3] == p[3]
                # if yes, continue
                found = true
                break
            end
        end
        if found
            continue
        end

        # if not, push it into
        push!(positions, item_to_handle)
        # for all connections
        index_from = size(positions, 1)

        # insert all connections to sites that are already inside the positions list and all other sites into the checklist
        for c in uc_connections
            # check if correct connections
            if c[1] != item_to_handle[3]
                continue
            end
            # search for the other element of the connection
            i_to = item_to_handle[2][1] + c[4][1]
            j_to = item_to_handle[2][2] + c[4][2]
            k_to = item_to_handle[2][3] + c[4][3]
            a_to = c[2]
            index_to = -1
            for (index,item_handled) in enumerate(positions)
                if item_handled[2] == [i_to, j_to, k_to] && item_handled[3] == a_to
                    # make a new connection
                    index_to = index
                    # break the loop
                    break
                end
            end
            # determine whether the element is already inside the list
            if index_to > 0
                # element is at index index_to an can be linked
                connection_new = [index_from; index_to; c[3]; (0, 0, 0)]
                # check if the connetion is already added
                if connection_new in connections
                    continue
                end
                # register as connection
                push!(connections, connection_new)
            else
                # element not in list yet, maybe should be added
                if item_to_handle[4] < bonddistance
                    # format for sites on checklist: ([pos_x, pos_y, pos_z], [uc_i, uc_j, uc_k], index_in_UC, bd_current)
                    push!(checklist, (
                        (i_to * uc_lattice_vectors[1]) .+ (j_to * uc_lattice_vectors[2]) .+ (k_to * uc_lattice_vectors[3]) .+ uc_basis[c[2]],
                        [i_to, j_to, k_to],
                        c[2],
                        item_to_handle[4]+1
                    ))
                end
            end
        end
    end

    # change the format of positions
    positions_TMP = positions

    # erase positions
    positions = Array{Float64, 1}[]
    positions_indices = Int64[]

    # insert the real positions
    for p in positions_TMP
        push!(positions, p[1])
        push!(positions_indices, p[3])
    end

    # insert missing connections (if (i to j) is present, insert (j to i))
    for c in connections
        c_proposed = [c[2]; c[1]; c[3]; (0, 0, 0)]
        if !(c_proposed in connections)
            push!(connections, c_proposed)
        else
            #println("connection already there")
        end
    end

    # save everything to a Lattice object
    lattice = Lattice(
            unitcell,
            Int64[],
            Array{Float64, 1}[],
            positions,
            positions_indices,
            connections,
            filename
        )
    # save the lattice object
    if save
        saveLattice(lattice)
    end

    # return the lattice object
    return lattice
end

"""
    getLatticeByBondDistance(unitcell::Unitcell, bonddistance::Int64 [; origin::Int64=1, load::Bool=false, save::Bool=false])

Function to construct a finite lattice (with open boundary conditions) as spreading around an origin site up to a certain bond distance. The lattice is
created from a passed `Unitcell` object.

Also, the newly created `Lattice` object can directly be saved. If this has been done before, passing a `load=true`
will allow to load the object instead of creating it again.

Note that this function works for both 2D and 3D unitcells.

# Examples

```julia-repl
julia> lattice = getLatticeByBondDistance(unitcell, 5)
LatticePhysics.Lattice(...)

julia> lattice = getLatticeByBondDistance(unitcell, 8, origin=2)
LatticePhysics.Lattice(...)
```
"""
function getLatticeByBondDistance(unitcell::Unitcell, bonddistance::Int64; origin::Int64=1, load=false, save=false)

    # check how many periodic dimensions the unitcell has
    N_dims = size(unitcell.lattice_vectors, 1)

    # depends on dimensions what to call
    if N_dims == 0
        # not possible
        println("cannot build a lattice without periodic lattice vectors!")
        return
    elseif N_dims == 1
        # not implemented yet
        println("building from bond distance not implemented yet for d=1")
    elseif N_dims == 2
        # just call the fitting routine with the first two entries
        return getLatticeByBondDistance2D(unitcell, bonddistance, origin=origin, load=load, save=save)
    elseif N_dims == 3
        # just call the fitting routine with the first two entries
        return getLatticeByBondDistance3D(unitcell, bonddistance, origin=origin, load=load, save=save)
    else
        # not implemented for any other lattice!
        println("building from bond distance not implemented for d=$(N_dims)")
        return
    end

end

export getLatticeByBondDistance

#-----------------------------------------------------------------------------------------------------------------------------
#
#   BUILDING FUNCTIONS FOR LATTICES THAT ARE BUILD BY BOND DISTANCE FROM AN ORIGIN SITE AND RESIDE IN A GIVEN SHAPE
#   These lattices have of course open boundaries
#
#   Parameters that have to be passed:
#   - Unitcell object from which lattice is build
#   - shape to determine the extent of the lattice from the origin site outwards
#     (a julia function that gives true / false if a site is inside the shape)
#   - the name of the shape for filename purposes
#   - OPTIONAL: origin site index
#
#-----------------------------------------------------------------------------------------------------------------------------

#-----------------------------------------------------------------------------------------------------------------------------
#
#   General building in shape for 2D and 3D and general lattices
#
#-----------------------------------------------------------------------------------------------------------------------------
function getLatticeInShape2D(unitcell::Unitcell, shape::Function, shapename::String="unknown"; origin::Int64=1, load::Bool=false, save::Bool=false)

    # generate the filename of the output
    if contains(unitcell.filename, FOLDER_UNITCELLS)
        filename = replace(unitcell.filename, FOLDER_UNITCELLS, FOLDER_LATTICES)
        filename = replace(filename, ".jld", "_in_shape_$(shapename).jld")
        filename = replace(filename, "_unitcell_", "_lattice_")
    else
        filename = "$(FOLDER_LATTICES)$(split(unitcell.filename, FOLDER_UNITCELLS[end])[end])"
        filename = replace(filename, ".jld", "_in_shape_$(shapename).jld")
        filename = replace(filename, "_unitcell_", "_lattice_")
    end

    # check if can be loaded
    if isfile(filename) && load
        return Lattice(filename)
    end

    # load the data from the unit cell
    uc_basis            = unitcell.basis
    uc_connections      = unitcell.connections
    uc_lattice_vectors  = unitcell.lattice_vectors

    # arrays for new positions and connections
    positions   = []
    connections = Array{Any,1}[]

    # checklist for all sites that are checked if added etc
    checklist   = []
    # format for sites on checklist: ([pos_x, pos_y], [uc_i, uc_j], index_in_UC, bd_current)

    # push the origin site to the checklist
    push!(
        checklist, (uc_basis[origin], [0,0], origin, 0)
    )

    # iterate while the checklist is not empty
    while size(checklist,1) > 0

        # get the item that was on the checklist for longest
        item_to_handle = shift!(checklist)

        # check if the item is already in the positions list
        found = false
        for p in positions
            if item_to_handle[2] == p[2] && item_to_handle[3] == p[3]
                # if yes, continue
                found = true
                break
            end
        end
        if found
            continue
        end

        # if not, push it into
        push!(positions, item_to_handle)
        # for all connections
        index_from = size(positions, 1)

        # insert all connections to sites that are already inside the positions list and all other sites into the checklist
        for c in uc_connections
            # check if correct connections
            if c[1] != item_to_handle[3]
                continue
            end
            # search for the other element of the connection
            i_to = item_to_handle[2][1] + c[4][1]
            j_to = item_to_handle[2][2] + c[4][2]
            a_to = c[2]
            index_to = -1
            for (index,item_handled) in enumerate(positions)
                if item_handled[2] == [i_to, j_to] && item_handled[3] == a_to
                    # make a new connection
                    index_to = index
                    # break the loop
                    break
                end
            end
            # determine whether the element is already inside the list
            if index_to > 0
                # element is at index index_to an can be linked
                connection_new = [index_from; index_to; c[3]; (0, 0)]
                # check if the connetion is already added
                if connection_new in connections
                    continue
                end
                # register as connection
                push!(connections, connection_new)
            else
                # element not in list yet, maybe should be added
                if shape((i_to * uc_lattice_vectors[1]) .+ (j_to * uc_lattice_vectors[2]) .+ uc_basis[c[2]] .- uc_basis[origin])
                    # format for sites on checklist: ([pos_x, pos_y], [uc_i, uc_j], index_in_UC, bd_current)
                    push!(checklist, (
                        (i_to * uc_lattice_vectors[1]) .+ (j_to * uc_lattice_vectors[2]) .+ uc_basis[c[2]],
                        [i_to, j_to],
                        c[2],
                        item_to_handle[4]+1
                    ))
                end
            end
        end

    end

    # change the format of positions
    positions_TMP = positions

    # erase positions
    positions = Array{Float64, 1}[]
    positions_indices = Int64[]

    # insert the real positions
    for p in positions_TMP
        push!(positions, p[1])
        push!(positions_indices, p[3])
    end

    # insert missing connections (if (i to j) is present, insert (j to i))
    for c in connections
        c_proposed = [c[2]; c[1]; c[3]; (0, 0)]
        if !(c_proposed in connections)
            push!(connections, c_proposed)
        else
            #println("connection already there")
        end
    end

    # save everything to a Lattice object
    lattice = Lattice(
            unitcell,
            Int64[],
            Array{Float64, 1}[],
            positions,
            positions_indices,
            connections,
            filename
        )
    # save the lattice object
    if save
        saveLattice(lattice)
    end

    # return the lattice object
    return lattice

end
function getLatticeInShape3D(unitcell::Unitcell, shape::Function, shapename::String="unknown"; origin::Int64=1, load::Bool=false, save::Bool=false)

    # generate the filename of the output
    if contains(unitcell.filename, FOLDER_UNITCELLS)
        filename = replace(unitcell.filename, FOLDER_UNITCELLS, FOLDER_LATTICES)
        filename = replace(filename, ".jld", "_in_shape_$(shapename).jld")
        filename = replace(filename, "_unitcell_", "_lattice_")
    else
        filename = "$(FOLDER_LATTICES)$(split(unitcell.filename, FOLDER_UNITCELLS[end])[end])"
        filename = replace(filename, ".jld", "_in_shape_$(shapename).jld")
        filename = replace(filename, "_unitcell_", "_lattice_")
    end

    # check if can be loaded
    if isfile(filename) && load
        return Lattice(filename)
    end

    # load the data from the unit cell
    uc_basis            = unitcell.basis
    uc_connections      = unitcell.connections
    uc_lattice_vectors  = unitcell.lattice_vectors

    # arrays for new positions and connections
    positions   = []
    connections = Array{Any,1}[]

    # checklist for all sites that are checked if added etc
    checklist   = []
    # format for sites on checklist: ([pos_x, pos_y, pos_z], [uc_i, uc_j, uc_k], index_in_UC, bd_current)

    # push the origin site to the checklist
    push!(
        checklist, (uc_basis[origin], [0,0,0], origin, 0)
    )

    # iterate while the checklist is not empty
    while size(checklist,1) > 0

        # get the item that was on the checklist for longest
        item_to_handle = shift!(checklist)

        # check if the item is already in the positions list
        found = false
        for p in positions
            if item_to_handle[2] == p[2] && item_to_handle[3] == p[3]
                # if yes, continue
                found = true
                break
            end
        end
        if found
            continue
        end

        # if not, push it into
        push!(positions, item_to_handle)
        # for all connections
        index_from = size(positions, 1)

        # insert all connections to sites that are already inside the positions list and all other sites into the checklist
        for c in uc_connections
            # check if correct connections
            if c[1] != item_to_handle[3]
                continue
            end
            # search for the other element of the connection
            i_to = item_to_handle[2][1] + c[4][1]
            j_to = item_to_handle[2][2] + c[4][2]
            k_to = item_to_handle[2][3] + c[4][3]
            a_to = c[2]
            index_to = -1
            for (index,item_handled) in enumerate(positions)
                if item_handled[2] == [i_to, j_to, k_to] && item_handled[3] == a_to
                    # make a new connection
                    index_to = index
                    # break the loop
                    break
                end
            end
            # determine whether the element is already inside the list
            if index_to > 0
                # element is at index index_to an can be linked
                connection_new = [index_from; index_to; c[3]; (0, 0, 0)]
                # check if the connetion is already added
                if connection_new in connections
                    continue
                end
                # register as connection
                push!(connections, connection_new)
            else
                # element not in list yet, maybe should be added
                if shape((i_to * uc_lattice_vectors[1]) .+ (j_to * uc_lattice_vectors[2]) .+ (k_to * uc_lattice_vectors[3]) .+ uc_basis[c[2]] .- uc_basis[origin])
                    # format for sites on checklist: ([pos_x, pos_y, pos_z], [uc_i, uc_j, uc_k], index_in_UC, bd_current)
                    push!(checklist, (
                        (i_to * uc_lattice_vectors[1]) .+ (j_to * uc_lattice_vectors[2]) .+ (k_to * uc_lattice_vectors[3]) .+ uc_basis[c[2]],
                        [i_to, j_to, k_to],
                        c[2],
                        item_to_handle[4]+1
                    ))
                end
            end
        end
    end

    # change the format of positions
    positions_TMP = positions

    # erase positions
    positions = Array{Float64, 1}[]
    positions_indices = Int64[]

    # insert the real positions
    for p in positions_TMP
        push!(positions, p[1])
        push!(positions_indices, p[3])
    end

    # insert missing connections (if (i to j) is present, insert (j to i))
    for c in connections
        c_proposed = [c[2]; c[1]; c[3]; (0, 0, 0)]
        if !(c_proposed in connections)
            push!(connections, c_proposed)
        else
            #println("connection already there")
        end
    end

    # save everything to a Lattice object
    lattice = Lattice(
            unitcell,
            Int64[],
            Array{Float64, 1}[],
            positions,
            positions_indices,
            connections,
            filename
        )
    # save the lattice object
    if save
        saveLattice(lattice)
    end

    # return the lattice object
    return lattice

end

"""
    getLatticeInShape(unitcell::Unitcell, shape::Function, [ shapename::String="unknown" ; origin::Int64=1, load::Bool=false, save::Bool=false])

Function to construct a finite lattice (with open boundary conditions)
as spreading around an origin site up to a certain shape which is given by the passed `Function` `shape`.
The lattice is created from a passed `Unitcell` object.

Also, the newly created `Lattice` object can directly be saved. If this has been done before, passing a `load=true`
will allow to load the object instead of creating it again.

Note that this function works for both 2D and 3D unitcells.

# Examples

```julia-repl
julia> lattice = getLatticeInShape(unitcell, shape_function)
LatticePhysics.Lattice(...)

julia> lattice = getLatticeInShape(unitcell, shape_function, "myshape")
LatticePhysics.Lattice(...)

julia> lattice = getLatticeInShape(unitcell, shape_function, "myshape", origin=2)
LatticePhysics.Lattice(...)
```
"""
function getLatticeInShape(unitcell::Unitcell, shape::Function, shapename::String="unknown"; origin::Int64=1, load::Bool=false, save::Bool=false)

    # check how many periodic dimensions the unitcell has
    N_dims = size(unitcell.lattice_vectors, 1)

    # depends on dimensions what to call
    if N_dims == 0
        # not possible
        println("cannot build a lattice without periodic lattice vectors!")
        return
    elseif N_dims == 1
        # not implemented yet
        println("building in shape not implemented yet for d=1")
    elseif N_dims == 2
        # just call the fitting routine with the first two entries
        return getLatticeInShape2D(unitcell, shape, shapename, origin=origin, load=load, save=save)
    elseif N_dims == 3
        # just call the fitting routine with the first two entries
        return getLatticeInShape3D(unitcell, shape, shapename, origin=origin, load=load, save=save)
    else
        # not implemented for any other lattice!
        println("building in shape not implemented for d=$(N_dims)")
        return
    end
end
export getLatticeInShape

#-----------------------------------------------------------------------------------------------------------------------------
#
#   Special cases for building in shape for 2D and 3D and general lattices
#
#-----------------------------------------------------------------------------------------------------------------------------

"""
    getLatticeInSphere(unitcell::Unitcell, radius::Float64 [; origin::Int64=1, load::Bool=false, save::Bool=false])

Function to construct a finite lattice (with open boundary conditions)
as spreading around an origin site up to a radius (forming a sphere / circle).
The lattice is created from a passed `Unitcell` object.

Also, the newly created `Lattice` object can directly be saved. If this has been done before, passing a `load=true`
will allow to load the object instead of creating it again.

Note that this function works for both 2D and 3D unitcells.

# Examples

```julia-repl
julia> lattice = getLatticeInSphere(unitcell, 3.0)
LatticePhysics.Lattice(...)

julia> lattice = getLatticeInSphere(unitcell, 3.0, origin=2)
LatticePhysics.Lattice(...)
```
"""
function getLatticeInSphere(unitcell::Unitcell, radius::Float64; origin::Int64=1, load::Bool=false, save::Bool=false)

    # check how many periodic dimensions the unitcell has
    N_dims = size(unitcell.lattice_vectors, 1)

    # depends on dimensions what to call
    if N_dims == 0
        # not possible
        println("cannot build a lattice without periodic lattice vectors!")
        return
    elseif N_dims == 1
        # not implemented yet
        println("building in a sphere not implemented yet for d=1")
        return
    elseif N_dims == 2
        # determine the shape name
        shapename = "sphere_$(radius)"
        # determine the shape function
        shape_circle(point) = sum(point.*point) < radius*radius
        # call the general shape function
        return getLatticeInShape2D(unitcell, shape_circle, shapename, origin=origin, load=load, save=save)
    elseif N_dims == 3
        # determine the shape name
        shapename = "sphere_$(radius)"
        # determine the shape function
        shape_sphere(point) = sum(point.*point) < radius*radius
        # call the general shape function
        return getLatticeInShape3D(unitcell, shape_sphere, shapename, origin=origin, load=load, save=save)
    else
        # not implemented for any other lattice!
        println("building in a sphere not implemented for d=$(N_dims)")
        return
    end
end
export getLatticeInSphere

"""
    getLatticeInBox(unitcell::Unitcell, extent_array::Array{Float64} [; origin::Int64=1, load::Bool=false, save::Bool=false])

Function to construct a finite lattice (with open boundary conditions)
as spreading around an origin site up to the boundaries of a given box / rectangle.
The lattice is created from a passed `Unitcell` object.

The box dimensions are given by the passed `extent_array` which contains the length in all elementary
directions as `Float64` numbers. The origin site is place exactly at the center of the box.

The newly created `Lattice` object can directly be saved. If this has been done before, passing a `load=true`
will allow to load the object instead of creating it again.

Note that this function works for both 2D and 3D unitcells.

# Examples

```julia-repl
julia> lattice = getLatticeInSphere(unitcell, 3.0)
LatticePhysics.Lattice(...)

julia> lattice = getLatticeInSphere(unitcell, 3.0, origin=2)
LatticePhysics.Lattice(...)
```
"""
function getLatticeInBox(unitcell::Unitcell, extent_array::Array{Float64}; origin::Int64=1, load=false, save=true)

    # check how many periodic dimensions the unitcell has
    N_dims = size(unitcell.lattice_vectors, 1)

    # depends on dimensions what to call
    if N_dims == 0
        # not possible
        println("cannot build a lattice without periodic lattice vectors!")
        return
    elseif N_dims == 1
        # not implemented yet
        println("building in box not implemented yet for d=1")
        return
    elseif N_dims == 2
        # get the dimensions
        length_x = extent_array[1]
        length_y = extent_array[2]
        # determine the shape name
        shapename = "box_$(length_x)x$(length_y)_around_$(origin)"
        # determine the shape function
        shape_box2d(point) = (abs(point[1])<=length_x/2.0) && (abs(point[2])<=length_y/2.0)
        # call the general shape function
        return getLatticeInShape2D(unitcell, shape_box2d, shapename, origin=origin, load=load, save=save)
    elseif N_dims == 3
        # get the dimensions
        length_x = extent_array[1]
        length_y = extent_array[2]
        length_z = extent_array[3]
        # determine the shape name
        shapename = "box_$(length_x)x$(length_y)x$(length_z)_around_$(origin)"
        # determine the shape function
        shape_box3d(point) = (abs(point[1])<=length_x/2.0) && (abs(point[2])<=length_y/2.0) && (abs(point[3])<=length_z/2.0)
        # call the general shape function
        return getLatticeInShape3D(unitcell, shape_box3d, shapename, origin=origin, load=load, save=save)
    else
        # not implemented for any other lattice!
        println("building in box not implemented for d=$(N_dims)")
        return
    end
end
export getLatticeInBox


text:
interval_markers(x::AbstractInterval{T,L,R}) where {T,L,R} = interval_markers(L, R)
interval_markers(L::Type{<:Bound}, R::Type{<:Bound}) = interval_marker.([L, R])
interval_marker(x::Type{Closed}) = :vline
interval_marker(x::Type{Open}) = :none

# TODO: Add support for plotting unbounded intervals
@recipe function f(xs::AbstractVector{<:AbstractInterval{T,L,R}}, ys) where {T, L <: Bounded, R <: Bounded}
    new_xs = T[]
    new_ys = []
    markers = Symbol[]
    for (x, y) in zip(xs, ys)
        # To cause line to not be connected, need to add a breaker point with
        # NaN in one of the coordinates. We put that in `new_ys` as that is probably a
        # float. where-as new_xs is potentially a DateTime etc that would not accept NaN
        append!(new_xs, [first(x), last(x), last(x)])
        append!(new_ys, [y, y, NaN])
        append!(markers, interval_markers(x))
        push!(markers, :none)
    end

    # Work around GR bug that shows :none as a marker
    # TODO: remove once https://github.com/jheinen/GR.jl/issues/295  is fixed
    markeralpha := [x == :none ? 0 : 1 for x in markers]

    markershape := markers
    seriestype  := :path  # always a path, even in a scatter plot
    new_xs, new_ys
end


text:
#=
The arithmetic sequence, 1487, 4817, 8147, in which each of the terms increases by 3330, is unusual in two ways: (i) each of the three terms are prime, and, (ii) each of the 4-digit numbers are permutations of one another.

There are no arithmetic sequences made up of three 1-, 2-, or 3-digit primes, exhibiting this property, but there is one other 4-digit increasing sequence.

What 12-digit number do you form by concatenating the three terms in this sequence?
=#
include("projecteulerutils.jl")
using Primes

function calc()
  b = c = 0
  for a in 1000:10000-3330-3330
    if a == 1487
      continue
    end
    b = a+3330
    c = b+3330
    if isprime(a) && isprime(b) && isprime(c) && ispermutation(digits(a),digits(b)) && ispermutation(digits(b), digits(c))
      return "$a$b$c"
    end
  end
end
@time println(calc())

max_stars_repo_name: daniel-beard/JuliaProjectEuler

max_issues_repo_name: daniel-beard/JuliaProjectEuler

max_forks_repo_name: daniel-beard/JuliaProjectEuler


text:
#
# This file is a part of MolecularGraph.jl
# Licensed under the MIT License http://opensource.org/licenses/MIT
#

export
    clear,
    fetchfile,
    fetchresource,
    pubchemsdf,
    molsfortest,
    resource_dir,
    pubchem_dir

resource_dir = joinpath(dirname(@__FILE__), "..", "_resources")
pubchem_dir = joinpath(resource_dir, "PubChem")

function initialize_resource_dir()
    if !isdir(resource_dir)
        mkdir(resource_dir)
        println("Resource directory created: $(resource_dir)")
    end
end

function clear()
    initialize_resource_dir()
    for p in readdir(resource_dir)
        rm(p, recursive=true)
    end
    println("Resource directory is now empty: $(resource_dir)")
end

function fetchfile(url, dest)
    if isfile(dest)
        println("file: $(basename(dest)) already exists")
        return
    end
    println("Downloading: $(url)")
    download(url, dest)
end

function fetchresource(url, filename)
    initialize_resource_dir()
    dest = joinpath(resource_dir, filename)
    fetchfile(url, dest)
end

function pubchemsdf(cid, name)
    if !isdir(pubchem_dir)
        initialize_resource_dir()
        mkdir(pubchem_dir)
        println("PubChem data directory created: $(pubchem_dir)")
    end
    dest = joinpath(resource_dir, "PubChem", "$(name).mol")
    if isfile(dest)
        println("file: $(name).mol already exists")
        return
    end
    url = "https://pubchem.ncbi.nlm.nih.gov/rest/pug/compound/cid/$(cid)/SDF"
    println("Downloading: $(url)")
    download(url, dest)
end

function molsfortest()
    pubchemfilelist = [
        ("Buckminsterfullerene", 123591)
    ]
    testdir = joinpath(resource_dir, "testmols")
    if !isdir(testdir)
        initialize_resource_dir()
        mkdir(testdir)
        println("Test data directory created: $(testdir)")
    end
    for (name, id) in pubchemfilelist
        dest = joinpath(testdir, "$(name).sdf")
        if isfile(dest)
            println("file: $(name).mol already exists")
            return
        end
        url = "https://pubchem.ncbi.nlm.nih.gov/rest/pug/compound/cid/$(id)/SDF"
        println("Downloading: $(url)")
        download(url, dest)
    end
    println("Finished: $(testdir)")
end