text:
@testset "fsm_loopback.jl" begin
    DEFAULT_CAPACITY = 65000
    TIMEOUT_DFLT = 1000
    NREPS = 64
    MAX_PAYLOAD_SIZE = 1460
    @testset "loop segments back into the same FSM" begin
        for i in 1:NREPS
            rx_offset = WrappingInt32(rand(UInt32))
            conn = fsm_in_established(rx_offset - 1, rx_offset - 1; cap=DEFAULT_CAPACITY)
            send_ack!(conn, rx_offset, rx_offset, 65000)
            d = rand(UInt8(1):UInt8(127), DEFAULT_CAPACITY) |> pointer |> unsafe_string
            sendoff = 0
            while sendoff < length(d)
                len = min(length(d) - sendoff, rand(UInt) % 8192)
                len == 0 && continue
                write!(conn, @view(d[sendoff+1:sendoff+len]))
                tick!(conn, 1)
                @test bytes_in_flight(conn) == len
                @test conn.segments_out |> !isempty
                n_segments = (len - 1) ÷ MAX_PAYLOAD_SIZE + 1
                bytes_remaining = len
                for i in 1:n_segments
                    expected_size = min(bytes_remaining, MAX_PAYLOAD_SIZE)
                    seg = expect_seg(conn; payload_size = expected_size)
                    bytes_remaining -= expected_size
                    send_seg!(conn, seg)
                    tick!(conn, 1)
                end
                for i in 1:n_segments
                    seg = expect_seg(conn; ack=true, payload_size=0)
                    send_seg!(conn, seg)
                    tick!(conn, 1)
                end
                expect_no_seg(conn)
                @test bytes_in_flight(conn) == 0
                sendoff += len
            end
            expect_data!(conn, d)
        end
    end
end

max_stars_repo_name: AquaIndigo/JLSponge

max_issues_repo_name: AquaIndigo/JLSponge

max_forks_repo_name: AquaIndigo/JLSponge


text:
# Copyright (c) Microsoft Corporation.
# Licensed under the MIT license.

module TestLoader
include("../shared/load.jl")
using ADPerfTest
using GMMData
using BAData
using HandData
using LSTMData

export get_gmm_test, get_ba_test, get_hand_test, get_lstm_test, get_test_for

module Box
end

loaded_modules = []

"""
    get_gmm_test(module_name)

Loads a `Test{GMMInput, GMMOutput}` object from the module `module_name`.

Callbacks of this object will be created in a world newer than the current one and as such will have to be
invoked via `Base.invokelatest`.
"""
get_gmm_test(module_name::AbstractString)::Test{GMMInput, GMMOutput} = get_test("gmm", module_name)

"""
    get_ba_test(module_name)

Loads a `Test{BAInput, BAOutput}` object from the module `module_name`.

Callbacks of this object will be created in a world newer than the current one and as such will have to be
invoked via `Base.invokelatest`.
"""
get_ba_test(module_name::AbstractString)::Test{BAInput, BAOutput} = get_test("ba", module_name)

"""
    get_hand_test(module_name)

Loads a `Test{HandInput, HandOutput}` object from the module `module_name`.

Callbacks of this object will be created in a world newer than the current one and as such will have to be
invoked via `Base.invokelatest`.
"""
get_hand_test(module_name::AbstractString)::Test{HandInput, HandOutput} = get_test("hand", module_name)

"""
    get_lstm_test(module_name)

Loads a `Test{LSTMInput, LSTMOutput}` object from the module `module_name`.

Callbacks of this object will be created in a world newer than the current one and as such will have to be
invoked via `Base.invokelatest`.
"""
get_lstm_test(module_name::AbstractString)::Test{LSTMInput, LSTMOutput} = get_test("lstm", module_name)

"""
    get_test_for(input, output, module_name)

Loads a `Test{Input, Output}` object appropriate for provided `input` and `output` from the module `module_name`.

Callbacks of this object will be created in a world newer than the current one and as such will have to be
invoked via `Base.invokelatest`.
"""
function get_test_for end
get_test_for(input::GMMInput, output::GMMOutput, module_name::AbstractString)::Test{GMMInput, GMMOutput} = get_gmm_test(module_name)
get_test_for(input::BAInput, output::BAOutput, module_name::AbstractString)::Test{BAInput, BAOutput} = get_ba_test(module_name)
get_test_for(input::HandInput, output::HandOutput, module_name::AbstractString)::Test{HandInput, HandOutput} = get_hand_test(module_name)
get_test_for(input::LSTMInput, output::LSTMOutput, module_name::AbstractString)::Test{LSTMInput, LSTMOutput} = get_lstm_test(module_name)

"""
    get_test(test_name, module_name)

Loads a `Test{Input, Output}` object appropriate for `test_name` objective from the module `module_name`.

Callbacks of this object will be created in a world newer than the current one and as such will have to be
invoked via `Base.invokelatest`.

Dynamically loads module `module_name` into an isolated environment if necessary and calls
`get_\$(test_name)_test()::Test{Input, Output}` that it is expected to export.
"""
function get_test(test_name::AbstractString, module_name::AbstractString)
    box_module_name = module_name * "Box"
    if !(module_name ∈ loaded_modules)
        # Creating a new submodule coorespondig to the module being loaded
        @eval(Box, $(Symbol(box_module_name)) = Module())
        # "using" the module being loaded in the new submodule
        @eval(Box, @eval($(Symbol(box_module_name)), using $(Symbol(module_name))))
        push!(loaded_modules, module_name)
    end
    get_test_function_name = "get_$(test_name)_test"
    # Checking that the submodule that uses the target module has a test factory function defined
    if @eval !isdefined($(Meta.parse("Box.$box_module_name")), Symbol($get_test_function_name))
        throw(ArgumentError("Module $module_name doesn't export a funtion $get_test_function_name."))
    end
    # Evaluating test factory function that we imported ("used") from the module $module_name
    # into submodule Box.$box_module_name
    @eval Base.invokelatest($(Meta.parse("Box.$box_module_name.$get_test_function_name")))
end

end

ext: jl

max_stars_repo_name: dsyme/ADBench

max_issues_repo_name: dsyme/ADBench

max_forks_repo_name: dsyme/ADBench


text:
module PDB

using Compat

export parse_entry, eachatom, eachmol, eachchain, name, models

import Base: convert, show, start, done, next

include("element.jl")
include("parser.jl")

end # module


text:
# # Structured concurrency

using Julio
using Test

# ## [Black box rule](@id black-box)
#
# In sequential programs, the side-effects of a function are "done" by the time
# the function returns [^closure]. However, many concurrent programming paradims
# do not let us assume such a simple but yet highly useful property. In [Notes
# on structured concurrency, or: Go statement considered harmful — njs
# blog](https://vorpus.org/blog/notes-on-structured-concurrency-or-go-statement-considered-harmful/)
# (See also: [Trio: Async concurrency for mere mortals - PyCon 2018 -
# YouTube](https://www.youtube.com/watch?v=oLkfnc_UMcE)), Nathaniel J. Smith
# called this property the *"black box rule"*.
#
# [^closure]: Programming constructs such as closures, coroutines, and "methods"
#     in class-based object oriented programming languages may be considered as
#     mechanisms for "resuming" the side-effects. However, they have visible
#     syntax (e.g., function call) for resuming the side-effect. Thus, these
#     constructs still follow the black box rule.
#
# ### Programming without the black box rule
#
# To demonstrate the pain comes with functions that does not follow the black
# box rule, let us consider the following simple function:

function channel_map_unstructured!(f, output, input; ntasks = Threads.nthreads())
    for _ in 1:ntasks
        Threads.@spawn for x in input
            y = f(x)
            put!(output, y)
        end
    end
end

# This function `channel_map_unstructured!` does not follow the black box rule
# because it does not wait for the spawned tasks; i.e., these "leaked" tasks are
# keep mutating `output` and `input` even after `channel_map_unstructured!`
# returns.  Programs using such functions like this are very hard to understand.

function test_channel_map_unstructured()
    #=
    Suppose we need to combine the result of two kinds of "computations" into
    one output channel (`source_channel` is defined below):
    =#
    output = Channel()
    channel_map_unstructured!(output, source_channel(1:100); ntasks = 10) do x
        sleep(0.01)
        2x
    end
    channel_map_unstructured!(output, source_channel(1:100); ntasks = 10) do x
        sleep(0.01)
        2x + 1
    end
    #=
    It is tricky to accumulate the result reliably (another disadvantage of
    unstructured concurrencty).  For this demonstration, we can "cheat" since we
    know the size of the output.
    =#
    results = []
    for _ in 1:200
        push!(results, take!(output))
    end
    close(output)
    #=
    Ideally, the program order is reflected in the result. If it were the case,
    we should see even numbers first and then odd numbers. However, since
    `channel_map_unstructured!` violates the black box rule and the "leaked"
    tasks are keep adding results to the `output` channel, we can't understand
    the program by looking at the invocations of `channel_map_unstructured!`.
    =#
    @test_broken all(iseven, results[1:end÷2])
    @test_broken all(isodd, results[end÷2+1:end])
    #=
    ---
    =#
end

# (The above example uses a simple utility function `source_channel` for
# generating the input source:)

function source_channel(xs)
    ch = Channel(Inf)
    for x in xs
        put!(ch, x)
    end
    close(ch)
    return ch
end

# ### Julio API enforces the black box rule
#
# Julio "enforces" the black box rule by providing API such that `Julio.spawn!`
# can be called only within a dynamical scope of `Julio.withtaskgroup`.  Thus,
# mechanically translating `channel_map_unstructured!` to use Julio API gives us
# a function that follows the black box rule.

function channel_map_structured!(f, output, input; ntasks = Threads.nthreads())
    Julio.withtaskgroup() do tg
        for _ in 1:ntasks
            Julio.spawn!(tg) do
                for x in input
                    y = f(x)
                    put!(output, y)
                end
            end
        end
    end
end

function test_channel_map_structured()
    send_endpoint, receive_endpoint = Julio.queue()
    try
        channel_map_structured!(send_endpoint, source_channel(1:100); ntasks = 10) do x
            sleep(0.01)
            2x
        end
        channel_map_structured!(send_endpoint, source_channel(1:100); ntasks = 10) do x
            sleep(0.01)
            2x + 1
        end
    finally
        close(send_endpoint)
    end
    results = collect(receive_endpoint)
    @test all(iseven, results[1:end÷2])
    @test all(isodd, results[end÷2+1:end])
end

# TODO: explain the nursery passing style

# ## Error handling
#
# It is possible to satisfy the black box rule using the `Base` API:

function channel_map_base!(f, output, input; ntasks = Threads.nthreads())
    @sync for _ in 1:ntasks
        Threads.@spawn for x in input
            y = f(x)
            put!(output, y)
        end
    end
end

# However, it is problematic when `f` throws. In the above example, the
# execution will not be finished until either all items in `input` are consumed
# or all tasks throw.

function error_on_10(x)
    x == 10 && error("error in one task")
    sleep(0.01)
    return x
end

function test_channel_map_base()
    output = Channel(Inf)
    try
        channel_map_base!(error_on_10, output, source_channel(1:100); ntasks = 10)
    catch
    end
    close(output)
    results = collect(output)
    push!(results, 10)
    sort!(results)
    @test results == 1:100
end

# !!! note
#     `Base.Experimental.@sync` can be used to throw an error as soon as the
#     first task throws. However, it then leaks unfinished tasks; i.e., we can't
#     assume the black box rule anymore.
#
# ### Manual concurrent error handling is hard
#
# In general, it is hard to implement robust error handling using the `Base`
# API. Even though it is possible to do so case-by-case basis, there is no
# simple mechanism that the users can rely on.  For example, we can introduce an
# intermediate channel in the above example. This intermediate channel will be
# closed on error and hence unblock all the tasks. However, this strategy
# results in a subtle code that obscures the core logic:

function channel_map_base2!(f, output, input; ntasks = Threads.nthreads())
    tmpch = Channel() do ch
        for y in ch
            put!(output, y)
        end
    end
    try
        @sync for _ in 1:ntasks
            Threads.@spawn try
                for x in input
                    y = f(x)
                    put!(tmpch, y)
                end
            catch
                close(tmpch)
                rethrow()
            end
        end
    finally
        close(tmpch)
    end
end

function test_channel_map_base2()
    output = Channel(Inf)
    try
        channel_map_base2!(error_on_10, output, source_channel(1:100); ntasks = 10)
    catch
    end
    close(output)
    results = collect(output)
    @test length(results) < 100
end

# ### Julio automates concurrent error handling
#
# When using Julio API, errors are automatically propagated.  In fact,
# `channel_map_structured!` defined above already have the desired property:

function test_channel_map_structured_error()
    send_endpoint, receive_endpoint = Julio.queue()
    try
        channel_map_structured!(
            error_on_10,
            send_endpoint,
            source_channel(1:100);
            ntasks = 10,
        )
    catch
    finally
        close(send_endpoint)
    end
    results = collect(receive_endpoint)
    @test length(results) < 100
end

# For more detailed controll on cancellation, see also:
# [Example: bounding search results](@ref ex-bounding-search)


text:
function returnsplot(df::DataFrame, xv::String, stock::String)

  df = copy(df)

  g  = plot(df,
            {:x => xv},
            Geom.bar)

  filename = strcat(stock, ".svg")
  img = SVG(filename, 6inch, 4inch)
  draw(img, g)
  finish(img)
  run(`open -a safari $filename`)
end


text:
abstract type AbstractLocalTurbulenceIntensityModel end

"""
    LocalTIModelNoLocalTI()

Don't calculate local turbulence intensity. Ambient TI will be used instead for all points

"""
struct LocalTIModelNoLocalTI{} <: AbstractLocalTurbulenceIntensityModel

end

"""
    LocalTIModelMaxTI(astar, bstar, k1, k2)

Calculate local turbulence intensity using the model presented in Niayifar and
Porte Agel (2015, 2016)

# Arguments
- `astar::Float`: wake spreading parameter from Bastankhah and Porte-Agel Gaussian wake model
- `bstar::Float`: wake spreading parameter from Bastankhah and Porte-Agel Gaussian wake model
- `k1::Float`: slope of k vs TI curve
- `k2::Float`: vertical offset of k vs TI curve
"""
struct LocalTIModelMaxTI{TF} <: AbstractLocalTurbulenceIntensityModel
    astar::TF
    bstar::TF
    k1::TF
    k2::TF
end
LocalTIModelMaxTI(x, y) = LocalTIModelMaxTI(x, y, 0.3837, 0.003678)
LocalTIModelMaxTI() = LocalTIModelMaxTI(2.32, 0.154, 0.3837, 0.003678)

"""
    calculate_local_ti(turbine_x, turbine_y, ambient_ti, rotor_diameter, hub_height, turbine_yaw, turbine_local_ti, sorted_turbine_index,
    turbine_inflow_velcities, turbine_ct, ti_model::LocalTIModelNoLocalTI; turbine_id=1, tol=1E-6)

Returns ambient turbulence intesity value whenever local turbulence intensity is requestesd

# Arguments
- `turbine_x::Array{Float,nTurbines}`: turbine wind direction locations in the wind direction
    reference frame
- `turbine_y::Array{Float,nTurbines}`: turbine cross wind locations in the wind direction
    reference frame
- `ambient_ti::Float`: ambient turbulence intensity
- `rotor_diameter::Array{Float,nTurbines}`: rotor diameters of all turbines
- `hub_height::Array{Float,nTurbines}`: hub heights of all turbines relative to the ground
- `turbine_yaw::Array{Float,nTurbines}`: yaw of all turbines for the current wind state in radians
- `turbine_local_ti::Array{Float,nTurbines}`: local turbulence intensity of all turbines for the current wind state`
- `sorted_turbine_index::Array{Float,nTurbines}`: turbine north-south locations in the
    global reference frame
- `turbine_inflow_velcities::Array{Float,nTurbines}`: effective inflow wind speed at each turbine for given state
- `turbine_ct::Array{Float,nTurbines}`: thrust coefficient of each turbine for the given state
- `ti_model::LocalTIModelNoLocalTI`: contains a struct defining the desired turbulence intensity model, no local TI in this case
- `turbine_id::Int`: index of wind turbine of interest. Provide 1 as default.
- `tol::Float`: How far upstream a turbine should be before being included in TI calculations
"""
function calculate_local_ti(turbine_x, turbine_y, ambient_ti, rotor_diameter, hub_height, turbine_yaw, turbine_local_ti, sorted_turbine_index,
                    turbine_inflow_velcities, turbine_ct, ti_model::LocalTIModelNoLocalTI; turbine_id=1, tol=1E-6)
    return ambient_ti
