.gitignore

@@ -3,5 +3,13 @@
 # committed for packages, but should be committed for applications that require a static
 # environment.
 Manifest.toml
+
+# VSCode settings
 .vscode/settings.json
-.DS_Store
\ No newline at end of file
+
+# Mac
+.DS_Store
+
+# Project-specific
+*.hdf5
+*.csv
\ No newline at end of file

dev/fix_scintillation_index.jl

@@ -3,122 +3,153 @@ using CSV
 include(joinpath("..", "examples", "repeated_wave_propagation.jl"))
 include(joinpath("..", "examples", "analyze_turbulent_channel_model.jl"))
 include(joinpath("..", "examples", "analyze_coupling_data.jl"))
+include(joinpath("..", "examples", "batch_analysis.jl"))
 include(joinpath("..", "examples", "from_dump.jl"))
 
 
-function analyze_batch_scintillation_index(dirname; dfcol=1)
-
+function plot_batch_scintillation_index(dirname; xcols=:, kwargs...)
     # gather all hdf5 files in the given directory `dirname`
     files = readdir(dirname, join=true)
-    csvidx = findfirst(contains("xvals.csv"), files)
-    h5files = filter(x -> splitext(x)[2] == ".hdf5", files)
-
-    Nf = length(h5files)
-
-    # gather xvals from corresponding csv if possible - otherwise 1:N
-    xval_df = if isnothing(csvidx)
-        DataFrame(:File => range(1, Nf))
+    xidx = findfirst(contains("xvals.csv"), files)
+    csvidx = findfirst(contains("scintillation_indices.csv"), files)
+    if isnothing(xidx)
+        throw(MissingException("File xvals.csv not found."))
+    end
+    if isnothing(csvidx)
+        throw(MissingException("File scintillation_indices.csv not found."))
+    end
+    # read files
+    xdf = CSV.read(files[xidx], DataFrame)
+    df = CSV.read(files[csvidx], DataFrame)
+    # gather all swept variables
+    xnames = names(xdf)
+    xnames_plt = if isone(length(xnames))
+        first(xnames)
     else
-        CSV.read(files[csvidx], DataFrame)
+        xnames[xcols]
     end
 
-    # choose column of xval if multiple parameters have been swept simultaneously
-    xname = names(xval_df)[dfcol]
-    xvals = xval_df[!, xname]
-
-    SI_est = zeros(Nf)
-    SI_pw = zeros(Nf)
-    SI_sw = zeros(Nf)
-    SI_gaussian = zeros(Nf)
-
-    # add columns for sampling requirements
-    req1 = falses(Nf)
-    req2 = falses(Nf)
-    req3 = falses(Nf)
-    req4 = falses(Nf)
+    plt = _plot_batch_scintillation_index(df, xnames, xnames_plt; kwargs...)
+    display(plt)
+    return
+end
 
-    p = Progress(Nf; showspeed=true)
 
-    for (k, file) ∈ enumerate(h5files)
-        SI_est[k] = scintillation_index(file)
+function _plot_batch_scintillation_index(df, xnames::AbstractVector{T}, xnames_plt::T; xscale="linear", yscale="linear", relative=false, reference="gaussian") where {T<:Union{Symbol,String}}
+    # get x-axis
+    x = df[!, xnames_plt]
+    yref = df[!, reference]
+    ye = df[!, "estimate"]
 
-        metadata = h5open(file, "r") do h5
-            Dict(attrs(h5))
-        end
+    if relative
+        plt = plot.scatter(x=x, y=100 * relrms(ye, yref, dims=()), name="Estimate", mode="markers", marker=Config(color="black"), showlegend=true)
 
-        Cn² = metadata["Cn2"]
-        λ = metadata["lambda"]
-        Dz = metadata["Dz"]
-        w0 = metadata["w0"]
+        # format y-axis
+        plt.layout.yaxis.title = "Rel. error w.r.t. reference / %"
+        plt.layout.yaxis.range = [0, 100]
+    else
+        plt = plot.scatter(x=x, y=yref, name="Reference ($reference)", mode="markers")
+        plt.scatter(x=x, y=ye, name="Estimate", mode="markers", marker=Config(color="black"))
         
