function convergence_parameters(p::Param , q0::Int)

    #Base number of nodes
    N0 = 2^q0 + 1

    #Create parameter list
    sim_param = []

    for q in 1:3

        #Number of nodes
        Nx = (2^(q-1))*(N0-1) + 1

        #Name of file
        filename    = string( p.fname , "_Nx_", Nx , "_Ny_", Nx , "_tmax_" , p.tmax)
        #Create simulation parameters
        s   = Param( p ,  xnodes = Nx , ynodes = Nx , fname = filename )

        #Add to list
        push!(sim_param , s)
    end

    return sim_param

end

function test_convergence( p::Param , q0::Int )

    #Create parameter list
    sim_param = convergence_parameters(p, q0)

    for x in sim_param
        solve_wave_equation_2D(x)
    end

end


function get_data_costum(p::Param , itt::Int , q0::Int, q::Int)

    #Get data at a given time
    ttt , xxx , yyy , zzz  , _ =  get_fields(p , itt)
    
    N0 = 2^q0 + 1

    #Obtain common nodes
    x_out   = [ xxx[ j + ( 2^(q-q0) - 1 )*(j-1)  ] for j in 1:N0  ]
    y_out   = [ yyy[ j + ( 2^(q-q0) - 1 )*(j-1)  ] for j in 1:N0  ]
    
    ψ = reshape(zzz , (p.xnodes , p.ynodes) )
    ψ_out  = zeros(N0,N0)


    for i in 1:N0 , j in 1:N0
        ψ_out[i,j]  = ψ[ i + ( 2^(q-q0) - 1 )*(i-1)  , j + ( 2^(q-q0) - 1 )*(j-1) ] 
    end
    
    return x_out , y_out , ψ_out
end



function explore_convergence(p::Param , q0::Int)

    sim_param = convergence_parameters(p , q0)

    #Max iteration
    max_iter , _ = get_t_max(sim_param[1])

    #Create Figure
    fig = Figure( resolution = (1000 , 950), font  = "Helvetica")

    ax1 = fig[4,1]      = Axis(fig  , width = 600 , height = 600 )
    ax2 = fig[4,2]      = Axis(fig  , width = 200 , height = 600 , xlabelvisible = :false , yaxisposition = :right)
    ax3 = fig[1:3,1]    = Axis(fig , width = 600 , height = 200 , yticksvisible = true , xaxisposition = :top)

    linkyaxes!(ax1, ax2)
    linkxaxes!(ax1, ax3)

    slider1 = Slider( fig[1,2] , range = 1:1:max_iter , height = 200/3 , width = 200)
    slider2 = Slider( fig[2,2] , range = 1:1:sim_param[1].xnodes  , height = 200/3 ,  width = 200)
    slider3 = Slider( fig[3,2] , range = 1:1:sim_param[1].ynodes  , height = 200/3 , width = 200)

    #Data for plotting

    x , y , _ = get_data_costum(sim_param[1] , 1 , q0 , q0 )

    #Get fields at a certain iteration that is controled by slider 1
    ψ1 =  lift(slider1.value) do v
        get_data_costum(sim_param[1] , v[] , q0 ,q0 )[3]