end

"""
    _k_star_func(ti_ust,k1, k2)

Calculate local wake spreading rate based on turbulence intensity using the model presented
in Niayifar and Porte Agel (2015, 2016)

# Arguments
- `ti_ust::Float`: upstream turbine local turbulence intensity
- `k1::Float`: slope of k vs TI curve
- `k2::Float`: vertical offset of k vs TI curve
"""
# compute wake spread parameter based on local turbulence intensity
function _k_star_func(ti_ust,k1,k2)
    # calculate wake spread parameter from Niayifar and Porte Agel (2015, 2016)

    k_star_ust = k1*ti_ust + k2

    return k_star_ust

end
_k_star_func(x) = _k_star_func(x, 0.3837, 0.003678)

"""
    _niayifar_added_ti_function(x, d_dst, d_ust, h_ust, h_dst, ct_ust, kstar_ust, delta_y,
        ti_amb, ti_ust, ti_dst, ti_area_ratio_in; s=700.0)

Main code for calculating the local turbulence intensity at a turbine using the method of
    Niayifar and Porte Agel (2015, 2016).

# Arguments
- `x::Float`: downstream distance from turbine to point of interest
- `d_dst::Float`: downstream turbine rotor diameter
- `d_ust::Float`: upstream turbine rotor diameter
- `h_ust::Float`: upstream turbine hub height
- `h_dst::Float`: downstream turbine hub height
- `ct_ust::Float`: upstream turbine thrust coefficient
- `kstar_ust::Float`: upstream turbine wake expansion rate
- `delta_y::Float`: cross wind separation from turbine to point of interest
- `ti_amb::Float`: ambient turbulence intensity
- `ti_ust::Float`: upstream turbine local turbulence intensity
- `ti_dst::Float`: downstream turbine local turbulence intensity
- `ti_area_ratio_in::Float`: current value of TI-area ratio for use in calculatin local TI
- `s::Float`: smooth max smootheness parameter
"""
function _niayifar_added_ti_function(x, d_dst, d_ust, h_ust, h_dst, ct_ust, kstar_ust, delta_y, ti_amb, ti_ust, ti_dst, ti_area_ratio_in; s=700.0)
    # Niayifar and Porte Agel 2015, 2016 using smooth max on area TI ratio

    # calculate axial induction based on the Ct value
    axial_induction_ust = _ct_to_axial_ind_func(ct_ust)

    # calculate BPA spread parameters Bastankhah and Porte Agel 2014
    beta = 0.5*((1.0 + sqrt(1.0 - ct_ust))/sqrt(1.0 - ct_ust))
    epsilon = 0.2*sqrt(beta)

    # calculate wake spread for TI calcs
    sigma = kstar_ust*x + d_ust*epsilon
    d_w = 4.0*sigma

    # calculate wake overlap ratio
    wake_overlap = overlap_area_func(0.0, h_dst, d_dst, delta_y, h_ust, d_w)

    # initialize the wake/rotor area overlap ratio
    ti_area_ratio = 0.0

    # only include turbines with area overlap in the softmax
    if wake_overlap > 0.0

        # Calculate the turbulence added to the inflow of the downstream turbine by the
        # wake of the upstream turbine
        ti_added = 0.73*(axial_induction_ust^0.8325)*(ti_ust^0.0325)*((x/d_ust)^(-0.32))

        rotor_area_dst = 0.25*pi*d_dst^2
        ti_area_ratio_tmp = ti_added*(wake_overlap/rotor_area_dst)

        # Run through the smooth max to get an approximation of the true max TI area ratio
        ti_area_ratio = smooth_max(ti_area_ratio_in, ti_area_ratio_tmp, s=s)

        # Calculate the total turbulence intensity at the downstream turbine based on
        # the result of the smooth max function
        ti_dst = sqrt(ti_amb^2 + ti_area_ratio^2)

    end

    return ti_dst, ti_area_ratio

end

"""
    calculate_local_ti(turbine_x, turbine_y, ambient_ti, rotor_diameter, hub_height, turbine_yaw, turbine_local_ti, sorted_turbine_index,
    turbine_inflow_velcities, turbine_ct, ti_model::LocalTIModelMaxTI; turbine_id=1, tol=1E-6)

Returns local turbulence intensity calculated using Niayifar and Porte Agel 2015, 2016 using smooth max on area TI ratio

# Arguments
- `turbine_x::Array{Float,nTurbines}`: turbine wind direction locations in the wind direction
    reference frame
- `turbine_y::Array{Float,nTurbines}`: turbine cross wind locations in the wind direction
    reference frame
- `ambient_ti::Float`: ambient turbulence intensity
- `rotor_diameter::Array{Float,nTurbines}`: rotor diameters of all turbines
- `hub_height::Array{Float,nTurbines}`: hub heights of all turbines relative to the ground
- `turbine_yaw::Array{Float,nTurbines}`: yaw of all turbines for the current wind state in radians
- `turbine_local_ti::Array{Float,nTurbines}`: local turbulence intensity of all turbines for the current wind state`
- `sorted_turbine_index::Array{Float,nTurbines}`: turbine north-south locations in the
    global reference frame
- `turbine_inflow_velcities::Array{Float,nTurbines}`: effective inflow wind speed at each turbine for given state
- `turbine_ct::Array{Float,nTurbines}`: thrust coefficient of each turbine for the given state
- `ti_model::LocalTIModelMaxTI`: contains a struct defining the desired turbulence intensity model, no local TI in this case
- `turbine_id::Int`: index of wind turbine of interest. Provide 1 as default.
- `tol::Float`: How far upstream a turbine should be before being included in TI calculations
"""
function calculate_local_ti(turbine_x, turbine_y, ambient_ti, rotor_diameter, hub_height, turbine_yaw, turbine_local_ti, sorted_turbine_index,
                    turbine_inflow_velcities, turbine_ct, ti_model::LocalTIModelMaxTI; turbine_id=1, tol=1E-6)

    # calculate local turbulence intensity at turbI

    # initialize the ri_dst and ti_area_ratio to 0.0 for current turbine
    ti_area_ratio = 0.0
    ti_dst = copy(ambient_ti)

    # extract downstream turbine information
    d_dst = rotor_diameter[turbine_id]
    h_dst = hub_height[turbine_id]

    # extract the number of turbines
    nturbines = length(rotor_diameter)

    # loop over upstream turbines
    for u=1:nturbines

        # get index of upstream turbine
        turb = sorted_turbine_index[u]

        # skip turbine's influence on itself
        if turb == turbine_id; continue; end

        # calculate downstream distance between wind turbines
        x = turbine_x[turbine_id] - turbine_x[turb]

        if x > tol

            # extract state and design info for current upstream turbine
            d_ust = rotor_diameter[turb]
            h_ust = hub_height[turb]
            yaw_ust = turbine_yaw[turb]
            ti_ust = turbine_local_ti[turb]
            ct_ust = turbine_ct[turb]

            # determine the far-wake onset location
            astar = ti_model.astar
            bstar = ti_model.bstar
            x0 = _gauss_yaw_potential_core(d_ust, yaw_ust, ct_ust, astar, ti_ust, bstar)

            # calculate the distance from the onset of far-wake
            deltax0 = x - x0

            # calculate wake spread rate for current upstream turbine
            kstar_ust = _k_star_func(ti_ust,ti_model.k1,ti_model.k2)

            # calculate horizontal and vertical spread standard deviations
            sigmay = sigmaz = _gauss_yaw_spread(x, x0, kstar_ust, d_ust, yaw_ust)

            # determine the initial wake angle at the onset of far wake
            theta0 = _bpa_theta_0(yaw_ust, ct_ust)

            # horizontal cross-wind wake displacement from hub
            wake_offset = _bpa_deflection(d_ust, ct_ust, yaw_ust, kstar_ust, kstar_ust, sigmay, sigmaz, theta0, x0)

            # cross wind distance from point location to upstream turbine wake center
            delta_y = turbine_y[turbine_id]  - (turbine_y[turb] + wake_offset)

            # save ti_area_ratio and ti_dst to new memory locations to avoid
            # aliasing during differentiation
            ti_area_ratio_tmp = deepcopy(ti_area_ratio)

            # update local turbulence intensity
            ti_dst, ti_area_ratio = _niayifar_added_ti_function(x, d_dst, d_ust, h_ust, h_dst, ct_ust, kstar_ust, delta_y, ambient_ti, ti_ust, ti_dst, ti_area_ratio_tmp)

        end

    end

    return ti_dst

end

"""
    GaussianTI(loc,turbine_x, turbine_y, rotor_diameter, hub_height, turbine_ct,
        sorted_turbine_index, ambient_ti; div_sigma=2.5, div_ti=1.2)

Calculate local turbulence intensity based on "On wake modeling, wind-farm gradients and AEP
    predictions at the Anholt wind farm" by Pena Diaz, Alfredo; Hansen, Kurt Schaldemose;
    Ott, Søren; van der Laan, Paul ??

# Arguments
- `loc::Array{Float,3}`: [x,y,z] location of point of interest in wind direction ref. frame
- `turbine_x::Array{Float,nTurbines}`: turbine wind direction locations in the wind direction
    reference frame
- `turbine_y::Array{Float,nTurbines}`: turbine cross wind locations in the wind direction
    reference frame
- `rotor_diameter::Array{Float,nTurbines}`: rotor diameters of all turbines
- `hub_height::Array{Float,nTurbines}`: hub heights of all turbines relative to the ground
- `turbine_ct::Array{Float,nTurbines}`: thrust coefficient of each turbine for the given state
- `sorted_turbine_index::Array{Float,nTurbines}`: turbine north-south locations in the
    global reference frame
- `ambient_ti::Float`: ambient turbulence intensity
- `div_sigma::Float`: ?
- `div_ti::Float`: ?
"""
function GaussianTI(loc,turbine_x, turbine_y, rotor_diameter, hub_height, turbine_ct, sorted_turbine_index, ambient_ti; div_sigma=2.5, div_ti=1.2)

    added_ti = 0.0
    e = 1.0*ambient_ti^0.1
    nturbines = length(turbine_x)

    for u=1:nturbines

        # get index of upstream turbine
        turb = sorted_turbine_index[u]

        # calculate downstream distance between wind turbines
        dx = loc[1] - turbine_x[turb]

        if dx > 1e-6
            dy = loc[2] - turbine_y[turb]
            dz = loc[3] - hub_height[turb]
            r = sqrt(dy^2 + dz^2)
            ct = turbine_ct[turb]

            kstar = 0.11*ct^1.07*ambient_ti^0.2
            epsilon = 0.23*ct^-0.25*ambient_ti^0.17
            d = 2.3*ct^-1.2
            f = 0.7*ct^-3.2*ambient_ti^-0.45

            dist = 0.5
            if r/rotor_diameter[turb] <= dist
                k1 = cos(pi/2.0*(r/rotor_diameter[turb]-dist))^2
                k2 = cos(pi/2.0*(r/rotor_diameter[turb]+dist))^2
            else
                k1 = 1.0
                k2 = 0.0
            end

            sigma = kstar*dx + epsilon*rotor_diameter[turb]
            if dz >= 0.0
                delta = 0.0
            else
                delta = ambient_ti*sin(pi*dz/hub_height[turb])^2
            end

            #2.5 for low TI 2.0 for high TI
            sigma = sigma/div_sigma

            new_ex = -2.0 #orig -2.0
            p1 = 1.0/(d + e*dx/rotor_diameter[turb] + f*(1.0+dx/rotor_diameter[turb])^new_ex)
            p2 = k1*exp(-(r-rotor_diameter[turb]/2.0)^2/(2.0*sigma^2)) + k2*exp(-(r+rotor_diameter[turb]/2.0)^2/(2.0*sigma^2))
            dI = p1*p2 - delta
            #1.2 for low TI 2.0 for high TI
            added_ti += dI/div_ti
        end
    end
    return ambient_ti + added_ti
end

# function GaussianTI_stanley(loc,windfarm,windfarmstate,ambient_ti)
#
#     added_ti = 0.0
#     e = 1.0*ambient_ti^0.1
#     nturbines = length(windfarm.turbine_x)
#
#     for u=1:nturbines
#
#         # get index of upstream turbine
#         turb = windfarmstate.sorted_turbine_index[u]
#
#         # calculate downstream distance between wind turbines
#         dx = loc[1] - windfarmstate.turbine_x[turb]
#
#         if dx > 1e-6
#             turbine_type = windfarm.turbine_definition_ids[turb]
#             rotor_diameter = windfarm.turbine_definitions[turbine_type].rotor_diameter[1]
#             hub_height = windfarm.turbine_definitions[turbine_type].hub_height[1]
#             dy = loc[2] - windfarmstate.turbine_y[turb]
#             dz = loc[3] - hub_height
#             r = sqrt(dy^2 + dz^2)
#
#             ct = windfarmstate.turbine_ct[turb]
#             kstar = 0.11*ct^1.07*ambient_ti^0.2
#             epsilon = 0.23*ct^-0.25*ambient_ti^0.17
#             d = 2.3*ct^-1.2
#             f = 0.7*ct^-3.2*ambient_ti^-0.45
#
#             # dist = 0.0
#             # if r/rotor_diameter <= dist
#             #     k1 = cos(pi/2.0*(r/rotor_diameter-dist))^2
#             #     k2 = cos(pi/2.0*(r/rotor_diameter+dist))^2
#             # else
#             #     k1 = 1.0
#             #     k2 = 0.0
#             # end
#             k1 = 1.0
#             k2 = 0.0
#
#             sigma = kstar*dx + epsilon*rotor_diameter
#             if dz >= 0.0
#                 delta = 0.0
#             else
#                 delta = ambient_ti*sin(pi*dz/hub_height)^2
#             end
#
#
#             sigma = sigma*0.5
#             p1 = 1.0/(d + e*dx/rotor_diameter + f*(1.0+dx/rotor_diameter)^-2.0)
#             # p2 = k1*exp(-(r-rotor_diameter/2.0)^2/(2.0*sigma^2)) + k2*exp(-(r+rotor_diameter/2.0)^2/(2.0*sigma^2))
#             p2 = exp(-(r/2.0)^2/(2.0*sigma^2))
#             dI = p1*p2 - delta
#             # if r < rotor_diameter*4.0/5.0
#             added_ti += dI
#             # end
#         end
#     end
#     return ambient_ti + added_ti
# end


text:
using Printf

problemName = "uniaxialExtension_paper"

generateBoolean   = true ;

outDir = "./output" ;

# ------------------------------------------------------------------------------
# geometria:  verificar coherencia .geo
ladoSecc = 0.05
largo    = 1

#boxGrid  = [ [ 0.1 ,          0.9 ], [ ladoSecc*.1,   ladoSecc*.5], [ ladoSecc*.1,   ladoSecc*.5 ] ]

boxGrid  = [ 0.1           0.9  ; ladoSecc*.1   ladoSecc*.5 ;  ladoSecc*.1   ladoSecc*.5 ] ;

# con este roi se asegura que solo el eje central de voxels queda en dentro
minsROI = [ 0.4,  ladoSecc*0.25,  ladoSecc*0.25 ]
maxsROI = [ 0.6,  ladoSecc*0.35,  ladoSecc*0.35 ]

voxelNums   = [ 100      ,  3        , 3        ]
# voxelNums   = [ 20      ,  3        , 3        ]' ;

pRef             = 2     # reference young modulus
nu               = 0.4
nLoadSteps       = 1
levelerror       = 0.0
plotParamsVector = [ 3 ]
tolDeltaP        = 1e-4   # converge criterion

realImagesBoolean = 0

tolDeltaP = 1e-5

# numero de parametros en los cuales se hace loop en las corridas
# numParsLoop = 3 ;

cpenal  = 1e5

#~ runsData         = cell ( numParsLoop, 1) ;

# identificacion: funcion / tension / Pini

#~ paramsVals{1} = [ 2 ]   ; paramsVals{2} = [ .1 ] ;
#~ paramsVals{3} = [ .4 ]' ; paramsVals{4} = [ 0 ]' ;

paramsVals = Vector{Vector{Float64}}(undef,4)

# intensity functions
paramsVals[1] = [ 2, 3 ]

