text:
#=##############################################################################
# DESCRIPTION
    Vahana vehicle geometry generation.

# AUTHORSHIP
  * Author    : Eduardo J. Alvarez
  * Email     : Edo.AlvarezR@gmail.com
  * Created   : Feb 2021
  * License   : MIT
=###############################################################################

#=

Geometry functions (e.g., `generategeometry_vahana`) are expected have the
following format:

    OUTPUTS
    * `system::vlm.WingSystem`:        System of all FLOWVLM objects.
    * `rotors::Array{vlm.Rotor, 1}`:   System of all FLOWVLM Rotor objects.
    * `tilting_systems::Tuple(vlm.WingSystem, ...)`:   Tuple of all FLOWVLM
                                        tilting objects, where
                                        `tilting_systems[i]` contains the i-th
                                        FLOWVLM system of lifting surfaces and
                                        rotors that tilt together.
    * `rotor_systems::Tuple(Array{vlm.Rotor,1}, ...)`:   Tuple of groups of
                                        Rotors that share a common RPM.
    * `vlm_system::vlm.WingSystem`:    System of all FLOWVLM objects to be
                                        solved through the VLM solver.
    * `wake_system::vlm.WingSystem`:   System of all FLOWVLM objects that will
                                        shed a VPM wake.
    * `fuselage::gt.Grid`:             Grid object of the fuselage.
    * `grounds::gt.Grid`:              Grid object of the ground.
    * `strn::String`                   String containing all vtk file names
                                        that will be generated if `system` is
                                        saved.

    NOTES
    * Maneuver functions are expected to return `(Vaircraft, angles, RPMs)`,
        where `Vaircraft(t)` is the translation velocity of both `system`
        and `fuselage`, `angles(t)[i]` is the tilt angle of
        `tilting_systems[i]`, `RPMs(t)[i]` is the RPM of
        `rotor_systems[i]`. `angles(t)[end]` is the tilt angle of the
        entire aircraft.
=#

