diff --git a/Project.toml b/Project.toml index 07f7365..3ab796b 100644 --- a/Project.toml +++ b/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"] diff --git a/examples/matter_oscillogram.jl b/examples/matter_oscillogram.jl index d0f41a0..e90e5a2 100644 --- a/examples/matter_oscillogram.jl +++ b/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 \ No newline at end of file diff --git a/test/examples_tests.jl b/test/examples_tests.jl new file mode 100644 index 0000000..3c33a49 --- /dev/null +++ b/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() \ No newline at end of file diff --git a/test/runtests.jl b/test/runtests.jl index 11b747b..1d93eff 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -15,3 +15,7 @@ end @testset "Matter" begin include("matter_tests.jl") end + +@testset "Examples" begin +include("examples_tests.jl") +end \ No newline at end of file