#Make comparison plot with theoretical model
function make_energy_plot(  p::Param ; 
                            theory_comparison = :true,
                            gamma_comparison  = :false,
                            gamma_list = nothing, 
                            gamma_legend = nothing, 
                            mscale = log10 ,
                            oname = "default_name",
                            xlims = nothing,
                            ylims = nothing ,
                            energy_data_every = 10)

    #Welcome information
    println("\n \n")        
    println("___________________________________________________________")        
    println(" ")
    println(" Making Energy Plot from file: /data" , p.folder  ,p.fname)
    println(" ")

    #Default values
    my_xlims = xlims
    my_ylims = ylims

    if( xlims === nothing )
        my_xlims = [0 , p.tmax]
    end
    if( ylims === nothing )
        my_ylims = [1E-4 , 5.0]
    end

    #Activatecairo
    set_my_cairo_theme()

    #Create figure
    fig = Figure(resolution = (650, 400))
    ax1 = fig[1, 1] = Axis(fig , yscale = mscale)

    if( theory_comparison)
    
        #Theoretical model
        x , y = get_simple_theoretical_model_data(p0 , my_xlims[2])
        lines!(ax1 , x , y , linewidth = 1 ,  linestyle = :dash , color = "red" , label = "Model")

    end

    if(gamma_comparison)

        for (i,g) in enumerate(gamma_list)
        
            mcolor = [ "green" , "blue", "red"]
            
            #Theoretical model
            x = 0:0.1:my_xlims[2]
            y = exp.((-2*g).*x)
            lines!(ax1 , x ,  y, linewidth = 1 ,  linestyle = :dash , color = mcolor[i] , label = gamma_legend[i] )#("γ = " * "$(g).") )
    
        end
    end

    #Actual data 
    x , y = get_energy(p , energy_data_every , p.tmax)
    lines!(ax1 , x , y/y[1] , linewidth = 1 ,  linestyle = :solid , color = "black", label = "Simulation")

    #Configure plot
    ax1.xlabel = "Time"
    ax1.ylabel = "Energy Density [ ϵ(t) / ϵ(0) ]"

    limits!(ax1 , my_xlims[1]  , my_xlims[2] ,   my_ylims[1]  ,   my_ylims[2]   )    

    #Legend
    axislegend(ax1 , position = :lb)
    
    fdata   = "./data"
    outname   = fdata * p.folder *  "$(oname).png"

    #outname = abspath(p.folder, p.fname * ".png")

    save(outname, fig, px_per_unit = 10)

    println(" Image saved to: " , outname )

    println("___________________________________________________________") 

end

#Make plot for several parameters
function make_energy_plot_params(   s  , legends ; 
                                    gamma_comparison  = :false,
                                    gamma_list = nothing, 
                                    gamma_legend = nothing, 
                                    mscale = log10 ,
                                    oname = "default_name",
                                    xlims = nothing,
                                    ylims = nothing ,
                                    energy_data_every = 10,
                                    plot_title = "")

    #Welcome information
    println("\n \n")        
    println("___________________________________________________________")        
    println(" ")
    println(" Making Energy Plot from files: ")
    
    for p in s
        println("/data" , p.folder  ,p.fname)
    end
    
    println(" ")

    #Default values
    my_xlims = xlims
    my_ylims = ylims

    if( xlims === nothing )
        my_xlims = [0 , s[1].tmax]
    end
    
    if( ylims === nothing )
        my_ylims = [1E-5 , 5.0]
    end

    #Activatecairo
    set_my_cairo_theme()
    
    #Create figure
    fig = Figure(resolution = (650, 400) )
    ax1 = fig[1, 1] = Axis(fig , yscale = mscale , title = plot_title )

    pallete = ["#111111", "#65ADC2", "#233B43", "#E84646", "#C29365", "#362C21", "#316675", "#168E7F", "#109B37"]
    #colormap( "Oranges" , 10) 
 
    for (i , p) in enumerate(s)

        x , y = get_energy(p , energy_data_every , p.tmax)
    
        lines!(ax1 , x , y./y[1] , color = pallete[4+i] , linewidth = 1, label = legends[i])

    end

    #Gamma comparison
    if(gamma_comparison)

        for (i,g) in enumerate(gamma_list)
                    
            #Theoretical model
            x = 0:0.1:3000
            y = exp.((-2*g).*x)
            lines!(ax1 , x ,  y, linewidth = 1 ,  linestyle = :dash , color = pallete[4+i]  )#("γ = " * "$(g).") )
    
        end
    end


    #Configure plot
    ax1.xlabel = "Time"
    ax1.ylabel = "Energy Density [ ϵ(t) / ϵ(0) ]"

    limits!(ax1 , my_xlims[1]  , my_xlims[2] ,   my_ylims[1]  ,   my_ylims[2]   )    

    #Legend
    axislegend(ax1 ,position = :rt)
    
    fdata   = "./data/"
    outname   = fdata *  "$(oname).png"
    save(outname, fig, px_per_unit = 10)

    println(" Image saved to: " , outname )

    println("___________________________________________________________") 

end
