Fix KM3NeT#12

Project.toml

@@ -25,6 +25,9 @@ julia = "^1.3"
 [extras]
 HDF5 = "f67ccb44-e63f-5c2f-98bd-6dc0ccc4ba2f"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
+Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
+LaTeXStrings = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f"
 
 [targets]
-test = ["Test", "HDF5"]
+test = ["Test", "HDF5", "Unitful", "Plots", "LaTeXStrings"]

examples/matter_oscillogram.jl

@@ -1,59 +1,64 @@
 #!/usr/bin/env julia
 
-println("Loading libraries")
 using Neurthino
 using Unitful
 using Plots
-using ProgressMeter
+using LaTeXStrings
 
+function example_matter_oscillogram()
+    # run the example_matter_oscillogram
 
-function main()
-    println("Initialising oscillation parameters")
+    # define oscillation parameters with 3 neutrino flavours
     osc = OscillationParameters(3);
+
+    # define values of the mixing angles, cpp phase and mass square differences
     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);
 
+    # compute the PMNS matrix and Hamiltonian before starting the loop
     U = PMNSMatrix(osc)
     H = Hamiltonian(osc)
 
+    # the zenith and energy domain to compute
     coszenith = range(-1; stop=0, length=100)
     zenith = acos.(coszenith)
     energies = 10 .^(range(0;stop=2,length=100))
-
+
+    # initialise data array
     data_points = zeros((3,3,length(zenith), length(energies)))
 
+    # loop over zenith directions
+    for (i,z) in enumerate(zenith)
+
+        # create a path at with zenith z and depth 2.5 km
+        path = Neurthino.prempath(z, 2.5, samples=50)
+
+        # compute the probabilities
+        prob = Neurthino.oscprob(osc, energies, path)
 
-    @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)
+        # save the results in data array
         for (k,l) in Iterators.product(fill(1:3, 2)...)
-            data_points[k,l,i,:] = map(x->x[k,l], P)
+            data_points[k,l,i,:] = prob[:,1,k,l]
         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()
-
-
-
-
 
+    # loop over the start and end flavours
+    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
+            # f = "nu_$(k)_nu_$(l).pdf"
+            # savefig(p, f)
+            # which is skipped for testing purposes
+        end
+    end
+end

test/examples_tests.jl

@@ -0,0 +1,9 @@
+# check if examples run without error or warning
+
+using Neurthino
+using Unitful
+using Plots
+using LaTeXStrings
+
+include("../examples/matter_oscillogram.jl")
+@test_nowarn example_matter_oscillogram()

test/runtests.jl

@@ -15,3 +15,7 @@ end
 @testset "Matter" begin
 include("matter_tests.jl")
 end
+
+@testset "Examples" begin
+include("examples_tests.jl")
+end