"""
    Generates the geometry of Vahana aircraft
"""
function generate_geometry_vahana(;
                                    # AIRCRAFT OPTIONS
                                    tiltrotor_file="VahanaTilt.csv", # Tilt-rotor
                                    stackedrotor_file="VahanaStacked.csv", # Stacked-rotor
                                    data_path=uns.def_data_path,
                                    xfoil=true,                     # Run XFOIL
                                    n_factor::Int=1,                # Refinement factor
                                    add_rotors=true,
                                    VehicleType=uns.VLMVehicle,     # Type of vehicle to generate
                                    # OUTPUT OPTIONS
                                    run_name="vahana",
                                    verbose=true, v_lvl=0)

    ############################################################################
    # PARAMETERS
    ############################################################################

    if verbose; println("\t"^(v_lvl)*"Defining parameters..."); end;

    # ------------ GEOMETRIC PARAMETERS ------------------------------------
    init_ori = 90.0               # Initial orientation of wings
    init_ori_prop = 0*90.0        # Initial orientation of rotors
    init_ori_stackedprop = 30.0   # Initial orientation of stacked rotors

    # Rotors
    # NOTE: The radii defined here must match the radius in rotor_file
    R_w = 0.75                    # (m) main wing rotor radius
    R_tw = R_w                    # (m) tandem wing rotor radius
    CW_w = true                   # Clockwise main wing rotor
    CW_tw = false                 # Clockwise tandem wing rotor
    np_w = 2                      # Number of propellers on wing on each side
    np_tw = 2                     # Number of propellers on tandem wing each side
    n_ccb = 7*n_factor            # Number of CCBlade elements
    blade_r = 1/20                # Geometric expansion between blade elements
    ReD07 = 1.5e6                 # Prop Reynolds number of relative velocity at r/R=0.7
    ReD = ReD07/0.7               # Reynolds at tip
    pitch = 0.0                   # (deg) pitch of propellers
    soD = 0.1                     # Tip-to-tip distance of rotors over diameter
    xoc_offset_main = 0.175       # Axial distance of rotors from LE over chord
    xoc_offset_tandem = 0.10      # Axial distance of rotors from LE over chord
    main_outtilt = 10             # Mounted out-tilt of tandem-wing rotors
    tandem_pitchtilt = 10         # Mounted pitch-tilt of tandem-wing rotors

    # Stacked rotors
    stacked = [np_w, np_w+1]      # Index of wing rotors that are stacked rotors
    stckd_xoc_offset = -1.40      # Axial distance of rotors from LE over chord
    stckd_zoc_offset = 0.10       # Distance of stacked rotors above wing over chord
    stckd_zoR_offset = -0.05      # Stacking distance between stacked rotors over R
    stckd_corotating = true       # Co-rotating stacked rotors
    stckd_phase = -10             # Initial phase difference of stacked rotors
    stckd_pitch = 7.5             # Pitch difference of stacked rotors

    # Wing
    b_w = 5.86                    # (m) span
    AR_w = 7.4                    # Aspect ratio
    tr_w = 1.0                    # Taper ratio
    twist_r_w = 7.5               # (deg) twist at root
    twist_t_w = twist_r_w         # (deg) twist at tip
    lambda_w = main_outtilt       # (deg) sweep
    gamma_w = 5.0                 # (deg) dihedral
    n_w = 2*12*n_factor           # Number of horseshoes per side of wing
    r_w = 2.0                     # Horseshoe expansion ratio
    md_w = 0.9                    # Ratio of length of middle section
    pivot_w = 1/4                 # Pivot point along chord of tilt-wing

    # Wing winglets
    b_wl = b_w/4                  # (m) span of winglet from top to bottom
    AR_wl = 3.0                   # Aspect ratio
    tr_wl = (b_wl/AR_wl)/(b_w/AR_w/tr_w)      # Taper ratio
    twist_r_wl = 2.5              # (deg) twist at root
    twist_t_wl = 0.0              # (deg) twist at tip
    lambda_wl = 40.0              # (deg) sweep
    gamma_wl = 15.0               # (deg) dihedral
    n_wl = 2*4*n_factor           # Number of horseshoes per side of wing
    r_wl = 2.0                    # Horseshoe expansion ratio

    # Tandem wing
    b_tw = b_w*1.0                # (m) span
    AR_tw = 9.5                   # Aspect ratio
    tr_tw = 1.0                   # Taper ratio
    twist_r_tw = 4.0              # (deg) twist at root
    twist_t_tw = twist_r_tw       # (deg) twist at tip
    lambda_tw = 0.0               # (deg) sweep
    gamma_tw = 0.0                # (deg) dihedral
    n_tw = 2*n_w                  # Number of horseshoes per side of wing
    r_tw = r_w                    # Horseshoe expansion ratio
    md_tw = 0.2                   # Ratio of length of middle section
    pivot_tw = pivot_w            # Pivot point along chord of tilt-wing

    # Fuselage
    l_f = 5.86                    # (m) length
    h_f = 2.81*2/3                # (m) height
    c_pos_f1 = 0                  # (m) span position of first chord
    c_pos_f2 = 0.2*h_f            # (m) span position of second chord
    c_pos_f3 = 0.7*h_f
    c_pos_f4 = h_f
    c_f1 = 0.66*l_f               # (m) length of first chord
    c_f2 = 0.75*l_f
    c_f3 = 0.75*l_f
    c_f4 = 0.50*l_f
    x_pos_f1 = 0.025*l_f          # (m) x-position of first chord
    x_pos_f2 = 0.0*l_f
    x_pos_f3 = 0.20*l_f
    x_pos_f4 = 0.50*l_f

    # ------------ ASSEMBLY PARAMETERS ------------------------------------
    # Position of wings on fuselage
    h_pos_w = 0.90*h_f            # (m) height position of wing
    h_pos_tw = 0.15*h_f           # (m) height position of tandem wing
    l_pos_w = 0.95*l_f-b_w/AR_w/tr_w   # (m) length position of wing
    l_pos_tw = 0.05*l_f           # (m) length position of tandem wing

    # Position of propellers on main wing
    # d_prop_w = d_factor*b_w/2/np_w# Distance between propellers on wing
    d_prop_w = (1+soD)*(2*R_w)
    y_pos_prop_w = Float64[b_w/2 - i*d_prop_w for i in 0:np_w-1] # y-positions

    # Position of propellers on tandem wing
    # d_prop_tw = d_factor*b_tw/2/np_tw# Distance between propellers on wing
    d_prop_tw = (1+soD)*(2*R_tw)
    y_pos_prop_tw = Float64[b_tw/2 - i*d_prop_tw for i in 0:np_tw-1] # y-positions

    ############################################################################
    # ASSEMBLY
    ############################################################################

    if verbose; println("\t"^(v_lvl)*"Generating geometry..."); end;

    # ------------ ROTORS ------------------------------------------------
    # Generates base propellers (one on each rotation orientation)
    if add_rotors
        tiltrotors = vlm.Rotor[]
        if verbose; println("\t"^(v_lvl+1)*"Generating first tilt-rotor..."); end;
        push!(tiltrotors, uns.generate_rotor(tiltrotor_file; pitch=pitch,
                                                n=n_ccb, blade_r=blade_r, CW=!CW_w, ReD=ReD,
                                                verbose=verbose, v_lvl=v_lvl+2, xfoil=xfoil,
                                                data_path=data_path, plot_disc=false))
        if verbose; println("\t"^(v_lvl+1)*"Generating second tilt-rotor..."); end;
        # push!(tiltrotors, generate_rotor(pitch; n=n_ccb, blade_r=blade_r, CW=CW_w, ReD=ReD,
        #                         verbose=verbose, xfoil=xfoil, rotor_file=tiltrotor_file))
        push!(tiltrotors, vlm.Rotor(!tiltrotors[1].CW, tiltrotors[1].r,
                                  tiltrotors[1].chord, tiltrotors[1].theta,
                                  tiltrotors[1].LE_x, tiltrotors[1].LE_z,
                                  tiltrotors[1].B, tiltrotors[1].airfoils))
        vlm.initialize(tiltrotors[2], tiltrotors[1].m)

        stackedrotors = vlm.Rotor[]
        if verbose; println("\t"^(v_lvl+1)*"Generating first stacked-rotor..."); end;
        push!(stackedrotors, uns.generate_rotor(stackedrotor_file; pitch=pitch,
                                                n=n_ccb, blade_r=blade_r, CW=!CW_w, ReD=ReD,
                                                verbose=verbose, v_lvl=v_lvl+2, xfoil=xfoil,
                                                data_path=data_path, plot_disc=false))
        if verbose; println("\t"^(v_lvl+1)*"Generating second stacked-rotor..."); end;
        # push!(stackedrotors, generate_rotor(pitch; n=n_ccb, blade_r=blade_r, CW=CW_w, ReD=ReD,
        #                         verbose=verbose, xfoil=xfoil, rotor_file=stackedrotor_file))
        push!(stackedrotors, vlm.Rotor(!stackedrotors[1].CW, stackedrotors[1].r,
                                  stackedrotors[1].chord, stackedrotors[1].theta,
                                  stackedrotors[1].LE_x, stackedrotors[1].LE_z,
                                  stackedrotors[1].B, stackedrotors[1].airfoils))
        vlm.initialize(stackedrotors[2], stackedrotors[1].m)
    end

    # ------------ MAIN WING ---------------------------------------------
    # Generates wing
    if verbose; println("\t"^(v_lvl+1)*"Generating main wing assembly..."); end;

    # wing = vlm.simpleWing(b_w, AR_w, tr_w, twist_r_w, lambda_w, gamma_w;
    #                                           twist_tip=twist_t_w, n=n_w, r=r_w)

    # Middle section
    # tip = b/ar
    # root = tip/tr
    # newtip = root + md*(tip-root) = tip/tr + md*tip*(1 - 1/tr) = tip * ( md + (1/tr)*(1-md) )
    # newtr = newtip/root = 1 + md*(tip/root - 1) = 1 + md*(tr - 1)
    # newar = b/newtip = ar / ( md + (1/tr)*(1-md) )
    pos_md_w = [0.0, md_w]
    clen_md_w = [1/tr_w, md_w + (1/tr_w)*(1-md_w)]
    twist_md_w = [twist_r_w, twist_r_w + md_w*(twist_t_w-twist_r_w)]
    wing_md = vlm.complexWing(b_w, AR_w, ceil(Int, md_w*n_w), pos_md_w, clen_md_w, twist_md_w,
                          lambda_w*ones(1), gamma_w*ones(1);
                          symmetric=true, chordalign=0.0,
                          _ign1=true)
    # Left section
    # pos_l_w = [-(1-md_w), -0]
    pos_l_w = [0, -(1-md_w)]   # Here I'm defining this wing section from right to left
    clen_l_w = [1, clen_md_w[end]]
    twist_l_w = [twist_t_w, twist_md_w[end]]
    wing_L = vlm.complexWing(b_w, AR_w, ceil(Int, (1-md_w)*n_w/2), pos_l_w, clen_l_w,
                            twist_l_w, -lambda_w*ones(1), -gamma_w*ones(1);
                            symmetric=false, chordalign=0.0,
                            _ign1=true)
    # Right section
    pos_r_w = [0, (1-md_w)]
    clen_r_w = [clen_md_w[end], 1]
    twist_r_w = [twist_md_w[end], twist_t_w]
    wing_R = vlm.complexWing(b_w, AR_w, ceil(Int, (1-md_w)*n_w/2), pos_r_w, clen_r_w,
                            twist_r_w, lambda_w*ones(1), gamma_w*ones(1);
                            symmetric=false, chordalign=0.0,
                            _ign1=true)

    O_w_R = (md_w*b_w/2)*[0, 1, 0]
    O_w_L = [1 0 0; 0 -1 0; 0 0 1]*O_w_R
    vlm.setcoordsystem(wing_R, O_w_R, Float64[1 0 0; 0 1 0; 0 0 1])
    vlm.setcoordsystem(wing_L, O_w_L, Float64[1 0 0; 0 1 0; 0 0 1])

    # Generates winglets
    winglet_R = vlm.simpleWing(b_wl, AR_wl, tr_wl, twist_r_wl, lambda_wl, gamma_wl;
                                          twist_tip=twist_t_wl, n=n_wl, r=r_wl)
    winglet_L = vlm.simpleWing(b_wl, AR_wl, tr_wl, twist_r_wl, lambda_wl, gamma_wl;
                                          twist_tip=twist_t_wl, n=n_wl, r=r_wl)

    O_wl_R = (b_w/2)*[0, 1, 0]
    O_wl_R += ((1-md_w)*b_w/2)*[tan(lambda_w*pi/180), 0, tan(gamma_w*pi/180)]
    O_wl_L = [1 0 0; 0 -1 0; 0 0 1]*O_wl_R
    Oaxis_wl_R = gt.rotation_matrix(0.0, 0.0, 90.0)
    Oaxis_wl_L = gt.rotation_matrix(0.0, 0.0, -90.0)
    vlm.setcoordsystem(winglet_R, O_wl_R, Oaxis_wl_R)
    vlm.setcoordsystem(winglet_L, O_wl_L, Oaxis_wl_L)

    ## Generates propellers on wing (from right to left)
    if add_rotors
        if verbose; println("\t"^(v_lvl+2)*"Generating main wing propellers..."); end;
        O_prop_w = [
                     (ypos - md_w*b_w/2)*[tan(lambda_w*pi/180), 0, -tan(gamma_w*pi/180)] +
                     ypos*[0, 1, 0] +
                     [-(i in stacked ? stckd_xoc_offset : xoc_offset_main)*AR_w/b_w, 0, 0] +
                     [0, 0, (i in stacked ? stckd_zoc_offset*tan(gamma_w*pi/180) : 0)*AR_w/b_w]
                     for (i, ypos) in enumerate(y_pos_prop_w)]
        props_w = vlm.Rotor[]
        props_w_stacked_up = vlm.Rotor[]
        props_w_stacked_low = vlm.Rotor[]
        for i in 1:2*np_w
            right = i<=np_w    # Indicates which side of the wing

            copy_prop = i in stacked ? stackedrotors[1+i%2] : tiltrotors[1+i%2]
            this_prop = deepcopy(copy_prop) # Alternates rotation orientation

            this_O = O_prop_w[ right ? i : np_w-(i-np_w-1)]   # Chooses position
            this_O = [1 0 0; 0 (-1)^!right 0; 0 0 1]*this_O   # Places it in correct side
            this_Oaxis = i in stacked ? gt.rotation_matrix(0.0, 90, 0.0) :
                                        gt.rotation_matrix((-1)^(i%2==0)*main_outtilt, 0.0, 0.0)

            vlm.setcoordsystem(this_prop, this_O, this_Oaxis; user=true)

            # Rotates props to be tip to tip
            vlm.rotate(this_prop, (-1)^(!CW_w) * (i in stacked ? init_ori_stackedprop : init_ori_prop))

            # Adds the original polars that don't get copied in deepcopy
            this_prop.airfoils = copy_prop.airfoils
            this_prop._polars = copy_prop._polars
            this_prop._polarroot = copy_prop._polarroot
            this_prop._polartip = copy_prop._polartip

            if !(i in stacked)
                push!(props_w, this_prop)

            else
                push!(props_w_stacked_up, this_prop)

                copy_prop = deepcopy(stackedrotors[ 1+(i+1*!stckd_corotating)%2 ])

                # Generate lower rotor of stack and add pitch
                if stckd_pitch != 0
                    this_prop = vlm.Rotor(copy_prop.CW, copy_prop.r,
                                              copy_prop.chord, copy_prop.theta .+ stckd_pitch,
                                              copy_prop.LE_x, copy_prop.LE_z,
                                              copy_prop.B, copy_prop.airfoils)
                    vlm.initialize(this_prop, copy_prop.m)
                else
                    this_prop = copy_prop
                end

                this_O += R_w*[0, 0, stckd_zoR_offset]
                vlm.setcoordsystem(this_prop, this_O, this_Oaxis; user=true)

                # Rotates props to be tip to tip
                vlm.rotate(this_prop, (-1)^(!CW_w) * (init_ori_stackedprop+stckd_phase))

                # Adds the original polars that don't get copied in deepcopy
                this_prop.airfoils = copy_prop.airfoils
                this_prop._polars = copy_prop._polars
                this_prop._polarroot = copy_prop._polarroot
                this_prop._polartip = copy_prop._polartip

                push!(props_w_stacked_low, this_prop)
            end
        end
    end

    # Assembles fixed section of the wing
    main_wing_fixed = vlm.WingSystem()
    vlm.addwing(main_wing_fixed, "WingM", wing_md)
    if add_rotors
        for (i, prop) in enumerate(props_w_stacked_up)
            vlm.addwing(main_wing_fixed, "StackedPropUp$i", prop)
        end
        for (i, prop) in enumerate(props_w_stacked_low)
            vlm.addwing(main_wing_fixed, "StackedPropLow$i", prop)
        end
    end

    main_wing_R = vlm.WingSystem()
    vlm.addwing(main_wing_R, "Tip", wing_R)
    vlm.addwing(main_wing_R, "Winglet", winglet_R)

    main_wing_L = vlm.WingSystem()
    vlm.addwing(main_wing_L, "Tip", wing_L)
    vlm.addwing(main_wing_L, "Winglet", winglet_L)

    # Assembles moving sections of the wing
    main_wing_moving = vlm.WingSystem()
    vlm.addwing(main_wing_moving, "WingR", main_wing_R)
    vlm.addwing(main_wing_moving, "WingL", main_wing_L)
    if add_rotors
        for (i, prop) in enumerate(props_w)
            vlm.addwing(main_wing_moving, "Prop$i", prop)
        end
    end

    # offset to align with pivot line
    x_off_w = pivot_w*b_w/AR_w

    # Aligns moving and fixed section with their pivot line
    O_off_w = [-x_off_w, 0.0, 0.0]

    vlm.setcoordsystem(wing_md, O_off_w, Array(1.0I, 3, 3))

    for vlmwing in main_wing_moving.wings
        if typeof(vlmwing)==vlm.Rotor
            vlm.setcoordsystem(vlmwing, vlmwing._wingsystem.O + O_off_w,
                                        vlmwing._wingsystem.Oaxis; user=false)
        else
            vlm.setcoordsystem(vlmwing, vlmwing.O + O_off_w, vlmwing.Oaxis)
        end
    end

    # Places tilting sections at the main wing tip
    O_mv = (md_w*b_w/2)*[tan(lambda_w*pi/180), 0, tan(gamma_w*pi/180)]

    # Initial rotation of moving sections
    Oaxis_wmv = gt.rotation_matrix(0.0, -init_ori, 0.0)
    vlm.setcoordsystem(main_wing_moving, O_mv, Oaxis_wmv)

    # Assembles main wing
    main_wing = vlm.WingSystem()
    vlm.addwing(main_wing, "Fixed", main_wing_fixed)
    vlm.addwing(main_wing, "Moving", main_wing_moving)

    # Position of main wing
    O_w = [l_pos_w + x_off_w, 0, h_pos_w]
    Oaxis_w = gt.rotation_matrix(0.0, 0.0, 0.0)
    vlm.setcoordsystem(main_wing, O_w, Oaxis_w)

    # ------------ TANDEM WING -------------------------------------------
    if verbose; println("\t"^(v_lvl+1)*"Generating tandem wing assembly..."); end;

    # Generates tandem wing
    twing = vlm.simpleWing(b_tw, AR_tw, tr_tw, twist_r_tw, lambda_tw,
                                gamma_tw; twist_tip=twist_t_tw, n=n_tw, r=r_tw)

    # Middle section
    pos_md_tw = [-md_tw, 0.0, md_tw]
    clen_md_tw = [md_tw + (1/tr_tw)*(1-md_tw), 1/tr_tw,
                                                   md_tw + (1/tr_tw)*(1-md_tw)]
    twist_md_tw = [twist_r_tw + md_tw*(twist_t_tw-twist_r_tw),
                        twist_r_tw, twist_r_tw + md_tw*(twist_t_tw-twist_r_tw)]
    twing_md = vlm.complexWing(b_tw, AR_tw, ceil(Int, md_tw*n_tw), pos_md_tw,
                                  clen_md_tw, twist_md_tw,
                                  lambda_tw*ones(2), gamma_tw*ones(2);
                                  symmetric=false, chordalign=0.0,
                                  _ign1=true)
    # Left section
    pos_l_tw = [-1, pos_md_tw[1]]
    clen_l_tw = [1, clen_md_tw[1]]
    twist_l_tw = [twist_t_tw, twist_md_tw[1]]
    twing_L = vlm.complexWing(b_tw, AR_tw, ceil(Int, (1-md_tw)*n_tw/2), pos_l_tw,
                                    clen_l_tw, twist_l_tw,
                                    lambda_tw*ones(1), gamma_tw*ones(1);
                                    symmetric=false, chordalign=0.0,
                                    _ign1=true)
    # Right section
    pos_r_tw = [pos_md_tw[end], 1]
    clen_r_tw = [clen_md_tw[end], 1]
    twist_r_tw = [twist_md_tw[end], twist_t_tw]
    twing_R = vlm.complexWing(b_tw, AR_tw, ceil(Int, (1-md_tw)*n_tw/2), pos_r_tw,
                                    clen_r_tw, twist_r_tw,
                                    lambda_tw*ones(1), gamma_tw*ones(1);
                                    symmetric=false, chordalign=0.0,
                                    _ign1=true)

    ## Generates propellers on tandem wing (from right to left)
    if add_rotors
        if verbose; println("\t"^(v_lvl+2)*"Generating tandem wing propellers..."); end;
        O_prop_tw = [ ypos*[tan(lambda_tw*pi/180), 1, tan(gamma_tw*pi/180)] +
                        [-xoc_offset_tandem*AR_tw/b_tw, 0, 0]
                                                  for ypos in y_pos_prop_tw]
        props_tw = vlm.Rotor[]
        for i in 1:2*np_tw
            right = i<=np_tw    # Indicates which side of the wing
            copy_prop = tiltrotors[1+(i+(CW_tw!=CW_w))%2]
            this_prop = deepcopy(copy_prop)       # Alternates rotation orientation
            this_O = O_prop_tw[ right ? i : np_tw-(i-np_tw-1)] # Chooses position
            this_O = [1 0 0; 0 (-1)^!right 0; 0 0 1]*this_O   # Places it in correct side
            this_Oaxis = gt.rotation_matrix(0, (-1)^(i%2==0)*(-1)^right*tandem_pitchtilt, 0)

            vlm.setcoordsystem(this_prop, this_O, this_Oaxis; user=true)

            # Rotates props to be tip to tip
            vlm.rotate(this_prop, (-1)^(!CW_tw) * init_ori_prop)

            # Adds the original polars that don't get copied in deepcopy
            this_prop.airfoils = copy_prop.airfoils
            this_prop._polars = copy_prop._polars
            this_prop._polarroot = copy_prop._polarroot
            this_prop._polartip = copy_prop._polartip

            push!(props_tw, this_prop)
        end
    end

    # Assembles moving sections of the wing
    tandem_wing_moving = vlm.WingSystem()
    vlm.addwing(tandem_wing_moving, "WingR", twing_R)
    vlm.addwing(tandem_wing_moving, "WingL", twing_L)
    if add_rotors
        for (i, prop) in enumerate(props_tw)
            vlm.addwing(tandem_wing_moving, "Prop$i", prop)
        end
    end

    # offset to align with pivot line
    x_off_tw = pivot_tw*b_tw/AR_tw

    # Aligns moving and fixed section with their pivot line
    O_off_tw = [-x_off_tw, 0.0, 0.0]

    vlm.setcoordsystem(twing_md, O_off_tw, Array(1.0I, 3, 3))

    for vlmwing in tandem_wing_moving.wings
        if typeof(vlmwing)==vlm.Rotor
            vlm.setcoordsystem(vlmwing, vlmwing._wingsystem.O + O_off_tw,
                                        vlmwing._wingsystem.Oaxis; user=false)
        else
            vlm.setcoordsystem(vlmwing, vlmwing.O + O_off_tw, vlmwing.Oaxis)
        end
    end

    # Initial rotation of moving sections
    Oaxis_twmv = gt.rotation_matrix(0.0, -init_ori, 0.0)
    vlm.setcoordsystem(tandem_wing_moving, zeros(3), Oaxis_twmv)

    # Assembles tandem wing
    tandem_wing = vlm.WingSystem()
    vlm.addwing(tandem_wing, "FixedWing", twing_md)
    vlm.addwing(tandem_wing, "Moving", tandem_wing_moving)

    # Position of tandem wing
    O_tw = [l_pos_tw + x_off_tw, 0, h_pos_tw]
    Oaxis_tw = gt.rotation_matrix(0.0, 0.0, 0.0)
    vlm.setcoordsystem(tandem_wing, O_tw, Oaxis_tw)

    # ------------ FUSELAGE ----------------------------------------------
    if verbose; println("\t"^(v_lvl+1)*"Generating fuselage..."); end;
    fuselage = generatefuselage_vahana(l_f; ncells=20)

    body = gt.MultiGrid(3)
    gt.addgrid(body, "Fuselage", fuselage)

    if add_rotors
        if verbose; println("\t"^(v_lvl+1)*"Generating pylons..."); end;

        pylon_pos = 0.70
        pylon_length = stckd_xoc_offset*AR_w/b_w * pylon_pos

        for i in stacked
            pylon = generatepylon(pylon_length; ncells=5)

            right = i<=np_w
            this_O = O_prop_w[ right ? i : np_w-(i-np_w-1)]   # Chooses position
            this_O = [1 0 0; 0 (-1)^!right 0; 0 0 1]*this_O   # Places it in correct side
            this_O += main_wing.O                             # Translates it with the wing
                                                              # Offsets it according to length
            this_O += [1/3*0.5*pylon_pos, 0, 3*R_w*stckd_zoR_offset]
            rotation = gt.rotation_matrix2(0, 0, 0)

            gt.lintransform!(pylon, rotation, this_O)

            gt.addgrid(body, "Pylon$i", pylon)
        end
    end

    # ------------ SYSTEM ------------------------------------------------
    # Creates system assembly
    system = vlm.WingSystem()
    vlm.addwing(system, "MainWing", main_wing)
    vlm.addwing(system, "TandemWing", tandem_wing)

    if add_rotors
        rotors = vcat(props_w, props_w_stacked_up, props_w_stacked_low, props_tw)
    end

    # ------------ GROUND SURFACE ----------------------------------------
    ground1 = generate_ground_vahana(l_f)          # Take-off pad
    ground2 = generate_ground_vahana(l_f)          # Landing pad

    grounds = (ground1, ground2)

    # ------------ OUTPUTS------------------------------------------------
    # Tilting systems
    tilting_systems = (main_wing_moving, tandem_wing_moving)

    # Rotors grouped by systems of the same RPM
    if add_rotors
        rotor_systems = (props_w, props_w_stacked_up, props_w_stacked_low, props_tw)
    else
        rotor_systems = ()
    end

    # System to solve through the VLM solver
    vlm_system_m = vlm.WingSystem()
    vlm.addwing(vlm_system_m, "middle", wing_md)
    vlm.addwing(vlm_system_m, "R", wing_R)
    vlm.addwing(vlm_system_m, "letR", winglet_R)
    vlm.addwing(vlm_system_m, "L", wing_L)
    vlm.addwing(vlm_system_m, "letL", winglet_L)
    vlm_system_t = vlm.WingSystem()
    vlm.addwing(vlm_system_t, "R", twing_R)
    vlm.addwing(vlm_system_t, "L", twing_L)

    vlm_system = vlm.WingSystem()
    vlm.addwing(vlm_system, "MWing", vlm_system_m)
    vlm.addwing(vlm_system, "TWing", vlm_system_t)

    # Wake-shedding system (`vlm_system`+`rotors`)
    wake_system = vlm.WingSystem()
    vlm.addwing(wake_system, "SolveVLM", vlm_system)
    if add_rotors
        if VehicleType==uns.VLMVehicle
            for (i, rotor) in enumerate(rotors)
                vlm.addwing(wake_system, "Rotor$i", rotor)
            end
        else
            # Mute colinear warnings. This is needed since the quasi-steady solver
            #   will probe induced velocities at the lifting line of the blade
            uns.vlm.VLMSolver._mute_warning(true)
        end
    end

    # Dummy grids that are rotated and translated along with the vehicle
    grids = [body]

    # FVS's Vehicle object
    vehicle = VehicleType(   system;
                                tilting_systems=tilting_systems,
                                rotor_systems=rotor_systems,
                                vlm_system=vlm_system,
                                # wake_system=vlm_system,
                                wake_system=wake_system,
                                grids=grids
                             )

    return vehicle, grounds