#            get_fields( sim_param[1] , v[] )[4]
        end

    ψ2 =  lift(slider1.value) do v
        get_data_costum(sim_param[2] , v[] , q0 ,q0+1 )[3]
        end

    ψ3 =  lift(slider1.value) do v
        get_data_costum(sim_param[3] , v[] , q0 ,q0+2 )[3]
        end

    #Get slices and differences for convergence analysis and scale them
    scale_factor = 100
    ψ1_x_slice = lift(slider2.value) do index
        scale_factor*(ψ1[][index,:] -ψ2[][index,:])
    end
    ψ2_x_slice = lift(slider2.value) do index
        scale_factor*2*(ψ2[][index,:] - ψ3[][index,:])
    end
    
    ψ3_x_slice = lift(slider2.value) do index
        scale_factor*4*(ψ2[][index,:] - ψ3[][index,:])
    end
    
    ψ4_x_slice = lift(slider2.value) do index
        scale_factor*16*(ψ2[][index,:] - ψ3[][index,:])
    end

    scale_factor = 100
    ψ1_y_slice = lift(slider3.value) do index
        scale_factor*(ψ1[][:,index] -ψ2[][:,index])
    end
    ψ2_y_slice = lift(slider3.value) do index
        scale_factor*2*(ψ2[][:,index] - ψ3[][:,index])
    end
    
    ψ3_y_slice = lift(slider3.value) do index
        scale_factor*4*(ψ2[][:,index] - ψ3[][:,index])
    end
    
    ψ4_y_slice = lift(slider3.value) do index
        scale_factor*16*(ψ2[][:,index] - ψ3[][:,index])
    end

    #Get vertical and horizontal lines to see cuts on histogram
    x_coord = lift(slider2.value) do index
        [ x[index] for i in 1:(2^q0+1) ]
    end
    
    y_coord = lift(slider3.value) do index
        [ y[index] for i in 1:(2^q0+1) ]
    end
    
    #y3 = lift(slider3.value) do index
    #    y[][:,index]
    #end

    #Plot histogram of lowest resolution and slices
    heatmap!(ax1 , x , y ,ψ1 ,colorrange=(-5,5) , colormap = :vik )
    lines!(ax1 , x_coord  , x )
    lines!(ax1 , x  , y_coord )

    #Plot x-slices with convergence analysis
    l1 = lines!(ax2 , ψ1_x_slice  , x,  linewidth  = 5.0  ,color = RGBA{Float32}(0.0f0, 0.8f0 , 0.9f0 , 1.0f0)) # color = ψ1_x_slice, colormap = :vik)
    l2 = lines!(ax2 , ψ2_x_slice  , x , color = :black , linestyle = :solid )
    l3 = lines!(ax2 , ψ3_x_slice  , x , color = :red   , linestyle = :solid)
    l4 = lines!(ax2 , ψ4_x_slice  , x , color = :green , linestyle = :solid)
    axislegend(ax1, [l1, l2 ,l3 ,l4], ["ψ1-ψ2", "2(ψ2-ψ3)" , "4(ψ2-ψ3)" , "16(ψ2-ψ3)"], position = :rb,
    orientation = :horizontal)

    #Plot y slices
    lines!(ax3 , y , ψ1_y_slice  ,  linewidth  = 5.0  ,color = RGBA{Float32}(0.0f0, 0.8f0 , 0.9f0 , 1.0f0)) # color = ψ1_x_slice, colormap = :vik)
    lines!(ax3 , y , ψ2_y_slice  ,  color = :black , linestyle = :solid )
    lines!(ax3 , y , ψ3_y_slice  ,  color = :red   , linestyle = :solid)
    lines!(ax3 , y , ψ4_y_slice  ,  color = :green , linestyle = :solid)    


    #Adjust limits of heatmap
    xlims!(ax1 , p.xmin ,p.xmax)
    ylims!(ax1 , p.ymin ,p.ymax)
    
    #Adjust limits of vertical
    xlims!(ax2 , -1 ,1)
    ylims!(ax3 , -1 ,1)
    
    #Do heatmap
    
    display(fig)

end


function see_energy_covergence(p::Param , q0::Int)

    sim_param = convergence_parameters(p , q0)

    #Max iteration
    max_iter , _ = get_t_max(sim_param[2])

    #Create Figure
    fig = Figure( resolution = (1000 , 950), font  = "Helvetica")
    ax = fig[1,1] = Axis(fig , title = "Test")
    
    ls = []
cols = [:red , :blue , :black]
    for (ix , x) in enumerate(sim_param)
        t , x = get_energy(x)
        push!(ls , lines!(ax , t , x  , color = cols[ix]) )    
    end

    axislegend(ax, ls, ["1" , "2" , "3"], position = :rb)

     #Adjust limits of heatmap
     ylims!(ax , 0,10)
    display(fig)

end

function plot_energy(p::Param , savetrue::Bool , name::String)

    #Max iteration
    max_iter , _ = get_t_max(p)

    #Create Figure
    fig = Figure( resolution = (800 , 450), font  = "Helvetica")
    ax = fig[1,1] = Axis(fig)
    
    t , x = get_energy(p)
    
    lines!(ax , t , x )  

    ax.ylabel = "Energy Density "
    ax.xlabel = "Time"
    
    ylims!(ax , 0,10)
    display(fig)

    if(savetrue)
        save( name , fig )
    end

end