# applied tension
paramsVals[2] = [ .1, .7 ]

# initial material parameter modulus
paramsVals[3] = [ .4, 3 ]

# error level
paramsVals[4] = [ 0 ]

# iniMet = 0; endMet = 0;  flagsPyOct = [1 1 ] ; % validacion
# iniMet = 2; endMet = 3; flagsPyOct = [] ;  % identif python
methodNums = [ 4, 5]
flagsPyOct = [] ; # identif octave

strPaper = "/home/jor/work/udelar/investigacion/papers/paper_IBHMI_repo/texs"

# --- genera tex con valores de corridas de ejemplo ---
paramsTexFile = strPaper * "/paramsTexs/" * problemName * "_Params.tex"

print(fileParams)

paramFID     = open( paramsTexFile, "w" ) ;
@printf( paramFID ,"\\def\\exUnoERef{%12.2e} \n", pRef )
@printf( paramFID ,"\\def\\exUnonu{%3.2f} \n"   , nu   )

@printf( paramFID, "\\def\\exUnoNumVoxY{%4i} \n", voxelNums[2] )
@printf( paramFID, "\\def\\exUnoNumVoxX{%4i} \n", voxelNums[1] )
@printf( paramFID, "\\def\\exUnoNumVoxZ{%4i} \n", voxelNums[3] )

@printf( paramFID, "\\def\\exUnoLength{%1.0f} \n", largo );
@printf( paramFID, "\\def\\exUnoLadoSecc{%5.2f} \n", ladoSecc );

@printf( paramFID, "\\def\\exUnoMinXROI{%5.3f} \n", minsROI[1] )
@printf( paramFID, "\\def\\exUnoMinYROI{%5.3f} \n", minsROI[2] )
@printf( paramFID, "\\def\\exUnoMinZROI{%5.3f} \n", minsROI[3] )

@printf( paramFID, "\\def\\exUnoMaxXROI{%5.3f} \n", maxsROI[1] )
@printf( paramFID, "\\def\\exUnoMaxYROI{%5.3f} \n", maxsROI[2] )
@printf( paramFID, "\\def\\exUnoMaxZROI{%5.3f} \n", maxsROI[3] )

@printf( paramFID, "\\def\\exUnoMinXFrame{%5.3f} \n", boxGrid[1,1] )
@printf( paramFID, "\\def\\exUnoMinYFrame{%5.3f} \n", boxGrid[2,1] )
@printf( paramFID, "\\def\\exUnoMinZFrame{%5.3f} \n", boxGrid[3,1] )

@printf( paramFID, "\\def\\exUnoMaxXFrame{%5.3f} \n", boxGrid[1,2] )
@printf( paramFID, "\\def\\exUnoMaxYFrame{%5.3f} \n", boxGrid[2,2] )
@printf( paramFID, "\\def\\exUnoMaxZFrame{%5.3f} \n", boxGrid[3,2] )

close( paramFID )

# run MatIPS
include("/home/jor/work/udelar/codigos/MatIPS_repo/src/julia/MatIPS.jl")