end

"""
    Returns a loft of the Vahana fuselage.
"""
function generatefuselage_vahana(fuselage_length; ncells=20)
    sections = [[] for i in 1:6]

    # --------------- SECTION #1------------------------------------------------
    Ps = [
            [0.5, 0.1+0.2],
            [0.0, 1.0],
            [-0.5, 0.1+0.2],
            [0, -0.20+0.2]
         ]

    CPs = [
            [0.4, 1.0],
            [-0.4, 1.0],
            [-0.1, -0.3+0.2],
            [0.1, -0.3+0.2],
         ]

    rhos = [
            0.6,
            0.6,
            0.2,
            0.2
         ]

    scaling = 1.0

    ss = [
            collect(range(0, 1, length=ncells)),
            collect(range(0, 1, length=ncells)),
            collect(range(0, 1, length=ncells)),
            collect(range(0, 1, length=ncells))
         ]

    points = gt.conic_cross_section(Ps, CPs, rhos, ss)
    points = scaling*points

    sections[1] = points

    # --------------- SECTION #2------------------------------------------------
    Ps = [
            [0.4, 0.4],
            [0.0, 1.0],
            [-0.4, 0.4],
            [0, -0.05]
         ]

    CPs = [
            [0.3, 1.0],
            [-0.3, 1.0],
            [-0.4, -0.08],
            [0.4, -0.08],
         ]

    rhos = [
            0.7,
            0.7,
            0.5,
            0.5
         ]

    scaling = 2.0

    ss = [
            collect(range(0, 1, length=ncells)),
            collect(range(0, 1, length=ncells)),
            collect(range(0, 1, length=ncells)),
            collect(range(0, 1, length=ncells))
         ]

    points = gt.conic_cross_section(Ps, CPs, rhos, ss)
    points = scaling*points

    sections[2] = points;

    # --------------- SECTION #3------------------------------------------------
    Ps = [
            [0.8, 1.0],
            [0.0, 2.5],
            [-0.8, 1.0],
            [0, 0.2]
         ]

    CPs = [
            [0.6, 2.5],
            [-0.6, 2.5],
            [-0.8, 0.2],
            [0.8, 0.2],
         ]

    rhos = [
            0.7,
            0.7,
            0.5,
            0.5
         ]

    scaling = 1.0

    ss = [
            collect(range(0, 1, length=ncells)),
            collect(range(0, 1, length=ncells)),
            collect(range(0, 1, length=ncells)),
            collect(range(0, 1, length=ncells))
         ]

    points = gt.conic_cross_section(Ps, CPs, rhos, ss)
    points = scaling*points

    sections[3] = points;

    # --------------- SECTION #4------------------------------------------------
    Ps = [
            [0.7, 2.0+0.25],
            [0.0, 3.0+0.25],
            [-0.7, 2.0+0.25],
            [0, 1.2]
         ]

    CPs = [
            [0.5, 3.0+0.25],
            [-0.5, 3.0+0.25],
            [-0.8, 1.2],
            [0.8, 1.2],
         ]

    rhos = [
            0.7,
            0.7,
            0.5,
            0.5
         ]

    scaling = 1.0

    ss = [
            collect(range(0, 1, length=ncells)),
            collect(range(0, 1, length=ncells)),
            collect(range(0, 1, length=ncells)),
            collect(range(0, 1, length=ncells))
         ]

    points = gt.conic_cross_section(Ps, CPs, rhos, ss)
    points = scaling*points

    sections[4] = points;

    # --------------- SECTION #5------------------------------------------------
    Ps = [
            [0.4, 2.75-0.2],
            [0.0, 3.25-0.2],
            [-0.4, 2.75-0.2],
            [0, 1.7]
         ]

    CPs = [
            [0.3, 3.25-0.2],
            [-0.3, 3.25-0.2],
            [-0.4, 1.7],
            [0.4, 1.7],
         ]

    rhos = [
            0.4,
            0.4,
            0.5,
            0.5
         ]

    scaling = 1.0

    ss = [
            collect(range(0, 1, length=ncells)),
            collect(range(0, 1, length=ncells)),
            collect(range(0, 1, length=ncells)),
            collect(range(0, 1, length=ncells))
         ]

    points = gt.conic_cross_section(Ps, CPs, rhos, ss)
    points = scaling*points

    sections[5] = points;

    # --------------- SECTION #6------------------------------------------------
    Ps = [
            [0.4-0.2, 3.25],
            [0.0, 3.5-0.1],
            [-0.4+0.2, 3.25],
            [0, 2.5+0.5]
         ]

    CPs = [
            [0.3-0.2, 3.5-0.1],
            [-0.3+0.2, 3.5-0.1],
            [-0.4, 2.5+0.5],
            [0.4, 2.5+0.5],
         ]

    rhos = [
            0.4,
            0.4,
            0.3,
            0.3
         ]

    scaling = 1.0

    ss = [
            collect(range(0, 1, length=ncells)),
            collect(range(0, 1, length=ncells)),
            collect(range(0, 1, length=ncells)),
            collect(range(0, 1, length=ncells))
         ]
    points = gt.conic_cross_section(Ps, CPs, rhos, ss)
    points = scaling*points

    sections[6] = points;

    # --------------- SECTION #0------------------------------------------------
    Ps = [
            [0.5, 0.1+0.2],
            [0.0, 1.0],
            [-0.5, 0.1+0.2],
            [0, -0.20+0.2]
         ]

    CPs = [
            [0.4, 1.0],
            [-0.4, 1.0],
            [-0.1, -0.3+0.2],
            [0.1, -0.3+0.2],
         ]

    rhos = [
            0.6,
            0.6,
            0.2,
            0.2
         ]

    scaling = 0.1

    ss = [
            collect(range(0, 1, length=ncells)),
            collect(range(0, 1, length=ncells)),
            collect(range(0, 1, length=ncells)),
            collect(range(0, 1, length=ncells))
         ]
    points = gt.conic_cross_section(Ps, CPs, rhos, ss)
    points = scaling*points

    section0 = points

    # --------------- LOFT -----------------------------------------------------
    # fuselage_length = 5.86              # (m) fuselage length
    bscale = fuselage_length            # Scales the entire geometry by this factor

    crosssections = [ # (l/bscale position , contour)
                        (0.0  , [[2.00*0.8, 0.25][k]*(point[k]) + [0, 0.3][k] for point in section0, k in 1:2]),
                        (0.02, [[0.6*0.8, 0.3][k]*(point[k]) + [0, 0.15][k] for point in sections[1], k in 1:2]),
                        (1/6*0.5, [[1.0*0.8, 0.75][k]*point[k] for point in sections[1], k in 1:2]),
    #                     (1/6*0.5, [[0.75, 0.75][k]*point[k] for point in sections[1], k in 1:2]),
                        (2/6, [[1.00*0.8, 0.75][k]*point[k] for point in sections[2], k in 1:2]),
                        (3/6, [[1.00*0.8, 0.75][k]*(point[k] - [0, 0.1*2][k]) for point in sections[3], k in 1:2]),
                        (4/6, [[0.75, 0.75][k]*(point[k] - [0, 0.4*2][k]) for point in sections[4], k in 1:2]),
                        (5/6, [[0.75, 0.75][k]*(point[k] - [0, 0.6][k]) for point in sections[5], k in 1:2]),
                        (6/6, [[0.75, 0.75][k]*(point[k] - [0, 1.0][k]) for point in sections[6], k in 1:2]),
                    ]

    # Dummy parameters
    b_pos = [sec[1] for sec in crosssections]
    chords = (1/bscale)*ones(size(b_pos))
    twists = zeros(size(b_pos))
    LE_x = zeros(size(b_pos))
    LE_z = zeros(size(b_pos))

    tilt_y = zeros(size(b_pos))
    # tilt_y[2] = 45

    grid_fuselage = gt.generate_loft(crosssections, bscale,
                                        b_pos, chords, twists, LE_x, LE_z;
                                        # MORE GEOMETRIC OPTIONS
                                        tilt_y=tilt_y,
                                        symmetric=false,
                                        loop_dim=0,
                                    )

    # Aligns the grid with the x-axis
    rotation = gt.rotation_matrix2(0, 0, 90)
    gt.lintransform!(grid_fuselage, rotation, zeros(3))

    dimsplit = 1
    triangrid_fuselage = gt.GridTriangleSurface(grid_fuselage, dimsplit)

    return triangrid_fuselage
end

"""
    Returns a loft of the pylon.
"""
function generatepylon(pylon_length; ncells=20)
    sections = [[] for i in 1:3]

    # --------------- SECTION #1------------------------------------------------
    Ps = [
            [0.5, 0.1+0.2],
            [0.0, 1.0],
            [-0.5, 0.1+0.2],
            [0, -0.20+0.2]
         ]

    CPs = [
            [0.4, 1.0],
            [-0.4, 1.0],
            [-0.1, -0.3+0.2],
            [0.1, -0.3+0.2],
         ]

    rhos = [
            0.6,
            0.6,
            0.2,
            0.2
         ]

    scaling = 1.0

    ss = [
            collect(range(0, 1, length=ncells)),
            collect(range(0, 1, length=ncells)),
            collect(range(0, 1, length=ncells)),
            collect(range(0, 1, length=ncells))
         ]

    points = gt.conic_cross_section(Ps, CPs, rhos, ss)
    points = scaling*points

    sections[1] = points

    # --------------- SECTION #2------------------------------------------------
    Ps = [
            [0.4, 0.4],
            [0.0, 1.0],
            [-0.4, 0.4],
            [0, -0.05]
         ]

    CPs = [
            [0.3, 1.0],
            [-0.3, 1.0],
            [-0.4, -0.08],
            [0.4, -0.08],
         ]

    rhos = [
            0.7,
            0.7,
            0.5,
            0.5
         ]

    scaling = 2.0

    ss = [
            collect(range(0, 1, length=ncells)),
            collect(range(0, 1, length=ncells)),
            collect(range(0, 1, length=ncells)),
            collect(range(0, 1, length=ncells))
         ]

    points = gt.conic_cross_section(Ps, CPs, rhos, ss)
    points = scaling*points

    sections[2] = points;

    # --------------- SECTION #3------------------------------------------------
    Ps = [
            [0.4, 2.75-0.2],
            [0.0, 3.25-0.2],
            [-0.4, 2.75-0.2],
            [0, 1.7]
         ]

    CPs = [
            [0.3, 3.25-0.2],
            [-0.3, 3.25-0.2],
            [-0.4, 1.7],
            [0.4, 1.7],
         ]

    rhos = [
            0.4,
            0.4,
            0.5,
            0.5
         ]

    scaling = 1.0

    ss = [
            collect(range(0, 1, length=ncells)),
            collect(range(0, 1, length=ncells)),
            collect(range(0, 1, length=ncells)),
            collect(range(0, 1, length=ncells))
         ]

    points = gt.conic_cross_section(Ps, CPs, rhos, ss)
    points = scaling*points

    points = [p .- [0, 1.7] for p in points]
    # points = [p .* [1, 0.75] for p in points]

    sections[3] = points;

    # --------------- SECTION #0------------------------------------------------
    Ps = [
            [0.5, 0.1+0.2],
            [0.0, 1.0],
            [-0.5, 0.1+0.2],
            [0, -0.20+0.2]
         ]

    CPs = [
            [0.4, 1.0],
            [-0.4, 1.0],
            [-0.1, -0.3+0.2],
            [0.1, -0.3+0.2],
         ]

    rhos = [
            0.6,
            0.6,
            0.2,
            0.2
         ]

    scaling = 0.1

    ss = [
            collect(range(0, 1, length=ncells)),
            collect(range(0, 1, length=ncells)),
            collect(range(0, 1, length=ncells)),
            collect(range(0, 1, length=ncells))
         ]
    points = gt.conic_cross_section(Ps, CPs, rhos, ss)
    points = scaling*points

    section0 = points

    # --------------- LOFT -----------------------------------------------------

    # --------------- LOFT -----------------------------------------------------
    # fuselage_length = 5.86              # (m) fuselage length
    bscale = pylon_length            # Scales the entire geometry by this factor

    crosssections = [ # (l/bscale position , contour)
                    (0.0  , [[2.00*0.8, 0.25][k]*(point[k]) + [0, 0.3][k] for point in section0, k in 1:2]*pylon_length/5.86),
                    (0.02, [[0.6*0.8, 0.3][k]*(point[k]) + [0, 0.15][k] for point in sections[1], k in 1:2]*pylon_length/5.86),
                    (1/3*0.25, [[1.0*0.8, 0.75*0.75][k]*point[k] for point in sections[1], k in 1:2]*pylon_length/5.86),
                    (2/3, [[1.00*0.8, 0.35*0.75][k]*point[k] for point in sections[2], k in 1:2]*pylon_length/5.86),
                    (3/3, [[1.00*0.8, 0.75*0.75][k]*(point[k] - [0, 0.1*2][k]) for point in sections[3], k in 1:2]*pylon_length/5.86),
                    ]

    # Dummy parameters
    b_pos = [sec[1] for sec in crosssections]
    chords = (1/bscale)*ones(size(b_pos))
    twists = zeros(size(b_pos))
    LE_x = zeros(size(b_pos))
    LE_z = zeros(size(b_pos))

    tilt_y = zeros(size(b_pos))
    # tilt_y[2] = 45

    grid_pylon = gt.generate_loft(crosssections, bscale,
                                        b_pos, chords, twists, LE_x, LE_z;
                                        # MORE GEOMETRIC OPTIONS
                                        tilt_y=tilt_y,
                                        symmetric=false,
                                        loop_dim=0,
                                    )

    # Aligns the grid with the x-axis
    rotation = gt.rotation_matrix2(0, 0, 90)
    gt.lintransform!(grid_pylon, rotation, zeros(3))

    dimsplit = 1
    triangrid_pylon = gt.GridTriangleSurface(grid_pylon, dimsplit)

    return triangrid_pylon