-        l0 = metadata["l0"]
-        L0 = metadata["L0"]
-
-        # the transverse grid resolution is inherently an inner scale of turbulence
-        dxm = maximum(metadata["dx"])
-        l0 = max(l0, dxm)
-
-        SI_pw[k] = scintillation_index_pw(2π / λ, Dz, Cn², l0, L0)
-        SI_sw[k] = scintillation_index_sw(2π / λ, Dz, Cn², l0, L0)
-        SI_gaussian[k] = scintillation_index_gaussian(2π / λ, Dz, Cn², w0, l0, L0)
-
-        req1[k], req2[k], req3[k], req4[k] = post_dump_sampling_check(file; silent=true)
-
-        next!(p, showvalues = [("file", file)])
+        # format y-axis
+        plt.layout.yaxis.title = "Scintillation index"
+        plt.layout.yaxis.range = [0, 1.1 * maximum.(select(df, Not(xnames)) |> eachcol) |> maximum]
     end
 
-    # dump results
-    df = DataFrame(xname => xvals, "pw" => SI_pw, "sw" => SI_sw, "gaussian" => SI_gaussian, "estimate" => SI_est, "req1" => req1, "req2" => req2, "req3" => req3, "req4" => req4)
-    open(joinpath(dirname, "scintillation_indices.csv"), "w") do fid
-        CSV.write(fid, df)
-    end
-
-    nothing
-end
-
-
-function plot_batch_scintillation_index(dirname; xscale="linear", yscale="linear")
-    # gather all hdf5 files in the given directory `dirname`
-    files = readdir(dirname, join=true)
-    csvidx = findfirst(contains("scintillation_indices.csv"), files)
-    # read file
-    df = CSV.read(files[csvidx], DataFrame)
-
-    xname = first(names(df))
-    xvals = df[!, xname]
-
-    # visualize results
-    plt = plot.scatter(x=xvals, y=df[!, "pw"], name="Plane wave")
-    plt.scatter(x=xvals, y=df[!, "sw"], name="Spherical wave")
-    plt.scatter(x=xvals, y=df[!, "gaussian"], name="Gaussian beam")
-    plt.scatter(x=xvals, y=df[!, "estimate"], name="Estimate", mode="markers", marker=Config(color="black"))
-    
-    plt.layout.xaxis.title = string(xname)
-    plt.layout.yaxis.title = "Scintillation index"
-    # set lower yaxis limit to 0
-    yl = plt.layout.yaxis.range
-    yl[1] = 0
-    plt.layout.yaxis.range = yl
+    plt.layout.xaxis.title = string(xnames_plt)
     # set x/y axis to logarithmic
     plt.layout.xaxis.type = xscale
     plt.layout.yaxis.type = yscale
 
-    display(plt)
+    return plt
+end
 
-    # # visualize relative errors
-    # plte = plot.scatter(x=xvals, y=100 * relrms(df[!, "estimate"], df[!, "pw"], dims=()), name="Plane wave")
-    # plte.scatter(x=xvals, y=100 * relrms(df[!, "estimate"], df[!, "sw"], dims=()), name="Spherical wave")
-    # plte.scatter(x=xvals, y=100 * relrms(df[!, "estimate"], df[!, "gaussian"], dims=()), name="Gaussian beam")
+function _plot_batch_scintillation_index(df, xnames::AbstractVector{T}, xnames_plt::AbstractVector{T}; xscale="linear", yscale="linear", relative=false, reference="gaussian") where {T<:Union{Symbol,String}}
+
+    # we do not support anything but two xnames in this function
+    @assert length(xnames_plt) == 2
+
+    # get z-data
+    z_df = unstack(df, xnames_plt..., "estimate")
+    # convert z-data to Matrix
+    z = Matrix(z_df[:, 2:end]) |> transpose
+
+    # detect dependent variables by checking whether the off-diagonals contain only NaNs
+    zv = vec(z) |> collect
+    offdiag = deleteat!(zv, diagind(z))
+
+    if all(ismissing, offdiag)
+        # we consider two dependent variables
+        ye = diag(z)
+        # gather corresponding reference values
+        zref_df = unstack(df, xnames_plt..., reference)
+        yref = diag(Matrix(zref_df[:, 2:end]))
+
+        # gather both x-axes
+        x1 = string.(z_df[!, xnames_plt[1]])
+        x2 = names(z_df)[2:end]
+
+        if relative
+            err = 100 * relrms(ye, yref, dims=())
+            plt = plot.scatter(x=x1, y=err, name="Estimate", mode="markers", marker=Config(color=err, cmin=0, cmax=100), showlegend=true)
+            plt.scatter(x=x2, y=fill(NaN, size(x2)), xaxis="x2", showlegend=false)
+
+            # format y-axis
+            plt.layout.yaxis.title = "Rel. error w.r.t. reference / %"
+            plt.layout.yaxis.range = [0, 100]
+        else
+            plt = plot.scatter(x=x1, y=yref, name="Reference ($reference)", mode="markers")
+            plt.scatter(x=x1, y=ye, name="Estimate", mode="markers", marker=Config(color="black"))
+            plt.scatter(x=x2, y=fill(NaN, size(x2)), xaxis="x2", showlegend=false)
+
+            # format y-axis
+            plt.layout.yaxis.title = "Scintillation index"
+            plt.layout.yaxis.range = [0, 1.1 * maximum.(select(df, Not(xnames)) |> eachcol) |> maximum]
+        end
 