text:
using Revise, LinearAlgebra, StaticArrays, PyCall  #FFTW, BenchmarkTools, LinearMaps, IterativeSolvers, Roots, GeometryPrimitives
using OptiMode
SHM3 = SHermitianCompact{3,Float64,6}
LinearAlgebra.ldiv!(c,A::LinearMaps.LinearMap,b) = mul!(c,A',b)
# MPB solve for reference
mp = pyimport("meep")
mpb = pyimport("meep.mpb")

w           = 1.7
t_core      = 0.7
edge_gap    = 0.5               # μm
n_core      = 2.4
n_subs      = 1.4
λ           = 1.55              # μm

nk          = 10
Δx       = 6.                    # μm
Δy       = 4.                    # μm
n_bands     = 1
res         = 16

n_guess = 0.9 * n_core
n_min = n_subs
n_max = n_core
t_subs = (Δy -t_core - edge_gap )/2.
c_subs_y = -Δy/2. + edge_gap/2. + t_subs/2.

# Set up MPB modesolver, use find-k to solve for one eigenmode `H` with prop. const. `k` at specified temporal freq. ω

k_pts = mp.interpolate(nk, [mp.Vector3(0.05, 0, 0), mp.Vector3(0.05*nk, 0, 0)] )
lat = mp.Lattice(size=mp.Vector3(Δx, Δy,0))
# core = mp.Block(size=mp.Vector3(w,t_core,10.0,),
#                     center=mp.Vector3(0,0,0,),
#                     material=mp.Medium(index=n_core),
#                    )
verts = [ mp.Vector3(-w/2., -t_core/2., -5.),mp.Vector3(w, 2*t_core, -5.), mp.Vector3(w, -t_core/2., -5.)  ]
core = mp.Prism(verts,
                   10.0,
                   axis=mp.Vector3(0.,0.,1.),
                   material=mp.Medium(index=n_core),
                  )
subs = mp.Block(size=mp.Vector3(Δx-edge_gap, t_subs , 10.0),
                center=mp.Vector3(0, c_subs_y, 0),
                material=mp.Medium(index=n_subs),
                )

ms = mpb.ModeSolver(geometry_lattice=lat,
                    geometry=[core,subs],
                    k_points=k_pts,
                    resolution=res,
                    num_bands=n_bands,
                    default_material=mp.vacuum)

ms.init_params(mp.NO_PARITY, false)
ε_mean_mpb = ms.get_epsilon()
nx_mpb = size(ε_mean_mpb)[1]
ny_mpb = size(ε_mean_mpb)[2]
dx_mpb = Δx / nx_mpb
dy_mpb = Δy / ny_mpb
x_mpb = (dx_mpb .* (0:(nx_mpb-1))) .- Δx/2. #(Δx/2. - dx_mpb)
y_mpb = (dy_mpb .* (0:(ny_mpb-1))) .- Δy/2. #(Δy/2. - dy_mpb)
k_mpb = ms.find_k(mp.NO_PARITY,             # parity (meep parity object)
                  1/λ,                    # ω at which to solve for k
                  1,                        # band_min (find k(ω) for bands
                  n_bands,                        # band_max  band_min:band_max)
                  mp.Vector3(0, 0, 1),      # k direction to search
                  1e-4,                     # fractional k error tolerance
                  n_guess/λ,              # kmag_guess, |k| estimate
                  n_min/λ,                # kmag_min (find k in range
                  n_max/λ,               # kmag_max  kmag_min:kmag_max)
)

neff_mpb = k_mpb * λ
ng_mpb = 1 / ms.compute_one_group_velocity_component(mp.Vector3(0, 0, 1), 1)
e_mpb = reshape(ms.get_efield(1),(nx_mpb,ny_mpb,3))
d_mpb = reshape(ms.get_dfield(1),(nx_mpb,ny_mpb,3))
h_mpb = reshape(ms.get_hfield(1),(nx_mpb,ny_mpb,3))
S_mpb = reshape(ms.get_poynting(1),(nx_mpb,ny_mpb,3))
U_mpb = reshape(ms.get_tot_pwr(1),(nx_mpb,ny_mpb))
# H_mpb = vec(reshape(ms.get_eigenvectors(1,1),(nx_mpb*ny_mpb,2)))
H_mpb_raw = reshape(ms.get_eigenvectors(1,1), (ny_mpb,nx_mpb,2))
H_mpb = zeros(ComplexF64,(3,nx_mpb,ny_mpb,1)) #Array{ComplexF64,4}(undef,(3,ny_mpb,ny_mpb,1))
for i=1:nx_mpb
    for j=1:ny_mpb
        H_mpb[1,i,j,1] = H_mpb_raw[j,i,1]
        H_mpb[2,i,j,1] = H_mpb_raw[j,i,2]
    end
end
SHM3 = SHermitianCompact{3,Float64,6}
z_mpb = [0.,]
nz_mpb = 1
ε⁻¹_mpb = [SHM3([ms.get_epsilon_inverse_tensor_point(mp.Vector3(x_mpb[i],y_mpb[j],z_mpb[k]))[a][b] for a=1:3,b=1:3]) for i=1:nx_mpb,j=1:ny_mpb,k=1:nz_mpb]
ε_mpb = [SHM3(inv(ε⁻¹_mpb[i,j,k])) for i=1:nx_mpb,j=1:ny_mpb,k=1:nz_mpb]
e_mpb = [SVector(e_mpb[i,j,:]...) for i=1:nx_mpb,j=1:ny_mpb]
d_mpb = [SVector(d_mpb[i,j,:]...) for i=1:nx_mpb,j=1:ny_mpb]
h_mpb = [SVector(h_mpb[i,j,:]...) for i=1:nx_mpb,j=1:ny_mpb]
S_mpb = [SVector(S_mpb[i,j,:]...) for i=1:nx_mpb,j=1:ny_mpb]

Nx = nx_mpb
Ny = ny_mpb
k = SVector(0.,0.,k_mpb[1])
g = MaxwellGrid(Δx,Δy,Nx,Ny)
ds = MaxwellData(k,g);

##

# Nz = 1
# Neigs = 1
# N = *(Nx,Ny,Nz)
# s = ridge_wg(w,t_core,edge_gap,n_core,n_subs,g)
# ε⁻¹ = εₛ⁻¹(s,g)
# ε = SHM3.(inv.(ε⁻¹))

# U_mpb_plot = heatmap(
#             x_mpb,
#             y_mpb,
#             transpose(U_mpb),
#             c=cgrad(:cherry),
#             aspect_ratio=:equal,
#             legend=false,
#             colorbar = true,
#             # clim = ( 0., max(ε_mpb...) ),
# )

# H1 = randn(ComplexF64,3,Neigs,Nx,Ny,Nz) #Array{ComplexF64}(undef,3,Neigs,Nx,Ny,Nz)
# Iₛₚ = CartesianIndices(size(H1)[3:5])
# HSA = [SVector{3,ComplexF64}(randn(ComplexF64,3)...) for p=1:Neigs,i=1:Nx,j=1:Ny,k=1:Nz];
# KpG = [KVec(SVector{3,Float64}(randn(Float64,3)...)) for p=1:1,i=1:Nx,j=1:Ny,k=1:Nz];
# ε⁻¹ = [SHM3(randn(Float64,6)) for p=1:1,i=1:Nx,j=1:Ny,k=1:Nz];
# P = plan_fft!(reinterpret(ComplexF64,HSA),(2:4));
# iP = plan_ifft!(reinterpret(ComplexF64,HSA),(2:4));
# HSV = vec(copy(HSA))

# HA = [HSA[p,i,j,k][d] for p=1:Neigs,i=1:Nx,j=1:Ny,k=1:Nz,d=1:3];

# 3D FFTW benchmarks for Nx=64,Ny=128,Nz=256 ##

# CPU #
# @btime fft!(H1,(3:5)); # 601.221 ms (28 allocations: 2.89 KiB)
# @btime fft!(H2,(1:3)); # 970.900 ms (28 allocations: 2.89 KiB)
# @btime fft!(reinterpret(ComplexF64,$HSA),(2:4)); #541.107 ms (29 allocations: 2.78 KiB)
# size(reinterpret(ComplexF64,H2M)) # (192, 128, 256, 5) --> mixes spatial and non-spatial dimension, non-spatial dims should go first

# non-Mutating Operators

# function t2c(Hin::AbstractArray{ComplexF64,4},ds::MaxwellData)::AbstractArray{ComplexF64,4}
#     Hout = Array{ComplexF64}(undef,(3,size(Hin)[2:end]...))
#     for i=1:ds.Nx,j=1:ds.Ny,k=1:ds.Nz
#         @inbounds scale = ds.kpG[i,j,k].mag
#         @inbounds Hout[1,i,j,k] = ( Hin[1,i,j,k] * ds.kpG[i,j,k].m[1] + Hin[2,i,j,k] * ds.kpG[i,j,k].n[1] ) * scale
#         @inbounds Hout[2,i,j,k] = ( Hin[1,i,j,k] * ds.kpG[i,j,k].m[2] + Hin[2,i,j,k] * ds.kpG[i,j,k].n[2] ) * scale
#         @inbounds Hout[3,i,j,k] = ( Hin[1,i,j,k] * ds.kpG[i,j,k].m[3] + Hin[2,i,j,k] * ds.kpG[i,j,k].n[3] ) * scale
#     end
#     return Hout
# end

# function c2t(Hin::AbstractArray{ComplexF64,4},ds::MaxwellData)::AbstractArray{ComplexF64,4}
#     Hout = Array{ComplexF64}(undef,(2,size(Hin)[2:end]...))
#     for i=1:ds.Nx,j=1:ds.Ny,k=1:ds.Nz
#         @inbounds Hout[1,i,j,k] =  Hin[1,i,j,k] * ds.kpG[i,j,k].m[1] + Hin[2,i,j,k] * ds.kpG[i,j,k].m[2] + Hin[3,i,j,k] * ds.kpG[i,j,k].m[3]
#         @inbounds Hout[2,i,j,k] =  Hin[1,i,j,k] * ds.kpG[i,j,k].n[1] + Hin[2,i,j,k] * ds.kpG[i,j,k].n[2] + Hin[3,i,j,k] * ds.kpG[i,j,k].n[3]
#     end
#     return Hout
# end

# function zcross_t2c(Hin::AbstractArray{ComplexF64,4},ds::MaxwellData)::AbstractArray{ComplexF64,4}
#     Hout = zeros(ComplexF64,(3,size(Hin)[2:end]...))
#     for i=1:ds.Nx,j=1:ds.Ny,k=1:ds.Nz
#         @inbounds Hout[1,i,j,k] = -Hin[1,i,j,k] * ds.kpG[i,j,k].m[2] - Hin[2,i,j,k] * ds.kpG[i,j,k].n[2]
#         @inbounds Hout[2,i,j,k] =  Hin[1,i,j,k] * ds.kpG[i,j,k].m[1] + Hin[2,i,j,k] * ds.kpG[i,j,k].n[1]
#     end
#     return Hout
# end

# function kcross_t2c(Hin::AbstractArray{ComplexF64,4},ds::MaxwellData)::AbstractArray{ComplexF64,4}
#     Hout = Array{ComplexF64}(undef,(3,size(Hin)[2:end]...))
#     for i=1:ds.Nx,j=1:ds.Ny,k=1:ds.Nz
#         @inbounds scale = -ds.kpG[i,j,k].mag
#         @inbounds Hout[1,i,j,k] = ( Hin[1,i,j,k] * ds.kpG[i,j,k].n[1] - Hin[2,i,j,k] * ds.kpG[i,j,k].m[1] ) * scale
#         @inbounds Hout[2,i,j,k] = ( Hin[1,i,j,k] * ds.kpG[i,j,k].n[2] - Hin[2,i,j,k] * ds.kpG[i,j,k].m[2] ) * scale
#         @inbounds Hout[3,i,j,k] = ( Hin[1,i,j,k] * ds.kpG[i,j,k].n[3] - Hin[2,i,j,k] * ds.kpG[i,j,k].m[3] ) * scale
#     end
#     return Hout
# end

# function kcross_c2t(Hin::AbstractArray{ComplexF64,4},ds::MaxwellData)::AbstractArray{ComplexF64,4}
#     Hout = Array{ComplexF64}(undef,(2,size(Hin)[2:end]...))
#     for i=1:ds.Nx,j=1:ds.Ny,k=1:ds.Nz
#         @inbounds scale = ds.kpG[i,j,k].mag
#         @inbounds at1 = Hin[1,i,j,k] * ds.kpG[i,j,k].m[1] + Hin[2,i,j,k] * ds.kpG[i,j,k].m[2] + Hin[3,i,j,k] * ds.kpG[i,j,k].m[3]
#         @inbounds at2 = Hin[1,i,j,k] * ds.kpG[i,j,k].n[1] + Hin[2,i,j,k] * ds.kpG[i,j,k].n[2] + Hin[3,i,j,k] * ds.kpG[i,j,k].n[3]
#         @inbounds Hout[1,i,j,k] =  -at2 * scale
#         @inbounds Hout[2,i,j,k] =  at1 * scale
#     end
#     return Hout
# end

# function kcrossinv_t2c(Hin::AbstractArray{ComplexF64,4},ds::MaxwellData)::AbstractArray{ComplexF64,4}
#     Hout = Array{ComplexF64}(undef,(3,size(Hin)[2:end]...))
#     for i=1:ds.Nx,j=1:ds.Ny,k=1:ds.Nz
#         @inbounds scale = 1 / ds.kpG[i,j,k].mag
#         @inbounds Hout[1,i,j,k] = ( Hin[1,i,j,k] * ds.kpG[i,j,k].n[1] - Hin[2,i,j,k] * ds.kpG[i,j,k].m[1] ) * scale
#         @inbounds Hout[2,i,j,k] = ( Hin[1,i,j,k] * ds.kpG[i,j,k].n[2] - Hin[2,i,j,k] * ds.kpG[i,j,k].m[2] ) * scale
#         @inbounds Hout[3,i,j,k] = ( Hin[1,i,j,k] * ds.kpG[i,j,k].n[3] - Hin[2,i,j,k] * ds.kpG[i,j,k].m[3] ) * scale
#     end
#     return Hout
# end

# function kcrossinv_c2t(Hin::AbstractArray{ComplexF64,4},ds::MaxwellData)::AbstractArray{ComplexF64,4}
#     Hout = Array{ComplexF64}(undef,(2,size(Hin)[2:end]...))
#     for i=1:ds.Nx,j=1:ds.Ny,k=1:ds.Nz
#         @inbounds scale = -1 / ds.kpG[i,j,k].mag
#         @inbounds at1 = Hin[1,i,j,k] * ds.kpG[i,j,k].m[1] + Hin[2,i,j,k] * ds.kpG[i,j,k].m[2] + Hin[3,i,j,k] * ds.kpG[i,j,k].m[3]
#         @inbounds at2 = Hin[1,i,j,k] * ds.kpG[i,j,k].n[1] + Hin[2,i,j,k] * ds.kpG[i,j,k].n[2] + Hin[3,i,j,k] * ds.kpG[i,j,k].n[3]
#         @inbounds Hout[1,i,j,k] =  -at2 * scale
#         @inbounds Hout[2,i,j,k] =  at1 * scale
#     end
#     return Hout
# end

# function ε⁻¹_dot(Hin::AbstractArray{ComplexF64,4},ε⁻¹::Array{SHermitianCompact{3,Float64,6},3},ds::MaxwellData)::AbstractArray{ComplexF64,4}
#     Hout = similar(Hin)
#     for i=1:ds.Nx,j=1:ds.Ny,k=1:ds.Nz
#         @inbounds Hout[1,i,j,k] =  ε⁻¹[i,j,k][1,1]*Hin[1,i,j,k] + ε⁻¹[i,j,k][2,1]*Hin[2,i,j,k] + ε⁻¹[i,j,k][3,1]*Hin[3,i,j,k]
#         @inbounds Hout[2,i,j,k] =  ε⁻¹[i,j,k][1,2]*Hin[1,i,j,k] + ε⁻¹[i,j,k][2,2]*Hin[2,i,j,k] + ε⁻¹[i,j,k][3,2]*Hin[3,i,j,k]
#         @inbounds Hout[3,i,j,k] =  ε⁻¹[i,j,k][1,3]*Hin[1,i,j,k] + ε⁻¹[i,j,k][2,3]*Hin[2,i,j,k] + ε⁻¹[i,j,k][3,3]*Hin[3,i,j,k]
#         # @inbounds Hout[1,i,j,k] =  ε⁻¹[i,j,k][1,1]*Hin[1,i,j,k] + ε⁻¹[i,j,k][1,2]*Hin[2,i,j,k] + ε⁻¹[i,j,k][1,3]*Hin[3,i,j,k]
#         # @inbounds Hout[2,i,j,k] =  ε⁻¹[i,j,k][2,1]*Hin[1,i,j,k] + ε⁻¹[i,j,k][2,2]*Hin[2,i,j,k] + ε⁻¹[i,j,k][2,3]*Hin[3,i,j,k]
#         # @inbounds Hout[3,i,j,k] =  ε⁻¹[i,j,k][3,1]*Hin[1,i,j,k] + ε⁻¹[i,j,k][3,2]*Hin[2,i,j,k] + ε⁻¹[i,j,k][3,3]*Hin[3,i,j,k]
#     end
#     return Hout
# end

# function ε_dot_approx(Hin::AbstractArray{ComplexF64,4},ε⁻¹::Array{SHermitianCompact{3,Float64,6},3},ds::MaxwellData)::AbstractArray{ComplexF64,4}
#     Hout = similar(Hin)
#     for i=1:ds.Nx,j=1:ds.Ny,k=1:ds.Nz
#         @inbounds ε_ave = 3 / tr(ε⁻¹[i,j,k])
#         @inbounds Hout[1,i,j,k] =  ε_ave * Hin[1,i,j,k]
#         @inbounds Hout[2,i,j,k] =  ε_ave * Hin[2,i,j,k]
#         @inbounds Hout[3,i,j,k] =  ε_ave * Hin[3,i,j,k]
#     end
#     return Hout
# end

# function M(Hin::AbstractArray{ComplexF64,4},ε⁻¹::Array{SHermitianCompact{3,Float64,6},3},ds::MaxwellData)::Array{ComplexF64,4}
#     d = ds.𝓕 * kcross_t2c(Hin,ds);
#     e = ε⁻¹_dot(d,ε⁻¹,ds); # (-1/(π)) .*
#     kcross_c2t(ds.𝓕⁻¹ * e,ds)
# end

# function M(Hin::AbstractArray{ComplexF64,1},ε⁻¹::Array{SHermitianCompact{3,Float64,6},3},ds::MaxwellData)::Array{ComplexF64,1}
#     HinA = reshape(Hin,(2,ds.Nx,ds.Ny,ds.Nz))
#     HoutA = M(HinA,ε⁻¹,ds)
#     return vec(HoutA)
# end

# M̂(ε⁻¹::Array{SHermitianCompact{3,Float64,6},3},ds::MaxwellData) = LinearMap{ComplexF64}(H::AbstractArray{ComplexF64,1} -> M(H,ε⁻¹,ds)::AbstractArray{ComplexF64,1},(2*ds.Nx*ds.Ny*ds.Nz),ishermitian=true,ismutating=false)

# function P(Hin::AbstractArray{ComplexF64,4},ε⁻¹::Array{SHermitianCompact{3,Float64,6},3},ds::MaxwellData)::Array{ComplexF64,4}
#     e = ds.𝓕⁻¹ * kcrossinv_t2c(Hin,ds);
#     d = ε_dot_approx(e,ε⁻¹,ds); # (-1/(π)) .*
#     kcrossinv_c2t(ds.𝓕 * d,ds)
# end

# function P(Hin::AbstractArray{ComplexF64,1},ε⁻¹::Array{SHermitianCompact{3,Float64,6},3},ds::MaxwellData)::Array{ComplexF64,1}
#     HinA = reshape(Hin,(2,ds.Nx,ds.Ny,ds.Nz))
#     HoutA = P(HinA,ε⁻¹,ds)
#     return vec(HoutA)
# end

# P̂(ε⁻¹::Array{SHermitianCompact{3,Float64,6},3},ds::MaxwellData) = LinearMap{ComplexF64}(H::AbstractArray{ComplexF64,1} -> P(H,ε⁻¹,ds)::AbstractArray{ComplexF64,1},(2*ds.Nx*ds.Ny*ds.Nz),ishermitian=true,ismutating=false)

# function Mₖ(Hin::AbstractArray{ComplexF64,4},ε⁻¹::Array{SHermitianCompact{3,Float64,6},3},ds::MaxwellData)::Array{ComplexF64,4}
#     d = ds.𝓕 * zcross_t2c(Hin,ds);
#     e = ε⁻¹_dot(d,ε⁻¹,ds); # (-1/(π)) .*
#     kcross_c2t(ds.𝓕⁻¹ * e,ds)
# end

# function Mₖ(Hin::AbstractArray{ComplexF64,1},ε⁻¹::Array{SHermitianCompact{3,Float64,6},3},ds::MaxwellData)::Array{ComplexF64,1}
#     HinA = reshape(Hin,(2,ds.Nx,ds.Ny,ds.Nz))
#     HoutA = Mₖ(HinA,ε⁻¹,ds)
#     return -vec(HoutA)
# end

# M̂ₖ(ε⁻¹::Array{SHermitianCompact{3,Float64,6},3},ds::MaxwellData) = LinearMap{ComplexF64}(H::AbstractArray{ComplexF64,1} -> Mₖ(H,ε⁻¹,ds)::AbstractArray{ComplexF64,1},(2*ds.Nx*ds.Ny*ds.Nz),ishermitian=true,ismutating=false)

# function solve_ω(k::SVector{3},ε⁻¹::Array{SHM3,3},g::MaxwellGrid;neigs=1,eigind=1,maxiter=3000,tol=1e-7)
#     ds = MaxwellData(k,g)
#     res = IterativeSolvers.lobpcg(M̂(ε⁻¹,ds),false,neigs;P=P̂(ε⁻¹,ds),maxiter,tol)
#     H =  res.X #[:,eigind]                       # eigenmode wavefn. magnetic fields in transverse pol. basis
#     ω =  √(real(res.λ[eigind]))                     # eigenmode temporal freq.,  neff = kz / ω, kz = k[3]
#     # ωₖ =   real( ( H' * M̂ₖ(ε⁻¹,ds) * H )[1]) / ω       # ωₖ/∂kz = group velocity = c / ng, c = 1 here
#     ωₖ =   real( ( H[:,eigind]' * M̂ₖ(ε⁻¹,ds) * H[:,eigind] )[1]) / ω       # ωₖ/∂kz = group velocity = c / ng, c = 1 here
#     return H, ω, ωₖ
# end

# function kpG(k::SVector{3,Float64},g::MaxwellGrid)::Array{KVec,3}
#     [KVec(k-gx-gy-gz) for gx=g.gx, gy=g.gy, gz=g.gz]
# end

# """
# modified solve_ω version for Newton solver, which wants (x -> f(x), f(x)/f'(x)) as input to solve f(x) = 0
# """
# function _solve_Δω²(k,ωₜ,ε⁻¹::Array{SHM3,3},ds::MaxwellData;neigs=1,eigind=1,maxiter=10000,tol=1e-8)
#     ds.k = SVector(0.,0.,k)
#     ds.kpG .= kpG(SVector(0.,0.,k),ds.grid)
#     # res = IterativeSolvers.lobpcg(M̂(ε⁻¹,ds),false,ds.H⃗;P=P̂(ε⁻¹,ds),maxiter,tol)
#     res = IterativeSolvers.lobpcg(M̂!(ε⁻¹,ds),false,ds.H⃗;P=P̂!(ε⁻¹,ds),maxiter,tol)
#     H =  res.X #[:,eigind]                      # eigenmode wavefn. magnetic fields in transverse pol. basis
#     ω² =  (real(res.λ[eigind]))                # eigenmode temporal freq.,  neff = kz / ωₖ, kz = k[3]
#     Δω² = ω² - ωₜ^2
#     ω²ₖ =   2 * real( ( H[:,eigind]' * M̂ₖ(ε⁻¹,ds) * H[:,eigind] )[1])       # ωₖ/∂kz = group velocity = c / ng, c = 1 here
#     ds.H⃗ .= H
#     ds.ω² = ω²
#     ds.ω²ₖ = ω²ₖ
#     return Δω² , Δω² / ω²ₖ
# end

# function solve_k(ω::Float64,k₀::Float64,ε⁻¹::Array{SHM3,3},g::MaxwellGrid;neigs=1,eigind=1,maxiter=10000,tol=1e-8)
#     ds = MaxwellData(k₀,g)
#     kz = Roots.find_zero(k -> _solve_Δω²(k,ω,ε⁻¹,ds;neigs,eigind,maxiter,tol), k₀, Roots.Newton())
#     ds.ω = √ds.ω²
#     ds.ωₖ = ds.ω²ₖ / ( 2 * ds.ω )
#     return kz, ds
# end

# function solve_k(ω::Float64,ε⁻¹::Array{SHM3,3},ds::MaxwellData;neigs=1,eigind=1,maxiter=10000,tol=1e-8)
#     kz = Roots.find_zero(k -> _solve_Δω²(k,ω,ε⁻¹,ds;neigs,eigind,maxiter,tol), ds.k[3], Roots.Newton())
#     ds.ω = √ds.ω²
#     ds.ωₖ = ds.ω²ₖ / ( 2 * ds.ω )
#     return kz
# end

# function compare_fields(f_mpb,f,xlim,ylim)
#     hm_f_mpb_real = [ heatmap(x_mpb,y_mpb,[real(f_mpb[i,j][ix]) for i=1:nx_mpb,j=1:ny_mpb]',aspect_ratio=:equal,c=cgrad(:RdBu),xlim=xlim,ylim=ylim) for ix=1:3]
#     hm_f_mpb_imag = [ heatmap(x_mpb,y_mpb,[imag(f_mpb[i,j][ix]) for i=1:nx_mpb,j=1:ny_mpb]',aspect_ratio=:equal,c=cgrad(:viridis),xlim=xlim,ylim=ylim) for ix=1:3]

#     hm_f_real = [ heatmap(x_mpb,y_mpb,[real(f[ix,i,j]) for i=1:nx_mpb,j=1:ny_mpb]',aspect_ratio=:equal,c=cgrad(:RdBu),xlim=xlim,ylim=ylim) for ix=1:3]
#     hm_f_imag = [ heatmap(x_mpb,y_mpb,[imag(f[ix,i,j]) for i=1:nx_mpb,j=1:ny_mpb]',aspect_ratio=:equal,c=cgrad(:viridis),xlim=xlim,ylim=ylim) for ix=1:3]

#     hm_f_ratio_real = [ heatmap(x_mpb,y_mpb,[real(f[ix,i,j])/real(f_mpb[i,j][ix]) for i=1:nx_mpb,j=1:ny_mpb]',aspect_ratio=:equal,c=cgrad(:RdBu),xlim=xlim,ylim=ylim) for ix=1:3]
#     hm_f_ratio_imag = [ heatmap(x_mpb,y_mpb,[imag(f[ix,i,j])/imag(f_mpb[i,j][ix]) for i=1:nx_mpb,j=1:ny_mpb]',aspect_ratio=:equal,c=cgrad(:viridis),xlim=xlim,ylim=ylim) for ix=1:3]

#     l = @layout [   a   b   c
#                     d   e   f
#                     g   h   i
#                     k   l   m
#                     n   o   p
#                     q   r   s    ]
#     plot(hm_f_mpb_real...,
#         hm_f_mpb_imag...,
#         hm_f_real...,
#         hm_f_imag...,
#         hm_f_ratio_real...,
#         hm_f_ratio_imag...,
#         layout=l,
#         size = (1300,1300),
#     )
# end

# # Mutating Operators

# # function t2c!(Hin::AbstractArray{ComplexF64,4},ds::MaxwellData)::AbstractArray{ComplexF64,4}
# #     for i=1:ds.Nx,j=1:ds.Ny,k=1:ds.Nz
# #         @inbounds scale = ds.kpG[i,j,k].mag
# #         @inbounds ds.e[1,i,j,k] = ( Hin[1,i,j,k] * ds.kpG[i,j,k].m[1] + Hin[2,i,j,k] * ds.kpG[i,j,k].n[1] ) * scale
# #         @inbounds ds.e[2,i,j,k] = ( Hin[1,i,j,k] * ds.kpG[i,j,k].m[2] + Hin[2,i,j,k] * ds.kpG[i,j,k].n[2] ) * scale
# #         @inbounds ds.e[3,i,j,k] = ( Hin[1,i,j,k] * ds.kpG[i,j,k].m[3] + Hin[2,i,j,k] * ds.kpG[i,j,k].n[3] ) * scale
# #     end
# #     return ds.e
# # end

# # function c2t!(Hin::AbstractArray{ComplexF64,4},ds::MaxwellData)::AbstractArray{ComplexF64,4}
# #     for i=1:ds.Nx,j=1:ds.Ny,k=1:ds.Nz
# #         @inbounds ds.e[1,i,j,k] =  Hin[1,i,j,k] * ds.kpG[i,j,k].m[1] + Hin[2,i,j,k] * ds.kpG[i,j,k].m[2] + Hin[3,i,j,k] * ds.kpG[i,j,k].m[3]
# #         @inbounds ds.e[2,i,j,k] =  Hin[1,i,j,k] * ds.kpG[i,j,k].n[1] + Hin[2,i,j,k] * ds.kpG[i,j,k].n[2] + Hin[3,i,j,k] * ds.kpG[i,j,k].n[3]
# #         # @inbounds ds.e[3,i,j,k] =  0.0
# #     end
# #     return ds.e
# # end

# # function zcross_t2c!(Hin::AbstractArray{ComplexF64,4},ds::MaxwellData)::AbstractArray{ComplexF64,4}
# #     for i=1:ds.Nx,j=1:ds.Ny,k=1:ds.Nz
# #         @inbounds ds.e[1,i,j,k] = -Hin[1,i,j,k] * ds.kpG[i,j,k].m[2] - Hin[2,i,j,k] * ds.kpG[i,j,k].n[2]
# #         @inbounds ds.e[2,i,j,k] =  Hin[1,i,j,k] * ds.kpG[i,j,k].m[1] + Hin[2,i,j,k] * ds.kpG[i,j,k].n[1]
# #     end
# #     return ds.e
# # end

# ###############################################################################
# function kcross_t2c!(ds::MaxwellData)::AbstractArray{ComplexF64,4}
#     @inbounds for i=1:ds.Nx,j=1:ds.Ny,k=1:ds.Nz
#         scale = -ds.kpG[i,j,k].mag #-ds.kpG[i,j,k].mag
#         ds.d[1,i,j,k] = ( ds.H[1,i,j,k] * ds.kpG[i,j,k].n[1] - ds.H[2,i,j,k] * ds.kpG[i,j,k].m[1] ) * scale
#         ds.d[2,i,j,k] = ( ds.H[1,i,j,k] * ds.kpG[i,j,k].n[2] - ds.H[2,i,j,k] * ds.kpG[i,j,k].m[2] ) * scale
#         ds.d[3,i,j,k] = ( ds.H[1,i,j,k] * ds.kpG[i,j,k].n[3] - ds.H[2,i,j,k] * ds.kpG[i,j,k].m[3] ) * scale
#     end
#     return ds.d
# end

# function kcross_c2t!(ds::MaxwellData)::AbstractArray{ComplexF64,4}
#     @inbounds for i=1:ds.Nx,j=1:ds.Ny,k=1:ds.Nz
#         scale = ds.kpG[i,j,k].mag
#         at1 = ds.e[1,i,j,k] * ds.kpG[i,j,k].m[1] + ds.e[2,i,j,k] * ds.kpG[i,j,k].m[2] + ds.e[3,i,j,k] * ds.kpG[i,j,k].m[3]
#         at2 = ds.e[1,i,j,k] * ds.kpG[i,j,k].n[1] + ds.e[2,i,j,k] * ds.kpG[i,j,k].n[2] + ds.e[3,i,j,k] * ds.kpG[i,j,k].n[3]
#         ds.H[1,i,j,k] =  -at2 * scale
#         ds.H[2,i,j,k] =  at1 * scale
#     end
#     return ds.H
# end

# function kcrossinv_t2c!(ds::MaxwellData)::AbstractArray{ComplexF64,4}
#     @inbounds for i=1:ds.Nx,j=1:ds.Ny,k=1:ds.Nz
#         scale = 1 / ds.kpG[i,j,k].mag
#         ds.e[1,i,j,k] = ( ds.H[1,i,j,k] * ds.kpG[i,j,k].n[1] - ds.H[2,i,j,k] * ds.kpG[i,j,k].m[1] ) * scale
#         ds.e[2,i,j,k] = ( ds.H[1,i,j,k] * ds.kpG[i,j,k].n[2] - ds.H[2,i,j,k] * ds.kpG[i,j,k].m[2] ) * scale
#         ds.e[3,i,j,k] = ( ds.H[1,i,j,k] * ds.kpG[i,j,k].n[3] - ds.H[2,i,j,k] * ds.kpG[i,j,k].m[3] ) * scale
#     end
#     return ds.e
# end

# function kcrossinv_c2t!(ds::MaxwellData)::AbstractArray{ComplexF64,4}
#     @inbounds for i=1:ds.Nx,j=1:ds.Ny,k=1:ds.Nz
#         scale = -1 / ds.kpG[i,j,k].mag
#         at1 = ds.d[1,i,j,k] * ds.kpG[i,j,k].m[1] + ds.d[2,i,j,k] * ds.kpG[i,j,k].m[2] + ds.d[3,i,j,k] * ds.kpG[i,j,k].m[3]
#         at2 = ds.d[1,i,j,k] * ds.kpG[i,j,k].n[1] + ds.d[2,i,j,k] * ds.kpG[i,j,k].n[2] + ds.d[3,i,j,k] * ds.kpG[i,j,k].n[3]
#         ds.H[1,i,j,k] =  -at2 * scale
#         ds.H[2,i,j,k] =  at1 * scale
#     end
#     return ds.H
# end

# # function kcrossinv_t2c!(ds::MaxwellData)::AbstractArray{ComplexF64,4}
# #     @inbounds for i=1:ds.Nx,j=1:ds.Ny,k=1:ds.Nz
# #         scale = 1 / ds.kpG[i,j,k].mag
# #         ds.einv[1,i,j,k] = ( ds.Hinv[1,i,j,k] * ds.kpG[i,j,k].n[1] - ds.Hinv[2,i,j,k] * ds.kpG[i,j,k].m[1] ) * scale
# #         ds.einv[2,i,j,k] = ( ds.Hinv[1,i,j,k] * ds.kpG[i,j,k].n[2] - ds.Hinv[2,i,j,k] * ds.kpG[i,j,k].m[2] ) * scale
# #         ds.einv[3,i,j,k] = ( ds.Hinv[1,i,j,k] * ds.kpG[i,j,k].n[3] - ds.Hinv[2,i,j,k] * ds.kpG[i,j,k].m[3] ) * scale
# #     end
# #     return ds.einv
# # end

# # function kcrossinv_c2t!(ds::MaxwellData)::AbstractArray{ComplexF64,4}
# #     @inbounds for i=1:ds.Nx,j=1:ds.Ny,k=1:ds.Nz
# #         scale = -1 / ds.kpG[i,j,k].mag
# #         at1 = ds.dinv[1,i,j,k] * ds.kpG[i,j,k].m[1] + ds.dinv[2,i,j,k] * ds.kpG[i,j,k].m[2] + ds.dinv[3,i,j,k] * ds.kpG[i,j,k].m[3]
# #         at2 = ds.dinv[1,i,j,k] * ds.kpG[i,j,k].n[1] + ds.dinv[2,i,j,k] * ds.kpG[i,j,k].n[2] + ds.dinv[3,i,j,k] * ds.kpG[i,j,k].n[3]
# #         ds.Hinv[1,i,j,k] =  -at2 * scale
# #         ds.Hinv[2,i,j,k] =  at1 * scale
# #     end
# #     return ds.Hinv
# # end

# function ε⁻¹_dot!(ε⁻¹::Array{SHermitianCompact{3,Float64,6},3},ds::MaxwellData)::AbstractArray{ComplexF64,4}
#     @inbounds for i=1:ds.Nx,j=1:ds.Ny,k=1:ds.Nz
#         ds.e[1,i,j,k] =  ε⁻¹[i,j,k][1,1]*ds.d[1,i,j,k] + ε⁻¹[i,j,k][2,1]*ds.d[2,i,j,k] + ε⁻¹[i,j,k][3,1]*ds.d[3,i,j,k]
#         ds.e[2,i,j,k] =  ε⁻¹[i,j,k][1,2]*ds.d[1,i,j,k] + ε⁻¹[i,j,k][2,2]*ds.d[2,i,j,k] + ε⁻¹[i,j,k][3,2]*ds.d[3,i,j,k]
#         ds.e[3,i,j,k] =  ε⁻¹[i,j,k][1,3]*ds.d[1,i,j,k] + ε⁻¹[i,j,k][2,3]*ds.d[2,i,j,k] + ε⁻¹[i,j,k][3,3]*ds.d[3,i,j,k]
#         # ds.e[1,i,j,k] =  ε⁻¹[i,j,k][1,1]*Hin[1,i,j,k] + ε⁻¹[i,j,k][1,2]*Hin[2,i,j,k] + ε⁻¹[i,j,k][1,3]*Hin[3,i,j,k]
#         # ds.e[2,i,j,k] =  ε⁻¹[i,j,k][2,1]*Hin[1,i,j,k] + ε⁻¹[i,j,k][2,2]*Hin[2,i,j,k] + ε⁻¹[i,j,k][2,3]*Hin[3,i,j,k]
#         # ds.e[3,i,j,k] =  ε⁻¹[i,j,k][3,1]*Hin[1,i,j,k] + ε⁻¹[i,j,k][3,2]*Hin[2,i,j,k] + ε⁻¹[i,j,k][3,3]*Hin[3,i,j,k]
#     end
#     return ds.e
# end

# # function ε_dot_approx!(ε⁻¹::Array{SHermitianCompact{3,Float64,6},3},ds::MaxwellData)::AbstractArray{ComplexF64,4}
# #     @inbounds for i=1:ds.Nx,j=1:ds.Ny,k=1:ds.Nz
# #         ε_ave = 3 / tr(ε⁻¹[i,j,k])
# #         ds.dinv[1,i,j,k] =  ε_ave * ds.einv[1,i,j,k]
# #         ds.dinv[2,i,j,k] =  ε_ave * ds.einv[2,i,j,k]
# #         ds.dinv[3,i,j,k] =  ε_ave * ds.einv[3,i,j,k]
# #     end
# #     return ds.dinv
# # end

# function ε_dot_approx!(ε⁻¹::Array{SHermitianCompact{3,Float64,6},3},ds::MaxwellData)::AbstractArray{ComplexF64,4}
#     @inbounds for i=1:ds.Nx,j=1:ds.Ny,k=1:ds.Nz
#         ε_ave = 3 / tr(ε⁻¹[i,j,k])
#         ds.d[1,i,j,k] =  ε_ave * ds.e[1,i,j,k]
#         ds.d[2,i,j,k] =  ε_ave * ds.e[2,i,j,k]
#         ds.d[3,i,j,k] =  ε_ave * ds.e[3,i,j,k]
#     end
#     return ds.d
# end

# function M!(ε⁻¹::Array{SHermitianCompact{3,Float64,6},3},ds::MaxwellData)::Array{ComplexF64,4}
#     kcross_t2c!(ds);
#     # ds.𝓕! * ds.d;
#     mul!(ds.d,ds.𝓕!,ds.d);
#     ε⁻¹_dot!(ε⁻¹,ds);
#     # ds.𝓕⁻¹! * ds.e;
#     mul!(ds.e,ds.𝓕⁻¹!,ds.e)
#     kcross_c2t!(ds)
# end

# function M!(Hout::AbstractArray{ComplexF64,1},Hin::AbstractArray{ComplexF64,1},ε⁻¹::Array{SHermitianCompact{3,Float64,6},3},ds::MaxwellData)::Array{ComplexF64,1}
#     # copyto!(ds.H,reshape(Hin,(2,ds.Nx,ds.Ny,ds.Nz)))
#     @inbounds ds.H .= reshape(Hin,(2,ds.Nx,ds.Ny,ds.Nz))
#     M!(ε⁻¹,ds);
#     # copyto!(Hout,vec(ds.H))
#     @inbounds Hout .= vec(ds.H)
# end

# function M̂!(ε⁻¹::Array{SHermitianCompact{3,Float64,6},3},ds::MaxwellData)
#     function f!(y::AbstractArray{ComplexF64,1},x::AbstractArray{ComplexF64,1})::AbstractArray{ComplexF64,1}
#         M!(y,x,ε⁻¹,ds)
#     end
#     return LinearMap{ComplexF64}(f!,(2*ds.Nx*ds.Ny*ds.Nz),ishermitian=true,ismutating=true)
# end

# # function P!(ε⁻¹::Array{SHermitianCompact{3,Float64,6},3},ds::MaxwellData)::Array{ComplexF64,4}
# #     kcrossinv_t2c!(ds);
# #     # ds.𝓕⁻¹! * ds.e;
# #     mul!(ds.einv,ds.𝓕⁻¹!,ds.einv)
# #     ε_dot_approx!(ε⁻¹,ds);
# #     # ds.𝓕! * ds.d;
# #     mul!(ds.dinv,ds.𝓕!,ds.dinv);
# #     kcrossinv_c2t!(ds)
# # end

# # function P!(Hout::AbstractArray{ComplexF64,1},Hin::AbstractArray{ComplexF64,1},ε⁻¹::Array{SHermitianCompact{3,Float64,6},3},ds::MaxwellData)::Array{ComplexF64,1}
# #     # copyto!(ds.H,reshape(Hin,(2,ds.Nx,ds.Ny,ds.Nz)))
# #     @inbounds ds.Hinv .= reshape(Hin,(2,ds.Nx,ds.Ny,ds.Nz))
# #     P!(ε⁻¹,ds);
# #     # copyto!(Hout,vec(ds.H))
# #     @inbounds Hout .= vec(ds.Hinv)
# # end

# function P!(ε⁻¹::Array{SHermitianCompact{3,Float64,6},3},ds::MaxwellData)::Array{ComplexF64,4}
#     kcrossinv_t2c!(ds);
#     # ds.𝓕⁻¹! * ds.e;
#     mul!(ds.e,ds.𝓕⁻¹!,ds.e)
#     ε_dot_approx!(ε⁻¹,ds);
#     # ds.𝓕! * ds.d;
#     mul!(ds.d,ds.𝓕!,ds.d);
#     kcrossinv_c2t!(ds)
# end

# function P!(Hout::AbstractArray{ComplexF64,1},Hin::AbstractArray{ComplexF64,1},ε⁻¹::Array{SHermitianCompact{3,Float64,6},3},ds::MaxwellData)::Array{ComplexF64,1}
#     # copyto!(ds.H,reshape(Hin,(2,ds.Nx,ds.Ny,ds.Nz)))
#     @inbounds ds.H .= reshape(Hin,(2,ds.Nx,ds.Ny,ds.Nz))
#     P!(ε⁻¹,ds);
#     # copyto!(Hout,vec(ds.H))
#     @inbounds Hout .= vec(ds.H)
# end

# function P̂!(ε⁻¹::Array{SHermitianCompact{3,Float64,6},3},ds::MaxwellData)
#     function fp!(y::AbstractArray{ComplexF64,1},x::AbstractArray{ComplexF64,1})::AbstractArray{ComplexF64,1}
#         P!(y,x,ε⁻¹,ds)
#     end
#     return LinearMap{ComplexF64}(fp!,(2*ds.Nx*ds.Ny*ds.Nz),ishermitian=true,ismutating=true)
# end

# ######################
# ######################3
# ########################
# ##

# eigind = 1;
# k₀ = 1.4 # k_mpb[1] # 1.47625...
# ds = MaxwellData(k₀,g)
# H, ω, ωₖ = solve_ω(SVector(0.,0.,k₀),ε⁻¹_mpb,g;neigs=1,eigind=1,maxiter=3000,tol=1e-7)

# ##

# # HA = reshape(ds.H,(2,ds.Nx,ds.Ny,ds.Nz))
# HA = M!(ε⁻¹_mpb,ds)
# # M!(ε⁻¹_mpb,ds)

# ##
# k₀ = 1.3 # k_mpb[1] # 1.47625...
# ds = MaxwellData(k₀,g)
# ω_mpb = k_mpb[1] / neff_mpb[1]
# @btime solve_k($ω_mpb,$ε⁻¹_mpb,$ds)
# # kz = solve_k(ω_mpb,ε⁻¹_mpb,ds)
# ##
# ω = ds.ω
# ωₖ = ds.ωₖ
# H = copy(reshape(ds.H⃗[:,1],(2,ds.Nx,ds.Ny,ds.Nz)))
# # HAout = similar(HAin)

# # HA = M(HAin,ε⁻¹_mpb,ds)
# e = ε⁻¹_dot(ds.𝓕 * kcross_t2c(HA,ds),ε⁻¹_mpb,ds)
# neff = kz / ω
# neff_err = neff - neff_mpb[1]
# neff_err_rel = neff_err / neff_mpb[1]
# println("k = $kz")
# println("ω = $ωₖ")
# println("ω_mpb = $ω_mpb")
# println("neff = $neff")
# println("neff_mpb = $(neff_mpb[1])")
# println("neff_err = $neff_err")
# println("neff_err_rel = $neff_err_rel")

# ng = 1 / ωₖ
# ng_err = ng - ng_mpb
# ng_err_rel = ng_err / ng_mpb
# println("ng = $ng")
# println("ng_mpb = $ng_mpb")
# println("ng_err = $ng_err")
# println("ng_err_rel = $ng_err_rel")

# compare_fields(e_mpb, e,(-3,3),(-2,2))

# # HAc = copy(HA)
# # MHAc = M(HA,ε⁻¹_mpb,ds)
# # # vec(MHAc ./ HAc)
# # # mean( vec(MHAc ./ HAc))
# # ##
# # HAout = similar(HAc)
# # MHAc2 = M!(HAout,HAc,ε⁻¹_mpb,ds)
# # ##
# # vec(MHAc2 ./ HAc)
# # ##
# # vec(HAout ./ HAc)
# # var( vec(MHAc ./ HAc))

# # abs2.( HAout - MHAc )

# ##

# function testfn2(kz,ε⁻¹_mpb,g)
#     ds = MaxwellData(kz,g)
#     res = IterativeSolvers.lobpcg(M̂!(ε⁻¹_mpb,ds),false,1;P=P̂!(ε⁻¹_mpb,ds),maxiter=1000,tol=1e-7)
# end

# ##

# @btime testfn2(1.4,$ε⁻¹_mpb,$g)

# ##
# k₀ = k_mpb[1] # 1.4 # k_mpb[1] # 1.47625...
# ds = MaxwellData(k₀,g)
# res = IterativeSolvers.lobpcg(M̂!(ε⁻¹_mpb,ds),false,1;P=P̂!(ε⁻¹_mpb,ds),maxiter=1000,tol=1e-7)
# # res = IterativeSolvers.lobpcg(M̂!(ε⁻¹_mpb,ds),false,1;maxiter=3000,tol=1e-7)

# H =  res.X #[:,eigind]                       # eigenmode wavefn. magnetic fields in transverse pol. basis
# ω =  √(abs(real(res.λ[eigind])))                     # eigenmode temporal freq.,  neff = kz / ω, kz = k[3]
# #ωₖ =   real( ( H[:,eigind]' * M̂ₖ(ε⁻¹,ds) * H[:,eigind] )[1]) / ω       # ωₖ/∂kz = group velocity = c / ng, c = 1 here
# ω_mpb = k_mpb[1] / neff_mpb[1]
# kz = k₀ #solve_k(ω_mpb,ε⁻¹_mpb,ds)
# HA = reshape(H,(2,ds.Nx,ds.Ny,ds.Nz))
# e = ε⁻¹_dot(ds.𝓕 * kcross_t2c(HA,ds),ε⁻¹_mpb,ds)
# neff = kz / ω
# neff_err = neff - neff_mpb[1]
# neff_err_rel = neff_err / neff_mpb[1]
# println("k = $kz")
# println("ω = $ωₖ")
# println("ω_mpb = $ω_mpb")
# println("neff = $neff")
# println("neff_mpb = $(neff_mpb[1])")
# println("neff_err = $neff_err")
# println("neff_err_rel = $neff_err_rel")

# ng = 1 / ωₖ
# ng_err = ng - ng_mpb
# ng_err_rel = ng_err / ng_mpb
# println("ng = $ng")
# println("ng_mpb = $ng_mpb")
# println("ng_err = $ng_err")
# println("ng_err_rel = $ng_err_rel")

# compare_fields(e_mpb, e,(-3,3),(-2,2))

# ##

# # @btime kcross_t2c!($Hout1,$HA,$ds);
# # @btime kcross_c2t!($Hout1,$HA,$ds);
# # @btime epsinv_dot!($Hout1,$HA,$ε⁻¹,$ds);
# # @btime $𝓕 * $Hout1;
# # @btime M!($Hout,$Hin,$ε⁻¹,$𝓕,$𝓕⁻¹,$ds,$Hw1,$Hw2); # 1.335 ms (0 allocations: 0 bytes)
# # @btime M!($Hout_v,$Hin_v,$ε⁻¹,$𝓕,$𝓕⁻¹,$ds,$Hw1,$Hw2); # 1.211 ms (4 allocations: 256 bytes)
# N = 3* ds.Nx * ds.Ny * ds.Nz
# f = (Hout_v,Hin_v) -> M!(Hout_v,Hin_v,ε⁻¹_mpb,𝓕,𝓕⁻¹,ds,Hw1,Hw2)
# Mop = LinearMap{ComplexF64}(f,N,ishermitian=true,ismutating=true)
# # @btime $Mop * $Hin_v  # 1.342 ms (6 allocations: 768.33 KiB
# res = IterativeSolvers.lobpcg(Mop,false,1;maxiter=3000,tol=1e-7)
# @show λ_sol = res.λ
# dsol_v = res.X[:,1]
# dsol = reshape(dsol_v,(3,size(ds.kpG)...)) # reinterpret(SVector{3,ComplexF64}, ... )[1,:,:,1]

# d =   copy(dsol)
# # H = permutedims(fftshift(reshape(vec(reinterpret(SVector{2,ComplexF64},permutedims(ms.get_eigenvectors(1,1),[2,1,3]))),(ny_mpb,nx_mpb))),[2,1])

# # function d_from_H2(H,ds::MaxwellData)
# #     kpG2 = [KVec(k-SVector(ggx, 0., 0.)-SVector(0., ggy, 0.)) for ggx=( twopi .* fftfreq(Nx,Nx/Δx)), ggy=(twopi  .* fftfreq(Ny,Ny/Δy)) ]
# #     d_recip = [ kcross_t2c(H[i,j],kpG2[i,j]) for i=1:ds.Nx, j=1:ds.Ny]
# #     temp =  (-1/(2π)) .* fft( reinterpret( ComplexF64, reshape( d_recip , (1,ds.Nx,ds.Ny) )), (2,3))
# #     return reshape(reinterpret(SVector{3,ComplexF64},temp),(ds.Nx,ds.Ny))
# # end

# # d = d_from_H2( H, ds )

# xlim =(-3,3)
# ylim =(-2,2)

# # xlim =(-0.5,0.5)
# # ylim =(-0.5,0.5)

# hm_d_mpb_real = [ heatmap(x_mpb,y_mpb,[real(d_mpb[i,j][ix]) for i=1:nx_mpb,j=1:ny_mpb]',aspect_ratio=:equal,c=cgrad(:RdBu),xlim=xlim,ylim=ylim) for ix=1:3]
# hm_d_mpb_imag = [ heatmap(x_mpb,y_mpb,[imag(d_mpb[i,j][ix]) for i=1:nx_mpb,j=1:ny_mpb]',aspect_ratio=:equal,c=cgrad(:viridis),xlim=xlim,ylim=ylim) for ix=1:3]

# hm_d_real = [ heatmap(x_mpb,y_mpb,[real(d[ix,i,j]) for i=1:nx_mpb,j=1:ny_mpb]',aspect_ratio=:equal,c=cgrad(:RdBu),xlim=xlim,ylim=ylim) for ix=1:3]
# hm_d_imag = [ heatmap(x_mpb,y_mpb,[imag(d[ix,i,j]) for i=1:nx_mpb,j=1:ny_mpb]',aspect_ratio=:equal,c=cgrad(:viridis),xlim=xlim,ylim=ylim) for ix=1:3]

# hm_d_ratio_real = [ heatmap(x_mpb,y_mpb,[real(d[ix,i,j])/real(d_mpb[i,j][ix]) for i=1:nx_mpb,j=1:ny_mpb]',aspect_ratio=:equal,c=cgrad(:RdBu),xlim=xlim,ylim=ylim) for ix=1:3]
# hm_d_ratio_imag = [ heatmap(x_mpb,y_mpb,[imag(d[ix,i,j])/imag(d_mpb[i,j][ix]) for i=1:nx_mpb,j=1:ny_mpb]',aspect_ratio=:equal,c=cgrad(:viridis),xlim=xlim,ylim=ylim) for ix=1:3]

# l = @layout [   a   b   c
#                 d   e   f
#                 g   h   i
#                 k   l   m
#                 n   o   p
#                 q   r   s    ]

# plot(
#     hm_d_mpb_real...,
#     hm_d_mpb_imag...,
#     hm_d_real...,
#     hm_d_imag...,
#     hm_d_ratio_real...,
#     hm_d_ratio_imag...,
#     layout=l,
#     size = (1300,1300),
# )

# ##

# function h_from_H!(h::AbstractArray{ComplexF64,4},H::AbstractArray{ComplexF64,4},ds::MaxwellData)::AbstractArray{ComplexF64,4}
#     t2c!(h,H,ds);
#     𝓕 * h
# end
# h_temp = copy(H_mpb)
# # h =   h_from_H!(h_temp, H_mpb, ds) .* (1/(2π))
# h =   h_from_H!(h_temp, Hsol, ds) .* (1/(2π))

# xlim =(-3,3)
# ylim =(-2,2)

# # xlim =(-0.5,0.5)
# # ylim =(-0.5,0.5)

# hm_h_mpb_real = [ heatmap(x_mpb,y_mpb,[real(h_mpb[i,j][ix]) for i=1:nx_mpb,j=1:ny_mpb]',aspect_ratio=:equal,c=cgrad(:RdBu),xlim=xlim,ylim=ylim) for ix=1:3]
# hm_h_mpb_imag = [ heatmap(x_mpb,y_mpb,[imag(h_mpb[i,j][ix]) for i=1:nx_mpb,j=1:ny_mpb]',aspect_ratio=:equal,c=cgrad(:viridis),xlim=xlim,ylim=ylim) for ix=1:3]

# hm_h_real = [ heatmap(x_mpb,y_mpb,[real(h[ix,i,j]) for i=1:nx_mpb,j=1:ny_mpb]',aspect_ratio=:equal,c=cgrad(:RdBu),xlim=xlim,ylim=ylim) for ix=1:3]
# hm_h_imag = [ heatmap(x_mpb,y_mpb,[imag(h[ix,i,j]) for i=1:nx_mpb,j=1:ny_mpb]',aspect_ratio=:equal,c=cgrad(:viridis),xlim=xlim,ylim=ylim) for ix=1:3]

# hm_h_ratio_real = [ heatmap(x_mpb,y_mpb,[real(h[ix,i,j])/real(h_mpb[i,j][ix]) for i=1:nx_mpb,j=1:ny_mpb]',aspect_ratio=:equal,c=cgrad(:RdBu),xlim=xlim,ylim=ylim) for ix=1:3]
# hm_h_ratio_imag = [ heatmap(x_mpb,y_mpb,[imag(h[ix,i,j])/imag(h_mpb[i,j][ix]) for i=1:nx_mpb,j=1:ny_mpb]',aspect_ratio=:equal,c=cgrad(:viridis),xlim=xlim,ylim=ylim) for ix=1:3]

# l = @layout [   a   b   c
#                 d   e   f
#                 g   h   i
#                 k   l   m
#                 n   o   p
#                 q   r   s    ]

# plot(
#     hm_h_mpb_real...,
#     hm_h_mpb_imag...,
#     hm_h_real...,
#     hm_h_imag...,
#     hm_h_ratio_real...,
#     hm_h_ratio_imag...,
#     layout=l,
#     size = (1300,1300),
# )

# ##

# # function d_from_H!(d::AbstractArray{ComplexF64,4},H::AbstractArray{ComplexF64,4},ds::MaxwellData)::AbstractArray{ComplexF64,4}
# #     kcross_t2c!(d,H,ds);
# #     𝓕 * d
# # end
# # d_temp = copy(H_mpb)
# # d =   d_from_H!(d_temp, H_mpb, ds) .* (-1/(π))

# xlim =(-3,3)
# ylim =(-2,2)

# hm_d_mpb_real = [ heatmap(x_mpb,y_mpb,[real(d_mpb[i,j][ix]) for i=1:nx_mpb,j=1:ny_mpb]',aspect_ratio=:equal,c=cgrad(:RdBu),xlim=xlim,ylim=ylim) for ix=1:3]
# hm_d_mpb_imag = [ heatmap(x_mpb,y_mpb,[imag(d_mpb[i,j][ix]) for i=1:nx_mpb,j=1:ny_mpb]',aspect_ratio=:equal,c=cgrad(:viridis),xlim=xlim,ylim=ylim) for ix=1:3]

# hm_d_real = [ heatmap(x_mpb,y_mpb,[real(d[ix,i,j]) for i=1:nx_mpb,j=1:ny_mpb]',aspect_ratio=:equal,c=cgrad(:RdBu),xlim=xlim,ylim=ylim) for ix=1:3]
# hm_d_imag = [ heatmap(x_mpb,y_mpb,[imag(d[ix,i,j]) for i=1:nx_mpb,j=1:ny_mpb]',aspect_ratio=:equal,c=cgrad(:viridis),xlim=xlim,ylim=ylim) for ix=1:3]

# hm_d_ratio_real = [ heatmap(x_mpb,y_mpb,[real(d[ix,i,j])/real(d_mpb[i,j][ix]) for i=1:nx_mpb,j=1:ny_mpb]',aspect_ratio=:equal,c=cgrad(:RdBu),xlim=xlim,ylim=ylim) for ix=1:3]
# hm_d_ratio_imag = [ heatmap(x_mpb,y_mpb,[imag(d[ix,i,j])/imag(d_mpb[i,j][ix]) for i=1:nx_mpb,j=1:ny_mpb]',aspect_ratio=:equal,c=cgrad(:viridis),xlim=xlim,ylim=ylim) for ix=1:3]

# l = @layout [   a   b   c
#                 d   e   f
#                 g   h   i
#                 k   l   m
#                 n   o   p
#                 q   r   s    ]

# plot(
#     hm_d_mpb_real...,
#     hm_d_mpb_imag...,
#     hm_d_real...,
#     hm_d_imag...,
#     hm_d_ratio_real...,
#     hm_d_ratio_imag...,
#     layout=l,
#     size = (1300,1300),
# )

# ##

# function e_from_H!(e::AbstractArray{ComplexF64,4},d::AbstractArray{ComplexF64,4},H::AbstractArray{ComplexF64,4},ε⁻¹::Array{SHermitianCompact{3,Float64,6},3},ds::MaxwellData)::AbstractArray{ComplexF64,4}
#     kcross_t2c!(d,H,ds);
#     𝓕 * d;
#     epsinv_dot!(e,(-1/(π)) .* d,ε⁻¹,ds);
# end

# d_temp = copy(H_mpb)
# e_temp = copy(H_mpb)
# e = e_from_H!(e_temp,d_temp,H_mpb,ε⁻¹_mpb,ds)

# xlim =(-3,3)
# ylim =(-2,2)

# # xlim =(-0.5,0.5)
# # ylim =(-0.5,0.5)

# hm_e_mpb_real = [ heatmap(x_mpb,y_mpb,[real(e_mpb[i,j][ix]) for i=1:nx_mpb,j=1:ny_mpb]',aspect_ratio=:equal,c=cgrad(:RdBu),xlim=xlim,ylim=ylim) for ix=1:3]
# hm_e_mpb_imag = [ heatmap(x_mpb,y_mpb,[imag(e_mpb[i,j][ix]) for i=1:nx_mpb,j=1:ny_mpb]',aspect_ratio=:equal,c=cgrad(:viridis),xlim=xlim,ylim=ylim) for ix=1:3]

# hm_e_real = [ heatmap(x_mpb,y_mpb,[real(e[ix,i,j]) for i=1:nx_mpb,j=1:ny_mpb]',aspect_ratio=:equal,c=cgrad(:RdBu),xlim=xlim,ylim=ylim) for ix=1:3]
# hm_e_imag = [ heatmap(x_mpb,y_mpb,[imag(e[ix,i,j]) for i=1:nx_mpb,j=1:ny_mpb]',aspect_ratio=:equal,c=cgrad(:viridis),xlim=xlim,ylim=ylim) for ix=1:3]

# hm_e_ratio_real = [ heatmap(x_mpb,y_mpb,[real(e[ix,i,j])/real(e_mpb[i,j][ix]) for i=1:nx_mpb,j=1:ny_mpb]',aspect_ratio=:equal,c=cgrad(:RdBu),xlim=xlim,ylim=ylim) for ix=1:3]
# hm_e_ratio_imag = [ heatmap(x_mpb,y_mpb,[imag(e[ix,i,j])/imag(e_mpb[i,j][ix]) for i=1:nx_mpb,j=1:ny_mpb]',aspect_ratio=:equal,c=cgrad(:viridis),xlim=xlim,ylim=ylim) for ix=1:3]

# l = @layout [   a   b   c
#                 d   e   f
#                 g   h   i
#                 k   l   m
#                 n   o   p
#                 q   r   s    ]
# plot(
#     hm_e_mpb_real...,
#     hm_e_mpb_imag...,
#     hm_e_real...,
#     hm_e_imag...,
#     hm_e_ratio_real...,
#     hm_e_ratio_imag...,
#     layout=l,
#     size = (1300,1300),
# )

# ##

# function d_from_d!(dout::AbstractArray{ComplexF64,4},e::AbstractArray{ComplexF64,4},H::AbstractArray{ComplexF64,4},din::AbstractArray{ComplexF64,4},ε⁻¹::Array{SHermitianCompact{3,Float64,6},3},ds::MaxwellData)::AbstractArray{ComplexF64,4}
#     epsinv_dot!(e,-din,ε⁻¹,ds);
#     𝓕⁻¹ * e;
#     kcross_c2t!(H,e,ds);
#     kcross_t2c!(dout,H,ds);
#     𝓕 * dout;
#     return dout
# end

# din = copy(d)
# H_temp = copy(H_mpb)
# e_temp = copy(H_mpb)
# dout_temp = copy(H_mpb)
# dout =  d_from_d!(dout_temp,e_temp,H_temp,din,ε⁻¹_mpb,ds)

# xlim =(-3,3)
# ylim =(-2,2)

# # xlim =(-0.5,0.5)
# # ylim =(-0.5,0.5)
# hm_d_mpb_real = [ heatmap(x_mpb,y_mpb,[real(d_mpb[i,j][ix]) for i=1:nx_mpb,j=1:ny_mpb]',aspect_ratio=:equal,c=cgrad(:RdBu),xlim=xlim,ylim=ylim) for ix=1:3]
# hm_d_mpb_imag = [ heatmap(x_mpb,y_mpb,[imag(d_mpb[i,j][ix]) for i=1:nx_mpb,j=1:ny_mpb]',aspect_ratio=:equal,c=cgrad(:viridis),xlim=xlim,ylim=ylim) for ix=1:3]

# hm_d_real = [ heatmap(x_mpb,y_mpb,[real(dout[ix,i,j]) for i=1:nx_mpb,j=1:ny_mpb]',aspect_ratio=:equal,c=cgrad(:RdBu),xlim=xlim,ylim=ylim) for ix=1:3]
# hm_d_imag = [ heatmap(x_mpb,y_mpb,[imag(dout[ix,i,j]) for i=1:nx_mpb,j=1:ny_mpb]',aspect_ratio=:equal,c=cgrad(:viridis),xlim=xlim,ylim=ylim) for ix=1:3]

# hm_d_ratio_real = [ heatmap(x_mpb,y_mpb,[real(dout[ix,i,j])/real(d_mpb[i,j][ix]) for i=1:nx_mpb,j=1:ny_mpb]',aspect_ratio=:equal,c=cgrad(:RdBu),xlim=xlim,ylim=ylim) for ix=1:3]
# hm_d_ratio_imag = [ heatmap(x_mpb,y_mpb,[imag(dout[ix,i,j])/imag(d_mpb[i,j][ix]) for i=1:nx_mpb,j=1:ny_mpb]',aspect_ratio=:equal,c=cgrad(:viridis),xlim=xlim,ylim=ylim) for ix=1:3]

# l = @layout [   a   b   c
#                 d   e   f
#                 g   h   i
#                 k   l   m
#                 n   o   p
#                 q   r   s    ]

# plot(
#     hm_d_mpb_real...,
#     hm_d_mpb_imag...,
#     hm_d_real...,
#     hm_d_imag...,
#     hm_d_ratio_real...,
#     hm_d_ratio_imag...,
#     layout=l,
#     size = (1300,1300),
# )

# ##
# # @btime kcross_t2c!($Hout1,$HA,$ds);
# # @btime kcross_c2t!($Hout1,$HA,$ds);
# # @btime epsinv_dot!($Hout1,$HA,$ε⁻¹,$ds);
# # @btime $𝓕 * $Hout1;
# # @btime M!($Hout,$Hin,$ε⁻¹,$𝓕,$𝓕⁻¹,$ds,$Hw1,$Hw2); # 1.335 ms (0 allocations: 0 bytes)
# # @btime M!($Hout_v,$Hin_v,$ε⁻¹,$𝓕,$𝓕⁻¹,$ds,$Hw1,$Hw2); # 1.211 ms (4 allocations: 256 bytes)
# N = 3* ds.Nx * ds.Ny * ds.Nz
# f = (Hout_v,Hin_v) -> M!(Hout_v,Hin_v,ε⁻¹,𝓕,𝓕⁻¹,ds,Hw1,Hw2)
# Mop = LinearMap{ComplexF64}(f,N,ismutating=true)
# # @btime $Mop * $Hin_v  # 1.342 ms (6 allocations: 768.33 KiB
# #res = IterativeSolvers.lobpcg(Mop,false,1;maxiter=3000,tol=1e-7)

# ##

# Hsol_v = res.X[:,1]
# Hsol = reshape(Hsol_v,(3,size(ds.kpG)...)) # reinterpret(SVector{3,ComplexF64}, ... )[1,:,:,1]
# dsol = similar(Hsol)
# kcross_t2c!(dsol,Hsol,ds);
# 𝓕 * dsol;

# heatmap([real(dsol[1,i,j]) for i=1:ds.Nx,j=1:ds.Ny])

# ##

# function testfn() #ε⁻¹,ds)
#     N = 3* ds.Nx * ds.Ny * ds.Nz
#     H1 = similar(HA);
#     H2 = similar(HA);
#     xout = similar(HA);
#     𝓕 = plan_fft!(HA,(2:4))
#     𝓕⁻¹ = inv(𝓕)
#     f = (xout,xin) -> M!(H1,H2,xout,xin,ε⁻¹,𝓕,𝓕⁻¹,ds)
#     Mop = LinearMap{ComplexF64}(f,f,N,N,ishermitian=true,ismutating=true)
#     res = IterativeSolvers.lobpcg(Mop,false,1;maxiter=3000,tol=1e-7)
#     # vals,vecs,info = eigsolve(x -> M!(H1,H2,H3,x,ε⁻¹,𝓕,𝓕⁻¹,ds),HA)
# end
# ##
# # function kcross_t2c!(Hout::Array{ComplexF64,4},Hin::Array{ComplexF64,4},ds::MaxwellData)::Array{ComplexF64,4}
# #     Hout[:,i,j,k] = [ Vector{ComplexF64}( ( Hin[1,i,j,k] * ds.kpG[i,j,k].n - Hin[2,i,j,k] * ds.kpG[i,j,k].m ) * ds.kpG[i,j,k].mag ) for i=1:ds.Nx,j=1:ds.Ny,k=1:ds.Nz ]
# # end

# # in place operations with staticarrays
# function kcross_t2c(v::SVector{3,ComplexF64},k::KVec)::SVector{3,ComplexF64}
#     return ( v[1] * k.n - v[2] * k.m ) * k.mag
# end

# function kcross_c2t(a::SVector{3,ComplexF64},k::KVec)::SVector{3,ComplexF64}
#     at1 = a ⋅ k.m
#     at2 = a ⋅ k.n
#     v0 = -at2 * k.mag
#     v1 = at1 * k.mag
#     return SVector(v0,v1,0.0)
# end

# function kcrossinv_t2c(v::SVector{2,ComplexF64},k::KVec)::SVector{3,ComplexF64}
#     return ( v[1] * k.n - v[2] * k.m ) * ( -1 / k.mag )
# end

# function kcrossinv_c2t(a::SVector{3,ComplexF64},k::KVec)::SVector{2,ComplexF64}
#     at1 = a ⋅ k.m
#     at2 = a ⋅ k.n
#     v0 = -at2 * (-1 / k.mag )
#     v1 = at1 * ( -1 / k.mag )
#     return SVector(v0,v1)
# end

# function myfft!(H::Array{SVector{3,ComplexF64}})::Nothing
#     P * reinterpret(ComplexF64,H);
#     return nothing
# end

# function myifft!(H::Array{SVector{3,ComplexF64}})::Nothing
#     iP * reinterpret(ComplexF64,H);
#     return nothing
# end

# function flatten(H::Array{SVector{3,ComplexF64}})
#     return vec(reinterpret(ComplexF64,H))
# end

# function unflatten(Hvec::Array{ComplexF64},shape::NTuple{4,Int64})::Array{SVector{3,ComplexF64}}
#     return reinterpret(SVector{3,ComplexF64},reshape(Hvec,(3*shape[1],shape[2:end]...)))
# end

# function M!(H::Array{SVector{3,ComplexF64}},eps::Array{SHM3},ds::MaxwellData)::Array{SVector{3,ComplexF64}}
#     H .= kcross_t2c.(H,ds.kpG);
#     ds.𝓕 * reinterpret(ComplexF64,H); # reinterpret(ComplexF64,reshape(H,(1,size(k)...)));
#     H .= eps .* -H;
#     ds.𝓕⁻¹ * reinterpret(ComplexF64,H); # reinterpret(ComplexF64,reshape(H,(1,size(k)...)));
#     H .= kcross_c2t.(conj.(H),ds.kpG)
# end

# function P!(H::Array{SVector{3,ComplexF64}},eps::Array{SHM3},ds::MaxwellData)::Array{SVector{3,ComplexF64}}
#     H .= kcrossinv_t2c.(H,ds.kpG);
#     ds.𝓕⁻¹ * reinterpret(ComplexF64,H); #reinterpret(ComplexF64,reshape(H,(1,size(k)...)));
#     H .= eps .* -H;
#     ds.𝓕 * reinterpret(ComplexF64,H);  # reinterpret(ComplexF64,reshape(H,(1,size(k)...)));
#     H .= kcrossinv_c2t.(conj.(H),ds.kpG)
# end

# function M!(Hv::Vector{SVector{3,ComplexF64}},eps::Array{SHM3},ds::MaxwellData)::Vector{SVector{3,ComplexF64}}
#     H = reshape(Hv,size(k))
#     H .= kcross_t2c.(H,ds.kpG);
#     myfft!(H);
#     H .= eps .* H;
#     myifft!(H);
#     H .= kcross_c2t.(H,ds.kpG)
#     return Hv
# end

# function M!(Hv::Vector{ComplexF64},eps::Array{SHM3},ds::MaxwellData)::Vector{ComplexF64}
#     HSv = copy(reinterpret(SVector{3,ComplexF64}, Hv))
#     M!(HSv,eps,ds.kpG);
#     Hv .= copy(reinterpret(ComplexF64, HSv))
# end

# function M!(Hv::Vector{ComplexF64},eps::Array{SHM3},ds::MaxwellData,Hw::Array{SVector{3,ComplexF64}})::Vector{ComplexF64}
#     copyto!(Hw, reinterpret(SVector{3,ComplexF64},reshape(Hv,(3,size(ds.kpG)...)))[1,:,:,:] )
#     M!(Hw,eps,ds.kpG);
#     copyto!(Hv, vec( reinterpret(ComplexF64,Hw) ) )
#     return Hv
# end

# function P!(Hv::Vector{ComplexF64},eps::Array{SHM3},ds::MaxwellData,Hw::Array{SVector{3,ComplexF64}})::Vector{ComplexF64}
#     copyto!(Hw, reinterpret(SVector{3,ComplexF64},reshape(Hv,(3,size(ds.kpG)...)))[1,:,:,:] )
#     P!(Hw,eps,ds.kpG);
#     copyto!(Hv, vec( reinterpret(ComplexF64,Hw) ) )
#     return Hv
# end

# Hv = copy(flatten(HSA));
# H₀ = copy(Hv);
# ##
# function M̂ₖ!(ε⁻¹::Array{SHM3},kpG::Array{KVec},Hw::Array{SVector{3,ComplexF64}})::LinearMaps.FunctionMap{ComplexF64}
#     N = 3 * length(kpG)
#     f = (Hout, Hin) -> ( M!(Hin,ε⁻¹,kpG,Hw), Hin )
#     return LinearMap{ComplexF64}(f,f,N,N,ishermitian=true,ismutating=true)
# end

# function P̂ₖ!(ε::Array{SHM3},ε⁻¹::Array{SHM3},kpG::Array{KVec},Hw::Array{SVector{3,ComplexF64}})::LinearMaps.FunctionMap{ComplexF64}
#     N = 3 * length(kpG)
#     f = (Hout, Hin) -> ( P!(Hin,ε,kpG,Hw), Hin )
#     fc = (Hout, Hin) -> ( M!(Hin,ε⁻¹,kpG,Hw), Hin )
#     return LinearMap{ComplexF64}(f,f,N,N,ishermitian=true,ismutating=true)
# end

# # f = x::Vector{ComplexF64}-> M!(x,ε⁻¹,KpG,Hw)
# # vals,vecs,info = eigsolve(f,H₀)

# function ftest()
#     Hw = copy(HSA)
#     Mop = M̂ₖ!(ε⁻¹,ds.kpG,Hw)
#     Pop = P̂ₖ!(ε,ε⁻¹,ds.kpG,Hw)
#     res = IterativeSolvers.lobpcg(Mop,false,1;maxiter=3000,tol=1e-7)
# end
# ftest()
# ##

# # function foo(x)
# #     x .= M!(x,ε⁻¹,ds.kpG,HSA)
# #     return x
# # end
# @show Hv[1:6]
# # foo(Hv);
# Mop * Hv;
# # M!(Hv,ε⁻¹,ds.kpG,HSA);
# @show Hv[1:6];
# ##
# M!(Hv,ε⁻¹,KpG,HSA)
# Hv[1:6]
# ##
# res = IterativeSolvers.lobpcg(Mop,false,1;P=Pop,maxiter=3000,tol=1e-7)

# ##
# # f = x::Vector{ComplexF64}-> M!(x,ε⁻¹,ds.kpG,Hw)

# # function foo(x::Vector{ComplexF64})::Vector{ComplexF64}
# #     x .= M!(x,ε⁻¹,ds.kpG,Hw)
# #     return x
# # end

# # vals,vecs,info = eigsolve(Mop,HSA)


text:
module RocketTupleWithOperatorTest

using Test
using Rocket

include("../test_helpers.jl")

@testset "operator: tuple_with_*()" begin

    println("Testing: operator tuple_with_*()")

    run_proxyshowcheck("TupleWithLeft", tuple_with_left(0))
    run_proxyshowcheck("TupleWithRight", tuple_with_right(0))

    run_testset([
        (
            source = from(1:3) |> tuple_with_left('0'),
            values = @ts([ ('0', 1), ('0', 2), ('0', 3), c ]),
            source_type = Tuple{Char, Int}
        ),
        (
            source = completed(Int) |> tuple_with_left('0'),
            values = @ts(c),
            source_type = Tuple{Char, Int}
        ),
        (
            source = never(Int) |> tuple_with_left('0'),
            values = @ts(),
            source_type = Tuple{Char, Int}
        ),
        (
            source = faulted("e") |> tuple_with_left('0'),
            values = @ts(e("e"))
        )
    ])

    run_testset([
        (
            source = from(1:3) |> tuple_with_right('0'),
            values = @ts([ (1, '0'), (2, '0'), (3, '0'), c ]),
            source_type = Tuple{Int, Char}
        ),
        (
            source = completed(Int) |> tuple_with_right('0'),
            values = @ts(c),
            source_type = Tuple{Int, Char}
        ),
        (
            source = never(Int) |> tuple_with_right('0'),
            values = @ts(),
            source_type = Tuple{Int, Char}
        ),
        (
            source = faulted("e") |> tuple_with_right('0'),
            values = @ts(e("e"))
        )
    ])

end

end


text:
### A Pluto.jl notebook ###
# v0.12.21

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ 7d550d46-7dc4-11eb-3df0-db7ad0594043
using PlutoUI

# ╔═╡ 8399ed34-7dc4-11eb-16ec-e572834e149d
@bind x Slider(1:10)

# ╔═╡ 88c83d42-7dc4-11eb-1b5e-a1ecad4b6ff1
x

# ╔═╡ Cell order:
# ╠═7d550d46-7dc4-11eb-3df0-db7ad0594043
# ╠═8399ed34-7dc4-11eb-16ec-e572834e149d
# ╠═88c83d42-7dc4-11eb-1b5e-a1ecad4b6ff1


text:
@testset "260.single-number-iii.jl" begin
    r1 = single_number_iii([1, 2, 1, 3, 2, 5])
    @test r1 == (3, 5) || r1 == (5, 3)
    r1 = single_number_iii([-1, 0])
    @test r1 == (-1, 0) || r1 == (0, -1)
    r1 = single_number_iii([1, 0])
    @test r1 == (1, 0) || r1 == (0, 1)
end


text:
module CommonSubexpressions

using MacroTools: @capture, postwalk, MacroTools
using Base.Iterators: drop

export @cse, cse, @binarize

struct Cache
    args_to_symbol::Dict{Symbol, Symbol}
    disqualified_symbols::Set{Symbol}
    setup::Vector{Expr}
end

Cache() = Cache(Dict{Symbol,Symbol}(), Set{Symbol}(), Vector{Expr}())

function add_element!(cache::Cache, name, expr::Expr)
    sym = gensym()
    cache.args_to_symbol[name] = sym
    push!(cache.setup, :($sym = $(expr)))
    sym
end

disqualify!(cache::Cache, x) = nothing
disqualify!(cache::Cache, s::Symbol) = push!(cache.disqualified_symbols, s)
disqualify!(cache::Cache, expr::Expr) = foreach(arg -> disqualify!(cache, arg), expr.args)

# fallback for non-Expr arguments
combine_subexprs!(setup, x; warn=true, mod=nothing) = x

const standard_expression_forms = Set{Symbol}(
    (:call,
     :block,
     :comprehension,
     :.,
     :(=>),
     :(:),
     :(&),
     :(&&),
     :(|),
     :(||),
     :tuple,
     :for,
     :ref,
     Symbol("'")))

const assignment_expression_forms = Set{Symbol}(
    (:(=),
     :(+=),
     :(-=),
     :(*=),
     :(/=)))

function combine_subexprs!(cache::Cache, expr::Expr;
        warn::Bool=true, mod::Union{Module, Nothing}=nothing)
    if expr.head == :macrocall
        if (mod === nothing)
            error("""
                `cse` cannot expand macro calls unless you explicitly pass in
                a `Module` in which to perform that expansion. You can pass
                `mod=@__MODULE__` to expand in the current module, or you can use
                the `@cse` macro which handles this automatically.""")
        end
        return combine_subexprs!(cache, macroexpand(mod, expr);
                                 warn=warn, mod=mod)
    elseif expr.head == :function
        # We can't continue CSE through a function definition, but we can
        # start over inside the body of the function:
        for i in 2:length(expr.args)
            expr.args[i] = combine_subexprs!(expr.args[i]; warn=warn, mod=mod)
        end
    elseif expr.head == :line
        # nothing
    elseif expr.head in assignment_expression_forms
        disqualify!(cache, expr.args[1])
        for i in 2:length(expr.args)
            expr.args[i] = combine_subexprs!(cache, expr.args[i]; warn=warn, mod=mod)
        end
    elseif expr.head == :generator
        for i in vcat(2:length(expr.args), 1)
            expr.args[i] = combine_subexprs!(cache, expr.args[i]; warn=warn, mod=mod)
        end
    elseif expr.head in standard_expression_forms
        for (i, child) in enumerate(expr.args)
            expr.args[i] = combine_subexprs!(cache, child; warn=warn, mod=mod)
        end
        if expr.head == :call
            for (i, child) in enumerate(expr.args)
                expr.args[i] = combine_subexprs!(cache, child, warn=warn, mod=mod)
            end
            if all(!isa(arg, Expr) && !(arg in cache.disqualified_symbols) for arg in drop(expr.args, 1))
                combined_args = Symbol(expr.args...)
                if !haskey(cache.args_to_symbol, combined_args)
                    sym = add_element!(cache, combined_args, expr)
                else
                    sym = cache.args_to_symbol[combined_args]
                end
                return sym
            else
            end
        end
    else
        warn && @warn("CommonSubexpressions can't yet handle expressions of this form: $(expr.head)")
    end
    return expr
end

combine_subexprs!(x; warn=true, mod=nothing) = x

function combine_subexprs!(expr::Expr; warn=true, mod=nothing)
    cache = Cache()
    expr = combine_subexprs!(cache, expr; warn=warn, mod=mod)
    Expr(:block, cache.setup..., expr)
end

function parse_cse_args(args)
    # Overly complicated way to look for `warn=true` or `warn=false`,
    # but should be easier to expand for other arguments later.
    params = Dict(:warn => true)
    for (i, arg) in enumerate(args)
        if @capture(arg, key_Symbol = val_Bool)
            if key in keys(params)
                params[key] = val
            else
                error("Unrecognized key: $key")
            end
        elseif i == 1 && arg isa Bool
            Base.depwarn("The `warn_enabled` positional argument is deprecated. Please use `warn=true` or `warn=false` instead", :cse_macro_kwargs)
        else
            error("Unrecognized argument: $arg. Expected `warn=true` or `warn=false`")

        end
    end
    params
end

"""
    @cse(expr; warn=true)

Perform naive common subexpression elimination under the assumption
that all functions called withing the body of the macro are pure,
meaning that they have no side effects. See [Readme.md](https://github.com/rdeits/CommonSubexpressions.jl/blob/master/Readme.md)
for more details.

This macro will recursively expand macro calls within the expression before
performing subexpression elimination. A useful macro to combine with this is
`@binarize`, which will turn n-ary function calls into nested binary calls and
can therefore provide more opportunities for subexpression elimination. Usage:

    @cse(@binarize(<your code here>))

If the macro encounters an expression which it does not know how to handle,
it will return that expression unmodified. If `warn=true`, then it
will also log a warning in that event.
"""
macro cse(expr, args...)
    params = parse_cse_args(args)
    result = combine_subexprs!(expr, warn=params[:warn], mod=__module__)
    esc(result)
end

Base.@deprecate cse(expr, warn_enabled::Bool) cse(expr, warn=warn_enabled)

function cse(expr; warn::Bool=true, mod::Union{Module, Nothing}=nothing)
    combine_subexprs!(copy(expr); warn=warn, mod=mod)
end

function _binarize(expr::Expr)
    if @capture(expr, f_(a_, b_, c_, args__))
        _binarize(:($f($f($a, $b), $c, $(args...))))
    else
        expr
    end
end

_binarize(x) = x

binarize(expr::Expr) = postwalk(_binarize, expr)
binarize(x) = x

"""
    @binarize(expr::Expr)

Convery all n-ary function calls within the given expression to nested binary
calls to the same function. That is, convert all calls of the form `f(a, b, c)`
to `f(f(a, b), c)` with as many layers of nesting as necessary. Operators like
`+` and `*` are handled just like any other function call, so

    @binarize a + b + c + d

will produce:

    ((a + b) + c) + d

This is intended to make subexpression elimination easier for long chained
function calls, such as https://github.com/rdeits/CommonSubexpressions.jl/issues/14
"""
macro binarize(expr::Expr)
    esc(binarize(expr))
end

end


text:
using Documenter, ScenTrees

#const ASSETS = readdir(joinpath(@__DIR__, "src", "assets"))

makedocs(
	sitename =  "ScenTrees.jl",
	authors = "Kipngeno Kirui",
	clean = true,
	doctest = false,
	format = Documenter.HTML(prettyurls = get(ENV, "CI", nothing) == "true"),
	pages = ["Home" => "index.md",
		"Tutorials" => Any["tutorial/tutorial1.md",
				    "tutorial/tutorial2.md",
				    "tutorial/tutorial3.md",
				    "tutorial/tutorial31.md",
				    "tutorial/tutorial4.md",
				    "tutorial/tutorial41.md",
				    "tutorial/tutorial5.md"
				]
		]
)

deploydocs(
	repo = "github.com/kirui93/ScenTrees.jl.git"
	#target = "build",
	#versions = ["stable" => "v^", "v#.#", "dev" => "master"]
)


text:
## This file autogenerated by BinaryProvider.write_deps_file().
## Do not edit.
##
## Include this file within your main top-level source, and call
## `check_deps()` from within your module's `__init__()` method

if isdefined((@static VERSION < v"0.7.0-DEV.484" ? current_module() : @__MODULE__), :Compat)
    import Compat.Libdl
elseif VERSION >= v"0.7.0-DEV.3382"
    import Libdl
end
const cclipper = joinpath(dirname(@__FILE__), "usr/lib/cclipper.so")
function check_deps()
    global cclipper
    if !isfile(cclipper)
        error("$(cclipper) does not exist, Please re-run Pkg.build(\"Clipper\"), and restart Julia.")
    end

    if Libdl.dlopen_e(cclipper) == C_NULL
        error("$(cclipper) cannot be opened, Please re-run Pkg.build(\"Clipper\"), and restart Julia.")
    end

end


text:
module DaemonMode

export serve, runfile, runargs, runexpr, sendExitCode

using Crayons
using Crayons.Box

using Sockets

const PORT = 3000

function add_packages(fname::AbstractString)
end

const first_time = Ref(true)
const token_runfile = "DaemonMode::runfile"
const token_runexpr = "DaemonMode::runexpr"
const token_exit = "DaemonMode::exit"
const token_end = "DaemonMode::end"

"""
    serve(port=3000, shared=false)

Run the daemon, running all files and expressions sended by the client function.

# Optionals

- port: port to listen (default=3000).
- shared: Share the environment between calls. If it is false (default) each run
  has its own environment, so the variables/functions are not shared.
- print_stack: Print the complete stack when there is an error. By default it is true.
"""
function serve(port=PORT, shared=missing; print_stack=true)
    server = Sockets.listen(Sockets.localhost, port)
    quit = false
    current = pwd()

    while !quit
        sock = accept(server)
        mode = readline(sock)

        if mode == token_runfile
            serverRunFile(sock, coalesce(shared, false), print_stack)
        elseif mode == token_runexpr
            serverRunExpr(sock, coalesce(shared, true), print_stack)
        elseif mode == token_exit
            println(sock, token_end)
            sleep(1)
            quit = true
        else
            println(sock, "Error, unrecognised mode, expected (\"runFile()\", \"runExpr()\" or \"exit()\", but received \"$mode\"")
            quit = true
        end

    end
    close(server)
end

function myshowerror(sock::IO, e)
    io = IOBuffer()
    showerror(io, e)
    msg = String(take!(io))
    line = first(split(msg, '\n'))
    posi = findfirst(':', line)
    print(sock, RED_FG, BOLD, "ERROR: ")
    print(sock, Crayon(foreground=:white, bold=false), line[1:posi])
    print(sock, Crayon(foreground=:red, bold=false), line[posi+1:end])
    print(sock, Crayon(reset=true))
end

function serverReplyError(sock, e)
    try
        myshowerror(sock, e)
        println(sock)
        println(sock, token_end)
    catch e
        if (e isa Base.IOError) && abs(e.code) == abs(Libc.EPIPE)
            # client disconnected early, ignore
        else
            rethrow()
        end
    end
end

function send_backtrace(sock, bt, fname)
    stacks = map(first, Base.process_backtrace(bt))
    fullname = ""

    if !isempty(stacks)
        println(sock)
        println(sock, "Stacktrace:")
    end

    if !isempty(fname)
        fullname = joinpath(pwd(), fname)
    end

    for (i, stack) in enumerate(stacks)
        file = String(stack.file)

        if occursin("string", file) && !isempty(fullname)
            file = fullname
        end

        if occursin("loading.jl", file)
            return
        end

        print(sock, " [$i] ")

        print(sock, BOLD, stack.func)
        print(sock, Crayon(reset=true), " at ")
        print(sock, BOLD, file, ":", stack.line)
        println(sock, Crayon(reset=true))
    end
end

function serverReplyError(sock, e, bt, fname)
    try
        myshowerror(sock, e)
        send_backtrace(sock, bt, fname)
        println(sock)
        println(sock, token_end)
    catch e
        if (e isa Base.IOError) && abs(e.code) == abs(Libc.EPIPE)
            # client disconnected early, ignore
        else
            rethrow()
        end
    end
end

function serverRun(run, sock, shared, print_stack, fname)
    redirect_stdout(sock) do
        redirect_stderr(sock) do

            try
                if shared
                    run(Main)
                else
                    m = Module()
                    add_include = Meta.parse("include(arg)=Base.include(@__MODULE__,arg)")
                    # Base.eval(m, "Base.exit(x)=return x")
                    Base.eval(m, add_include)
                    run(m)
                end
                println(sock, token_end)

            catch e
                if print_stack
                    serverReplyError(sock, e, catch_backtrace(), fname)
                else
                    serverReplyError(sock, e)
                end
            end

        end
    end

end

"""
serverRunFile(sock, shared)

Run the source code of the filename push through the socket.

# Parameters

- sock: socket in which is going to receive the dir, the filename, and args to run.
- shared: Share the environment between calls. If it is false (default) each run
has its own environment, so the variables/functions are not shared.
"""
function serverRunFile(sock, shared, print_stack)
    fname = ""

    try
        dir = readline(sock)
        fname = readline(sock)
        args_str = readline(sock)

        args = split(args_str, " ")
        append!(ARGS, args)

        first_time[] = true

        cd(dir) do
            content = read(fname, String)
            serverRun(sock, shared, print_stack, fname) do mod
                include_string(mod, content)
            end
        end
    catch e
        serverReplyError(sock, e)
    end

    empty!(ARGS)
end

"""
serverRunExpr(sock, shared)

Run the source code of the filename push through the socket.

# Parameters

- sock: socket in which is going to receive the code to run.
- shared: Share the environment between calls. If it is false (default) each run
has its own environment, so the variables/functions are not shared.
"""
function serverRunExpr(sock, shared, print_stack)

    try
        dir = readline(sock)
        expr = readuntil(sock, token_end) # Read until token_end to handle multi-line expressions
        parsedExpr = Meta.parse(strip(expr))

        cd(dir) do
            serverRun(sock, shared, print_stack, "") do mod
                Base.eval(mod, parsedExpr)
            end
        end
    catch e
        serverReplyError(sock, e)
    end

    empty!(ARGS)
end

"""
    runexpr(expr::AbstractString)

Ask the server to run julia code in a string pass as parameters.

# Parameters

- expr: Julia code to run in the server.

# Optionals

- port: Port (default=3000).
- output: stream in which it is shown the output of the run.
"""
function runexpr(expr::AbstractString ; output = stdout, port = PORT)
    try
        sock = Sockets.connect(port)
        println(sock, token_runexpr)
        println(sock, pwd())
        println(sock, expr)
        println(sock, token_end)

        line = readline(sock)
        while (line != token_end)
            println(output, line)
            line = readline(sock)
        end

    catch e
        println(stderr, "Error, cannot connect with server. Is it running?")
    end
end

"""
    runfile(fname::AbstractString)

Ask the server to run a specific filename.

# Parameters

- fname: Filename to run.

# Optionals

- args: List of arguments (array of String, default=[]).
- port: Port (default=3000)
- output: stream in which it is shown the output of the run.
"""
function runfile(fname::AbstractString; args=String[], port = PORT, output=stdout)
    dir = dirname(fname)

    if isempty(dir)
        fcompletename = joinpath(dir, fname)
    else
        fcompletename = fname
    end
    try
        sock = Sockets.connect(port)
        println(sock, token_runfile)
        println(sock, pwd())
        println(sock, fcompletename)
        println(sock, join(args, " "))
        line = readline(sock)
        token_size = length(token_end)

        while (length(line) < token_size || !occursin(token_end, line))
            println(output, line)
            line = readline(sock)
        end

        if length(line) > token_size
            end_line = replace(line, token_end => "")
            print(output, end_line)
        end
    catch e
        println(stderr, "Error, cannot connect with server. Is it running?")
        exit(1)
    end
    return
end

"""
    sendExitCode(port)

send the exit code, it closes the server.

# Optionals

- port: port to connect (default=3000).
"""
function sendExitCode(port=PORT)
    sock = Sockets.connect(port)
    println(sock, token_exit)
end

"""
    runargs(port=PORT)

Ask the server to run all files in ARGS.

# Optionals

- port: Port of the server to ask (default=3000).
"""
function runargs(port=PORT)
    if isempty(ARGS)
        println(stderr, "Error: missing filename")
        exit(1)
    elseif !isfile(ARGS[1])
        println(stderr, "Error: file '$(ARGS[1])' doest not exist")
        exit(1)
    end

    runfile(ARGS[1], args=ARGS[2:end], port=port)
end

end # module