diff --git a/.gitignore b/.gitignore
index abdcd06..e144696 100644
--- a/.gitignore
+++ b/.gitignore
@@ -12,3 +12,4 @@ _research
 *.synctex.gz(buzy)
 *.out
 ~$*
+\build
\ No newline at end of file
diff --git a/Project.toml b/Project.toml
index 07f7365..76a5d12 100644
--- a/Project.toml
+++ b/Project.toml
@@ -27,4 +27,4 @@ HDF5 = "f67ccb44-e63f-5c2f-98bd-6dc0ccc4ba2f"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Test", "HDF5"]
+test = ["Test", "HDF5"]
\ No newline at end of file
diff --git a/docs/Manifest.toml b/docs/Manifest.toml
deleted file mode 100644
index 0e4b168..0000000
--- a/docs/Manifest.toml
+++ /dev/null
@@ -1,92 +0,0 @@
-# This file is machine-generated - editing it directly is not advised
-
-[[Base64]]
-uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"
-
-[[Dates]]
-deps = ["Printf"]
-uuid = "ade2ca70-3891-5945-98fb-dc099432e06a"
-
-[[Distributed]]
-deps = ["Random", "Serialization", "Sockets"]
-uuid = "8ba89e20-285c-5b6f-9357-94700520ee1b"
-
-[[DocStringExtensions]]
-deps = ["LibGit2", "Markdown", "Pkg", "Test"]
-git-tree-sha1 = "88bb0edb352b16608036faadcc071adda068582a"
-uuid = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
-version = "0.8.1"
-
-[[Documenter]]
-deps = ["Base64", "Dates", "DocStringExtensions", "InteractiveUtils", "JSON", "LibGit2", "Logging", "Markdown", "REPL", "Test", "Unicode"]
-git-tree-sha1 = "0be9bf63e854a2408c2ecd3c600d68d4d87d8a73"
-uuid = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
-version = "0.24.2"
-
-[[InteractiveUtils]]
-deps = ["Markdown"]
-uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
-
-[[JSON]]
-deps = ["Dates", "Mmap", "Parsers", "Unicode"]
-git-tree-sha1 = "b34d7cef7b337321e97d22242c3c2b91f476748e"
-uuid = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
-version = "0.21.0"
-
-[[LibGit2]]
-uuid = "76f85450-5226-5b5a-8eaa-529ad045b433"
-
-[[Libdl]]
-uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb"
-
-[[Logging]]
-uuid = "56ddb016-857b-54e1-b83d-db4d58db5568"
-
-[[Markdown]]
-deps = ["Base64"]
-uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"
-
-[[Mmap]]
-uuid = "a63ad114-7e13-5084-954f-fe012c677804"
-
-[[Parsers]]
-deps = ["Dates", "Test"]
-git-tree-sha1 = "0139ba59ce9bc680e2925aec5b7db79065d60556"
-uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0"
-version = "0.3.10"
-
-[[Pkg]]
-deps = ["Dates", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "UUIDs"]
-uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
-
-[[Printf]]
-deps = ["Unicode"]
-uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7"
-
-[[REPL]]
-deps = ["InteractiveUtils", "Markdown", "Sockets"]
-uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb"
-
-[[Random]]
-deps = ["Serialization"]
-uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
-
-[[SHA]]
-uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce"
-
-[[Serialization]]
-uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"
-
-[[Sockets]]
-uuid = "6462fe0b-24de-5631-8697-dd941f90decc"
-
-[[Test]]
-deps = ["Distributed", "InteractiveUtils", "Logging", "Random"]
-uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
-
-[[UUIDs]]
-deps = ["Random", "SHA"]
-uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"
-
-[[Unicode]]
-uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5"
diff --git a/docs/Project.toml b/docs/Project.toml
index dfa65cd..0068ca7 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -1,2 +1,6 @@
 [deps]
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
+LaTeXStrings = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f"
+Neurthino = "ea3506b1-7c4e-45d8-9afb-6e8d18a7dde4"
+Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
+Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
diff --git a/docs/make.jl b/docs/make.jl
index b1ebbe1..9ba9fa6 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -9,6 +9,7 @@ makedocs(;
     ),
     pages=[
         "Introduction" => "index.md",
+        "Basic usage" => "basic_usage.md",
         "API" => "api.md",
     ],
     repo="https://github.com/KM3NeT/Neurthino.jl/blob/{commit}{path}#L{line}",