-    # plte.layout.yaxis.range = [0, 200]
-    # plte.layout.xaxis.title = string(xname)
-    # plte.layout.yaxis.title = "Rel. error / %"
-    
-    # # set x/y axis to logarithmic
-    # plte.layout.xaxis.type = xscale
-    # plte.layout.yaxis.type = yscale
+        # format primary xaxis
+        plt.layout.xaxis.title = string(xnames_plt[1])
+        # format second x-axis
+        plt.layout.xaxis2.position = 1.0
+        plt.layout.xaxis2.anchor = "free"
+        plt.layout.xaxis2.overlaying = "x"
+        plt.layout.xaxis2.side = "top"
+        plt.layout.xaxis2.title = string(xnames_plt[2])
+        # legend
+        plt.layout.legend.orientation = "h"
+        plt.layout.legend.xanchor = "center"
+        plt.layout.legend.x = 0.5
+        plt.layout.legend.y = 1.4
+
+        return plt
+    else
+        # get x-axis, should be categorical!
+        x = string.(z_df[!, xnames_plt[1]])
+        # get y-axis, should be categorical!
+        y = names(z_df)[2:end]
+
+        if relative
+            # get z-data
+            zref_df = unstack(df, xnames_plt..., reference)
+            # convert z-data to Matrix
+            zref = Matrix(zref_df[:, 2:end]) |> transpose
+
+            plt = plot.heatmap(x=x, y=y, z=100 * relrms(z, zref, dims=()), name="Estimate", zmin=0, zmax=100, colorscale="Reds")
+            
+            plt.layout.title = "Rel. error w.r.t. reference / %"
+        else
+            plt = plot.heatmap(x=x, y=y, z=z, name="Estimate", zmin=0, zmax=1.1 * maximum.(select(df, Not(xnames)) |> eachcol) |> maximum)
+            
+            plt.layout.title = "Scintillation index"
+        end
 
-    # display(plte)
+        plt.layout.xaxis.title = string(xnames_plt[1])
+        plt.layout.yaxis.title = string(xnames_plt[2])
 
-    return
+        # set x/y axis to logarithmic if needed
+        plt.layout.xaxis.type = xscale
+        plt.layout.yaxis.type = yscale
+
+        return plt
+    end
 end
\ No newline at end of file

examples/batch_analysis.jl

 include("analyze_turbulent_channel_model.jl")
 include("analyze_coupling_data.jl")
+include("repeated_wave_propagation.jl")
 
-function batch_repeated_wave_propagation(args=Dict(:Nz => 4:16, :dz => 1e3 ./ (3:15)); iters=100, λ=1.55e-6, showprogress=true, samplingcheck=false, kwargs...)
+using CSV
 
-    # make sure, that all value arrays of the argument dict are of the same size
-    valuelen = length.(values(args))
-    Nv = minimum(valuelen)
-    @assert isone(length(unique(valuelen)))
 