end

"""
    Generates ground surface.
"""
function generate_ground_vahana(bscale; verify_spline=false, spl_s=0.001,
                                spl_k="automatic", visualize=false)

    ## Ellipse
    fx, fy = 1.7, 1.05
    lx = bscale*3.5
    ly = bscale*3
    cx = -lx*3/4
    cy = 0.0
    cz = -0.06*bscale
    a, b = lx, ly
    perimeter = [ a*b/sqrt((b*cos(tht))^2 + (a*sin(tht))^2
                    ) * [fx*cos(tht), fy*sin(tht), 0] + [cx, cy, cz] for
                                            tht in range(0, 2*pi, length=179)]

    # Fluid domain
    NDIVSx = 50              # Cells in the parametric x-direction
    NDIVSy = 50              # Cells in the parametric y-direction
    NDIVSz = 0               # Cells in the geometric z-direction

    z_off = cz               # z-offset for ground position

    if visualize
        opt_args = [(:save_path, "temps/"), (:file_name, "vahanaground"),
                    (:paraview, true)]
    else
        opt_args = []
    end

    ground = generate_perimetergrid(perimeter, NDIVSx, NDIVSy, NDIVSz;
                                       z_min=z_off, z_max=z_off,
                                        verify_spline=verify_spline, spl_s=spl_s,
                                        spl_k=spl_k,
                                        opt_args...
                                        )

    return ground
end

# -------------- END OF VAHANA -------------------------------------------------

################################################################################
# UTILITIES
################################################################################

"""
`generate_perimetergrid(perimeter::Array{Array{T, 1}, 1},
                                  NDIVSx, NDIVSy, NDIVSz;
                                  z_min::Real=0, z_max::Real=0,
                                  # SPLINE OPTIONS
                                  verify_spline::Bool=true,
                                  spl_s=0.001, spl_k="automatic",
                                  # FILE OPTIONS
                                  save_path=nothing, file_name="perimeter",
                                  paraview=true
                                )`

  Generates the perimeter grid with `perimeter` the array of points of the
  contour (must be a closed contour), and `NDIVS_` the number of cells in
  each parametric dimension (give it `NDIVSz=0` for a flat surface, otherwise
  it'll generate a volumetric grid between `z_min` and `z_max`).
"""
function generate_perimetergrid(perimeter::Array{Array{T, 1}, 1},
                                  NDIVSx, NDIVSy, NDIVSz;
                                  z_min::Real=0, z_max::Real=0,
                                  # SPLINE OPTIONS
                                  verify_spline::Bool=true,
                                  spl_s=0.001, spl_k="automatic",
                                  # FILE OPTIONS
                                  save_path=nothing, file_name="perimeter",
                                  paraview=true
                                ) where{T<:Real}

  # Error cases
  multidiscrtype = Array{Tuple{Float64,Int64,Float64,Bool},1}
  if typeof(NDIVSx)==Int64
    nz = NDIVSz
  elseif typeof(NDIVSz)==multidiscrtype
    nz = 0
    for sec in NDIVSz
      nz += sec[2]
    end
  else
    error("Expected `NDIVSz` to be type $(Int64) or $MultiDiscrType,"*
            " got $(typeof(NDIVSz)).")
  end

  # --------- REPARAMETERIZES THE PERIMETER ---------------------------
  org_x = [p[1] for p in perimeter]
  org_y = [p[2] for p in perimeter]
  # Separate upper and lower sides to make the contour injective in x
  upper, lower = gt.splitcontour(org_x, org_y)

  # # Parameterize both sides independently
  # fun_upper = gt.parameterize(upper[1], upper[2], zeros(upper[1]); inj_var=1,
  #                                                     s=spl_s, kspl=spl_k)
  # fun_lower = gt.parameterize(lower[1], lower[2], zeros(lower[1]); inj_var=1,
  #                                                     s=spl_s, kspl=spl_k)
  # # Discretizes both sides
  # if NDIVSx==multidiscrtype
  #   new_upper = gt.multidiscretize(fun_upper, 0, 1, NDIVSx)
  #   new_lower = gt.multidiscretize(fun_lower, 0, 1, NDIVSx)
  # else
  #   new_upper = gt.discretize(fun_upper, 0, 1, NDIVSx, 1.0)
  #   new_lower = gt.discretize(fun_lower, 0, 1, NDIVSx, 1.0)
  # end

  splt_up = Int(ceil((size(upper[1],1)/2)))
  splt_low = Int(ceil((size(lower[1],1)/2)))

  # Splits the perimeter into four faces
  upper1 = [[x for x in upper[1][1:splt_up]], [y for y in upper[2][1:splt_up]]]
  upper2 = [[x for x in upper[1][splt_up:end]], [y for y in upper[2][splt_up:end]]]
  lower1 = [[x for x in lower[1][1:splt_low]], [y for y in lower[2][1:splt_low]]]
  lower2 = [[x for x in lower[1][splt_low:end]], [y for y in lower[2][splt_low:end]]]

  # Parameterize both sides independently
  fun_upper1 = gt.parameterize(upper1[1], upper1[2], zeros(size(upper1[1])); inj_var=1,
                                                      s=spl_s, kspl=spl_k)
  fun_upper2 = gt.parameterize(upper2[1], upper2[2], zeros(size(upper2[1])); inj_var=1,
                                                      s=spl_s, kspl=spl_k)
  fun_lower1 = gt.parameterize(lower1[1], lower1[2], zeros(size(lower1[1])); inj_var=1,
                                                      s=spl_s, kspl=spl_k)
  fun_lower2 = gt.parameterize(lower2[1], lower2[2], zeros(size(lower2[1])); inj_var=1,
                                                      s=spl_s, kspl=spl_k)
  # Discretizes both sides
  if NDIVSx==multidiscrtype
    new_upper1 = gt.multidiscretize(fun_upper1, 0, 1, NDIVSx)
    new_upper2 = gt.multidiscretize(fun_upper2, 0, 1, NDIVSy)
    new_lower1 = gt.multidiscretize(fun_lower1, 0, 1, NDIVSy)
    new_lower2 = gt.multidiscretize(fun_lower2, 0, 1, NDIVSx)
  else
    new_upper1 = gt.discretize(fun_upper1, 0, 1, NDIVSx, 1.0)
    new_upper2 = gt.discretize(fun_upper2, 0, 1, NDIVSy, 1.0)
    new_lower1 = gt.discretize(fun_lower1, 0, 1, NDIVSy, 1.0)
    new_lower2 = gt.discretize(fun_lower2, 0, 1, NDIVSx, 1.0)
  end

  # ----------------- SPLINE VERIFICATION --------------------------------------

  if verify_spline
  #   new_points = vcat(reverse(new_upper), new_lower)
    # new_x = [p[1] for p in new_points]
    # new_y = [p[2] for p in new_points]
    plt.plot(org_x, org_y, "--^k", label="Original", alpha=0.5)
    # plt.plot(new_x, new_y, ":.b", label="Parameterized")

    plt.plot(upper1[1], upper1[2], "*b", label="upper1", alpha=0.75)
    plt.plot(upper2[1], upper2[2], "^b", label="upper2", alpha=0.75)
    plt.plot(lower1[1], lower1[2], "*r", label="lower1", alpha=0.75)
    plt.plot(lower2[1], lower2[2], "^r", label="lower2", alpha=0.75)

    new_points = vcat(new_upper1, new_upper2)
    new_x = [p[1] for p in new_points]
    new_y = [p[2] for p in new_points]
    plt.plot(new_x, new_y, ":.b", label="Parameterized Upper", alpha=0.75)
    new_points = vcat(new_lower1, new_lower2)
    new_x = [p[1] for p in new_points]
    new_y = [p[2] for p in new_points]
    plt.plot(new_x, new_y, ":.r", label="Parameterized Lower", alpha=0.75)
    # plt.xlabel(plt.L"x")
    # plt.ylabel(plt.L"y")
    plt.legend(loc="best")
    plt.grid(true, color="0.8", linestyle="--")
  end

  # --------- GRIDS THE INSIDE OF THE PERIMETER ---------------------
  # Parametric grid
  P_min = zeros(3)
  P_max = [1, 1, 1*(nz!=0)]
  param_grid = gt.Grid(P_min, P_max, [NDIVSx, NDIVSy, NDIVSz])

  # function my_space_transform(X, ind)
  #     i = ind[1]                      # Arc length point
  #     w = X[2]                        # Weight
  #     z = z_min + X[3]*(z_max-z_min)  # z-position
  #
  #     Y = new_lower[i] + w*(new_upper[i]-new_lower[i])
  #     Y[3] = z
  #
  #     return Y
  # end

  rev_new_upper1 = reverse(new_upper1) # Left edge
  rev_new_upper2 = reverse(new_upper2) # Upper edge
  rev_new_lower1 = reverse(new_lower1) # Lower edge
  rev_new_lower2 = reverse(new_lower2) # Right edge

  function my_space_transform(X, ind)
      beta = 1e-2*pi/2                 # Warping parameter

      z = z_min + X[3]*(z_max-z_min)   # z-position

      xw = X[2]                        # x weight
      yw = X[1]                        # y weight

      # Linearly interpolate the edges
      x = new_upper1[ind[1]]*xw + new_lower2[ind[1]]*(1-xw) # Left-right edges
      y = rev_new_upper2[ind[2]]*yw + rev_new_lower1[ind[2]]*(1-yw) # Upper-lower edges
      # Y = x + y

      # A clever way of merging the two interpolations
      # First weight
      yw = abs(X[1]-0.5)/0.5
      yw = tan(yw*pi/2 - beta)/tan(pi/2 - beta)
      # Second weight
      xw = abs(X[2]-0.5)/0.5
      xw = tan(xw*pi/2 - beta)/tan(pi/2 - beta)
      # Merged weight
      w = ( (1-xw) + yw ) /2
      # Merged point
      Y = x*(1-w) + y*w

      Y[3] = z

      return Y
  end

  # Applies the space transformation to the parametric grid
  gt.transform3!(param_grid, my_space_transform)

  if save_path!=nothing
    gt.save(param_grid, file_name; path=save_path)

    if paraview
      strn = file_name*".vtk"
      run(`paraview --data=$save_path/$strn`)
    end

  end

  return param_grid::gt.Grid
end
# -------------- END OF COMMONS ------------------------------------------------


text:
module TestTikZ

using Test
using Catlab.Graphics.TikZ

# Pretty-print
##############

spprint(expr::Expression) = sprint(pprint, expr)

# Node statement
@test spprint(Node("f")) == "\\node (f) {};"
@test spprint(Node("f"; props=[Property("circle")])) ==
  "\\node[circle] (f) {};"
@test spprint(Node("f"; props=[Property("circle"),Property("fill","red")])) ==
  "\\node[circle,fill=red] (f) {};"
@test spprint(Node("f"; coord=Coordinate("0","1"))) ==
  "\\node (f) at (0,1) {};"
@test spprint(Node("f"; coord=Coordinate("1","0"), content="fun")) ==
  "\\node (f) at (1,0) {fun};"

# Edge statement
@test spprint(Edge("f","g")) == "\\draw (f) to (g);"
@test spprint(Edge("f","g"; props=[Property("red")])) == "\\draw[red] (f) to (g);"
@test spprint(Edge("f","g"; node=EdgeNode())) == "\\draw (f) to node (g);"
@test spprint(Edge("f","g"; node=EdgeNode(;props=[Property("circle")],content="e"))) ==
  "\\draw (f) to node[circle] {e} (g);"
@test spprint(Edge("f","g"; op=PathOperation("to"; props=[Property("out","0")]))) ==
  "\\draw (f) to[out=0] (g);"

# Picture statement
@test spprint(
Picture(
  Node("f"),
  Node("g"),
  Edge("f","g")
)) == """
\\begin{tikzpicture}
  \\node (f) {};
  \\node (g) {};
  \\draw (f) to (g);
\\end{tikzpicture}"""

@test spprint(Picture(Node("f");props=[Property("node distance","1cm")])) == """
\\begin{tikzpicture}[node distance=1cm]
  \\node (f) {};
\\end{tikzpicture}"""

@test spprint( # Nested pictures
Picture(
  Node("outer1"; content=
    Picture(
      Node("inner1"),
      Node("inner2")
    )),
  Node("outer2"; content=
    Picture(
      Node("inner3"),
      Node("inner4")
    ));
  props=[Property("remember picture")]
)) == """
\\begin{tikzpicture}[remember picture]
  \\node (outer1) {
    \\begin{tikzpicture}
      \\node (inner1) {};
      \\node (inner2) {};
    \\end{tikzpicture}
  };
  \\node (outer2) {
    \\begin{tikzpicture}
      \\node (inner3) {};
      \\node (inner4) {};
    \\end{tikzpicture}
  };
\\end{tikzpicture}"""

# Scope statement
@test spprint(
Picture(
  Scope(Node("f"); props=[Property("red")]),
  Scope(Node("g"); props=[Property("blue")]),
  Edge("f","g")
)) == """
\\begin{tikzpicture}
  \\begin{scope}[red]
    \\node (f) {};
  \\end{scope}
  \\begin{scope}[blue]
    \\node (g) {};
  \\end{scope}
  \\draw (f) to (g);
\\end{tikzpicture}"""

# Graph, GraphScope, GraphNode, GraphEdge statements
@test spprint(
Graph(
  GraphNode("f";content="fun"),
  GraphEdge("f","g")
)) == """
\\graph{
  f/\"fun\";
  f -> g;
};"""

@test spprint(
Graph(
  GraphNode("f";props=[Property("circle")]),
  GraphEdge("f","g";props=[Property("edge label","X")])
)) == """
\\graph{
  f [circle];
  f ->[edge label=X] g;
};"""

