Skip to content

Commit

Permalink
add height of max ionization to ionization plot
Browse files Browse the repository at this point in the history
  • Loading branch information
@egavazzi
egavazzi committed Jan 25, 2025
1 parent 092ddb1 commit 7b8cffd
Showing 1 changed file with 85 additions and 77 deletions.
162 changes: 85 additions & 77 deletions scripts/plotting/Plot_ionization_rate.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -10,7 +10,7 @@ set_theme!(fontsize = 20)
## Directory to plot, absolute path
full_path_to_directory = joinpath(REVONTULI_MOUNT,
"mnt/data/etienne/Julia/AURORA.jl/data/Visions2/" *
"InvertedV_480s_fixed-secondaries")
"Alfven_train_474s")

# Read the ionization data
ionization_file = joinpath(full_path_to_directory, "Qzt_all_L.mat")
Expand All @@ -19,10 +19,88 @@ QO2i = data["QO2i"] .* 1e4 # Q was calculated with Ie in (/cm²) when it sould h
QOi = data["QOi"] .* 1e4
QN2i = data["QN2i"] .* 1e4
Qtot = QN2i + QO2i + QOi
h_atm = vec(data["h_atm"])
h_atm = vec(data["h_atm"]) ./ 1e3
t = vec(data["t"])