-    dir = timeprefix(fmt="YYYY-mm-dd")
+function outer_product_across_dicts(dicts...)
+    if isempty(dicts)
+        return []
+    end
+
+    # Make sure the value arrays of each dict are of the same size
+    for dict ∈ dicts
+        valuelen = length.(values(dict))
+        
+        if !isone(length(unique(valuelen)))
+            @info "Length of values:" Dict(k => length(v) for (k, v) ∈ dict)
+            throw(ArgumentError("All values of an argument must have the same length."))
+        end
+    end
+
+    # Extract keys and vectors from each dictionary
+    dict_keys_vectors = [collect(pairs(dict)) for dict ∈ dicts]
+
+    # Calculate the number of combinations
+    dims = [length(first(kv)[2]) for kv ∈ dict_keys_vectors]
+
+    # Initialize the result array
+    result = Array{NamedTuple, length(dims)}(undef, dims...)
+
+    # Generate all combinations of indices
+    indices_combinations = Iterators.product(map(x -> 1:x, dims)...)
+
+    # Populate the result array, use an OrderedDict to maintain the order!
+    for indices ∈ indices_combinations
+        current_dict = OrderedDict{Symbol, Any}()
+        for (i, kv_pairs) ∈ enumerate(dict_keys_vectors)
+            for (key, vector) ∈ kv_pairs
+                current_dict[key] = vector[indices[i]]
+            end
+        end
+        result[indices...] = NamedTuple(current_dict)
+    end
+
+    return result
+end
+
+
+function batch_repeated_wave_propagation(arg1=Dict(:Nz => 4:16, :dz => 1e3 ./ (3:15)), additional...; showprogress=true, samplingcheck=false, kwargs...)
+
+    args = (arg1, additional...)
 
-    for (k, v) ∈ pairs(args)
-        vmin, vmax = extrema(v)
-        dir *= "$(k)_$(round(vmin,sigdigits=3))_$(round(vmax,sigdigits=3))"
+    dir = timeprefix(; fmt="YYYY-mm-dd", postfix="")
+
+    # Build name of target directory
+    for arg ∈ args
+        for (k, v) ∈ pairs(arg)
+            vmin, vmax = if eltype(v) <: Integer
+                extrema(v)
+            else
+                round.(extrema(v), sigdigits=3)
+            end
+            
+            dir *= "_$(k)_$(vmin)_$(vmax)"
+        end
     end
 
     if isdir(dir)
         if !isempty(readdir(dir))
-            @warn "Directory $dir already exists! Do you want to write into this directory anyway? (y/N)"
+            @warn "Directory $dir already exists and is not empty! Do you want to write into this directory anyway? (y/N)"
             response = lowercase(strip(readline()))
             if response == "y"
                 println("Continuing...")