@test spprint(
Graph(
  GraphScope(
    GraphEdge("a","b"); props=[Property("same layer")]
  ),
  GraphScope(
    GraphEdge("c","d"); props=[Property("same layer")]
  )
)) == """
\\graph{
  {[same layer]
    a -> b;
  };
  {[same layer]
    c -> d;
  };
};"""

# Matrix statement
@test spprint(
MatrixNode(
  [ [[Node("f")]]              [[Edge("g1","g2")]];
    [[Node("h1"),Node("h2")]]  [[Node("i1"),Node("i2")]] ];
  props=[Property("draw","red")]
)) == """
\\matrix[draw=red]{
  \\node (f) {}; &
  \\draw (g1) to (g2); \\\\
  \\node (h1) {};
  \\node (h2) {}; &
  \\node (i1) {};
  \\node (i2) {}; \\\\
};"""

end

ext: jl


text:
#=
Code for simulating, visualizing and manipulating additive and
multiplicative functionals.

@authors: Shunsuke Hori
=#
using QuantEcon
using PyPlot
using Distributions

"""
This type transforms an additive (multipilcative)
functional into a QuantEcon linear state space system.
"""
struct AMF_LSS_VAR{TF<:AbstractFloat, TI<:Integer}
    A::Array{TF, 2}
    B::Array{TF, 2}
    D::Array{TF, 2}
    F::Array{TF, 2}
    nu::Array{TF, 2}
    nx::TI
    nk::TI
    nm::TI
    lss::LSS
end

function AMF_LSS_VAR(A::Array, B::Array,
                D::Union{RowVector, Array},
                F::Union{Void, Array}=nothing;
                nu::Union{Void, Array}=nothing)

    if typeof(B) <: Vector
        B = reshape(B, length(B), 1)
    end
    # Unpack required elements
    nx, nk = size(B)

    # checking the dimension of D (extended from the scalar case)
    if ndims(D) > 1
        nm = size(D, 1)
        if typeof(D) <: RowVector
            D = convert(Matrix, D)
        end
    else
        nm = 1
        D = reshape(D, 1, length(D))
    end

    # Set F
    if F==nothing
        F = zeros(nk, 1)
    elseif ndims(F) == 1
        F = reshape(F, length(F), 1)
    end

    # Set nu
    if nu==nothing
        nu = zeros(nm, 1)
    elseif ndims(nu) == 1
        nu = reshape(nu, length(nu), 1)
    else
        throw(ArgumentError("nu must be column vector!"))
    end

    if size(nu, 1) != size(D, 1)
        error("The size of nu is inconsistent with D!")
    end

    # Construct BIG state space representation
    lss = construct_ss(A, B, D, F, nu, nx, nk, nm)

    return AMF_LSS_VAR(A, B, D, F, nu, nx, nk, nm, lss)
end

AMF_LSS_VAR(A::Array, B::Array, D::Union{RowVector, Array}) =
    AMF_LSS_VAR(A, B, D, nothing, nu=nothing)
AMF_LSS_VAR(A::Array, B::Array, D::Union{RowVector, Array}, F::Real, nu::Real) =
    AMF_LSS_VAR(A, B, D, [F], nu=[nu])

"""
This creates the state space representation that can be passed
into the quantecon LSS class.
"""
function construct_ss(A::Array, B::Array, D::Union{RowVector, Array}, F::Array,
                      nu, nx::TI, nk::TI, nm::TI) where TI <: Integer

    H, g = additive_decomp(A, B, D, F, nx)

    # Auxiliary blocks with 0's and 1's to fill out the lss matrices
    nx0c = zeros(nx, 1)
    nx0r = zeros(1, nx)
    nx1 = ones(1, nx)
    nk0 = zeros(1, nk)
    ny0c = zeros(nm, 1)
    ny0r = zeros(1, nm)
    ny1m = eye(nm)
    ny0m = zeros(nm, nm)
    nyx0m = zeros(D)

    # Build A matrix for LSS
    # Order of states is: [1, t, xt, yt, mt]
    A1 = hcat(1, 0, nx0r, ny0r, ny0r)          # Transition for 1
    A2 = hcat(1, 1, nx0r, ny0r, ny0r)          # Transition for t
    A3 = hcat(nx0c, nx0c, A, nyx0m', nyx0m')   # Transition for x_{t+1}
    A4 = hcat(nu, ny0c, D, ny1m, ny0m)         # Transition for y_{t+1}
    A5 = hcat(ny0c, ny0c, nyx0m, ny0m, ny1m)   # Transition for m_{t+1}
    Abar = vcat(A1, A2, A3, A4, A5)

    # Build B matrix for LSS
    Bbar = vcat(nk0, nk0, B, F, H)

    # Build G matrix for LSS
    # Order of observation is: [xt, yt, mt, st, tt]
    G1 = hcat(nx0c, nx0c, eye(nx), nyx0m', nyx0m')    # Selector for x_{t}
    G2 = hcat(ny0c, ny0c, nyx0m, ny1m, ny0m)          # Selector for y_{t}
    G3 = hcat(ny0c, ny0c, nyx0m, ny0m, ny1m)          # Selector for martingale
    G4 = hcat(ny0c, ny0c, -g, ny0m, ny0m)             # Selector for stationary
    G5 = hcat(ny0c, nu, nyx0m, ny0m, ny0m)            # Selector for trend
    Gbar = vcat(G1, G2, G3, G4, G5)

    # Build LSS type
    x0 = hcat(1, 0, nx0r, ny0r, ny0r)
    S0 = zeros(length(x0), length(x0))
    lss = LSS(Abar, Bbar, Gbar, x0, S0)

    return lss
end

"""
Return values for the martingale decomposition
    - nu        : unconditional mean difference in Y
    - H         : coefficient for the (linear) martingale component (kappa_a)
    - g         : coefficient for the stationary component g(x)
    - Y_0       : it should be the function of X_0 (for now set it to 0.0)
"""
function additive_decomp(A::Array, B::Array, D::Array, F::Union{Array, Real},
                          nx::Integer)
    I = eye(nx)
    A_res = \(I-A, I)
    g = D * A_res
    H = F + D * A_res * B

    return H, g
end