diff --git a/docs/src/basic_usage.md b/docs/src/basic_usage.md
new file mode 100644
index 0000000..394df0c
--- /dev/null
+++ b/docs/src/basic_usage.md
@@ -0,0 +1,142 @@
+# Basic usage
+
+## Setting up the Oscillation Parameters
+We start by creating an `OscillationParameters` struct in vacuum with 3 neurtrino flavours:
+```@example 1
+using Neurthino
+
+osc = OscillationParameters(3);
+nothing # hide
+```
+
+The mixing angles, CP-phase and mass squared differences are all initialised to be zero. We can set them individually by:
+```@example 1
+mixingangle!(osc, 1 => 2, 0.58);
+mixingangle!(osc, 1 => 3, 0.15);
+mixingangle!(osc, 2 => 3, 0.737);
+
+cpphase!(osc, 1 => 3, 5.34);
+
+masssquareddiff!(osc, 2 => 3, -2.457e-3);
+masssquareddiff!(osc, 1 => 2, -7.5e-5);
+nothing # hide
+```
+One could also use `setθ!`, `setδ!` and `setΔm²!`.
+
+We can then use `oscprob`, or similar `Pνν`, to compute the oscillation probabilities for all 9 flavour changes. For an single energy of 1.5 GeV and baseline of 10000 km (we could also pass in a vector), we get:
+```@example 1
+prob = oscprob(osc, 1.5, 10000)
+```
+The output is an `AxisArray` which provides intuitive indexing, e.g.
+for `prob(νμ→ντ)` at the given energy and baseline:
+```@example 1
+prob[Energy=1, Baseline=1, InitFlav=Muon, FinalFlav=Tau]
+```
+The probabilities are calculated based on the transition matrix
+(the so-called PMNS-Matrix) between flavour and mass eigenstates,
+as well as the Hamiltonian in the mass eigenbasis. We can compute PMNS matrix:
+```@example 1
+U = PMNSMatrix(osc)
+```
+and the Hamiltonian:
+```@example 1
+H = Hamiltonian(osc)
+```
+once, and avoid repeated calculations during a loop. Then, we can use `oscprob(U, H, 1, 10000)` to get the same results as above.
+
+!!! note "Oscillations in vacuum"
+    Up to this point, the `osc`, `U` and `H` are all defined in vacuum. Using multiple dispatch, the same function `oscprob` will be used to get the oscillation probabilities along a path in matter.
+
+## Neutrino oscillations in water
+`Neurthino.jl` comes with the PREM model, a layered Earth model, which allows us to create a neutrino path based on just the depth and zenith angle. Instead of moving through a vacuum, we want to model a Water-Cherenkov detector in sea or ice. For this example, we place the detector at 2km depth:
+```@example 1
+waterdepth_km = 2;
+nothing # hide
+```
+The zenith angles for up going neutrinos  and
+subsequently the neutrino paths are generated first:
+
+We first set the desired zenith angles (cos(θ)ϵ[-1,0]) for up going neutrinos and with an energy range of [0, 2] GeV:
+```@example 1
+coszenith = range(-1; stop=0, length=100);
+zenith = acos.(coszenith);
+
+energies = 10 .^(range(0;stop=2,length=100));
+nothing # hide
+```
+
+
+With the zenith angles defined, we will create a loop that computes the oscillation probabilities for each angle.
+
+First we pre-allocate the data array
+```@example 1
+data_points = zeros((3,3,length(zenith), length(energies)));
+nothing # hide
+```
+
+We can pass the the matrices `U` and `H` directly to avoid repeating their computation. Although `U` and `H` were defined in vacuum, the call to `oscprob` will automatically convert them using the densities defined by `path`. Let's start the loop over zeniths:
+```@example 1
+for (i,z) in enumerate(zenith)
+
+    # create a path with zenith z and waterdepth_km
+    path = Neurthino.prempath(z, waterdepth_km, samples=50);
+
+    # compute the probabilities for all flavours
+    prob = oscprob(U, H, energies, path);
+
+    # save the results in data array
+    for (k,l) in Iterators.product(fill(1:3, 2)...)
+        data_points[k,l,i,:] = prob[:,1,k,l];
+    end
+end
+```
+Here, we used `Neurthino.prempath` to create a `path`, which defines the depths and the densities according to the PREM Earth model. The layered model is defined by `Neurthino.PREM`. The first 3km is water.
+
+Then we call `oscprob` to compute the oscillation probabilities of all neutrinos on that path. Note that we can use both `oscprob(osc, energies, path)` and `oscprob(U, H, energies, path)` (latter is prefered during loops to avoid recalculating `U` and `H` in each iteration). Lastly, we save the probabilities into our pre-allocated array.
+
+Now we have our data, let's visualize the oscillogram! First we import the packages we will use:
+```@example 1
+using Unitful
+using Plots
+using LaTeXStrings
+```
+We again loop over all flavours and create an heatmap for each combination:
+```@example 1
+for (k, startFlav) in enumerate([L"\nu_e", L"\nu_\mu", L"\nu_\tau"])
+    for (l, endFlav) in enumerate([L"\nu_e", L"\nu_\mu", L"\nu_\tau"])
+
+        # create an oscillogram for every start- and endFlav combination
+        p = heatmap(energies*u"GeV", coszenith, data_points[k,l,:,:], xscale=:log10,
+            xlabel="Energy", ylabel=L"\cos(\theta_z)", title=latexstring(startFlav, L"\rightarrow", endFlav),
+            clim=(0,1));
+
+        # save each plot
+        filename = "water_nu_$(k)_nu_$(l).svg"
+        savefig(p, filename)
+
+    end
+    # instead of showing the output for each flavour, we break here
+    break
+end
+```
+
+For the electron neutrino oscillations, we then get:
+![water_oscillogram](water_nu_1_nu_1.svg)
+
+![water_oscillogram](water_nu_1_nu_2.svg)
+
+![water_oscillogram](water_nu_1_nu_3.svg)
+
+
+## Vectorize operations
+The for-loops can be avoided by vectorizing the calculations. We first define all the paths based on all zenith angles and the waterdepth:
+```@example 1
+paths = Neurthino.prempath(zenith, waterdepth_km, samples=100);
+nothing # hide
+```
+
+Then we compute all probabilites at once:
+```@example 1
+prob = oscprob(U, H, energies, paths)
+```
+The returned array `prob` is again of type `AxisArray` with an axis `:Energy` and now `:Path` for the path index (instead of the `:Baseline` axis).
\ No newline at end of file
diff --git a/examples/matter_oscillogram.jl b/examples/matter_oscillogram.jl
deleted file mode 100644
index d0f41a0..0000000
--- a/examples/matter_oscillogram.jl
+++ /dev/null
@@ -1,59 +0,0 @@
-#!/usr/bin/env julia
-
-println("Loading libraries")
-using Neurthino
-using Unitful
-using Plots
-using ProgressMeter
-
-
-function main()
-    println("Initialising oscillation parameters")
-    osc = OscillationParameters(3);
-    mixingangle!(osc, 1 => 2, 0.58);
-    mixingangle!(osc, 1 => 3, 0.15);
-    mixingangle!(osc, 2 => 3, 0.737);
-    cpphase!(osc, 1 => 3, 5.34);
-    masssquareddiff!(osc, 2 => 3, -2.457e-3);
-    masssquareddiff!(osc, 1 => 3, -2.457e-3-7.5e-5);
-    masssquareddiff!(osc, 1 => 2, -7.5e-5);
-
-    U = PMNSMatrix(osc)
-    H = Hamiltonian(osc)
-
-    coszenith = range(-1; stop=0, length=100)
-    zenith = acos.(coszenith)
-    energies = 10 .^(range(0;stop=2,length=100))
-    
-    data_points = zeros((3,3,length(zenith), length(energies)))
-
-
-    @showprogress 1 "Calculating transition probabilities for different zenith angles" for (i,z) in enumerate(zenith)
-        sec, dens = Neurthino.prempath(z, 2.5, samples=50)
-        P = Neurthino.mattertransprob(osc, energies, dens, sec)
-        for (k,l) in Iterators.product(fill(1:3, 2)...)
-            data_points[k,l,i,:] = map(x->x[k,l], P)
-        end
-    end
-    
-    plot_arr = Array{Plots.Plot{Plots.GRBackend},1}()
-    for (k,l) in Iterators.product(fill(1:3, 2)...)
-        p = heatmap(data_points[k,l,:,:],  yticks=(1:20:100, coszenith[1:10:100]), xticks=(1:40:200, energies[1:10:50]), clim=(0,1))
-        push!(plot_arr, p)
-    end
-    files = Vector{AbstractString}()
-    @showprogress 1 "Creating plots" for (idx, p) in enumerate(plot_arr)
-        f = "$idx.pdf"
-        push!(files, f)
-        savefig(p, f)
-    end
-end
-
-main()
-
-
-
-
-
-
-
diff --git a/test/runtests.jl b/test/runtests.jl
index 11b747b..ef05b88 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -14,4 +14,4 @@ end
 
 @testset "Matter" begin
 include("matter_tests.jl")
-end
+end
\ No newline at end of file