@@ -30,31 +81,21 @@ function batch_repeated_wave_propagation(args=Dict(:Nz => 4:16, :dz => 1e3 ./ (3
         mkdir(dir)
     end
 
-    p = Progress(Nv; showspeed=true, enabled=showprogress)
-    
-    if :iters ∈ keys(args)
-        # assume no other parameter is swept
-        (length(args) > 1) && @warn "Other parameters besides iters are ignored..."
-        for i ∈ 1:Nv
-            # round! iters must take integer values
-            _iters = round(Int, args[:iters][i])
-            # we have to sort the outputs, why we prepend a number to the filenames | simply provide a namedtuple (; iters) to the kwarg meta
-            repeated_wave_propagation(_iters; λ, dir, fnprefix=lpad(i, ceil(Int, log10(Nv)), "0") * "_", showprogress=false, samplingcheck, meta=(; iters=_iters), kwargs...)
-            sleep(0.1)
-            next!(p, showvalues = [(string(k), v) for (k, v) ∈ pairs(xargs)])
-        end
-    else
-        for i ∈ 1:Nv
-            xargs = (; (k => v[i] for (k, v) ∈ pairs(args))...)
-            # we have to sort the outputs, why we prepend a number to the filenames
-            repeated_wave_propagation(iters; xargs..., λ, dir, fnprefix=lpad(i, ceil(Int, log10(Nv)), "0") * "_", showprogress=false, samplingcheck, meta=xargs, kwargs...)
-            sleep(0.1)
-            next!(p, showvalues = [(string(k), v) for (k, v) ∈ pairs(xargs)])
-        end
+    # Set up the array of NamedTuples for the iteration
+    itr = outer_product_across_dicts(args...)
+    Nv = length(itr)    
+
+    p = Progress(length(itr); showspeed=true, enabled=showprogress)
+
+    for (i, xarg) ∈ enumerate(itr)
+        # we have to sort the outputs, why we prepend a number to the filenames
+        repeated_wave_propagation(; xarg..., dir, fnprefix=lpad(i, ceil(Int, log10(Nv + 1)), "0") * "_", showprogress=false, samplingcheck, meta=xarg, kwargs...)
+        sleep(0.1)
+        next!(p, showvalues=[(string(k), v) for (k, v) ∈ pairs(xarg)])
     end
 
     # Save to CSV
-    df = DataFrame(args)
+    df = DataFrame(itr |> vec)
     open(joinpath(dir, "xvals.csv"), "w") do fid
         CSV.write(fid, df)
     end
@@ -62,6 +103,71 @@ function batch_repeated_wave_propagation(args=Dict(:Nz => 4:16, :dz => 1e3 ./ (3
     nothing
 end
 
+function batch_scintillation_index(dirname)
+
+    # gather all hdf5 files in the given directory `dirname`
+    files = readdir(dirname, join=true)
+    csvidx = findfirst(contains("xvals.csv"), files)
+    h5files = filter(x -> splitext(x)[2] == ".hdf5", files)
+
+    Nf = length(h5files)
+
+    # gather xvals from corresponding csv if possible - otherwise 1:N
+    df = if isnothing(csvidx)
+        DataFrame(:File => range(1, Nf))
+    else
+        CSV.read(files[csvidx], DataFrame)
+    end
+
+    SI_est = zeros(Nf)
+    SI_pw = zeros(Nf)
+    SI_sw = zeros(Nf)
+    SI_gaussian = zeros(Nf)
+
+    p = Progress(Nf; showspeed=true)
+
+    for (k, file) ∈ enumerate(h5files)
+        SI_est[k] = scintillation_index(file)
+
+        metadata = attributes_from_dump(file)
+
+        # make sure the order the files are read matches the xvals in the DataFrame
+        xname = first(names(df))
+        @assert df[k, xname] == metadata[string(xname)]
+
+        Cn² = metadata["Cn2"]
+        λ = metadata["lambda"]
+        Dz = metadata["Dz"]
+        w0 = metadata["w0"]
+
+        l0 = metadata["l0"]
+        L0 = metadata["L0"]
+
+        # the transverse grid resolution is inherently an inner scale of turbulence
+        dxm = maximum(metadata["dx"])
+        l0 = max(l0, dxm)
+
+        SI_pw[k] = scintillation_index_pw(2π / λ, Dz, Cn², l0, L0)
+        SI_sw[k] = scintillation_index_sw(2π / λ, Dz, Cn², l0, L0)
+        SI_gaussian[k] = scintillation_index_gaussian(2π / λ, Dz, Cn², w0, l0, L0)
+
+        next!(p, showvalues=[("file", file)])
+    end
+
+    # dump results
+    df[!, "pw"] = SI_pw
+    df[!, "sw"] = SI_sw
+    df[!, "gaussian"] = SI_gaussian
+    df[!, "estimate"] = SI_est
+
+    open(joinpath(dirname, "scintillation_indices.csv"), "w") do fid
+        CSV.write(fid, df)
+    end
+
+    nothing
+end
+
+
 function batch_coupling(; basedir=pwd())
     fns = readdir(basedir, join=true)
     filter!(contains(r"\d+\.hdf5"), fns)

examples/from_dump.jl

@@ -29,10 +29,13 @@ function plot_from_dump(fn; iter=1, mode=1, collimate=false)
     Plot(tr, lt)
 end
 
-function attributes_from_dump(fn)
-    h5open(fn, "r") do fid
+function attributes_from_dump(fn; disp=false)
+    meta = h5open(fn, "r") do fid
         OrderedDict(attrs(fid))
     end
+    disp && show(IOContext(stdout, :limit => false), "text/plain", meta)
+
+    return meta
 end
 
 function post_dump_sampling_check(fn; kwargs...)

examples/repeated_wave_propagation.jl

@@ -21,9 +21,9 @@ function repeated_wave_propagation(utx::AbstractArray{Complex{T}}, midx, domain:
     # estimate Fried parameters for partial propagations
     r0 = partial_r0_pw(2π / λ, domain.z, Cn²)
 
+    σr²eff = rytov_pw(2π / λ, domain.Dz, channel.Cn²)
     # compute Rytov variance per step
-    σr² = rytov_pw(2π / λ, domain.dz, channel.Cn²)
-    (σr² >= 0.1) && @warn "Insufficient number of planes, not within weak fluctuation regime."
+    σr² = rytov_pw.(2π / λ, diff(domain.z), channel.Cn²)
 
     # estimated beam waist radius at receiver
     zr = 2π / λ * w0^2 / 2
@@ -74,7 +74,7 @@ function repeated_wave_propagation(utx::AbstractArray{Complex{T}}, midx, domain:
             attrs(fid)["modes"] = midx
             attrs(fid)["iters"] = iters
             attrs(fid)["nchunks"] = nchunks
-            attrs(fid)["chunksize"] = chunksize
+            attrs(fid)["chunksize"] = min(chunksize, iters)
             attrs(fid)["dx"] = collect(domain.dx)
             attrs(fid)["dy"] = collect(domain.dy)
             attrs(fid)["dz"] = domain.dz
@@ -84,6 +84,8 @@ function repeated_wave_propagation(utx::AbstractArray{Complex{T}}, midx, domain:
             attrs(fid)["D_r0"] = D_r0
             attrs(fid)["r0eff"] = r0eff
             attrs(fid)["r0"] = r0
+            attrs(fid)["sigmaR2eff"] = σr²eff
+            attrs(fid)["sigmaR2"] = σr²
             attrs(fid)["Cn2"] = channel.Cn²
             attrs(fid)["l0"] = channel.l0
             attrs(fid)["L0"] = channel.L0
@@ -111,7 +113,7 @@ function repeated_wave_propagation(utx::AbstractArray{Complex{T}}, midx, domain:
     nothing
 end
 
-function repeated_wave_propagation(iters, T::DataType=Float64; Nx::Int=128, λ=1.55e-6, w0=1e-2, dx1=1e-3, dxn=dx1, dz=5e1, Nz::Int=21, m=0, n=0, kwargs...)
+function repeated_wave_propagation(T::DataType=Float64; iters=100, Nx::Int=128, λ=1.55e-6, w0=1e-2, dx1=1e-3, dxn=dx1, dz=5e1, Nz::Int=21, m=0, n=0, kwargs...)
     beam = Beam{T}(m, n, λ, w0)
 
     domain = if dx1 == dxn
@@ -126,7 +128,7 @@ function repeated_wave_propagation(iters, T::DataType=Float64; Nx::Int=128, λ=1
     repeated_wave_propagation(utx, midx, domain, iters; beam.λ, beam.w0, kwargs...)
 end
 
-function repeated_wave_propagation_multi(iters, T::DataType=Float64; Nx::Int=128, λ=1.55e-6, w0=1e-2, dx1=1e-3, dxn=dx1, dz=5e1, Nz::Int=21, ng=3, kwargs...)
+function repeated_wave_propagation_multi(T::DataType=Float64; iters=100, Nx::Int=128, λ=1.55e-6, w0=1e-2, dx1=1e-3, dxn=dx1, dz=5e1, Nz::Int=21, ng=3, kwargs...)
     # simulation domain
     domain = if dx1 == dxn
         Domain{T}(dx1, dz, Nx, Nz)

examples/sdm_channel_stats.jl

@@ -167,6 +167,9 @@ function sdm_statistics(fn; ng=nothing, basedir=pwd(), Nopt=10, statfun=mdl, tar
 
     if showplot
         x, y = estimate_ccdf(stat)
-        plot.scatter(; x, y)
+        plt = plot.scatter(; x, y)
+        plt.layout.yaxis.type = "log"
+
+        display(plt)
     end
 end
\ No newline at end of file

src/estimate_statistics.jl

@@ -229,7 +229,7 @@ end
 """
 function scv(u::AbstractArray{T,N}; dims=:) where {T,N}
     I = abs2.(u)
-    var(I; dims, corrected=false) ./ mean(I; dims) .^ 2
+    mean(abs2, I; dims) ./ mean(I; dims) .^ 2 .- 1 # it might be advantageous to not use `var` to avoid computing the mean twice
 end
 
 

src/helpers.jl

@@ -18,7 +18,7 @@ containertype(::T) where {T} = T.name.wrapper
 
 
 # generate current time prefix
-timeprefix(t::DateTime=now(); fmt="YYYY-mm-dd-HH-MM-SS") = "$(Dates.format(t, fmt))_"
+timeprefix(t::DateTime=now(); fmt="YYYY-mm-dd-HH-MM-SS", postfix="_") = "$(Dates.format(t, fmt))" * postfix
 
 function cartgrid(M, dx, N, dy)
     x = fftfreq(M, M * dx)