"""
Return values for the multiplicative decomposition (Example 5.4.4.)
    - nu_tilde  : eigenvalue
    - H         : vector for the Jensen term
"""
function multiplicative_decomp(A::Array, B::Array, D::Array, F::Union{Array, Real},
                                nu::Union{Array, Real}, nx::Integer)
    H, g = additive_decomp(A, B, D, F, nx)
    nu_tilde = nu + 0.5*diag(H*H')

    return H, g, nu_tilde
end

function loglikelihood_path(amf::AMF_LSS_VAR, x::Array, y::Array)
    A, B, D, F = amf.A, amf.B, amf.D, amf.F
    k, T = size(y)
    FF = F*F'
    FFinv = inv(FF)
    temp = y[:, 2:end]-y[:, 1:end-1] - D*x[:, 1:end-1]
    obs =  temp .* FFinv .* temp
    obssum = cumsum(obs)
    scalar = (log(det(FF)) + k*log(2*pi))*collect(1:T)

    return -0.5*(obssum + scalar)
end

function loglikelihood(amf::AMF_LSS_VAR, x::Array, y::Array)
    llh = loglikelihood_path(amf, x, y)

    return llh[end]
end

"""
Plots for the additive decomposition
"""
function plot_additive(amf::AMF_LSS_VAR, T::Integer;
                       npaths::Integer=25, show_trend::Bool=true)

    # Pull out right sizes so we know how to increment
    nx, nk, nm = amf.nx, amf.nk, amf.nm

    # Allocate space (nm is the number of additive functionals - we want npaths for each)
    mpath = Array{Real}(nm*npaths, T)
    mbounds = Array{Real}(nm*2, T)
    spath = Array{Real}(nm*npaths, T)
    sbounds = Array{Real}(nm*2, T)
    tpath = Array{Real}(nm*npaths, T)
    ypath = Array{Real}(nm*npaths, T)

    # Simulate for as long as we wanted
    moment_generator = moment_sequence(amf.lss)
    state = start(moment_generator)
    # Pull out population moments
    for t in 1:T
        tmoms, state = next(moment_generator, state)
        ymeans = tmoms[2]
        yvar = tmoms[4]

        # Lower and upper bounds - for each additive functional
        for ii in 1:nm
            li, ui = (ii-1)*2+1, ii*2
            if sqrt(yvar[nx+nm+ii, nx+nm+ii]) != 0.0
                madd_dist = Normal(ymeans[nx+nm+ii], sqrt(yvar[nx+nm+ii, nx+nm+ii]))
                mbounds[li, t] = quantile(madd_dist, 0.01)
                mbounds[ui, t] = quantile(madd_dist, 0.99)
            elseif sqrt(yvar[nx+nm+ii, nx+nm+ii]) == 0.0
                mbounds[li, t] = ymeans[nx+nm+ii]
                mbounds[ui, t] = ymeans[nx+nm+ii]
            else
                error("standard error is negative")
            end

            if sqrt(yvar[nx+2*nm+ii, nx+2*nm+ii]) != 0.0
                sadd_dist = Normal(ymeans[nx+2*nm+ii], sqrt(yvar[nx+2*nm+ii, nx+2*nm+ii]))
                sbounds[li, t] = quantile(sadd_dist, 0.01)
                sbounds[ui, t] = quantile(sadd_dist, 0.99)
            elseif sqrt(yvar[nx+2*nm+ii, nx+2*nm+ii]) == 0.0
                sbounds[li, t] = ymeans[nx+2*nm+ii]
                sbounds[ui, t] = ymeans[nx+2*nm+ii]
            else
                error("standard error is negative")
            end

        end
    end

    # Pull out paths
    for n in 1:npaths
        x, y = simulate(amf.lss,T)
        for ii in 0:nm-1
            ypath[npaths*ii+n, :] = y[nx+ii+1, :]
            mpath[npaths*ii+n, :] = y[nx+nm + ii+1, :]
            spath[npaths*ii+n, :] = y[nx+2*nm + ii+1, :]
            tpath[npaths*ii+n, :] = y[nx+3*nm + ii+1, :]
        end
    end

    add_figs = Array{Any}(nm)

    for ii in 0:nm-1
        li, ui = npaths*(ii), npaths*(ii+1)
        LI, UI = 2*(ii), 2*(ii+1)
        add_figs[ii+1] =
            plot_given_paths(T, ypath[li+1:ui, :], mpath[li+1:ui, :], spath[li+1:ui, :],
                             tpath[li+1:ui, :], mbounds[LI+1:UI, :], sbounds[LI+1:UI, :],
                             show_trend=show_trend)

        add_figs[ii+1][:suptitle]( L"Additive decomposition of $y_{$(ii+1)}$", fontsize=14 )
    end

    return add_figs
end

"""
Plots for the multiplicative decomposition
"""
function plot_multiplicative(amf::AMF_LSS_VAR, T::Integer,
                            npaths::Integer=25, show_trend::Bool=true)
    # Pull out right sizes so we know how to increment
    nx, nk, nm = amf.nx, amf.nk, amf.nm
    # Matrices for the multiplicative decomposition
    H, g, nu_tilde = multiplicative_decomp(A, B, D, F, nu, nx)

    # Allocate space (nm is the number of functionals - we want npaths for each)
    mpath_mult = Array{Real}(nm*npaths, T)
    mbounds_mult = Array{Real}(nm*2, T)
    spath_mult = Array{Real}(nm*npaths, T)
    sbounds_mult = Array{Real}(nm*2, T)
    tpath_mult = Array{Real}(nm*npaths, T)
    ypath_mult = Array{Real}(nm*npaths, T)

    # Simulate for as long as we wanted
    moment_generator = moment_sequence(amf.lss)
    state = start(moment_generator)
    # Pull out population moments
    for t in 1:T
        tmoms, state = next(moment_generator, state)
        ymeans = tmoms[2]
        yvar = tmoms[4]

        # Lower and upper bounds - for each multiplicative functional
        for ii in 1:nm
            li, ui = (ii-1)*2+1, ii*2
            if yvar[nx+nm+ii, nx+nm+ii] != 0.0
                Mdist = LogNormal(ymeans[nx+nm+ii]- t*0.5*diag(H * H')[ii],
                                sqrt(yvar[nx+nm+ii, nx+nm+ii]))
                mbounds_mult[li, t] = quantile(Mdist, 0.01)
                mbounds_mult[ui, t] = quantile(Mdist, 0.99)
            elseif yvar[nx+nm+ii, nx+nm+ii] == 0.0
                mbounds_mult[li, t] = exp.(ymeans[nx+nm+ii]- t*0.5*diag(H * H')[ii])
                mbounds_mult[ui, t] = exp.(ymeans[nx+nm+ii]- t*0.5*diag(H * H')[ii])
            else
                error("standard error is negative")
            end
            if yvar[nx+2*nm+ii, nx+2*nm+ii] != 0.0
                Sdist = LogNormal(-ymeans[nx+2*nm+ii],
                                sqrt(yvar[nx+2*nm+ii, nx+2*nm+ii]))
                sbounds_mult[li, t] = quantile(Sdist, 0.01)
                sbounds_mult[ui, t] = quantile(Sdist, 0.99)
            elseif yvar[nx+2*nm+ii, nx+2*nm+ii] == 0.0
                sbounds_mult[li, t] = exp.(-ymeans[nx+2*nm+ii])
                sbounds_mult[ui, t] = exp.(-ymeans[nx+2*nm+ii])
            else
                error("standard error is negative")
            end
        end
    end

    # Pull out paths
    for n in 1:npaths
        x, y = simulate(amf.lss,T)
        for ii in 0:nm-1
            ypath_mult[npaths*ii+n, :] = exp.(y[nx+ii+1, :])
            mpath_mult[npaths*ii+n, :] =
                exp.(y[nx+nm + ii+1, :] - collect(1:T)*0.5*diag(H * H')[ii+1])
            spath_mult[npaths*ii+n, :] = 1./exp.(-y[nx+2*nm + ii+1, :])
            tpath_mult[npaths*ii+n, :] =
                exp.(y[nx+3*nm + ii+1, :] + collect(1:T)*0.5*diag(H * H')[ii+1])
        end
    end

    mult_figs = Array{Any}(nm)

    for ii in 0:nm-1
        li, ui = npaths*(ii), npaths*(ii+1)
        LI, UI = 2*(ii), 2*(ii+1)
        mult_figs[ii+1] =
            plot_given_paths(T, ypath_mult[li+1:ui, :], mpath_mult[li+1:ui, :],
                             spath_mult[li+1:ui, :], tpath_mult[li+1:ui, :],
                             mbounds_mult[LI+1:UI, :], sbounds_mult[LI+1:UI, :],
                             horline = 1.0, show_trend=show_trend)
        mult_figs[ii+1][:suptitle]( L"Multiplicative decomposition of $y_{$(ii+1)}$",
                                    fontsize=14)
    end

    return mult_figs
end

function plot_martingales(amf::AMF_LSS_VAR, T::Integer, npaths::Integer=25)

    # Pull out right sizes so we know how to increment
    nx, nk, nm = amf.nx, amf.nk, amf.nm
    # Matrices for the multiplicative decomposition
    H, g, nu_tilde = multiplicative_decomp(amf.A, amf.B, amf.D, amf.F, amf.nu, amf.nx)

    # Allocate space (nm is the number of functionals - we want npaths for each)
    mpath_mult = Array{Real}(nm*npaths, T)
    mbounds_mult = Array{Real}(nm*2, T)

    # Simulate for as long as we wanted
    moment_generator = moment_sequence(amf.lss)
    state = start(moment_generator)
    # Pull out population moments
    for t in 1:T
        tmoms, state = next(moment_generator, state)
        ymeans = tmoms[2]
        yvar = tmoms[4]

        # Lower and upper bounds - for each functional
        for ii in 1:nm
            li, ui = (ii-1)*2+1, ii*2
            if yvar[nx+nm+ii, nx+nm+ii] != 0.0
                Mdist = LogNormal(ymeans[nx+nm+ii]-t*(.5)*diag(H*H')[ii],
                            sqrt(yvar[nx+nm+ii, nx+nm+ii]))
                mbounds_mult[li, t] = quantile(Mdist, 0.01)
                mbounds_mult[ui, t] = quantile(Mdist, 0.99)
            elseif yvar[nx+nm+ii, nx+nm+ii] == 0.0
                mbounds_mult[li, t] = ymeans[nx+nm+ii]-t*(.5)*diag(H*H')[ii]
                mbounds_mult[ui, t] = ymeans[nx+nm+ii]-t*(.5)*diag(H*H')[ii]
            else
                error("standard error is negative")
            end
        end
    end

    # Pull out paths
    for n in 1:npaths
        x, y = simulate(amf.lss, T)
        for ii in 0:nm-1
            mpath_mult[npaths*ii+n, :] =
                exp.(y[nx+nm + ii+1, :] - (1:T)*0.5*diag(H*H')[ii+1])
        end
    end

    mart_figs = Array{Any}(nm)

    for ii in 0:nm-1
        li, ui = npaths*(ii), npaths*(ii+1)
        LI, UI = 2*(ii), 2*(ii+1)
        mart_figs[ii+1] = plot_martingale_paths(T, mpath_mult[li+1:ui, :],
                                                    mbounds_mult[LI+1:UI, :], horline=1)
        mart_figs[ii+1][:suptitle](L"Martingale components for many paths of $y_{ii+1}$",
                                   fontsize=14)
    end

    return mart_figs
end

function plot_given_paths(T::Integer,
                        ypath::Array, mpath::Array, spath::Array, tpath::Array,
                        mbounds::Array, sbounds::Array; horline::Real=0.0,
                        show_trend::Bool = true)

    # Allocate space
    trange = 1:T

    # Create figure
    fig, ax = subplots(2, 2, sharey=true, figsize=(15, 8))

    # Plot all paths together
    ax[1, 1][:plot](trange, ypath[1, :], label=L"$y_t$", color="k")
    ax[1, 1][:plot](trange, mpath[1, :], label=L"$m_t$", color="m")
    ax[1, 1][:plot](trange, spath[1, :], label=L"$s_t$", color="g")
    if show_trend == true
        ax[1, 1][:plot](trange, tpath[1, :], label=L"t_t", color="r")
    end
    ax[1, 1][:axhline](horline, color="k", linestyle = "-.")
    ax[1, 1][:set_title]("One Path of All Variables")
    ax[1, 1][:legend](loc="top left")

    # Plot Martingale Component
    ax[1, 2][:plot](trange, mpath[1, :], "m")
    ax[1, 2][:plot](trange, mpath', alpha=0.45, color="m")
    ub = mbounds[2, :]
    lb = mbounds[1, :]
    ax[1, 2][:fill_between](trange, lb, ub, alpha=0.25, color="m")
    ax[1, 2][:set_title]("Martingale Components for Many Paths")
    ax[1, 2][:axhline](horline, color="k", linestyle = "-.")

    # Plot Stationary Component
    ax[2, 1][:plot](spath[1, :], color="g")
    ax[2, 1][:plot](spath', alpha=0.25, color="g")
    ub = sbounds[2, :]
    lb = sbounds[1, :]
    ax[2, 1][:fill_between](trange, lb, ub, alpha=0.25, color="g")
    ax[2, 1][:axhline](horline, color="k", linestyle = "-.")
    ax[2, 1][:set_title]("Stationary Components for Many Paths")

    # Plot Trend Component
    if show_trend == true
        ax[2, 2][:plot](tpath', color="r")
    end
    ax[2, 2][:set_title]("Trend Components for Many Paths")
    ax[2, 2][:axhline](horline, color="k", linestyle = "-.")

    return fig
end

function plot_martingale_paths(T::Integer,
                        mpath::Array, mbounds::Array;
                        horline::Real=1,
                        show_trend::Bool = false)
    # Allocate space
    trange = 1:T

    # Create figure
    fig, ax = subplots(1, 1, figsize=(10, 6))

    # Plot Martingale Component
    ub = mbounds[2, :]
    lb = mbounds[1, :]
    ax[:fill_between](trange, lb, ub, color="#ffccff")
    ax[:axhline](horline, color="k", linestyle = "-.")
    ax[:plot](trange, mpath', linewidth=0.25, color="#4c4c4c")

    return fig
end


text:
# julia
using ArgParse
using BioStructures
# input parameters function
function parse_commandline()

    s = ArgParseSettings(description = "read a pdb file from PDB and convert it to fasta without the need to download it choosing the model and chain")
    @add_arg_table s begin
        "--pdb"
            help = "input pdb ids"
        "--model"
            help = "model to make the fasta file for"
        "--chain"
            help = "chain in model to make the fasta file for"
        "--out"
            help = "output fasta file"
    end
    return parse_args(s)
end

function main()
    parsed_args = parse_commandline()
    println(parsed_args)
# main
    downloadpdb(parsed_args["pdb"]) do fp
        s = read(fp, PDB)
        open(parsed_args["out"],"a") do io
            println(io,">",parsed_args["pdb"],"_",parsed_args["chain"],"\n",LongAminoAcidSeq(s[parse(Int, parsed_args["model"])][parsed_args["chain"]], standardselector, gaps=false))
        end
    end
end

main()


text:
# Custom, hardcoded streamers

abstract type CustomROOTStruct end

# KM3NeT

struct KM3NETDAQHit <: CustomROOTStruct
    dom_id::Int32
    channel_id::UInt8
    tdc::Int32
    tot::UInt8
end
function readtype(io::IO, T::Type{KM3NETDAQHit})
    T(readtype(io, Int32), read(io, UInt8), read(io, Int32), read(io, UInt8))
end
# Experimental implementation for maximum performance (using reinterpret)
primitive type DAQHit 80 end
function Base.getproperty(hit::DAQHit, s::Symbol)
    r = Ref(hit)
    GC.@preserve r begin
        if s === :dom_id
            return ntoh(unsafe_load(Ptr{Int32}(Base.unsafe_convert(Ptr{Cvoid}, r))))
        elseif s === :channel_id
            return unsafe_load(Ptr{UInt8}(Base.unsafe_convert(Ptr{Cvoid}, r)+4))
        elseif s === :tdc
            return unsafe_load(Ptr{UInt32}(Base.unsafe_convert(Ptr{Cvoid}, r)+5))
        elseif s === :tot
            return unsafe_load(Ptr{UInt8}(Base.unsafe_convert(Ptr{Cvoid}, r)+9))
        end
    end
    error("unknown field $s of type $(typeof(hit))")
end
Base.show(io::IO, h::DAQHit) = print(io, "DAQHit(", h.dom_id, ',', h.channel_id, ',', h.tdc, ',', h.tot, ')')

struct KM3NETDAQTriggeredHit
    dom_id::Int32
    channel_id::UInt8
    tdc::Int32
    tot::UInt8
    trigger_mask::UInt64
end
function readtype(io::IO, T::Type{KM3NETDAQTriggeredHit})
    dom_id = readtype(io, Int32)
    channel_id = read(io, UInt8)
    tdc = read(io, Int32)
    tot = read(io, UInt8)
    skip(io, 6)
    trigger_mask = readtype(io, UInt64)
    T(dom_id, channel_id, tdc, tot, trigger_mask)
end

struct KM3NETDAQEventHeader
    detector_id::Int32
    run::Int32
    frame_index::Int32
    UTC_seconds::UInt32
    UTC_16nanosecondcycles::UInt32
    trigger_counter::UInt64
    trigger_mask::UInt64
    overlays::UInt32
end

function readtype(io::IO, T::Type{KM3NETDAQEventHeader})
    skip(io, 18)
    detector_id = readtype(io, Int32)
    run = readtype(io, Int32)
    frame_index = readtype(io, Int32)
    skip(io, 6)
    UTC_seconds = readtype(io, UInt32)
    UTC_16nanosecondcycles = readtype(io, UInt32)
    skip(io, 6)
    trigger_counter = readtype(io, UInt64)
    skip(io, 6)
    trigger_mask = readtype(io, UInt64)
    overlays = readtype(io, UInt32)
    T(detector_id, run, frame_index, UTC_seconds, UTC_16nanosecondcycles, trigger_counter, trigger_mask, overlays)
end


text:
using Test
using LogicCircuits
import .LogicCircuits.IO: zoo_cnf, zoo_dnf

@testset "CNF file parser tests" begin

   circuit = zoo_cnf("8.cnf")

   @test circuit isa UnstLogicalΔ
   @test circuit[end] isa ⋀Node
   @test circuit[end - 1] isa ⋁Node
   @test num_variables(circuit) == 227
   @test num_nodes(circuit) == 1168

end

@testset "DNF file parser tests" begin

   circuit = zoo_dnf("8.dnf")

   @test circuit isa UnstLogicalΔ
   @test circuit[end] isa ⋁Node
   @test circuit[end - 1] isa ⋀Node
   @test num_variables(circuit) == 227
   @test num_nodes(circuit) == 1168

end


text:
function reroot(root::T, new_root::String)::T where T<:AbstractNode

    new_tree = deepcopy(root)
    root_node = find_by_name(new_tree, new_root)

    mother = get_mother(root_node)

    recursive_invert(mother, root_node)

    root_node.root = true
    new_tree.root = false

    set_binary!(root_node)

    return root_node
end

function recursive_invert(old_mother::T, old_daughter::T)::T where T
    if old_mother.root == true
        # arrived at the root
        od = remove_child!(old_mother, old_daughter)
        add_child!(od, old_mother)
        return od
    end
        od1 = recursive_invert(get_mother(old_mother), old_mother)
        od = remove_child!(od1, old_daughter)
        add_child!(od, od1)
        return od
end


text:
using HetaSimulator, Plots

p = load_platform("./time-events"; rm_out = true);
model = models(p)[:nameless]

# for s1 : dose

res = sim(model; tspan = (0, 1200), observables = [:s1])
plot(res)

res = sim(
    model;
    tspan = (0, 1200),
    observables = [:comp1, :s1],
    events_active = [:dose_sw1=>false]
)
plot(res)

# for s2 : injection

res = sim(model;
    tspan = (0, 69600),
    saveat = collect(48000.:720.:69600.),
    observables = [:s2],
    events_active = [:dose_sw1=>false]
    )
plot(res)

res = sim(model;
    tspan = (0, 69_600),
    saveat = collect(48000.:720.:69600.),
    observables = [:comp2, :s2],
    events_active = [:infusion_on2=>false, :infusion_off2=>false]
)
plot(res)

max_stars_repo_name: hetalang/qsp-benchmark

max_issues_repo_name: hetalang/qsp-benchmark

max_forks_repo_name: hetalang/qsp-benchmark


text:
using DataFrames
using Gadfly
using Colors

include("theory.jl")
include("transitions.jl")

tr_chains_d1(θ, N_H, N_E, p_H, p_E) = tr_chains(θ, N_H, N_E, p_H, p_E, 1)
tr_chains_d2(θ, N_H, N_E, p_H, p_E) = tr_chains(θ, N_H, N_E, p_H, p_E, 2)
tr_chains_d10(θ, N_H, N_E, p_H, p_E) = tr_chains(θ, N_H, N_E, p_H, p_E, 10)
tr_chains_d50(θ, N_H, N_E, p_H, p_E) = tr_chains(θ, N_H, N_E, p_H, p_E, 50)
tr_chains_d100(θ, N_H, N_E, p_H, p_E) = tr_chains(θ, N_H, N_E, p_H, p_E, 100)
conj_chains_d1(λ_E, λ_H, p_H, p_E) = conj_chains(λ_E, λ_H, p_H, p_E, 1)
theory_chains_d1(λ_E, λ_H, p_H, p_E) = theory_chains_upper(λ_E, λ_H, p_H, p_E, 1)
conj_chains_d50(λ_E, λ_H, p_H, p_E) = conj_chains(λ_E, λ_H, p_H, p_E, 50)
theory_chains_d50(λ_E, λ_H, p_H, p_E) = theory_chains_upper(λ_E, λ_H, p_H, p_E, 50)

function run_sim(θ, p_E, p_H, T, transition)
    N_H = 0
    N_E = 0
    hist = DataFrame(time = Array{Float64}(0),
                     N_E = Array{Float64}(0),
                     N_H = Array{Float64}(0))
    for i in 1:T
        (N_H, N_E) = transition(θ, N_H, N_E, p_H, p_E)
        push!(hist, (i, N_E, N_H))
    end
    return hist #mean(hist[convert(3*T/4):T,:N_H])
end

#Runs markov chain, not the Jump Process.
function run_sim(λ_E, λ_H, p_E, p_H, T, transition)
    θ= λ_H/(λ_H + λ_E)
    if transition in [theory_prioE_lower, theory_prioE_upper,
                      conj_prioE, theory_chains_d1,
                      conj_chains_d1, theory_chains_d50,
                      conj_chains_d50]
        return transition(λ_E, λ_H, p_E, p_H)
    else
        h = run_sim(θ, p_E, p_H, T, transition)
        h[:w_H] = h[:N_H]/λ_H
        h[:w_E] = h[:N_E]/λ_E
        return mean(h[Int(3*T/4):T, :w_H])
    end
end

function run_sim_E(λ_E, λ_H, p_E, p_H, T, transition)
    θ= λ_H/(λ_H + λ_E)
    h = run_sim(θ, p_E, p_H, T, transition)
    h[:w_H] = h[:N_H]/λ_H
    h[:w_E] = h[:N_E]/λ_E
    return mean(h[Int(3*T/4):T, :w_E])
end

function run_multiple_sim(λ_E_range, λ_H_range, p_E_range, p_H_range, transitions_list; w_E = false)
    df = DataFrame(λ_H = Array{Float64}(0),
                   λ_E = Array{Float64}(0),
                   θ = Array{Float64}(0),
                   p_H = Array{Float64}(0),
                   p_E = Array{Float64}(0))
    for transition in transitions_list
        df[Symbol(transition)] = Array{Float64}(0)
    end
    v = zeros(5 + length(transitions_list))
    for λ_H in λ_H_range
        v[1] = λ_H
        for λ_E in λ_E_range
            v[2] = λ_E
            θ = λ_H/(λ_H + λ_E)
            v[3] = θ
            for p_H in p_H_range
                v[4] = p_H
                for p_E in p_E_range
                    v[5] = p_E
                    for (k, transition) in enumerate(transitions_list)
                        v[5 + k] = run_sim(λ_E, λ_H, p_E, p_H, T, transition)
                        if w_E
                            v[5 + k] = run_sim_E(λ_E, λ_H, p_E, p_H, T, transition)
                        end
                    end
                    push!(df, v)
                end
            end
        end
    end
    return df
end


text:
using SnoopCompile

botconfig = BotConfig(
  "AcuteML";
  os = ["linux", "windows", "macos"],
  version = [v"1.4.2", v"1.3.1"],
  blacklist = [],
  exhaustive = false,
)

println("Benchmarking the inference time of `using AcuteML`")
snoop_bench(
  botconfig,
  :(using AcuteML),
)

println("Benchmarking the inference time of `using AcuteML` and runtests")
snoop_bench(
  botconfig,
  "$(@__DIR__)/example_script.jl",
)


text:
## port methods of simplify module of SymPy
## from http://docs.sympy.org/0.7.2/modules/simplify/simplify.html

## simple methods (x, args) -> y (y coercion happens via PyCall)
simplify_sympy_meths = (:collect, :rcollect, :separatevars,
                        :radsimp, :ratsimp, :trigsimp, :besselsimp,
                        :powsimp, :combsimp, :hypersimp,
                        :fraction,
                        :simplify, :nsimplify,
                        :posify, :powdenest, :sqrtdenest,
                        :logcombine, :hyperexpand)

expand_sympy_meths = (:expand_trig,
                      :expand_power_base, :expand_power_exp,
                      :expand_log,
                      :expand_func,
                      :hyperexpand
                      )

## special case cse due to output
## If this pattern is common, will need to collect and use metaprogramming.
@doc """
`$(cse)`: a SymPy function.
The SymPy documentation can be found through: http://docs.sympy.org/latest/search.html?q=cse
""" ->
cse{T<:SymbolicObject}(ex::T, args...; kwargs...) = sympy_meth(:cse, ex, args...; kwargs...)
cse{T<:SymbolicObject}(ex::Vector{T}, args...; kwargs...) = sympy_meth(:cse, map(project,ex), args...; kwargs...)

function cse{T<:SymbolicObject, N}(ex::Array{T, N}, args...; kwargs...)
    a,b = cse(ex[:], args...; kwargs...)
    bb = convert(Array{Sym},  reshape(b, size(ex)))
    a, bb
end
export(cse)

## didn't do traversal tools, EPath tools


text:
using GBQ
using Base.Test

# write your own tests here
@test 1 == 1
@test gbq_query("SELECT * FROM `nyc-tlc.green.trips_2014` LIMIT 5;")


text:
module TimeSteppers

export
    AdamsBashforthTimeStepper,
    time_step!,
    compute_w_from_continuity!,
    tendencies

using CUDA
using GPUifyLoops: @launch, @loop, @unroll

import Oceananigans: TimeStepper

using Oceananigans.Architectures: @hascuda, device
using Oceananigans.Architectures
using Oceananigans.Grids
using Oceananigans.Fields
using Oceananigans.Operators
using Oceananigans.Coriolis
using Oceananigans.Buoyancy
using Oceananigans.SurfaceWaves
using Oceananigans.BoundaryConditions
using Oceananigans.Solvers
using Oceananigans.Models
using Oceananigans.Utils

using Oceananigans.TurbulenceClosures:
    calculate_diffusivities!, ∂ⱼ_2ν_Σ₁ⱼ, ∂ⱼ_2ν_Σ₂ⱼ, ∂ⱼ_2ν_Σ₃ⱼ, ∇_κ_∇c

"""
    AbstractTimeStepper

Abstract supertype for time steppers.
"""
abstract type AbstractTimeStepper end

"""
    TimeStepper(name, args...)

Returns a timestepper with name `name`, instantiated with `args...`.

Example
=======

julia> stepper = TimeStepper(:AdamsBashforth, Float64, CPU(), grid, tracernames)
"""
function TimeStepper(name::Symbol, args...)
    fullname = Symbol(name, :TimeStepper)
    return eval(Expr(:call, fullname, args...))
end

# Fallbacks
TimeStepper(stepper::AbstractTimeStepper, args...) = stepper

include("generic_time_stepping.jl")
include("velocity_and_tracer_tendencies.jl")
include("time_stepping_kernels.jl")
include("adams_bashforth.jl")

end # module


text:
#####
##### BaseCase methods
#####

function free_convection_windspeed(surf::SurfaceBase, grid, state, gm::GridMeanVariables, param_set, ::BaseCase)
    prog_gm = center_prog_grid_mean(state)
    aux_tc = center_aux_turbconv(state)
    zi = get_inversion(grid, state, param_set, surf.Ri_bulk_crit)
    wstar = get_wstar(surf.bflux, zi) # yair here zi in TRMM should be adjusted
    surf.windspeed = sqrt(surf.windspeed * surf.windspeed + (1.2 * wstar) * (1.2 * wstar))
    return
end

#####
##### Default SurfaceBase behavior
#####

free_convection_windspeed(surf::SurfaceBase, grid, state, gm::GridMeanVariables, param_set) =
    free_convection_windspeed(surf, grid, state, gm, param_set, BaseCase())

#####
##### SurfaceNone
#####

function update(surf::SurfaceBase{SurfaceNone}, grid, state, gm::GridMeanVariables, param_set)
    # JH: assigning small fluxed so that simulation won"t crash when computing mixing length
    kc_surf = kc_surface(grid)
    surf.windspeed = 0.0001
    surf.zrough = 1e-4
    surf.bflux = 1e-4
    surf.ustar = compute_friction_velocity(param_set, surf.windspeed, surf.bflux, surf.zrough, grid.zc[kc_surf])
    return
end
free_convection_windspeed(surf::SurfaceBase{SurfaceNone}, grid, state, gm::GridMeanVariables, param_set) = nothing

function update(surf::SurfaceBase{SurfaceFixedFlux}, grid, state, gm::GridMeanVariables, param_set)
    kc_surf = kc_surface(grid)
    kf_surf = kf_surface(grid)
    zb = grid.zc[kc_surf]
    aux_gm = center_aux_grid_mean(state)
    prog_gm = center_prog_grid_mean(state)
    p0_f_surf = face_ref_state(state).p0[kf_surf]
    ρ0_f_surf = face_ref_state(state).ρ0[kf_surf]
    α0_f_surf = face_ref_state(state).α0[kf_surf]
    u_gm_surf = prog_gm.u[kc_surf]
    v_gm_surf = prog_gm.v[kc_surf]
    q_tot_gm_surf = prog_gm.q_tot[kc_surf]
    θ_liq_ice_gm_surf = prog_gm.θ_liq_ice[kc_surf]
    T_gm_surf = aux_gm.T[kc_surf]

    ts = thermo_state_pθq(param_set, p0_f_surf, surf.Tsurface, surf.qsurface)
    rho_tflux = surf.shf / TD.cp_m(ts)

    surf.windspeed = sqrt(u_gm_surf^2 + v_gm_surf^2)
    surf.rho_qtflux = surf.lhf / TD.latent_heat_vapor(param_set, surf.Tsurface)

    ts = thermo_state_pθq(param_set, p0_f_surf, θ_liq_ice_gm_surf, q_tot_gm_surf)
    Π = TD.exner(ts)
    surf.rho_hflux = rho_tflux / Π
    surf.bflux = buoyancy_flux(param_set, surf.shf, surf.lhf, surf.Tsurface, surf.qsurface, α0_f_surf, ts)

    if !surf.ustar_fixed
        # Correction to windspeed for free convective cases (Beljaars, QJRMS (1994), 121, pp. 255-270)
        # Value 1.2 is empirical, but should be O(1)
        if surf.windspeed < 0.1  # Limit here is heuristic
            if surf.bflux > 0.0
                free_convection_windspeed(surf, grid, state, gm, param_set)
            else
                print("WARNING: Low windspeed + stable conditions, need to check ustar computation")
                print("surf.bflux ==>", surf.bflux)
                print("surf.shf ==>", surf.shf)
                print("surf.lhf ==>", surf.lhf)
                print("u_gm_surf ==>", u_gm_surf)
                print("v_gm_surf ==>", v_gm_surf)
                print("q_tot_gm_surf ==>", q_tot_gm_surf)
                print("α0_f_surf ==>", α0_f_surf)
                error("Terminating execution.")
            end
        end

        surf.ustar = compute_friction_velocity(param_set, surf.windspeed, surf.bflux, surf.zrough, zb)
    end

    von_karman_const = CPSGS.von_karman_const(param_set)
    surf.obukhov_length = obukhov_length(param_set, surf.ustar, surf.bflux)
    surf.rho_uflux = -ρ0_f_surf * surf.ustar^2 / surf.windspeed * u_gm_surf
    surf.rho_vflux = -ρ0_f_surf * surf.ustar^2 / surf.windspeed * v_gm_surf
    return
end

function update(surf::SurfaceBase{SurfaceFixedCoeffs}, grid, state, gm::GridMeanVariables, param_set)
    kc_surf = kc_surface(grid)
    kf_surf = kf_surface(grid)
    p0_f_surf = face_ref_state(state).p0[kf_surf]
    ρ0_f_surf = face_ref_state(state).ρ0[kf_surf]
    α0_f_surf = face_ref_state(state).α0[kf_surf]
    aux_gm = center_aux_grid_mean(state)
    prog_gm = center_prog_grid_mean(state)
    u_gm_surf = prog_gm.u[kc_surf]
    v_gm_surf = prog_gm.v[kc_surf]
    T_gm_surf = aux_gm.T[kc_surf]
    q_tot_gm_surf = prog_gm.q_tot[kc_surf]
    θ_liq_ice_gm_surf = prog_gm.θ_liq_ice[kc_surf]

    ts = thermo_state_pθq(param_set, p0_f_surf, surf.Tsurface, surf.qsurface)
    surf.windspeed = max(sqrt(u_gm_surf^2 + v_gm_surf^2), 0.01)
    cp_ = TD.cp_m(ts)
    lv = TD.latent_heat_vapor(ts)

    surf.rho_qtflux = -surf.cq * surf.windspeed * (q_tot_gm_surf - surf.qsurface) * ρ0_f_surf
    surf.lhf = lv * surf.rho_qtflux

    ts = thermo_state_pθq(param_set, p0_f_surf, θ_liq_ice_gm_surf, q_tot_gm_surf)
    Π = TD.exner(ts)
    surf.rho_hflux = -surf.ch * surf.windspeed * (θ_liq_ice_gm_surf - surf.Tsurface / Π) * ρ0_f_surf
    surf.shf = cp_ * surf.rho_hflux

    surf.bflux = buoyancy_flux(param_set, surf.shf, surf.lhf, surf.Tsurface, surf.qsurface, α0_f_surf, ts)

    surf.ustar = friction_velocity_given_windspeed(surf.cm, surf.windspeed)
    surf.obukhov_length = obukhov_length(param_set, surf.ustar, surf.bflux)

    surf.rho_uflux = -ρ0_f_surf * surf.ustar^2 / surf.windspeed * u_gm_surf
    surf.rho_vflux = -ρ0_f_surf * surf.ustar^2 / surf.windspeed * v_gm_surf
    return
end

function update(surf::SurfaceBase{SurfaceMoninObukhov}, grid, state, gm::GridMeanVariables, param_set)
    g = CPP.grav(param_set)
    kc_surf = kc_surface(grid)
    kf_surf = kf_surface(grid)
    zb = grid.zc[kc_surf]
    p0_f_surf = face_ref_state(state).p0[kf_surf]
    p0_c_surf = center_ref_state(state).p0[kc_surf]
    ρ0_f_surf = face_ref_state(state).ρ0[kf_surf]
    α0_f_surf = face_ref_state(state).α0[kf_surf]
    u_gm_surf = prog_gm.u[kc_surf]
    v_gm_surf = prog_gm.v[kc_surf]
    q_tot_gm_surf = prog_gm.q_tot[kc_surf]
    θ_liq_ice_gm_surf = prog_gm.θ_liq_ice[kc_surf]
    T_gm_surf = aux_gm.T[kc_surf]
    Pg = surf.ref_params.Pg

    pvg = TD.saturation_vapor_pressure(param_set, TD.PhaseEquil, surf.Tsurface)
    surf.qsurface = TD.q_vap_saturation_from_density(param_set, surf.Tsurface, ρ0_f_surf, pvg)
    lv = TD.latent_heat_vapor(param_set, T_gm_surf)

    phase_part = TD.PhasePartition(surf.qsurface, 0.0, 0.0)
    h_star = TD.liquid_ice_pottemp_given_pressure(param_set, surf.Tsurface, Pg, phase_part)

    ts_g = thermo_state_pθq(param_set, Pg, h_star, surf.qsurface)
    θ_virt_g = TD.virtual_pottemp(ts_g)

    ts_b = thermo_state_pθq(param_set, p0_c_surf, θ_liq_ice_gm_surf, surf.qsurface)
    θ_virt_b = TD.virtual_pottemp(ts_b)

    surf.windspeed = sqrt(u_gm_surf^2 + v_gm_surf^2)
    Nb2 = g / θ_virt_g * (θ_virt_b - θ_virt_g) / zb
    Ri = Nb2 * zb * zb / (surf.windspeed * surf.windspeed)

    surf.cm, surf.ch, surf.obukhov_length = exchange_coefficients_byun(param_set, Ri, zb, surf.zrough)
    surf.rho_uflux = -surf.cm * surf.windspeed * u_gm_surf * ρ0_f_surf
    surf.rho_vflux = -surf.cm * surf.windspeed * v_gm_surf * ρ0_f_surf

    surf.rho_hflux = -surf.ch * surf.windspeed * (θ_liq_ice_gm_surf - h_star) * ρ0_f_surf
    surf.rho_qtflux = -surf.ch * surf.windspeed * (q_tot_gm_surf - surf.qsurface) * ρ0_f_surf
    surf.lhf = lv * surf.rho_qtflux

    ts = thermo_state_pθq(param_set, p0_f_surf, surf.Tsurface, surf.qsurface)
    surf.shf = TD.cp_m(ts) * surf.rho_hflux

    surf.bflux = buoyancy_flux(param_set, surf.shf, surf.lhf, surf.Tsurface, surf.qsurface, α0_f_surf, ts)
    surf.ustar = friction_velocity_given_windspeed(surf.cm, surf.windspeed)

    return
end

function update(surf::SurfaceBase{SurfaceMoninObukhovDry}, grid, state, gm::GridMeanVariables, param_set)
    g = CPP.grav(param_set)
    kc_surf = kc_surface(grid)
    kf_surf = kf_surface(grid)

    zb = grid.zc[kc_surf]
    p0_f_surf = face_ref_state(state).p0[kf_surf]
    p0_c_surf = center_ref_state(state).p0[kc_surf]
    ρ0_f_surf = face_ref_state(state).ρ0[kf_surf]
    α0_f_surf = face_ref_state(state).α0[kf_surf]
    prog_gm = center_prog_grid_mean(state)
    aux_gm = center_aux_grid_mean(state)
    u_gm_surf = prog_gm.u[kc_surf]
    v_gm_surf = prog_gm.v[kc_surf]
    q_tot_gm_surf = prog_gm.q_tot[kc_surf]
    θ_liq_ice_gm_surf = prog_gm.θ_liq_ice[kc_surf]
    T_gm_surf = aux_gm.T[kc_surf]
    Pg = surf.ref_params.Pg

    pvg = TD.saturation_vapor_pressure(param_set, TD.PhaseEquil, surf.Tsurface)
    surf.qsurface = TD.q_vap_saturation_from_density(param_set, surf.Tsurface, ρ0_f_surf, pvg)
    lv = TD.latent_heat_vapor(param_set, T_gm_surf)

    phase_part = TD.PhasePartition(surf.qsurface, 0.0, 0.0)
    h_star = TD.liquid_ice_pottemp_given_pressure(param_set, surf.Tsurface, Pg, phase_part)

    ts_g = thermo_state_pθq(param_set, Pg, h_star, surf.qsurface)
    θ_virt_g = TD.virtual_pottemp(ts_g)

    ts_b = thermo_state_pθq(param_set, p0_c_surf, θ_liq_ice_gm_surf, surf.qsurface)
    θ_virt_b = TD.virtual_pottemp(ts_b)

    surf.windspeed = sqrt(u_gm_surf^2 + v_gm_surf^2)
    Nb2 = g / θ_virt_g * (θ_virt_b - θ_virt_g) / zb
    Ri = Nb2 * zb * zb / (surf.windspeed * surf.windspeed)

    surf.cm, surf.ch, surf.obukhov_length = exchange_coefficients_byun(param_set, Ri, zb, surf.zrough)
    surf.rho_uflux = -surf.cm * surf.windspeed * u_gm_surf * ρ0_f_surf
    surf.rho_vflux = -surf.cm * surf.windspeed * v_gm_surf * ρ0_f_surf

    surf.rho_hflux = -surf.ch * surf.windspeed * (θ_liq_ice_gm_surf - h_star) * ρ0_f_surf
    surf.rho_qtflux = -surf.ch * surf.windspeed * (q_tot_gm_surf - surf.qsurface) * ρ0_f_surf
    surf.lhf = lv * 0.0

    ts = thermo_state_pθq(param_set, p0_f_surf, surf.Tsurface, surf.qsurface)
    surf.shf = TD.cp_m(ts) * surf.rho_hflux

    surf.bflux = buoyancy_flux(param_set, surf.shf, surf.lhf, surf.Tsurface, surf.qsurface, α0_f_surf, ts)
    surf.ustar = friction_velocity_given_windspeed(surf.cm, surf.windspeed)

    return
end

function update(surf::SurfaceBase{SurfaceSullivanPatton}, grid, state, gm::GridMeanVariables, param_set)
    g = CPP.grav(param_set)
    R_d = CPP.R_d(param_set)
    kc_surf = kc_surface(grid)
    kf_surf = kf_surface(grid)

    zb = grid.zc[kc_surf]

    p0_f_surf = face_ref_state(state).p0[kf_surf]
    p0_c_surf = center_ref_state(state).p0[kc_surf]
    ρ0_f_surf = face_ref_state(state).ρ0[kf_surf]
    α0_c_surf = center_ref_state(state).α0[kc_surf]
    u_gm_surf = prog_gm.u[kc_surf]
    v_gm_surf = prog_gm.v[kc_surf]
    q_tot_gm_surf = prog_gm.q_tot[kc_surf]
    θ_liq_ice_gm_surf = prog_gm.θ_liq_ice[kc_surf]
    T_gm_surf = aux_gm.T[kc_surf]
    Pg = surf.ref_params.Pg

    ts = thermo_state_pθq(param_set, p0_f_surf, surf.Tsurface, surf.qsurface)
    lv = TD.latent_heat_vapor(param_set, T_gm_surf)
    T0 = p0_c_surf * α0_c_surf / R_d # TODO: can we use a thermo state here?

    θ_flux = 0.24
    phase_part = TD.PhasePartition(q_tot_gm_surf, 0.0, 0.0)
    Π = TD.exner_given_pressure(param_set, p0_f_surf, phase_part)
    surf.bflux = g * θ_flux * Π / T0

    pvg = TD.saturation_vapor_pressure(param_set, TD.PhaseEquil, surf.Tsurface)
    surf.qsurface = TD.q_vap_saturation_from_density(param_set, surf.Tsurface, ρ0_f_surf, pvg)
    h_star = TD.liquid_ice_pottemp_given_pressure(param_set, surf.Tsurface, Pg, phase_part)

    ts_g = thermo_state_pθq(param_set, Pg, h_star, surf.qsurface)
    θ_virt_g = TD.virtual_pottemp(ts_g)

    ts_b = thermo_state_pθq(param_set, p0_c_surf, θ_liq_ice_gm_surf, surf.qsurface)
    θ_virt_b = TD.virtual_pottemp(ts_b)

    surf.windspeed = sqrt(u_gm_surf^2 + v_gm_surf^2)
    Nb2 = g / θ_virt_g * (θ_virt_b - θ_virt_g) / zb
    Ri = Nb2 * zb * zb / (surf.windspeed * surf.windspeed)

    surf.cm, surf.ch, surf.obukhov_length = exchange_coefficients_byun(param_set, Ri, zb, surf.zrough)
    surf.rho_uflux = -surf.cm * surf.windspeed * u_gm_surf * ρ0_f_surf
    surf.rho_vflux = -surf.cm * surf.windspeed * v_gm_surf * ρ0_f_surf

    surf.rho_hflux = -surf.ch * surf.windspeed * (θ_liq_ice_gm_surf - h_star) * ρ0_f_surf
    surf.rho_qtflux = -surf.ch * surf.windspeed * (q_tot_gm_surf - surf.qsurface) * ρ0_f_surf
    surf.lhf = lv * surf.rho_qtflux
    surf.shf = TD.cp_m(ts) * surf.rho_hflux

    surf.ustar = friction_velocity_given_windspeed(surf.cm, surf.windspeed)
    return
end


text:
module CatcherModule

import ..GridWorlds as GW
import Random
import ReinforcementLearningBase as RLBase

#####
##### game logic
#####

const NUM_OBJECTS = 2
const AGENT = 1
const GEM = 2
const NUM_ACTIONS = 3

mutable struct Catcher{R, RNG} <: GW.AbstractGridWorld
    tile_map::BitArray{3}
    agent_position::CartesianIndex{2}
    reward::R
    rng::RNG
    done::Bool
    gem_position::CartesianIndex{2}
    gem_reward::R
    terminal_penalty::R
end

function Catcher(; R = Float32, height = 8, width = 8, rng = Random.GLOBAL_RNG)
    tile_map = falses(NUM_OBJECTS, height, width)

    gem_position = CartesianIndex(1, 1)
    tile_map[GEM, gem_position] = true

    agent_position = CartesianIndex(height, 1)
    tile_map[AGENT, agent_position] = true

    reward = zero(R)
    done = false
    gem_reward = one(R)
    terminal_penalty = -one(R)

    env = Catcher(tile_map, agent_position, reward, rng, done, gem_position, gem_reward, terminal_penalty)

    GW.reset!(env)

    return env
end

function GW.reset!(env::Catcher)
    tile_map = env.tile_map
    rng = env.rng

    _, height, width = size(tile_map)

    tile_map[AGENT, env.agent_position] = false
    tile_map[GEM, env.gem_position] = false

    new_gem_position = CartesianIndex(1, rand(rng, 1:width))
    env.gem_position = new_gem_position
    tile_map[GEM, new_gem_position] = true

    new_agent_position = CartesianIndex(height, rand(rng, 1:width))
    env.agent_position = new_agent_position
    tile_map[AGENT, new_agent_position] = true

    env.reward = zero(env.reward)
    env.done = false

    return nothing
end

function GW.act!(env::Catcher, action)
    @assert action in Base.OneTo(NUM_ACTIONS) "Invalid action $(action). Action must be in Base.OneTo($(NUM_ACTIONS))"

    tile_map = env.tile_map
    rng = env.rng
    _, height, width = size(tile_map)
    agent_position = env.agent_position
    gem_position = env.gem_position

    if action == 1
        if agent_position[2] == 1
            new_agent_position = agent_position
        else
            new_agent_position = GW.move_left(agent_position)
        end
    elseif action == 2
        if agent_position[2] == width
            new_agent_position = agent_position
        else
            new_agent_position = GW.move_right(agent_position)
        end
    else
        new_agent_position = agent_position
    end

    tile_map[AGENT, agent_position] = false
    env.agent_position = new_agent_position
    tile_map[AGENT, new_agent_position] = true

    if gem_position[1] == height
        new_gem_position = CartesianIndex(1, rand(rng, 1 : width))
    else
        new_gem_position = GW.move_down(gem_position)
    end

    tile_map[GEM, gem_position] = false
    env.gem_position = new_gem_position
    tile_map[GEM, new_gem_position] = true

    if new_gem_position[1] == height
        if new_gem_position[2] == new_agent_position[2]
            env.done = false
            env.reward = env.gem_reward
        else
            env.done = true
            env.reward = env.terminal_penalty
        end
    else
        env.done = false
        env.reward = zero(env.reward)
    end

    return nothing
end

#####
##### miscellaneous
#####

GW.get_height(env::Catcher) = size(env.tile_map, 2)
GW.get_width(env::Catcher) = size(env.tile_map, 3)

function GW.get_pretty_tile_map(env::Catcher, position::CartesianIndex{2})
    characters = ('☻', '♦', '⋅')

    object = findfirst(@view env.tile_map[:, position])
    if isnothing(object)
        return characters[end]
    else
        return characters[object]
    end
end

GW.get_object_names(env::Catcher) = (:AGENT, :GEM)
GW.get_action_names(env::Catcher) = (:MOVE_LEFT, :MOVE_RIGHT, :NO_MOVE)

function Base.show(io::IO, ::MIME"text/plain", env::Catcher)
    str = "tile_map:\n"
    str = str * GW.get_pretty_tile_map(env)
    str = str * "\nreward: $(env.reward)\ndone: $(env.done)"
    str = str * "\naction_names: $(GW.get_action_names(env))"
    str = str * "\nobject_names: $(GW.get_object_names(env))"
    print(io, str)
    return nothing
end

GW.get_action_keys(env::Catcher) = ('a', 'd', 's')

#####
##### RLBase API
#####

RLBase.StateStyle(env::GW.RLBaseEnv{E}) where {E <: Catcher} = RLBase.InternalState{Any}()
RLBase.state_space(env::GW.RLBaseEnv{E}, ::RLBase.InternalState) where {E <: Catcher} = nothing
RLBase.state(env::GW.RLBaseEnv{E}, ::RLBase.InternalState) where {E <: Catcher} = env.env.tile_map

RLBase.reset!(env::GW.RLBaseEnv{E}) where {E <: Catcher} = GW.reset!(env.env)

RLBase.action_space(env::GW.RLBaseEnv{E}) where {E <: Catcher} = Base.OneTo(NUM_ACTIONS)
(env::GW.RLBaseEnv{E})(action) where {E <: Catcher} = GW.act!(env.env, action)

RLBase.reward(env::GW.RLBaseEnv{E}) where {E <: Catcher} = env.env.reward
RLBase.is_terminated(env::GW.RLBaseEnv{E}) where {E <: Catcher} = env.env.done

end # module


text:
# This file is a part of ChangesOfVariables.jl, licensed under the MIT License (MIT).

import ForwardDiff

torv_and_back(V::AbstractVector{<:Real}) = V, identity
torv_and_back(x::Real) = [x], V -> V[1]
torv_and_back(x::Complex) = [real(x), imag(x)], V -> Complex(V[1], V[2])
torv_and_back(x::NTuple{N}) where N = [x...], V -> ntuple(i -> V[i], Val(N))

function torv_and_back(x::Ref)
    xval = x[]
    V, to_xval = torv_and_back(xval)
    back_to_ref(V) = Ref(to_xval(V))
    return (V, back_to_ref)
end

torv_and_back(A::AbstractArray{<:Real}) = vec(A), V -> reshape(V, size(A))

function torv_and_back(A::AbstractArray{Complex{T}, N}) where {T<:Real, N}
    RA = cat(real.(A), imag.(A), dims = N+1)
    V, to_array = torv_and_back(RA)
    function back_to_complex(V)
        RA = to_array(V)
        Complex.(view(RA, map(_ -> :, size(A))..., 1), view(RA, map(_ -> :, size(A))..., 2))
    end
    return (V, back_to_complex)
end

function getjacobian(f, x)
    V, to_x = torv_and_back(x)
    vf(V) = torv_and_back(f(to_x(V)))[1]
    ForwardDiff.jacobian(vf, V)
end

foo(x) = inv(exp(-x) + 1)