## Plot ionization rate as heatmap + top incoming flux
#region
# Read the Ie_top flux
Ietop_file = find_Ietop_file(full_path_to_directory)
data = matread(Ietop_file)
Ietop = data["Ie_total"]
θ_Visions_lims = data["theta_lims"]
t_top = data["t_top"]; t_top = [t_top; t_top[end] + diff(t_top)[end]] .- t_top[1]
Ie_file = joinpath(full_path_to_directory, "IeFlickering-01.mat")
Egrid = matopen(Ie_file) do file; read(file, "E"); end # extract Egrid from the simulation
dEgrid = diff(Egrid); dEgrid = [dEgrid; dEgrid[end] + diff(dEgrid)[end]] # calculate dE
# Plot e- precipitation
fig = Figure(size = (850, 950))
angle_cone = [90 180] # angle for the cone of precipitation to plot, 180° is field-aligned down
ax_Ietop = Axis(fig[1, 1], yscale = log10, ylabel = " Energy (eV)",
title = string(180 - angle_cone[2], " - ", 180 - angle_cone[1], "° DOWN"),
yminorticksvisible = true, yminorticks = IntervalsBetween(9),
xticklabelsvisible = false, xminorticksvisible = true,
xticksmirrored = true, yticksmirrored = true,
limits = ((0, 1), nothing))
idx_θ = vec(angle_cone[1] .<= abs.(acosd.(mu_avg(θ_Visions_lims))) .<= angle_cone[2])
BeamW = beam_weight([angle_cone[1], angle_cone[2]])
data_heatmap = dropdims(sum(Ietop[idx_θ, :, :]; dims=1); dims=1) ./ BeamW ./ dEgrid' .* Egrid'
hm = Makie.heatmap!(t_top, Egrid, data_heatmap; colorrange = (1e6, maximum(data_heatmap)),
colorscale = log10, colormap = :inferno)
custom_formatter(values) = map(v -> rich("10", superscript("$(round(Int64, v))")), values)
Colorbar(fig[1, 2], hm; label = "IeE (eV/m²/s/eV/ster)")
# Plot ionization rate as heatmap
ax_Q = Axis(fig[2, 1], xlabel = "t (s)", ylabel = "altitude (km)", aspect = AxisAspect(1),
xminorticksvisible = true, yminorticksvisible = true, xticksmirrored = true,
yticksmirrored = true)
hm = heatmap!(t, h_atm, Qtot'; colorrange = (1e6, maximum(Qtot)), colorscale=log10, rasterize = true)
cb = Colorbar(fig[2, 2], hm; label = "Ionization rate (/m³/s)")
Qtot_contour = log10.(Qtot) |> (x -> replace(x, -Inf => 0))
ct = contour!(t, h_atm, Qtot_contour'; levels = 6:10, color = :black, labels = true)
colsize!(fig.layout, 1, Aspect(2, 1.0))
rowsize!(fig.layout, 1, Relative(1/5))
linkxaxes!(ax_Ietop, ax_Q)
xlims!(ax_Q, 0, 3)
ylims!(ax_Q, 100, 500)
# Add line to indicate height of max ionization
height_Qmax = [h_atm[i_max[1]] for i_max in vec(findmax(Qtot, dims=1)[2])]
# height_Qmax[vec(maximum(Qtot, dims=1)) .< 1e6] .= NaN # remove when values of Q are very low
height_Qmax[vec(maximum(Qtot, dims=1)) .< maximum(Qtot) / 10] .= NaN # remove when values of Q are less than 10% of the max
lines!(ax_Q, t, height_Qmax; color = :red, linestyle = :dash, linewidth = 2)
# display(fig, backend = GLMakie)
display(fig)

## Save the figure
savefile = joinpath(full_path_to_directory, "Ionization_rate_zt.png")
save(savefile, fig; backend = CairoMakie)
println("Saved $savefile")
savefile = joinpath(full_path_to_directory, "Ionization_rate_zt.svg")
save(savefile, fig; backend = CairoMakie)
println("Saved $savefile")
savefile = joinpath(full_path_to_directory, "Ionization_rate_zt.pdf")
save(savefile, fig; backend = CairoMakie)
println("Saved $savefile")
savefile = joinpath(full_path_to_directory, "Ionization_rate_zt.eps")
save(savefile, fig; backend = CairoMakie)
println("Saved $savefile")
#endregion










## ====================================================================================== ##
## ====================================================================================== ##



## Plot ionization rate as an animated lineplot
#region
fig = Figure(size=(800, 1200))
Expand All @@ -44,11 +122,11 @@ ax1 = Axis(fig[1, 1], title = time,
xscale = log10, xticks = LogTicks(2:15))
Q_max = maximum(maximum.([QO2i, QOi, QN2i]))
xlims!(ax1, 1e3, Q_max * 10)
ylims!(h_atm[1] / 1e3 - 30, h_atm[end] / 1e3 + 30)
l_QO2i = lines!(QO2i_slice, h_atm / 1e3; linestyle =:dash)
l_QOi = lines!(QOi_slice, h_atm / 1e3; linestyle =:dash)
l_QN2i = lines!(QN2i_slice, h_atm / 1e3; linestyle =:dash)
l_tot = lines!(Qtot_slice, h_atm / 1e3; color = :black)
ylims!(h_atm[1] - 30, h_atm[end] + 30)
l_QO2i = lines!(QO2i_slice, h_atm; linestyle =:dash)
l_QOi = lines!(QOi_slice, h_atm; linestyle =:dash)
l_QN2i = lines!(QN2i_slice, h_atm; linestyle =:dash)
l_tot = lines!(Qtot_slice, h_atm; color = :black)
Legend(fig[1, 2], [l_QO2i, l_QOi, l_QN2i, l_tot], ["O2", "O", "N2", "tot"],
"Ionization rates")
# make button to run/stop the animation
Expand Down Expand Up @@ -91,73 +169,3 @@ record(fig, video_file, 1:length(t); framerate = 60, px_per_unit = 2) do i_t
end
println("Saved $video_file")
#endregion





## Plot ionization rate as heatmap + top incoming flux
#region
# File with Ie_top flux datafilename
# Ietop_file = joinpath(full_path_to_directory, "Ie_incoming_475s.mat")
incoming_files = filter(file -> startswith(file, "Ie_incoming_"), readdir(full_path_to_directory))
if length(incoming_files) > 1
error("More than one file contains incoming flux. This is not normal")
else
global Ietop_file = joinpath(full_path_to_directory, incoming_files[1])
end

# Read the Ie_top flux
data = matread(Ietop_file)
Ietop = data["Ie_total"]
θ_Visions_lims = data["theta_lims"]
t_top = data["t_top"]; t_top = [t_top; t_top[end] + diff(t_top)[end]] .- t_top[1]
Ie_file = joinpath(full_path_to_directory, "IeFlickering-01.mat")
Egrid = matopen(Ie_file) do file; read(file, "E"); end # extract Egrid from the simulation
dEgrid = diff(Egrid); dEgrid = [dEgrid; dEgrid[end] + diff(dEgrid)[end]] # calculate dE
# Plot e- precipitation
fig = Figure(size = (850, 950))
angle_cone = [90 180] # angle for the cone of precipitation to plot, 180° is field-aligned down
ax_Ietop = Axis(fig[1, 1], yscale = log10, ylabel = " Energy (eV)",
title = string(180 - angle_cone[2], " - ", 180 - angle_cone[1], "° DOWN"),
yminorticksvisible = true, yminorticks = IntervalsBetween(9),
xticklabelsvisible = false, xminorticksvisible = true,
xticksmirrored = true, yticksmirrored = true,
limits = ((0, 1), nothing))
idx_θ = vec(angle_cone[1] .<= abs.(acosd.(mu_avg(θ_Visions_lims))) .<= angle_cone[2])
BeamW = beam_weight([angle_cone[1], angle_cone[2]])
data_heatmap = dropdims(sum(Ietop[idx_θ, :, :]; dims=1); dims=1) ./ BeamW ./ dEgrid' .* Egrid'
hm = Makie.heatmap!(t_top, Egrid, data_heatmap; colorrange = (1e6, maximum(data_heatmap)),
colorscale = log10, colormap = :inferno)
custom_formatter(values) = map(v -> rich("10", superscript("$(round(Int64, v))")), values)
Colorbar(fig[1, 2], hm; label = "IeE (eV/m²/s/eV/ster)")
# Plot ionization rate as heatmap
ax_Q = Axis(fig[2, 1], xlabel = "t (s)", ylabel = "altitude (km)", aspect = AxisAspect(1),
xminorticksvisible = true, yminorticksvisible = true, xticksmirrored = true,
yticksmirrored = true)
hm = heatmap!(t, h_atm / 1e3, Qtot'; colorrange = (1e6, maximum(Qtot)), colorscale=log10, rasterize = true)
cb = Colorbar(fig[2, 2], hm; label = "Ionization rate (/m³/s)")
Qtot_contour = log10.(Qtot) |> (x -> replace(x, -Inf => 0))
ct = contour!(t, h_atm / 1e3, Qtot_contour'; levels = 6:10, color = :black, labels = true)
colsize!(fig.layout, 1, Aspect(2, 1.0))
rowsize!(fig.layout, 1, Relative(1/5))
linkxaxes!(ax_Ietop, ax_Q)
xlims!(ax_Q, 0, 1)
ylims!(ax_Q, 100, 500)
# display(fig, backend = GLMakie)
display(fig)

## Save the figure
savefile = joinpath(full_path_to_directory, "Ionization_rate_zt.png")
save(savefile, fig; backend = CairoMakie)
println("Saved $savefile")
savefile = joinpath(full_path_to_directory, "Ionization_rate_zt.svg")
save(savefile, fig; backend = CairoMakie)
println("Saved $savefile")
savefile = joinpath(full_path_to_directory, "Ionization_rate_zt.pdf")
save(savefile, fig; backend = CairoMakie)
println("Saved $savefile")
savefile = joinpath(full_path_to_directory, "Ionization_rate_zt.eps")
save(savefile, fig; backend = CairoMakie)
println("Saved $savefile")
#endregion

0 comments on commit 7b8cffd

Please sign in to comment.