text:
#------- SET UP DESIGN PARAMETERS --------#

struct plyproperties
    names#::Array{String,1}
    plies::Array{Composites.material,1}
end

struct geometry
    t::Array{Float64,1}
    s::Array{Float64,1}
    x::Array{Float64,1} #x-location of airfoil section (positive x = downstream)
    y::Array{Float64,1} #y-location of airfoil section (positive y = spanwise)
    z::Array{Float64,1} #z-location of airfoil section (positive z = up)
    chord::Array{Float64,1}
    normalchord::Array{Float64,1} #normal chord length
    twist::Array{Float64,1} #airfoil twist in degrees
    sweep::Array{Float64,1}
    dihedral::Array{Float64,1}
    airfoilthickness::Array{Float64,1}
    xaf::Array{Float64,2}
    yaf::Array{Float64,2}
end

struct compstructure
    material::Array{Array{Composites.material,1},1}
    laminate::Array{Composites.laminate,1}
    A::Array{Array{Float64,2},1}
    B::Array{Array{Float64,2},1}
    D::Array{Array{Float64,2},1}
    bucklingstrain::Array{Float64,2}
    precompinput::Array{PreComp.input,1}
    precompoutput::Array{PreComp.output,1}
end

struct designparameters
    # configuration assumptions
    p_n_akima::Int64
    p_radii::Array{Float64,1}
    p_webloc::Array{Float64,1}

    # operating point assumptions
    altitude::Float64 #m

    # material properties
    plyprops::plyproperties
    p_usedmaterials
    p_orientation::Array{Float64,1}

    # Propeller Design
    p_lam_t::Array{Float64,1}

    p_chord::Array{Float64,1}
    p_sweep_d::Array{Float64,1}
    p_twist_d::Array{Float64,1}
    p_pitch_d::Float64
    p_dihedral_d::Array{Float64,1}
    p_airfoilthickness::Array{Float64,1}
    p_aerocenter::Array{Float64,1}
    p_x_offset::Array{Float64,1}
    velocity::Float64
    p_Rtip::Float64
    blades::Float64
    numprops::Int64
    machtip::Float64
    p_rpm::Float64
    kv::Float64
    i0::Float64
    batt_mass::Float64

    # Wing Design
    w_n_linear::Int64
    w_webloc::Array{Float64,1}
    w_lam_t::Array{Float64,1}
    w_chord::Array{Float64,1}
    w_aoa_d::Float64
    w_twist_d::Array{Float64,1}
    w_halfspan::Float64
    w_nondim_halfspany::Array{Float64,1}
    w_sweep_d::Array{Float64,1}
    w_dihedral_d::Array{Float64,1}
    w_usedmaterials #array of strings???
    w_orientation::Array{Float64,1}
    w_x_offset::Array{Float64,1}
    w_aerocenter::Array{Float64,1}

    n_CP::Float64
    payload::Float64
    max_mass::Float64
    gravity::Float64
    h_set::Float64
    specif_energy::Float64
    required_range::Float64
    required_takeoff_dist::Float64
    altitude_accel::Float64
    altitude_climb::Float64
    altitude_cruise::Float64
    required_noise::Float64
    required_cruise::Float64

end

function designparameters(;
    p_n_akima = 20,
    p_radii::Array{Float64,1}=Float64[0.0,0.2,.4,.6,.8,1.0],
    p_webloc::Array{Float64,1}=[1.0],
    altitude::Float64=18288.0,
    plyprops::plyproperties=plyproperties(),
    p_usedmaterials = ["highmodulus_uni","highmodulus_weave","taylor_foam"],
    p_orientation=[0.0,-45],
    # structural design variables
    p_lam_t = ones(length(p_radii)*length(p_usedmaterials)),
    # geometric design variables
    p_chord = Float64[0.1, 0.1, 0.1, 0.1, 0.1],
    p_sweep_d = Float64[0.0,0.0,0.0,0.0,0.0],
    p_twist_d = Float64[30.0, 30.0, 30.0, 30.0],
    p_pitch_d = 0.0,
    p_dihedral_d = fill(0.0,5),
    p_airfoilthickness = Float64[0.147538, 0.132307, 0.127404, 0.130604, 0.245836, 0.0629818],
    p_aerocenter = ones(p_radii)*0.25,
    p_x_offset = zeros(p_radii),
    # operating point design variables
    velocity = 15.40662113981876,
    # propulsion design variables
    p_Rtip = 2.25,
    blades = 2,
    numprops = 4,
    machtip = 0.8,
    p_rpm = 1020.45,
    kv = 4.271513933940295,
    i0 = 0.1,
    batt_mass = 1E8,

    w_n_linear = 20,
    w_halfspan = 10.0,
    w_webloc::Array{Float64,1}=[0.40],
    w_usedmaterials = ["highmodulus_uni","highmodulus_weave","taylor_foam"],
    w_nondim_halfspany = linspace(0,w_halfspan,length(w_chord)),
    w_lam_t = ones(length(w_nondim_halfspany)*length(w_usedmaterials)),
    w_chord = [0.5,0.5],
    w_aoa_d = 0.0,
    w_twist_d = [0.0,0.0],
    w_sweep_d = [0.0,0.0],
    w_dihedral_d = [0.0,0.0],
    w_orientation = [0.0,45.0,45.0], #uni, weave, webweave
    w_x_offset = ones(w_nondim_halfspany)*0.25,
    w_aerocenter = ones(w_nondim_halfspany)*0.275,

    n_CP = 60.0,
    payload = 684.0,
    max_mass = 1500.0,
    gravity = 9.81,
    h_set = 20.0, #meters
    specif_energy = 300.0, #wh/kg
    required_range = 150000.0, #km
    required_takeoff_dist = 395.0, #meters
    altitude_accel = 0.0,
    altitude_climb = 0.0,
    altitude_cruise = 3000.0,
    required_noise = 68.0,
    required_cruise = 28.0,

    )

    parameters = designparameters(p_n_akima,p_radii,p_webloc,altitude,plyprops,p_usedmaterials,p_orientation,p_lam_t,
    p_chord,p_sweep_d,p_twist_d,p_pitch_d,p_dihedral_d,p_airfoilthickness,p_aerocenter,p_x_offset,velocity,p_Rtip,
    blades,numprops,machtip,p_rpm,kv,i0,batt_mass,w_n_linear,w_webloc,w_lam_t,w_chord,w_aoa_d,w_twist_d,w_halfspan,w_nondim_halfspany,w_sweep_d,w_dihedral_d,w_usedmaterials,
    w_orientation,w_x_offset,w_aerocenter,n_CP,payload,max_mass,gravity,h_set,
    specif_energy,required_range,required_takeoff_dist,altitude_accel,altitude_climb,
    altitude_cruise,required_noise,required_cruise)
    return parameters
end

struct propwingdesign
    #Design Variables at this scope
    batt_mass::Float64
    w_aoa::Float64
    p_rpm::Float64
    p_pitch::Float64
    velocity::Float64

    #Design Parameters at this scope
    constraints
    af
    p_xaf
    p_yaf
    p_xafstrain
    p_yafstrain
    w_xaf
    w_yaf
    w_xafstrain
    w_yafstrain
    spline_3D_freestream
    spline_3D_blown

    # Design variables unique to the system function
    p_lam_tin
    p_orientation::Array{Float64,1}
    p_chord_pts::Array{Float64,1}
    p_twist_pts::Array{Float64,1}
    p_Rtip::Float64
    kv::Float64
    i0::Float64
    w_chord_pts::Array{Float64,1}
    w_twist_pts::Array{Float64,1}
    w_halfspan::Float64
    w_sweep_pts::Array{Float64,1}
    w_dihedral_pts::Array{Float64,1}
    w_lam_tin
    w_orientation::Array{Float64,1}
    w_airfoilthickness::Array{Float64,1}

    #Design variables unique to the system function
    etaprop_center::Array{Float64,1}
    p_n_akima::Int64
    p_radii::Array{Float64,1}
    altitude::Float64
    numprops::Int64
    p_x_offset::Array{Float64,1}
    p_aerocenter::Array{Float64,1}
    p_usedmaterials
    p_webloc::Array{Float64,1}
    plyprops
    blades::Float64
    w_n_linear::Int64
    w_x_offset::Array{Float64,1}
    w_aerocenter::Array{Float64,1}
    w_nondim_halfspany::Array{Float64,1}
    w_usedmaterials
    w_webloc::Array{Float64,1}
    prop_tilt::Float64

    #Control parameters
    detailed::Bool
    VTKfilename
    savefinalvtk::Bool
    verification::Bool
    plots::Bool
    n_CP::Int64

end

function propwingdesign(;
    #Design Variables at this scope
    batt_mass = 0.0,
    w_aoa = 0.0,
    p_rpm = 0.0,
    p_pitch = 0.0,
    velocity = 0.0,

    #Design Parameters at this scope
    constraints = [],
    af = [],
    p_xaf = [],
    p_yaf = [],
    p_xafstrain = [],
    p_yafstrain = [],
    w_xaf = [],
    w_yaf = [],
    w_xafstrain = [],
    w_yafstrain = [],
    spline_3D_freestream = [],
    spline_3D_blown = [],

    # Design variables unique to the system function
    p_lam_tin= [],
    p_orientation = zeros(2),
    p_chord_pts = zeros(2),
    p_twist_pts = zeros(2),
    p_Rtip =0.0,
    kv =0.0,
    i0 =0.0,
    w_chord_pts = zeros(2),
    w_twist_pts = zeros(2),
    w_halfspan =0.0,
    w_sweep_pts = zeros(2),
    w_dihedral_pts = zeros(2),
    w_lam_tin= [],
    w_orientation = zeros(2),
    w_airfoilthickness = zeros(2),

    #Design variables unique to the system function
    etaprop_center = zeros(2),
    p_n_akima = 10 ,
    p_radii = zeros(2),
    altitude =0.0,
    numprops = 10 ,
    p_x_offset = zeros(2),
    p_aerocenter = zeros(2),
    p_usedmaterials = [],
    p_webloc = zeros(2),
    plyprops = [],
    blades = 10 ,
    w_n_linear = 10 ,
    w_x_offset = zeros(2),
    w_aerocenter = zeros(2),
    w_nondim_halfspany = zeros(2),
    w_usedmaterials = [],
    w_webloc = zeros(2),
    prop_tilt = 0.0,

    #Control parameters
    detailed = false,
    VTKfilename = "notdefined.vtk",
    savefinalvtk = false,
    verification = false,
    plots = false,
    n_CP = 60,
    )

    return propwingdesign(batt_mass, w_aoa, p_rpm, p_pitch, velocity,
        constraints, af, p_xaf, p_yaf, p_xafstrain, p_yafstrain, w_xaf, w_yaf,
        w_xafstrain, w_yafstrain, spline_3D_freestream, spline_3D_blown,
        p_lam_tin, p_orientation, p_chord_pts, p_twist_pts, p_Rtip, kv, i0,
        w_chord_pts, w_twist_pts, w_halfspan, w_sweep_pts, w_dihedral_pts,
         w_lam_tin, w_orientation, w_airfoilthickness, etaprop_center,
        p_n_akima, p_radii, altitude, numprops, p_x_offset, p_aerocenter,
        p_usedmaterials, p_webloc, plyprops, blades, w_n_linear, w_x_offset,
        w_aerocenter, w_nondim_halfspany, w_usedmaterials, w_webloc, prop_tilt,
        detailed, VTKfilename, savefinalvtk, verification, plots, n_CP)
end

struct propwing_out
    #Propulsion
    ETA_total::Float64
    ETA_m::Float64
    ETA_p::Float64
    # Propeller
    CT::Float64
    CQ::Float64
    Np::Array{Float64,1}
    Tp::Array{Float64,1}
    p_alphas::Array{Float64,1}
    uprop::Array{Float64,1}
    vprop::Array{Float64,1}
    total_thrust::Float64
    advance_ratio::Float64
    torque::Float64
    db_test::Array{Float64,1}

    # Motor
    V_m::Float64
    I_m::Float64
    # Masses
    propwing_mass::Float64
    mass_prop::Float64
    mass_ESC::Float64
    mass_motor::Float64
    mass_wing::Float64

    # Wing
    panels#::Array{Float64,1}
    Lift::Float64
    Drag::Float64
    viscousDrag::Float64
    a::Float64
    Fp#::Array{Float64,1}
    cllocal::Array{Float64,1}
    cl::Array{Float64,1}
    Vinfeff::Array{Float64,1}
    alphaeff::Array{Float64,1}

    #Structures
    p_c_stress #::Array{Float64,2}
    w_c_stress #::Array{Float64,2}
    p_c_buckling #::Array{Float64,2}
    w_c_buckling #::Array{Float64,2}

    fail::Bool
end

function propwing_out(;
    #Propulsion
    ETA_total = 0.0,
    ETA_m = 0.0,
    ETA_p = 0.0,
    # Propeller
    CT = 0.0,
    CQ = 0.0,
    Np = zeros(2),
    Tp = zeros(2),
    p_alphas = zeros(2),
    uprop = zeros(2),
    vprop = zeros(2),
    total_thrust = 0.0,
    advance_ratio = 0.0,
    torque = 0.0,
    db_test = zeros(2),

    # Motor
    V_m = 0.0,
    I_m = 0.0,
    # Masses
    propwing_mass = 0.0,
    mass_prop = 0.0,
    mass_ESC = 0.0,
    mass_motor = 0.0,
    mass_wing = 0.0,

    # Wing
    panels = zeros(2),
    Lift = 0.0,
    Drag = 0.0,
    viscousDrag = 0.0,
    a = 0.0,
    Fp = [],# = zeros(2)
    cllocal = zeros(2),
    cl = zeros(2),
    Vinfeff = zeros(2),
    alphaeff = zeros(2),

    #Structures
    p_c_stress = zeros(2),
    w_c_stress = zeros(2),
    p_c_buckling = zeros(2),
    w_c_buckling = zeros(2),
    fail = false,
)

    return propwing_out(ETA_total,ETA_m,ETA_p,CT,CQ,Np,Tp,p_alphas,uprop,vprop,total_thrust,
    advance_ratio,torque,db_test,V_m,I_m,propwing_mass,mass_prop,mass_ESC,
    mass_motor,mass_wing,panels,Lift,Drag,viscousDrag,a,Fp,cllocal,cl,Vinfeff,alphaeff,
    p_c_stress,w_c_stress,p_c_buckling,w_c_buckling, fail)

end

struct system_output
    range_true::Float64
    dist_accel::Float64
    dist_climb::Float64

    energy_accel::Float64
    energy_climb::Float64
    energy_cruise::Float64

    t_accel::Float64
    t_transition::Float64
    t_climb::Float64
    t_cruise::Float64

    total_mass::Float64
    total_weight::Float64

    gama::Float64
    radius::Float64
    xtr::Float64
    ytr::Float64
    H_left::Float64
    x_left::Float64

    optconvals::Dict{Symbol,Array{Float64,1}}

    payload::Float64
    max_mass::Float64
    gravity::Float64
    h_set::Float64
    specif_energy::Float64
    required_takeoff_dist::Float64
    altitude_accel::Float64
    altitude_climb::Float64
    altitude_cruise::Float64
    required_noise::Float64
    etaprop_center::Array{Float64,1}
    batt_mass_accel::Float64
    batt_mass_climb::Float64
    batt_mass_cruise::Float64
    Fl::Float64
end

function system_output(;
    range_true = 0.0,
    dist_accel = 0.0,
    dist_climb = 0.0,

    energy_accel = 0.0,
    energy_climb = 0.0,
    energy_cruise = 0.0,

    t_accel = 0.0,
    t_transition = 0.0,
    t_climb = 0.0,
    t_cruise = 0.0,

    total_mass = 0.0,
    total_weight = 0.0,

    gama = 0.0,
    radius = 0.0,
    xtr = 0.0,
    ytr = 0.0,
    H_left = 0.0,
    x_left = 0.0,

    optconvals = Dict{Symbol,Array{Float64,1}}(),

    payload = 0.0,
    max_mass = 0.0,
    gravity = 0.0,
    h_set = 0.0,
    specif_energy = 0.0,
    required_takeoff_dist = 0.0,
    altitude_accel = 0.0,
    altitude_climb = 0.0,
    altitude_cruise = 0.0,
    required_noise = 0.0,
    etaprop_center = [0.0,0.0],
    batt_mass_accel = 0.0,
    batt_mass_climb = 0.0,
    batt_mass_cruise = 0.0,
    Fl = 0.0,
    )

    return system_output(range_true,dist_accel,dist_climb, energy_accel,
    energy_climb,energy_climb, t_accel,t_transition,t_climb,t_cruise,
    total_mass,total_weight, gama,radius,xtr,ytr,H_left,x_left, optconvals,
    payload,max_mass,gravity,h_set,specif_energy,required_takeoff_dist,
    altitude_accel,altitude_climb,altitude_cruise,required_noise,
    etaprop_center,batt_mass_accel,batt_mass_climb,batt_mass_cruise,Fl)

end

# This file sets up default values and bounds for design variables
function getoptxbounds!(paramdict,designparams)

    ##############################################################################
    #----------------------- Prop Design Variables ------------------------#
    ##############################################################################

    lb_scale = 1.0 #1E-6
    ub_scale = 1E3
    # Skin Fabric Ply Thickness Design Variables
    p_nstations = length(designparams.p_radii)
    matnames = designparams.plyprops.names
    p_nplies = length(designparams.p_usedmaterials)

    name = :p_lam_t
    default = zeros(designparams.p_lam_t)
    lowerbound = zeros(designparams.p_lam_t)
    upperbound = zeros(designparams.p_lam_t)
    scaling = zeros(designparams.p_lam_t)
    for i = 1:p_nplies
        idx = find(matnames -> matnames == designparams.p_usedmaterials[i],matnames)
        material = designparams.plyprops.plies[idx[1]]

        default1 = designparams.p_lam_t[(i-1)*p_nstations+1:(i)*p_nstations]
        lowerbound1 = fill(designparams.plyprops.plies[idx[1]].t*lb_scale,p_nstations)
        upperbound1 = fill(designparams.plyprops.plies[idx[1]].t*ub_scale,p_nstations)
        scaling1 = fill(0.01/designparams.plyprops.plies[idx[1]].t,p_nstations)

        default[(i-1)*p_nstations+1:(i)*p_nstations] = default1
        lowerbound[(i-1)*p_nstations+1:(i)*p_nstations] = lowerbound1
        upperbound[(i-1)*p_nstations+1:(i)*p_nstations] = upperbound1
        scaling[(i-1)*p_nstations+1:(i)*p_nstations] = scaling1

    end
    OptParams.addparam!(paramdict,name,default,lowerbound,upperbound,scaling)

    name = :p_orientation
    default = designparams.p_orientation
    lowerbound = fill(-90.0,length(default))
    upperbound = fill(90.0,length(default))
    scaling = fill(1.0./90.0,length(default))
    OptParams.addparam!(paramdict,name,default,lowerbound,upperbound,scaling)

    name = :p_chord
    default = designparams.p_chord
    lowerbound = fill(0.01,p_nstations)
    lowerbound[end] = 0.001 #tip, works well for splines
    upperbound = fill(0.8,p_nstations)
    scaling = fill(1E2,p_nstations)
    OptParams.addparam!(paramdict,name,default,lowerbound,upperbound,scaling)

    name = :p_sweep_d
    default = designparams.p_sweep_d
    lowerbound = fill(0.0,p_nstations-1)
    upperbound = fill(45.0,p_nstations-1)
    scaling = fill(pi/180.0,p_nstations-1)
    OptParams.addparam!(paramdict,name,default,lowerbound,upperbound,scaling)

    OptParams.addparam!(paramdict,:p_pitch_d_accel, designparams.p_pitch_d*10.0, -25.0, 45.0, 1E1)
    OptParams.addparam!(paramdict,:p_pitch_d_climb, designparams.p_pitch_d*10.0, -25.0, 45.0, 1E0)
    OptParams.addparam!(paramdict,:p_pitch_d_cruise, designparams.p_pitch_d*-2, -5.0, 80.0, 1E2)

    name = :p_twist_d
    default = designparams.p_twist_d
    lowerbound = fill(-0.0,p_nstations)
    upperbound = fill(45.0,p_nstations)
    scaling = fill(0.1,p_nstations)
    OptParams.addparam!(paramdict,name,default,lowerbound,upperbound,scaling)

    name = :p_dihedral_d
    default = designparams.p_dihedral_d
    lowerbound = fill(0.0,p_nstations-1)
    upperbound = fill(10.0,p_nstations-1)
    scaling = fill(pi/180.0,p_nstations-1)
    OptParams.addparam!(paramdict,name,default,lowerbound,upperbound,scaling)

    name = :p_airfoilthickness
    default = designparams.p_airfoilthickness
    lowerbound = fill(0.050,p_nstations)
    upperbound = fill(0.250,p_nstations)
    scaling = fill(10.0,p_nstations)
    OptParams.addparam!(paramdict,name,default,lowerbound,upperbound,scaling)

    OptParams.addparam!(paramdict,:p_Rtip,designparams.p_Rtip,0.1,5.0,1E-1)

    ##############################################################################
    #----------------------- Wing Design Variables ------------------------#
    ##############################################################################

    # Skin Fabric Ply Thickness Design Variables
    w_nstations = length(designparams.w_chord)

    name = :w_lam_t
    default = zeros(designparams.w_lam_t)
    lowerbound = zeros(designparams.w_lam_t)
    upperbound = zeros(designparams.w_lam_t)
    scaling = zeros(designparams.w_lam_t)
    for i = 1:length(designparams.w_usedmaterials)
        idx = find(matnames -> matnames == designparams.w_usedmaterials[i],matnames)
        material = designparams.plyprops.plies[idx[1]]

        default1 = designparams.w_lam_t[(i-1)*w_nstations+1:(i)*w_nstations]
        lowerbound1 = fill(designparams.plyprops.plies[idx[1]].t*lb_scale,w_nstations)
        upperbound1 = fill(designparams.plyprops.plies[idx[1]].t*ub_scale,w_nstations)
        scaling1 = fill(0.01/designparams.plyprops.plies[idx[1]].t,w_nstations)

        default[(i-1)*w_nstations+1:(i)*w_nstations] = default1
        lowerbound[(i-1)*w_nstations+1:(i)*w_nstations] = lowerbound1
        upperbound[(i-1)*w_nstations+1:(i)*w_nstations] = upperbound1
        scaling[(i-1)*w_nstations+1:(i)*w_nstations] = scaling1

    end
    OptParams.addparam!(paramdict,name,default,lowerbound,upperbound,scaling)

    name = :w_orientation
    default = designparams.w_orientation
    lowerbound = fill(-90.0,length(default))
    upperbound = fill(90.0,length(default))
    scaling = fill(1.0./90.0,length(default))
    OptParams.addparam!(paramdict,name,default,lowerbound,upperbound,scaling)

    name = :w_chord
    default = designparams.w_chord
    lowerbound = fill(0.1,w_nstations)
    upperbound = fill(5.0,w_nstations)
    scaling = fill(100.0/designparams.w_halfspan,w_nstations)
    OptParams.addparam!(paramdict,name,default,lowerbound,upperbound,scaling)

    OptParams.addparam!(paramdict,:w_aoa_d_accel,designparams.w_aoa_d*1.1,2.0,40.0,1E-2)
    OptParams.addparam!(paramdict,:w_aoa_d_climb,designparams.w_aoa_d*0.5,5.0,40.0,1E0)
    OptParams.addparam!(paramdict,:w_aoa_d_cruise,designparams.w_aoa_d*0.5,0.20,20.0,1E0)

    name = :w_twist_d
    default = designparams.w_twist_d
    lowerbound = fill(-10.0,w_nstations)
    upperbound = fill(20.0,w_nstations)
    scaling = fill(10.0,w_nstations)
    OptParams.addparam!(paramdict,name,default,lowerbound,upperbound,scaling)

    name = :w_halfspan
    default = designparams.w_halfspan
    lowerbound = 5.699
    upperbound = 12.0
    scaling = 1.0/designparams.w_halfspan
    OptParams.addparam!(paramdict,name,default,lowerbound,upperbound,scaling)

    name = :w_sweep_d
    default = designparams.w_sweep_d
    lowerbound = fill(0.0,w_nstations)
    upperbound = fill(20.0,w_nstations)
    scaling = fill(10.0,w_nstations)
    OptParams.addparam!(paramdict,name,default,lowerbound,upperbound,scaling)

    name = :w_dihedral_d
    default = designparams.w_dihedral_d
    lowerbound = fill(0.0,w_nstations)
    upperbound = fill(20.0,w_nstations)
    scaling = fill(10.0,w_nstations)
    OptParams.addparam!(paramdict,name,default,lowerbound,upperbound,scaling)

    name = :w_airfoilthickness
    default = designparams.p_airfoilthickness
    lowerbound = fill(0.050,p_nstations)
    upperbound = fill(0.250,p_nstations)
    scaling = fill(10.0,p_nstations)
    OptParams.addparam!(paramdict,name,default,lowerbound,upperbound,scaling)

    ##############################################################################
    #-----------------  Operating Point Design Variables ------------------------#
    ##############################################################################

    OptParams.addparam!(paramdict,:velocity_accel,designparams.velocity/1.4,5.0,28.0,1E1)
    OptParams.addparam!(paramdict,:velocity_climb,designparams.velocity*1.0,5.0,60.0,1E-2)
    OptParams.addparam!(paramdict,:velocity_cruise,designparams.velocity*1.0,20.0,100.0,1E1)

    ##############################################################################
    #-------------------------  Propulsion Design Variables ---------------------#
    ##############################################################################

    OptParams.addparam!(paramdict,:p_rpm_accel,designparams.p_rpm*4.5,200.0,40000.0,1E-1)
    OptParams.addparam!(paramdict,:p_rpm_climb,designparams.p_rpm*4.6,200.0,40000.0,1E-2)
    OptParams.addparam!(paramdict,:p_rpm_cruise,designparams.p_rpm*2.5,100.0,10000.0,1E0)
    OptParams.addparam!(paramdict,:kv,designparams.kv,5.0,5000.0,1E-0)
    OptParams.addparam!(paramdict,:i0,designparams.i0,0.1,6.0,1E0)
    OptParams.addparam!(paramdict,:batt_mass_accel,designparams.batt_mass,0.001,100.0,1E0)
    OptParams.addparam!(paramdict,:batt_mass_climb,designparams.batt_mass,0.01,100.0,1E-1)
    OptParams.addparam!(paramdict,:batt_mass_cruise,designparams.batt_mass,5.0,1000.0,1E-1)
    OptParams.addparam!(paramdict,:number_blades,2.0,2.0,5.0,1E-0)
    OptParams.addparam!(paramdict,:radius,100.0,2.0,2000.0,1E-1)
    OptParams.addparam!(paramdict,:gama_d,30.0,0.1,89.9,1E-1)
end

#------ SET UP OPTIMIZATION PARAMETERS -------#

mutable struct controlparams
    objective::Symbol
    printfreq::Int
    savevtk::Int
    printobjective::Bool
    printdetailedoutput::Bool
    printconstraintviolations::Bool
end

function controlparams(;
    objective::Symbol = :mass,
    printfreq::Int = 1,
    savevtk::Int = 3,
    printobjective::Bool = true,
    printdetailedoutput::Bool = true,
    printconstraintviolations::Bool = true,
    )
    controlparams(objective,printfreq,savevtk,printobjective,printdetailedoutput,
    printconstraintviolations)
end

# Define design variables
function getoptparams()
    designvar =  [
    :p_lam_t, #1-12
    # :p_orientation,
    :p_chord, #13-18
    :p_twist_d, #19-24
    :p_pitch_d_accel, #25
    :p_pitch_d_climb, #25
    :p_pitch_d_cruise, #25

    :velocity_accel,
    :velocity_climb,
    :velocity_cruise,
    :p_Rtip,
    :p_rpm_accel,
    :p_rpm_climb,
    :p_rpm_cruise,
    # :kv, #28
    # :i0, #29

    # :w_orientation,
    :w_chord, #30-32
    :w_aoa_d_accel,   #multi
    :w_aoa_d_climb,   #multi
    :w_aoa_d_cruise,   #multi
    # :w_twist_d, #33-35
    # :w_halfspan,
    # :w_sweep_d,
    # :w_dihedral_d,
    # :w_lam_t, #36-end
    :batt_mass_accel,
    :batt_mass_climb,
    :batt_mass_cruise,
    ]

    # Define general constraints
    constraints = [
    :p_materialfailure,
    :p_localbuckling,
    # :w_materialfailure,
    # :w_localbuckling,

    :totalmass,
    :machtip,
    :Lift,
    :Drag,
    :noise,
    # :p_chordlamthick,
    # :range,
    :takeoff_dist,
    :stall,
    :capacity,
    :p_alphas,
    ]

    # define run controls
    ctlparams = controlparams(
    objective= :range, #:mass,#:range #:ETA_total_cruise,
    printfreq = 50,
    savevtk = 500,
    printobjective=true,
    printdetailedoutput=true,
    printconstraintviolations=true,
    )

    return designvar, constraints, ctlparams
end

#------- MATERIAL PROPERTIES -------#
function plyproperties()
    names = ["highmodulus_uni",
    "standard_uni",
    "MR60H",
    "T3900_uni",
    "T700_uni",
    "ELT5500",
    "UDCarbon",
    "highmodulus_weave",
    "standard_weave",
    "T3900_weave",
    "T700_weave",
    "Gelcoat",
    "Triax",
    "Saertex",
    "taylor_foam",
    "SNL_foam"]

    e1  =zeros(16)
    e2  =zeros(16)
    g12 =zeros(16)
    anu =zeros(16)
    rho =zeros(16)
    xt  =zeros(16)
    xc  =zeros(16)
    yt  =zeros(16)
    yc  =zeros(16)
    s   =zeros(16)
    plythickness =zeros(16)

    # "highmodulus_uni"
    e1[1]  = 175.0e9
    e2[1]  = 8.0e9
    g12[1] = 5.0e9
    anu[1] = 0.30
    rho[1] = 1600.0
    xt[1]  = min(1.0,1355.656/1682.011)*1000e6 #mean -> A-basis
    xc[1]  = min(1.0,1103.943/1396.504)*850e6 #mean -> A-basis
    yt[1]  = min(1.0,39.226/52.975)*40e6 #mean -> A-basis
    yc[1]  = min(1.0,235.434/282.439)*200e6 #mean -> A-basis
    s[1]   = min(1.0,142.411/159.516)*60e6 #mean -> A-basis
    plythickness[1] = 0.152e-3
    # "standard_uni"
    e1[2]  = 135.0e9
    e2[2]  = 10.0e9
    g12[2] = 5.0e9
    anu[2] = 0.30
    rho[2] = 1600.0
    xt[2]  = min(1.0,1355.656/1682.011)*1500e6 #mean -> A-basis
    xc[2]  = min(1.0,1103.943/1396.504)*1200e6 #mean -> A-basis
    yt[2]  = min(1.0,39.226/52.975)*50e6 #mean -> A-basis
    yc[2]  = min(1.0,235.434/282.439)*250e6 #mean -> A-basis
    s[2]   = min(1.0,142.411/159.516)*70e6 #mean -> A-basis
    plythickness[2] = 0.152e-3
    # "MR60H"
    e1[3]  = (165e9+150e9)/2.0
    e2[3]  = 8.56e9
    g12[3] = 4.39e9
    anu[3] = 0.326
    rho[3] = 1810.0
    xt[3]  = min(1.0,1355.656/1682.011)*3190e6 #mean -> A-basis
    xc[3]  = min(1.0,1103.943/1396.504)*1440e6 #mean -> A-basis
    yt[3]  = min(1.0,39.226/52.975)*82.0e6 #mean -> A-basis
    yc[3]  = min(1.0,235.434/282.439)*200.0e6 #mean -> A-basis
    s[3]   = min(1.0,142.411/159.516)*141e6 #mean -> A-basis
    plythickness[3] = 0.152e-3
    # "T3900_uni"
    e1[4]  = (148e9+131e9)/2.0
    e2[4]  = (9.7e9+9.7e9)/2.0
    g12[4] = 4.83e9
    anu[4] = 0.33
    rho[4] = 1573.0
    xt[4]  = min(1.0,1355.656/1682.011)*2830e6 #CTD mean -> A-basis
    xc[4]  = min(1.0,1103.943/1396.504)*1772e6 #CTD mean -> A-basis
    yt[4]  = min(1.0,39.226/52.975)*56.9e6 #CTD mean -> A-basis
    yc[4]  = min(1.0,235.434/282.439)*303e6 #CTD mean -> A-basis
    s[4]   = min(1.0,142.411/159.516)*89.6e6 #CTD mean -> A-basis
    plythickness[4] = 0.191e-3
    # "T700_uni"
    e1[5]=120.8e9
    e2[5]=11.57e9
    g12[5]=5.219e9
    anu[5]=0.350
    rho[5]=1525.0
    xt[5]=1356e6
    xc[5]=1104e6
    yt[5]=39.23e6
    yc[5]=235.4e6
    s[5]=142.4e6
    plythickness[5] = 0.152e-3
    # "ELT5500"
    e1[6]=41.8e9
    e2[6]=14.00e9
    g12[6]=2.630e9
    anu[6]=0.280
    rho[6]=1920.0
    xt[6]=972.0e6
    xc[6]=702.0e6
    yt[6]=100.0e6 #made up
    yc[6]=100.0e6 #made up
    s[6]=100.0e6 #made up
    plythickness[6] = 0.47e-3
    # "UDCarbon"
    e1[7]=114.5e9
    e2[7]=8.39e9
    g12[7]=5.990e9
    anu[7]=0.270
    rho[7]=1220.0
    xt[7]=1546.0e6
    xc[7]=1047.0e6
    yt[7]=100.0e6 #made up
    yc[7]=100.0e6 #made up
    s[7]=100.0e6 #made up
    plythickness[7] = 0.47e-3

    # FABRICS
    # "highmodulus_weave"
    e1[8]  = 85.0e9
    e2[8]  = 85.0e9
    g12[8] = 5.0e9
    anu[8] = 0.10
    rho[8] = 1600.0
    xt[8]  = min(1.0,1355.656/1682.011)*350e6 #mean -> A-basis
    xc[8]  = min(1.0,1103.943/1396.504)*150e6 #mean -> A-basis
    yt[8]  = min(1.0,39.226/52.975)*350e6 #mean -> A-basis
    yc[8]  = min(1.0,235.434/282.439)*150e6 #mean -> A-basis
    s[8]   = min(1.0,142.411/159.516)*35e6 #mean -> A-basis
    plythickness[8] = 0.218e-3
    # "standard_weave"
    e1[9]  = 70.0e9
    e2[9]  = 70.0e9
    g12[9] = 5.0e9
    anu[9] = 0.10
    rho[9] = 1600.0
    xt[9]  = min(1.0,1355.656/1682.011)*600e6 #mean -> A-basis
    xc[9]  = min(1.0,1103.943/1396.504)*570e6 #mean -> A-basis
    yt[9]  = min(1.0,39.226/52.975)*600e6 #mean -> A-basis
    yc[9]  = min(1.0,235.434/282.439)*570e6 #mean -> A-basis
    s[9]   = min(1.0,142.411/159.516)*90e6 #mean -> A-basis
    plythickness[9] = 0.218e-3
    # "T3900_weave"
    e1[10]  = (70.3e9+71.0e9)/2.0
    e2[10]  = (68.9e9+67.6e9)/2.0
    g12[10] = 4.6e9
    anu[10] = 0.032
    rho[10] = 1551.0
    xt[10]  = min(1.0,701.302/803.236)*1055e6 #CTD mean -> A-basis
    xc[10]  = min(1.0,549.748/749.955)*676e6 #CTD mean -> A-basis
    yt[10]  = min(1.0,557.575/722.602)*945e6 #CTD mean -> A-basis
    yc[10]  = min(1.0,604.067/741.866)*614e6 #CTD mean -> A-basis
    s[10]   = min(1.0,138.440/154.888)*79.3e6 #CTD mean -> A-basis
    plythickness[10] = 0.218e-3
    # "T700_weave"
    e1[11]=55.82e9
    e2[11]=52.10e9
    g12[11]=4.295e9
    anu[11]=0.085
    rho[11]=1501.0
    xt[11]=701.32e6
    yt[11]=557.59e6
    xc[11]=549.77e6
    yc[11]=604.08e6
    s[11]=138.44e6
    plythickness[11] = 0.218e-3
    # "Gelcoat"
    e1[12]=3.44e9
    e2[12]=3.44e9
    g12[12]=1.38e9
    anu[12]=0.3
    rho[12]=1235.0
    xt[12]=100.0e6 #made up
    xc[12]=100.0e6 #made up
    yt[12]=100.0e6 #made up
    yc[12]=100.0e6 #made up
    s[12]=100.0e6 #made up
    plythickness[12] = 0.05e-3
    # "Triax"
    e1[13]=27.7e9
    e2[13]=13.65e9
    g12[13]=7.2e9
    anu[13]=0.39
    rho[13]=1850.0
    xt[13]=700.0e6
    xc[13]=100.0e6 #made up
    yt[13]=700.0e6
    yc[13]=100.0e6 #made up
    s[13]=100.0e6 #made up
    plythickness[13] = 0.94e-3
    # "Saertex"
    e1[14]=13.6e9
    e2[14]=13.3e9
    g12[14]=11.8e9
    anu[14]=0.49
    rho[14]=1780.0
    xt[14]=144.0e6
    xc[14]=213.0e6
    yt[14]=144.0e6
    yc[14]=213.0e6
    s[14]=100.0e6 #made up
    plythickness[14] = 1.0e-3

    #FOAMS
    # "taylor_foam"
    e1[15]=48.0e6
    e2[15]=48.0e6
    g12[15]=28.0e6
    anu[15]=0.3
    rho[15]=75.0
    xt[15]=100.0e6 #made up
    yt[15]=100.0e6 #made up
    xc[15]=100.0e6 #made up
    yc[15]=100.0e6 #made up
    s[15]=100.0e6 #made up
    plythickness[15]=1.0E-3
    # "SNL_foam"
    e1[16]=256.0e6
    e2[16]=256.0e6
    g12[16]=22.0e6
    anu[16]=0.3
    rho[16]=200.0
    xt[16]=100.0e6 #made up
    yt[16]=100.0e6 #made up
    xc[16]=100.0e6 #made up
    yc[16]=100.0e6 #made up
    s[16]=100.0e6 #made up
    plythickness[16]=1.0E-3

    return plyproperties(names,Composites.material(e1,e2,g12,anu,rho,xt,xc,yt,yc,s,plythickness))
end


text:
"""
figures.jl
"""

WDIR = joinpath(@__DIR__, "../../")
include(joinpath(WDIR, "src", "spin", "spin.jl"))

using LaTeXStrings
using Printf
import Plots
using Statistics

# paths
const EXPERIMENT_NAME = "figures"
const SAVE_PATH = joinpath(SPIN_OUT_PATH, EXPERIMENT_NAME)
const F3C_DATA_FILE_PATH = joinpath(SAVE_PATH, "f3c.h5")
const F3D_DATA_FILE_PATH = joinpath(SAVE_PATH, "f3d.h5")

# types
@enum PulseType begin
    analytic = 1
    qoc = 2
    s2 = 3
    s4 = 4
    d2 = 5
    d3 = 6
    s2b = 7
    s4b = 8
    d2b = 9
    d3b = 10
    corpse = 11
    d1 = 12
    sut8 = 13
    d1b = 14
    sut8b = 15
    d1bb = 16
    d1bbb = 17
end

const GT_LIST = [zpiby2, ypiby2, xpiby2]

# common dict keys
const SAVEFP_KEY = 1
const SAVET_KEY = 2
const DATAFP_KEY = 3
const COLOR_KEY = 4
const ACORDS_KEY = 5
const MARKER_KEY = 6
const LCORDS_KEY = 7
const DATA2FP_KEY = 8
const INVAL = 99999

### FIGURE 1 ###

const F1_DATA = Dict(
    zpiby2 => Dict(
        qoc => joinpath(SPIN_OUT_PATH, "spin15/00209_spin15.h5"),
        analytic => joinpath(SPIN_OUT_PATH, "spin14/00000_spin14.h5"),
    ),
    ypiby2 => Dict(
        qoc => joinpath(SPIN_OUT_PATH, "spin15/00205_spin15.h5"),
        analytic => joinpath(SPIN_OUT_PATH, "spin14/00003_spin14.h5"),
    ),
    xpiby2 => Dict(
        qoc => joinpath(SPIN_OUT_PATH, "spin15/00262_spin15.h5"),
        analytic => joinpath(SPIN_OUT_PATH, "spin14/00004_spin14.h5"),
    ),
)
const F1_PT_LIST = [analytic, qoc]

function gen_1a()
    gate_types = [Integer(gate_type) for gate_type in GT_LIST]
    pulse_types = [Integer(pulse_type) for pulse_type in F1_PT_LIST]
    save_file_paths = Array{String, 1}([])
    for (i, gate_type) in enumerate(GT_LIST)
        for (j, pulse_type) in enumerate(F1_PT_LIST)
            save_file_path = F1_DATA[gate_type][pulse_type]
            push!(save_file_paths, save_file_path)
        end
    end

    data_file_path = generate_file_path("h5", "f1a", SAVE_PATH)
    h5open(data_file_path, "w") do data_file
        write(data_file, "save_file_paths", save_file_paths)
        write(data_file, "gate_types", gate_types)
        write(data_file, "pulse_types", pulse_types)
    end
    println("Saved f1a data to $(data_file_path)")
end

const F1B_SAMPLE_COUNT = Integer(5e2)
function gen_1b()
    gate_types = GT_LIST
    gate_types_integer = [Integer(gt) for gt in gate_types]
    pulse_types = F1_PT_LIST
    pulse_types_integer = [Integer(pt) for pt in pulse_types]
    amps_fit = Array(range(0, stop=MAX_CONTROL_NORM_0, length=F1B_SAMPLE_COUNT))
    t1s_fit =  map(amp_t1_spline_cubic, amps_fit)
    amps_data = -1 .* map(fbfq_amp, FBFQ_ARRAY)
    t1s_data = T1_ARRAY
    t1s_data_err = T1_ARRAY_ERR
    avg_amps = Array{Float64, 1}([])
    avg_amps_t1 = Array{Float64, 1}([])

    for (i, gate_type) in enumerate(gate_types)
        for (j, pulse_type) in enumerate(pulse_types)
            save_file_path = F1_DATA[gate_type][pulse_type]
            (controls, controls_dt_inv, evolution_time) = grab_controls(save_file_path)
            avg_amp = mean(map(abs, controls))
            avg_amp_t1 = amp_t1_spline_cubic(avg_amp)
            push!(avg_amps, avg_amp)
            push!(avg_amps_t1, avg_amp_t1)
        end
    end

    data_file_path = generate_file_path("h5", "f1b", SAVE_PATH)
    h5open(data_file_path, "w") do data_file
        write(data_file, "amps_fit", amps_fit)
        write(data_file, "t1s_fit", t1s_fit)
        write(data_file, "amps_data", amps_data)
        write(data_file, "t1s_data", t1s_data)
        write(data_file, "t1s_data_err", t1s_data_err)
        write(data_file, "avg_amps", avg_amps)
        write(data_file, "avg_amps_t1", avg_amps_t1)
        write(data_file, "pulse_types", pulse_types_integer)
        write(data_file, "gate_types", gate_types_integer)
    end
    println("Saved f1b data to $(data_file_path)")
end

const F1C_GATE_COUNT = Integer(2e3)
const F1C_AVG_COUNT = Integer(1e3)
function gen_1c(;use_previous=true)
    gate_types_integer = [Integer(gt) for gt in GT_LIST]
    gate_type_count = size(GT_LIST)[1]
    pulse_types_integer = [Integer(pt) for pt in F1_PT_LIST]
    pulse_type_count = size(F1_PT_LIST)[1]
    gate_errors = zeros(gate_type_count, pulse_type_count, F1C_AVG_COUNT, F1C_GATE_COUNT + 1)
    save_file_paths = Array{String, 2}(undef, gate_type_count, pulse_type_count)

    # check for previous computation
    if use_previous
        data_file_path_old = latest_file_path("h5", "f1c", SAVE_PATH)
    else
        data_file_path_old = nothing
    end
    if isnothing(data_file_path_old)
        (save_file_paths_old = gate_errors_old = gate_type_count_old = pulse_type_count_old
         = avg_count_old = nothing)
    else
        (save_file_paths_old, gate_errors_old,
         ) = h5open(data_file_path_old, "r") do data_file_old
             save_file_paths_old = read(data_file_old, "save_file_paths")
             gate_errors_old = read(data_file_old, "gate_errors")
             gate_type_count_old = size(gate_errors_old)[1]
             pulse_type_count_old = size(gate_errors_old)[2]
             avg_count_old = size(gate_errors_old)[3]
             gate_count_old = size(gate_errors_old)[4] - 1
             if gate_count_old != F1C_GATE_COUNT
                 save_file_paths_old = gate_errors_old = nothing
             end
             return (save_file_paths_old, gate_errors_old, gate_type_count_old,
                     pulse_type_count_old, avg_count_old)
         end
    end

    for (i, gate_type) in enumerate(GT_LIST)
        println("gt[$(i)]: $(gate_type)")
        for (j, pulse_type) in enumerate(F1_PT_LIST)
            print("pt[$(j)]: $(pulse_type) ")
            save_file_path = F1_DATA[gate_type][pulse_type]
            save_file_paths[i, j] = save_file_path
            save_file_path_old = ((isnothing(save_file_paths_old) || i > gate_type_count_old
                                   || j > pulse_type_count_old)
                                  ? nothing : save_file_paths_old[i, j])
            if save_file_path == ""
                println("")
                continue
            end
            if pulse_type == analytic
                if gate_type == zpiby2
                    dynamics_type = zpiby2t1
                elseif gate_type == ypiby2
                    dynamics_type = ypiby2t1
                elseif gate_type == xpiby2
                    dynamics_type = xpiby2t1
                end
            else
                dynamics_type = lindbladt1
            end
            for k = 1:F1C_AVG_COUNT
                # skip redundant computation
                if (!isnothing(save_file_path_old) && save_file_path == save_file_path_old
                    && k <= avg_count_old)
                    gate_errors[i, j, k, :] = gate_errors_old[i, j, k, :]
                else
                    res = run_sim_prop(
                        F1C_GATE_COUNT, gate_type; dynamics_type=dynamics_type,
                        save_file_path=save_file_path, save=false, seed=k
                    )
                    gate_errors[i, j, k, :] = 1 .- res["fidelities"]
                end
            end
            println("")
        end
    end

    data_file_path = generate_file_path("h5", "f1c", SAVE_PATH)
    h5open(data_file_path, "w") do data_file
        write(data_file, "pulse_types", pulse_types_integer)
        write(data_file, "gate_types", gate_types_integer)
        write(data_file, "gate_errors", gate_errors)
        write(data_file, "save_file_paths", save_file_paths)
    end
    println("Saved f1c data to $(data_file_path)")
end

### FIGURE 2 ###
const F2A_DATA_ZPIBY2 = Dict(
    analytic => joinpath(SPIN_OUT_PATH, "spin14/00000_spin14.h5"),
    s2 => joinpath(SPIN_OUT_PATH, "spin12/00808_spin12.h5"),
    sut8 => joinpath(SPIN_OUT_PATH, "spin30/lanthanum/00062_spin30.h5"),
    d1 => joinpath(SPIN_OUT_PATH, "spin11/00539_spin11.h5"),
    d2 => joinpath(SPIN_OUT_PATH, "spin11/00547_spin11.h5"),
)
const F2A_PT_LIST = [analytic, s2, sut8, d1, d2]
function gen_2a(;gate_type=zpiby2)
    if gate_type == zpiby2
        data = F2A_DATA_ZPIBY2
    end
    pulse_types_integer = [Integer(pulse_type) for pulse_type in F2A_PT_LIST]
    pulse_type_count = size(F2A_PT_LIST)[1]
    save_file_paths = Array{String, 1}(undef, pulse_type_count)
    for (i, pulse_type) in enumerate(F2A_PT_LIST)
        save_file_paths[i] = data[pulse_type]
    end

    data_file_path = generate_file_path("h5", "f2a", SAVE_PATH)
    h5open(data_file_path, "w") do data_file
        write(data_file, "save_file_paths", save_file_paths)
        write(data_file, "pulse_types", pulse_types_integer)
        write(data_file, "gate_type", Integer(gate_type))
    end
    println("Saved f2a data to $(data_file_path)")
end

const F2B_DATA_ZPIBY2 = Dict(
    analytic => joinpath(SPIN_OUT_PATH, "spin14/00000_spin14.h5"),
    d1 => joinpath(SPIN_OUT_PATH, "spin11/00468_spin11.h5"),
    d1b => joinpath(SPIN_OUT_PATH, "spin11/00539_spin11.h5"),
    d1bb => joinpath(SPIN_OUT_PATH, "spin11/00556_spin11.h5"),
    d1bbb => joinpath(SPIN_OUT_PATH, "spin11/00585_spin11.h5"),
)
const F2B_TRIAL_COUNT = Integer(1e2)
const F2B_FQ_DEV = 3e-2
const F2B_PT_LIST = [analytic, d1, d1b, d1bb, d1bbb]
const F2B_AVG_COUNT = 1000
function gen_2b(;use_previous=true, gate_type=zpiby2)
    @assert iseven(F2B_TRIAL_COUNT)
    if gate_type == zpiby2
        data = F2B_DATA_ZPIBY2
        dynamics_type_analytic = zpiby2nodis
    end
    pulse_types_integer = [Integer(pt) for pt in F2B_PT_LIST]
    pulse_type_count = size(F2B_PT_LIST)[1]
    fq_devs = Array(range(-F2B_FQ_DEV, stop=F2B_FQ_DEV, length=2 * F2B_TRIAL_COUNT))
    push!(fq_devs, 0)
    negi_h0s = [(FQ + FQ * fq_dev) * NEGI_H0_ISO for fq_dev in fq_devs]
    gate_errors = ones(pulse_type_count, 2 * F2B_TRIAL_COUNT + 1, F2B_AVG_COUNT)
    save_file_paths = Array{String, 1}(undef, pulse_type_count)

    # check for previous computation
    if use_previous
        data_file_path_old = latest_file_path("h5", "f2b", SAVE_PATH)
    else
        data_file_path_old = nothing
    end
    if isnothing(data_file_path_old)
        save_file_paths_old = nothing
    else
        (save_file_paths_old, gate_errors_old, avg_count_old, pulse_type_count_old
         ) = h5open(data_file_path_old, "r") do data_file_old
             save_file_paths_old = read(data_file_old, "save_file_paths")
             gate_errors_old = read(data_file_old, "gate_errors")
             pulse_type_count_old = size(gate_errors_old)[1]
             avg_count_old = read(data_file_old, "avg_count")
             fq_devs_old = read(data_file_old, "fq_devs")
             if fq_devs_old != fq_devs
                 save_file_paths_old = nothing
             end
            return (save_file_paths_old, gate_errors_old, avg_count_old, pulse_type_count_old)
        end
    end

    for (i, pulse_type) in enumerate(F2B_PT_LIST)
        println("pt[$(i)]: $(pulse_type)")
        save_file_path = data[pulse_type]
        save_file_paths[i] = save_file_path
        save_file_path_old = ((isnothing(save_file_paths_old) || i > pulse_type_count_old)
                              ? nothing : save_file_paths_old[i])
        dynamics_type = pulse_type == analytic ? dynamics_type_analytic : schroed
        for j = 1:2 * F2B_TRIAL_COUNT + 1
            negi_h0 = negi_h0s[j]
            for k = 1:F2B_AVG_COUNT
                # skip redundant computation
                if (!isnothing(save_file_path_old) && save_file_path == save_file_path_old
                    && k <= avg_count_old)
                    gate_errors[i, j, k] = gate_errors_old[i, j, k]
                else
                    res = run_sim_prop(
                        1, gate_type; save_file_path=save_file_path,
                        dynamics_type=dynamics_type, save=false, negi_h0=negi_h0, seed=k
                    )
                    gate_errors[i, j, k] = 1 - res["fidelities"][end]
                end
            end
        end
    end

    data_file_path = generate_file_path("h5", "f2b", SAVE_PATH)
    h5open(data_file_path, "w") do data_file
        write(data_file, "pulse_types", pulse_types_integer)
        write(data_file, "gate_errors", gate_errors)
        write(data_file, "save_file_paths", save_file_paths)
        write(data_file, "avg_count", F2B_AVG_COUNT)
        write(data_file, "fq_devs", fq_devs)
    end
    println("Saved f2b data to $(data_file_path)")
end

const F2C_GATE_TIMES_XPIBY2 = [50., 56.8, 60, 70, 80, 90, 100, 110, 120, 130, 140, 150, 160] # 13
const F2C_DATA_XPIBY2 = Dict(
    analytic => [joinpath(SPIN_OUT_PATH, "spin14/$(lpad(index, 5, '0'))_spin14.h5") for index in [
        INVAL, 4, INVAL, INVAL, INVAL, INVAL, INVAL, INVAL, INVAL, INVAL, INVAL, INVAL, INVAL
    ]],
    s2 => [joinpath(SPIN_OUT_PATH, "spin12/$(lpad(index, 5, '0'))_spin12.h5") for index in [
        629, 627, 638, 655, 654, 666, 658, 665, 661, 667, 668, 681, 685
    ]],
    sut8 => [joinpath(SPIN_OUT_PATH, "spin23/$(lpad(index, 5, '0'))_spin23.h5") for index in [
        38, 94, 103, 113, 132, 135, 125, 140, 143, 145, 158, 154, 162
    ]],
    d1 => [joinpath(SPIN_OUT_PATH, "spin11/$(lpad(index, 5, '0'))_spin11.h5") for index in [
        505, 504, 506, 508, 509, 510, 512, 514, 516, 517, 519, 521, 524
    ]],
    d2 => [joinpath(SPIN_OUT_PATH, "spin11/$(lpad(index, 5, '0'))_spin11.h5") for index in [
        520, 528, 522, 523, 525, 526, 527, 530, 531, 533, 535, 536, 537
    ]],
)
const F2C_GATE_TIMES_ZPIBY2 = Array(18:2:72) # 28
const F2C_DATA_ZPIBY2 = Dict(
    analytic => [joinpath(SPIN_OUT_PATH, "spin14/$(lpad(index, 5, '0'))_spin14.h5") for index in [
        0, INVAL, INVAL, INVAL, INVAL, INVAL, INVAL, INVAL, INVAL, INVAL,
        INVAL, INVAL, INVAL, INVAL, INVAL, INVAL, INVAL, INVAL, INVAL, INVAL,
        INVAL, INVAL, INVAL, INVAL, INVAL, INVAL, INVAL, INVAL,
    ]],
    s2 => [joinpath(SPIN_OUT_PATH, "spin12/$(lpad(index, 5, '0'))_spin12.h5") for index in [
        779, 781, 782, 783, 784, 785, 786, 787, 788, 808,
        790, 793, 791, 792, 794, 795, 796, 797, 798, 799,
        800, 801, 802, 803, 805, 804, 806, 807,
    ]],
    sut8 => [joinpath(SPIN_OUT_PATH, "spin30/lanthanum/$(lpad(index, 5, '0'))_spin30.h5") for index in [
        53, 54, 55, 56, 57, 58, 59, 60, 61, 62,
        63, 66, 65, 69, 67, 70, 68, 64, 81, 77,
        71, 73, 74, 78, 75, 76, 79, 82
    ]],
    d1 => [joinpath(SPIN_OUT_PATH, "spin11/$(lpad(index, 5, '0'))_spin11.h5") for index in [
        468, 471, 472, 474, 475, 476, 477, 478, 480, 539,
        541, 542, 543, 544, 545, 551, 553, 554, 556, 569,
        570, 571, 572, 573, 574, 576, 567, 585
    ]],
    d2 => [joinpath(SPIN_OUT_PATH, "spin11/$(lpad(index, 5, '0'))_spin11.h5") for index in [
        487, 488, 489, 491, 492, 493, 494, 498, 499, 547,
        549, 552, 558, 560, 561, 562, 563, 564, 565, 575,
        577, 588, 597, 590, 591, 592, 600, 598
    ]],
)
const F2C_PT_LIST = [analytic, s2, sut8, d1, d2]
const F2C_AVG_COUNT = 1000
const F2C_SIGMA = 1e-2
const F2C_S1_NEGI_H0 = (FQ + FQ * F2C_SIGMA) * NEGI_H0_ISO
const F2C_S2_NEGI_H0 = (FQ - FQ * F2C_SIGMA) * NEGI_H0_ISO
function gen_2c(;use_previous=true, gate_type=zpiby2)
    if gate_type == zpiby2
        gate_times = F2C_GATE_TIMES_ZPIBY2
        data = F2C_DATA_ZPIBY2
        dynamics_type_analytic = zpiby2nodis
    elseif gate_type == xpiby2
        gate_times = F2C_GATE_TIMES_XPIBY2
        data = F2C_DATA_XPIBY2
        dynamics_type_analytic = xpiby2nodis
    end
    pulse_type_count = size(F2C_PT_LIST)[1]
    pulse_types_integer = [Integer(pulse_type) for pulse_type in F2C_PT_LIST]
    gate_time_count = size(gate_times)[1]
    gate_errors = ones(pulse_type_count, gate_time_count, 2 * F2C_AVG_COUNT)
    save_file_paths = Array{String, 2}(undef, pulse_type_count, gate_time_count)

    # check for previous computation
    if use_previous
        data_file_path_old = latest_file_path("h5", "f2c", SAVE_PATH)
    else
        data_file_path_old = nothing
    end
    if isnothing(data_file_path_old)
        save_file_paths_old = nothing
    else
        (save_file_paths_old, gate_errors_old, avg_count_old
         ) = h5open(data_file_path_old, "r") do data_file_old
            save_file_paths_old = read(data_file_old, "save_file_paths")
            gate_errors_old = read(data_file_old, "gate_errors")
            avg_count_old = read(data_file_old, "avg_count")
            return (save_file_paths_old, gate_errors_old, avg_count_old)
         end
        pulse_type_count_old = size(gate_errors_old)[1]
        gate_time_count_old = size(gate_errors_old)[2]
    end

    for (i, pulse_type) in enumerate(F2C_PT_LIST)
        println("pt[$(i)]: $(pulse_type)")
        for (j, gate_time) in enumerate(gate_times)
            print("gt[$(j)]: $(gate_time) ")
            save_file_path = data[pulse_type][j]
            save_file_paths[i, j] = save_file_path
            save_file_path_old = ((isnothing(save_file_paths_old) ||
                                   i > pulse_type_count_old || j > gate_time_count_old)
                                  ? nothing : save_file_paths_old[i, j])
            dynamics_type = pulse_type == analytic ? dynamics_type_analytic : schroed
            if !isnothing(findfirst("$(INVAL)", save_file_path))
                println("INVAL")
                continue
            end
            for k = 1:F2C_AVG_COUNT
                # skip redundant computation
                if (!isnothing(save_file_path_old) && save_file_path == save_file_path_old
                    && k <= avg_count_old)
                    gate_errors[i, j, (k - 1) * 2 + 1] = gate_errors_old[i, j, (k - 1) * 2 + 1]
                    gate_errors[i, j, (k - 1) * 2 + 2] = gate_errors_old[i, j, (k - 1) * 2 + 2]
                else
                    res1 = run_sim_prop(
                        1, gate_type; save_file_path=save_file_path, dynamics_type=dynamics_type,
                        negi_h0=F2C_S1_NEGI_H0, save=false, seed=k
                    )
                    res2 = run_sim_prop(
                        1, gate_type; save_file_path=save_file_path, dynamics_type=dynamics_type,
                        negi_h0=F2C_S2_NEGI_H0, save=false, seed=k
                    )
                    ge1 = 1 - res1["fidelities"][end]
                    ge2 = 1 - res2["fidelities"][end]
                    gate_errors[i, j, (k - 1) * 2 + 1] = ge1
                    gate_errors[i, j, (k - 1) * 2 + 2] = ge2
                end
            end
            println(" $(mean(gate_errors[i, j, :]))")
        end
    end

    data_file_path = generate_file_path("h5", "f2c", SAVE_PATH)
    h5open(data_file_path, "w") do data_file
        write(data_file, "gate_errors", gate_errors)
        write(data_file, "gate_times", gate_times)
        write(data_file, "pulse_types", pulse_types_integer)
        write(data_file, "save_file_paths", save_file_paths)
        write(data_file, "avg_count", F2C_AVG_COUNT)
    end
    print("Saved f2c data to $(data_file_path)")
end

### FIGURE 3 ###
const F3_DATA = Dict(
    analytic => joinpath(SPIN_OUT_PATH, "spin14/00004_spin14.h5"),
    s2 => joinpath(SPIN_OUT_PATH, "spin18/00059_spin18.h5"),
    sut8 => joinpath(SPIN_OUT_PATH, "spin25/lanthanum/00000_spin25.h5"),
    d1 => joinpath(SPIN_OUT_PATH, "spin17/00083_spin17.h5"),
    d2 => joinpath(SPIN_OUT_PATH, "spin17/00085_spin17.h5"),
)

const F3A_PT_LIST = [analytic, s2, sut8, d1, d2]
function gen_3a()
    pulse_types_integer = [Integer(pulse_type) for pulse_type in F3A_PT_LIST]
    pulse_type_count = size(F3A_PT_LIST)[1]
    save_file_paths = Array{String, 1}(undef, pulse_type_count)
    for (i, pulse_type) in enumerate(F3A_PT_LIST)
        save_file_paths[i] = F3_DATA[pulse_type]
    end

    data_file_path = generate_file_path("h5", "f3a", SAVE_PATH)
    h5open(data_file_path, "w") do data_file
        write(data_file, "save_file_paths", save_file_paths)
        write(data_file, "pulse_types", pulse_types_integer)
    end
    println("Saved f3a data to $(data_file_path)")
end

const F3B_PT_LIST = [analytic, s2, sut8, d1, d2]
const F3B_AVG_COUNT = 1000
const F3B_GATE_COUNT = 200
function gen_3b(;use_previous=true)
    gate_type = xpiby2
    pulse_types_integer = [Integer(pulse_type) for pulse_type in F3B_PT_LIST]
    pulse_type_count = size(F3B_PT_LIST)[1]
    save_file_paths = Array{String, 1}(undef, pulse_type_count)
    gate_errors = ones(pulse_type_count, F3B_GATE_COUNT + 1, F3B_AVG_COUNT)
    states = ones(pulse_type_count, F3B_GATE_COUNT + 1, F3B_AVG_COUNT, HDIM_ISO)

    # check for previous computation
    if use_previous
        data_file_path_old = latest_file_path("h5", "f3b", SAVE_PATH)
    else
        data_file_path_old = nothing
    end
    if !isnothing(data_file_path_old)
        (save_file_paths_old, gate_errors_old, avg_count_old, states_old
         ) = h5open(data_file_path_old, "r") do data_file_old
             gate_count_old = read(data_file_old, "gate_count")
             pulse_types_integer_old = read(data_file_old, "pulse_types")
             if gate_count_old == F3B_GATE_COUNT && pulse_types_integer_old == pulse_types_integer
                 save_file_paths_old = read(data_file_old, "save_file_paths")
                 gate_errors_old = read(data_file_old, "gate_errors")
                 avg_count_old = read(data_file_old, "avg_count")
                 states_old = read(data_file_old, "states")
             else
                 save_file_paths_old = gate_errors_old = avg_count_old = states_old = nothing
             end
             return (save_file_paths_old, gate_errors_old, avg_count_old, states_old)
         end
    else
        save_file_paths_old = gate_errors_old = avg_count_old = states_old = nothing
    end

    for (i, pulse_type) in enumerate(F3B_PT_LIST)
        print("pt[$(i)]: $(pulse_type) ")
        save_file_path = F3_DATA[pulse_type]
        save_file_paths[i] = save_file_path
        if isempty(save_file_path)
            println("")
            continue
        end
        save_file_path_old = isnothing(save_file_paths_old) ? nothing : save_file_paths_old[i]
        save_file_path_sim = pulse_type == analytic ? nothing : save_file_path
        dynamics_type = pulse_type == analytic ? xpiby2da : schroedda
        for j = 1:F3B_AVG_COUNT
            # avoid redundant computation
            if (!isnothing(save_file_path_old) && save_file_path == save_file_path_old
                && j <= avg_count_old)
                gate_errors[i, :, j] = gate_errors_old[i, :, j]
                states[i, :, j, :] = states_old[i, :, j, :]
            else
                res = run_sim_prop(F3B_GATE_COUNT, gate_type; dynamics_type=dynamics_type,
                                   save_file_path=save_file_path_sim, seed=j, save=false)
                gate_errors[i, :, j] = 1 .- res["fidelities"]
                states[i, :, j, :] = res["states"]
            end
        end
        println("")
    end

    data_file_path = generate_file_path("h5", "f3b", SAVE_PATH)
    h5open(data_file_path, "w") do data_file
        write(data_file, "save_file_paths", save_file_paths)
        write(data_file, "pulse_types", pulse_types_integer)
        write(data_file, "gate_errors", gate_errors)
        write(data_file, "states", states)
        write(data_file, "avg_count", F3B_AVG_COUNT)
        write(data_file, "gate_count", F3B_GATE_COUNT)
    end
    println("Saved f3b data to $(data_file_path)")
end

"""
gen_3c - generate an average of the flux noise used in fig3b,
see rbqoc/src/spin/spin.jl/run_sim_deqjl for noise construction
"""
function gen_3c()
    noise_dt_inv = DT_NOISE_INV
    noise_dist = STD_NORMAL
    noise_amp = NAMP_PREFACTOR
    evolution_time = 56.8 * 5e3
    noise_count = Int(ceil(evolution_time * noise_dt_inv)) + 1

    # generate noise
    noise = pink_noise_from_white(noise_count, noise_dt_inv, noise_dist; seed=0)
    for seed = 1:9
        noise .+= pink_noise_from_white(noise_count, noise_dt_inv, noise_dist; seed=seed)
    end
    # average
    noise ./= 10

    # take fft
    noise_fft = Array{Complex{Float64}, 1}(noise)
    fft!(noise_fft)
    # normalize and take absolute value
    for i = 1:size(noise_fft)[1]
        noise_fft[i] = abs(noise_fft[i]) / noise_count
    end
    noise_fft = Array{Float64, 1}(noise_fft)

    # heh, linearity
    noise .*= NAMP_PREFACTOR
    noise_fft .*= NAMP_PREFACTOR

    times = Array(0:(noise_count - 1)) / noise_dt_inv
    freqs = Array(fftfreq(noise_count, noise_dt_inv))

    h5open(F3C_DATA_FILE_PATH, "w") do data_file
        write(data_file, "noise", noise)
        write(data_file, "times", times)
        write(data_file, "noise_fft", noise_fft)
        write(data_file, "freqs", freqs)
    end
end

const F3D_SAVE_STEP = 1e-1
const F3D_PTS = [analytic, s2, s4, d2, d3]
"""
gen_3d - used for generating expectation values of an evolution
to plot on bloch sphere
"""
function gen_3d()
    pulse_types_int = [Integer(pulse_type) for pulse_type in F3D_PTS]
    point_count = Integer(ceil(56.80 * F3D_SAVE_STEP^(-1)))
    pointss = zeros(size(F3D_PTS)[1], point_count, 3)
    for (i, pulse_type) in enumerate(F3D_PTS)
        data = F3_DATA[pulse_type]
        if !(DATAFP_KEY in keys(data))
            if pulse_type == analytic
                data_file_path = run_sim_fine_deqjl(
                    dynamics_type=xpiby2nodis, dt=1e-3,
                    save_step=F3D_SAVE_STEP)
            else
                data_file_path = run_sim_fine_deqjl(
                    save_file_path=data[SAVEFP_KEY],
                    save_type=data[SAVET_KEY],
                    dt=1e-3, save_step=F3D_SAVE_STEP,
                )
            end
        else
            data_file_path = data[DATAFP_KEY]
        end

        (states_) = h5open(data_file_path, "r") do data_file
            states_ = read(data_file, "states")
            return (states_)
        end
        for j = 1:point_count
            state = get_vec_uniso(states_[j, :])
            state_dag = state'
            pointss[i, j, 1] = Real(state_dag * SIGMAX * state)
            pointss[i, j, 2] = Real(state_dag * SIGMAY * state)
            pointss[i, j, 3] = Real(state_dag * SIGMAZ * state)
        end
    end
    h5open(F3D_DATA_FILE_PATH, "w") do data_file
        write(data_file, "pointss", pointss)
        write(data_file, "pulse_types", pulse_types_int)
    end
end

const F4_DATA = Dict(
    corpse => joinpath(SPIN_OUT_PATH, "spin14/00124_spin14.h5"),
    s2 => nothing,
    s4 => nothing,
    d1 => joinpath(SPIN_OUT_PATH, "spin19/00025_spin19.h5"),
    d2 => joinpath(SPIN_OUT_PATH, "spin19/00026_spin19.h5"),
)

const F4A_PTS = [corpse, s2, s4, d1, d2]
function gen_4a()
    pulse_types_integer = [Integer(pulse_type) for pulse_type in F4A_PTS]
    pulse_type_count = size(F4A_PTS)[1]
    save_file_paths = Array{String, 1}(undef, pulse_type_count)
    for (i, pulse_type) in enumerate(F4A_PTS)
        save_file_path = F4_DATA[pulse_type]
        save_file_paths[i] = isnothing(save_file_path) ? "" : save_file_path
    end

    data_file_path = generate_file_path("h5", "f4a", SAVE_PATH)
    h5open(data_file_path, "w") do data_file
        write(data_file, "save_file_paths", save_file_paths)
        write(data_file, "pulse_types", pulse_types_integer)
    end
    println("Saved f4a data to $(data_file_path)")
end

const F4B_PTS = [corpse, s2, s4, d1, d2]
const F4B_TRIAL_COUNT = Integer(1e2)
const F4B_FQ_DEV = 2e-2
const F4B_AVG_COUNT = 10
const F4B_DT_INV = Int64(1e4) # go down to 1e5 if want ge < 10^-10
function gen_4b(;use_previous=true)
    gate_type = xpi
    pulse_types_integer = [Integer(pt) for pt in F4B_PTS]
    pulse_type_count = size(F4B_PTS)[1]
    fq_devs = Array(range(-F4B_FQ_DEV, stop=F4B_FQ_DEV, length=2 * F4B_TRIAL_COUNT))
    push!(fq_devs, 0)
    negi_h0s = [fq_dev * NEGI_H0_ISO for fq_dev in fq_devs]
    gate_errors = ones(pulse_type_count, 2 * F4B_TRIAL_COUNT + 1, F4B_AVG_COUNT)
    save_file_paths = Array{String, 1}(undef, pulse_type_count)

    # check for previous computation
    if use_previous
        data_file_path_old = latest_file_path("h5", "f4b", SAVE_PATH)
    else
        data_file_path_old = nothing
    end
    if isnothing(data_file_path_old)
        save_file_paths_old = nothing
    else
        (save_file_paths_old, gate_errors_old, avg_count_old, pulse_type_count_old
         ) = h5open(data_file_path_old, "r") do data_file_old
             save_file_paths_old = read(data_file_old, "save_file_paths")
             gate_errors_old = read(data_file_old, "gate_errors")
             pulse_type_count_old = size(gate_errors_old)[1]
             avg_count_old = read(data_file_old, "avg_count")
             fq_devs_old = read(data_file_old, "fq_devs")
             if fq_devs_old != fq_devs
                 save_file_paths_old = nothing
             end
            return (save_file_paths_old, gate_errors_old, avg_count_old, pulse_type_count_old)
        end
    end

    for (i, pulse_type) in enumerate(F4B_PTS)
        println("pt[$(i)]: $(pulse_type)")
        save_file_path = F4_DATA[pulse_type]
        save_file_paths[i] = isnothing(save_file_path) ? "" : save_file_path
        if isnothing(save_file_path)
            continue
        end
        save_file_path_old = ((isnothing(save_file_paths_old) || i > pulse_type_count_old)
                              ? nothing : save_file_paths_old[i])
        save_file_path_sim = pulse_type == corpse ? nothing : save_file_path
        if pulse_type == corpse
            dynamics_type = xpicorpse
        else
            dynamics_type = schroed
        end

        for j = 1:2 * F4B_TRIAL_COUNT + 1
            negi_h0 = negi_h0s[j]
            for k = 1:F4B_AVG_COUNT
                # skip redundant computation
                if (!isnothing(save_file_path_old) && save_file_path == save_file_path_old
                    && k <= avg_count_old)
                    gate_errors[i, j, k] = gate_errors_old[i, j, k]
                else
                    res = run_sim_deqjl(
                        1, gate_type; dynamics_type=dynamics_type,
                        save_file_path=save_file_path_sim,
                        dt_inv=F4B_DT_INV, save=false, negi_h0=negi_h0, seed=k
                    )
                    gate_errors[i, j, k] = 1 - res["fidelities"][end]
                end
            end
        end
    end

    data_file_path = generate_file_path("h5", "f4b", SAVE_PATH)
    h5open(data_file_path, "w") do data_file
        write(data_file, "pulse_types", pulse_types_integer)
        write(data_file, "gate_errors", gate_errors)
        write(data_file, "save_file_paths", save_file_paths)
        write(data_file, "avg_count", F4B_AVG_COUNT)
        write(data_file, "fq_devs", fq_devs)
    end
    println("Saved f4b data to $(data_file_path)")
end

"""
generate all required figure data
"""
function gen_all(;use_previous=true)
    gen_1a()
    gen_1b()
    gen_1c(;use_previous=use_previous)
    gen_2a()
    gen_2b(;use_previous=use_previous)
    gen_2c(;use_previous=use_previous)
    gen_3a()
    gen_3b(;use_previous=use_previous)
end

# If you run `julia figures.jl` all of the figure data will be
# generated from scratch.
if String(@__FILE__) == String(joinpath(pwd(), PROGRAM_FILE))
    gen_all(;use_previous=false)
end

ext: jl

max_stars_repo_name: SchusterLab/rbqoc

max_issues_repo_name: SchusterLab/rbqoc

max_forks_repo_name: SchusterLab/rbqoc


text:
include("rice2010.jl")
using .Rice2010

function getmarginal_rice_models(;emissionyear=2005,datafile=joinpath(@__DIR__, "..", "data", "RICE_2010_base_000.xlsm"))

    RICE = getrice()
    run(RICE)

    mm = MarginalModel(RICE)
    m1 = mm.base
    m2 = mm.marginal

    add_comp!(m2, adder, :marginalemission, before=:co2cycle)

    time = dimension(m1, :time)
    addem = zeros(length(time))
    addem[time[emissionyear]] = 1.0

    set_param!(m2,:marginalemission,:add,addem)
    connect_param!(m2,:marginalemission,:input,:emissions,:E)
    connect_param!(m2, :co2cycle,:E,:marginalemission,:output)

    run(m1)
    run(m2)

    return m1, m2
end

# This function returns a matrix of marginal damages per one t of carbon emission in the
# emissionyear parameter year.
function getmarginaldamages_rice(;emissionyear=2005,datafile=joinpath(@__DIR__, "..", "data", "RICE_2010_base_000.xlsm"))
    m1, m2 = getmarginal_rice_models(emissionyear=emissionyear, datafile=datafile)

    run(m1)
    run(m2)

    damage1 = m1[:damages,:DAMAGES]
    damage2 = m2[:damages,:DAMAGES]

    # Calculate the marginal damage between run 1 and 2 for each
    # year/region
    marginaldamage = (damage2.-damage1) .* 10^12 / 10^9 / 10

    return marginaldamage
end

# Old Way

#= function getmarginal_rice_models(;emissionyear=2005,datafile=joinpath(@__DIR__, "..", "data", "RICE_2010_base_000.xlsm"))
    m1 = getrice(datafile=datafile)

    m2 = getrice(datafile=datafile)
    add_comp!(m2, adder, :marginalemission, before=:co2cycle)
    addem = zeros(60)
    addem[getindexfromyear_rice_2010(emissionyear)] = 1.0
    setparameter(m2,:marginalemission,:add,addem)
    connectparameter(m2,:marginalemission,:input,:emissions,:E)
    connectparameter(m2, :co2cycle,:E,:marginalemission,:output)

    run(m1)
    run(m2)

    return m1, m2
end

# This function returns a matrix of marginal damages per one t of carbon emission in the
# emissionyear parameter year.
function getmarginaldamages_rice(;emissionyear=2005,datafile=joinpath(@__DIR__, "..", "data", "RICE_2010_base_000.xlsm"))
    m1, m2 = getmarginal_rice_models(emissionyear=emissionyear, datafile=datafile)

    run(m1)
    run(m2)

    damage1 = m1[:damages,:DAMAGES]
    damage2 = m2[:damages,:DAMAGES]

    # Calculate the marginal damage between run 1 and 2 for each
    # year/region
    marginaldamage = (damage2.-damage1) .* 10^12 / 10^9 / 10

    return marginaldamage
end =#


text:
@parameters t

# New simplification rules
exponent_rules = [
    @rule(exp(log(~x)) => ~x),
    @rule(log(exp(~x)) => ~x),
    @acrule(exp(~x + ~y) => exp(~x) * exp(~y)),
    @acrule(exp(~x * ~y) => exp(~y)^~x),
    @acrule(log(~x) + log(~y) => log(~x * ~y)),
    @rule(log((~x)^(~a)) => ~a * log(~x))
]
rw_exp = RestartedChain(exponent_rules)
rw_chain = RestartedChain([SymbolicUtils.default_simplifier(), rw_exp])
rewriter = Postwalk(rw_chain)

# Constitutive relations
constitutive_relations(n::AbstractNode) = equations(n)
function constitutive_relations(n::EqualEffort)
    if all(freeports(n)) # all ports are empty
        return Equation[]
    end

    N = numports(n)
    @variables E[1:N](t) F[1:N](t)

    flow_constraint = [0 ~ sum(collect(F))]
    effort_constraints = [0 ~ E[1] - e for e in E[2:end]]
    return vcat(flow_constraint, effort_constraints)
end
function constitutive_relations(n::EqualFlow)
    if all(freeports(n)) # all ports are empty
        return Equation[]
    end

    N = numports(n)
    @variables E[1:N](t) F[1:N](t)

    W = weights(n)
    weighted_e = W .* collect(E)
    weighted_f = W .* collect(F)

    effort_constraint = [0 ~ sum(weighted_e)]
    flow_constraints = [0 ~ weighted_f[1] - f for f in weighted_f[2:end]]
    return vcat(effort_constraint, flow_constraints)
end
function constitutive_relations(m::BondGraph)
    return simplify.(full_equations(ODESystem(m)); expand = true, rewriter = rewriter)
end
function constitutive_relations(bgn::BondGraphNode)
    return constitutive_relations(bgn.bondgraph)
end

@connector function MTKPort(; name)
    vars = @variables E(t) F(t) [connect = Flow]
    ODESystem(Equation[], t, vars, []; name = name)
end

function get_connection_eq(b::Bond,subsystems)
    (s,d) = b
    src_port = getproperty(subsystems[s.node], Symbol("p$(s.index)"))
    dst_port = getproperty(subsystems[d.node], Symbol("p$(d.index)"))
    connect(src_port, dst_port)
end

function ModelingToolkit.ODESystem(n::AbstractNode)
    N = numports(n)
    ps = [MTKPort(name = Symbol("p$i")) for i in 1:N]

    @variables E[1:N](t) F[1:N](t)
    e_sub_rules = Dict(E[i] => ps[i].E for i in 1:N)
    f_sub_rules = Dict(F[i] => ps[i].F for i in 1:N)
    sub_rules = merge(e_sub_rules, f_sub_rules)
    eqs = Equation[substitute(eq, sub_rules) for eq in constitutive_relations(n)]

    sys = ODESystem(eqs, t, states(n), parameters(n);
        name = n.name, defaults = defaults(n))
    return compose(sys, ps...)
end
function ModelingToolkit.ODESystem(m::BondGraph; simplify_eqs = true)
    (subsystems, connections) = get_subsys_and_connections(m)
    compose_bg_model(subsystems, connections, m.name, simplify_eqs)
end
function ModelingToolkit.ODESystem(bgn::BondGraphNode; simplify_eqs = false)
    N = numports(bgn)
    ps = [MTKPort(name = Symbol("p$i")) for i in 1:N]

    (subsystems, connections) = get_subsys_and_connections(bgn.bondgraph)
    es = [subsystems[c].p1.E for c in exposed(bgn)]
    fs = [subsystems[c].p1.F for c in exposed(bgn)]
    port_eqs = [
        [0 ~ p.E - E for (E,p) in zip(es,ps)];
        [0 ~ p.F + F for (F,p) in zip(fs,ps)]
    ]
    eqs = [connections; port_eqs]
    sys = compose_bg_model(subsystems, eqs, bgn.name, simplify_eqs)
    compose(sys, ps...)
end

function get_subsys_and_connections(bg::BondGraph)
    # Collect constitutive relations from components
    subsystems = OrderedDict(n => ODESystem(n) for n in nodes(bg))
    # Add constraints from bonds
    connections = [get_connection_eq(b,subsystems) for b in bonds(bg)]
    return subsystems, connections
end

function compose_bg_model(subsystems, eqs, name, simplify_eqs)
    @named _model = ODESystem(eqs, t; name = name)
    model = compose(_model, collect(values(subsystems)))

    if simplify_eqs
        return structural_simplify(model)
    else
        return model
    end
end

function simulate(m::BondGraph, tspan; u0 = [], pmap = [], probtype::Symbol = :default, kwargs...)
    sys = ODESystem(m)
    flag_ODE = !any([isequal(eq.lhs, 0) for eq in ModelingToolkit.equations(sys)])
    if probtype == :ODE || flag_ODE
        prob = ODEProblem(sys, u0, tspan, pmap)
    else
        prob = ODAEProblem(sys, u0, tspan, pmap)
    end
    return solve(prob; kwargs...)
end

max_stars_repo_name: mic-pan/BondGraphs

max_issues_repo_name: mic-pan/BondGraphs

max_forks_repo_name: mic-pan/BondGraphs


text:
#****************************************************************************
#                     Molecular Dynamics Potentials (MDP)
#                            CESMIX-MIT Project
#
#  Contributing authors: Ngoc-Cuong Nguyen (cuongng@mit.edu, exapde@gmail.com)
#****************************************************************************

__precompile__()

module Gencode

using Revise, SymPy

#export syminit, gencode, compilecode
export syminit, gencode, compile, runcode, tring2cmd

#include("initializeapp.jl");
include("syminit.jl");
include("varsassign.jl");
include("symsassign.jl");
include("genpotnum.jl");
include("getpotstr.jl");
include("cppfiles.jl");
include("contains.jl");
include("gendensity.jl");
include("gensingle.jl");
include("genpair.jl");
include("gentriplet.jl");
include("genquadruplet.jl");
include("genbo2.jl");
include("genbo3.jl");
include("genpotential.jl");
include("nopotential.jl");
include("string2cmd.jl");
include("compile.jl");
include("runcode.jl");

include("derivassign.jl");
include("genpotentialgrad.jl");
include("nopotentialgrad.jl");
include("gensinglegrad.jl");
include("genpairgrad.jl");
include("gentripletgrad.jl");
include("genquadrupletgrad.jl");
include("gendensitygrad.jl");
include("genbo2grad.jl");
include("genbo3grad.jl");

function gencode(app)

print("generate code...\n");
if !isdir("app")
    mkdir("app");
end

potential = getfield(Main, Symbol("Main"))

xij, xik, xil, xi, xj, xk, xl, qi, qj, qk, ql, ti, tj, tk, tl, ai, aj, ak, al = syminit(app);
eta = [SymPy.symbols("eta$i") for i=1:app.nceta];
kappa = [SymPy.symbols("kappa$i") for i=1:app.nckappa];

foldername = app.sourcepath * "C++/Potentials";
#foldername = app.sourcepath * "Julia/Gencode/app";

potform = 1;
filename = "Singlea";
mu = [SymPy.symbols("mu$i") for i=1:app.ncmu1a];
if isdefined(potential, Symbol(filename))
    u = potential.Singlea(xi, qi, ti, mu, eta, kappa);
    genpotential(foldername, filename, u, xij, xik, xil, xi, xj, xk, xl, qi, qj, qk, ql, ti, tj, tk, tl, ai, aj, ak, al, mu, eta, kappa, potform);
    genpotentialgrad(foldername, filename * "Gradient", u, xij, xik, xil, xi, xj, xk, xl, qi, qj, qk, ql, ti, tj, tk, tl, ai, aj, ak, al, mu, eta, kappa, potform);
else
    nopotential(foldername, filename, [], xij, xik, xil, xi, xj, xk, xl, qi, qj, qk, ql, ti, tj, tk, tl, ai, aj, ak, al, mu, eta, kappa, potform);
    genpotentialgrad(foldername, filename * "Gradient", [], xij, xik, xil, xi, xj, xk, xl, qi, qj, qk, ql, ti, tj, tk, tl, ai, aj, ak, al, mu, eta, kappa, potform);
end

filename = "Singleb";
mu = [SymPy.symbols("mu$i") for i=1:app.ncmu1b];
if isdefined(potential, Symbol(filename))
    u = potential.Singleb(xi, qi, ti, mu, eta, kappa);
    genpotential(foldername, filename, u, xij, xik, xil, xi, xj, xk, xl, qi, qj, qk, ql, ti, tj, tk, tl, ai, aj, ak, al, mu, eta, kappa, potform);
    genpotentialgrad(foldername, filename * "Gradient", u, xij, xik, xil, xi, xj, xk, xl, qi, qj, qk, ql, ti, tj, tk, tl, ai, aj, ak, al, mu, eta, kappa, potform);
else
    nopotential(foldername, filename, [], xij, xik, xil, xi, xj, xk, xl, qi, qj, qk, ql, ti, tj, tk, tl, ai, aj, ak, al, mu, eta, kappa, potform);
    genpotentialgrad(foldername, filename * "Gradient", [], xij, xik, xil, xi, xj, xk, xl, qi, qj, qk, ql, ti, tj, tk, tl, ai, aj, ak, al, mu, eta, kappa, potform);
end

potform = 2;
filename = "Paira";
mu = [SymPy.symbols("mu$i") for i=1:app.ncmu2a];
if isdefined(potential, Symbol(filename))
    u = potential.Paira(xij, qi, qj, ti, tj, mu, eta, kappa);
    genpotential(foldername, filename, u, xij, xik, xil, xi, xj, xk, xl, qi, qj, qk, ql, ti, tj, tk, tl, ai, aj, ak, al, mu, eta, kappa, potform);
    genpotentialgrad(foldername, filename * "Gradient", u, xij, xik, xil, xi, xj, xk, xl, qi, qj, qk, ql, ti, tj, tk, tl, ai, aj, ak, al, mu, eta, kappa, potform);
else
    nopotential(foldername, filename, [], xij, xik, xil, xi, xj, xk, xl, qi, qj, qk, ql, ti, tj, tk, tl, ai, aj, ak, al, mu, eta, kappa, potform);
    genpotentialgrad(foldername, filename * "Gradient", [], xij, xik, xil, xi, xj, xk, xl, qi, qj, qk, ql, ti, tj, tk, tl, ai, aj, ak, al, mu, eta, kappa, potform);
end

filename = "Pairb";
mu = [SymPy.symbols("mu$i") for i=1:app.ncmu2b];
if isdefined(potential, Symbol(filename))
    u = potential.Pairb(xij, qi, qj, ti, tj, mu, eta, kappa);
    genpotential(foldername, filename, u, xij, xik, xil, xi, xj, xk, xl, qi, qj, qk, ql, ti, tj, tk, tl, ai, aj, ak, al, mu, eta, kappa, potform);
    genpotentialgrad(foldername, filename * "Gradient", u, xij, xik, xil, xi, xj, xk, xl, qi, qj, qk, ql, ti, tj, tk, tl, ai, aj, ak, al, mu, eta, kappa, potform);
else
    nopotential(foldername, filename, [], xij, xik, xil, xi, xj, xk, xl, qi, qj, qk, ql, ti, tj, tk, tl, ai, aj, ak, al, mu, eta, kappa, potform);
    genpotentialgrad(foldername, filename * "Gradient", [], xij, xik, xil, xi, xj, xk, xl, qi, qj, qk, ql, ti, tj, tk, tl, ai, aj, ak, al, mu, eta, kappa, potform);
end

filename = "Pairc";
mu = [SymPy.symbols("mu$i") for i=1:app.ncmu2c];
rho = [SymPy.symbols("rho$i") for i=1:1];
if isdefined(potential, Symbol(filename))
    u = potential.Pairc(xij, qi, qj, ti, tj, rho, mu, eta, kappa);
    genbo2(foldername, filename, u, xij, qi, qj, ti, tj, ai, aj, rho, mu, eta, kappa, app.natomtype);
    genbo2grad(foldername, filename * "Gradient", u, xij, qi, qj, ti, tj, ai, aj, rho, mu, eta, kappa, app.natomtype);
    #function genbo2grad(foldername, filename, u, xij, qi, qj, ti, tj, ai, aj, rho, mu, eta, kappa, natomtype)
else
    ifile = 0;
    gen = 0;
    stropu, strcpu, strgpu = genpair(foldername, filename, [], xij, qi, qj, ti, tj, ai, aj, mu, eta, kappa, gen, ifile);
    tmpopu, tmpcpu, tmpgpu = gendensity(foldername, "PaircDensity", [], mu, mu, eta, kappa, gen, ifile);
    stropu = stropu  * "\n" * tmpopu * "\n";
    strcpu = strcpu  * "\n" * tmpcpu * "\n";
    strgpu = strgpu  * "\n" * tmpgpu * "\n";
    cppfiles(foldername, filename, stropu, strcpu, strgpu, 1);
    genbo2grad(foldername, filename * "Gradient", [], xij, qi, qj, ti, tj, ai, aj, rho, mu, eta, kappa, app.natomtype);
end

potform = 3;
filename = "Tripleta";
mu = [SymPy.symbols("mu$i") for i=1:app.ncmu3a];
if isdefined(potential, Symbol(filename))
    u = potential.Tripleta(xij, xik, qi, qj, qk, ti, tj, tk, mu, eta, kappa);
    genpotential(foldername, filename, u, xij, xik, xil, xi, xj, xk, xl, qi, qj, qk, ql, ti, tj, tk, tl, ai, aj, ak, al, mu, eta, kappa, potform);
    genpotentialgrad(foldername, filename * "Gradient", u, xij, xik, xil, xi, xj, xk, xl, qi, qj, qk, ql, ti, tj, tk, tl, ai, aj, ak, al, mu, eta, kappa, potform);
else
    nopotential(foldername, filename, [], xij, xik, xil, xi, xj, xk, xl, qi, qj, qk, ql, ti, tj, tk, tl, ai, aj, ak, al, mu, eta, kappa, potform);
    genpotentialgrad(foldername, filename * "Gradient", [], xij, xik, xil, xi, xj, xk, xl, qi, qj, qk, ql, ti, tj, tk, tl, ai, aj, ak, al, mu, eta, kappa, potform);
end

filename = "Tripletb";
mu = [SymPy.symbols("mu$i") for i=1:app.ncmu3b];
if isdefined(potential, Symbol(filename))
    u = potential.Tripletb(xij, xik, qi, qj, qk, ti, tj, tk, mu, eta, kappa);
    genpotential(foldername, filename, u, xij, xik, xil, xi, xj, xk, xl, qi, qj, qk, ql, ti, tj, tk, tl, ai, aj, ak, al, mu, eta, kappa, potform);
    genpotentialgrad(foldername, filename * "Gradient", u, xij, xik, xil, xi, xj, xk, xl, qi, qj, qk, ql, ti, tj, tk, tl, ai, aj, ak, al, mu, eta, kappa, potform);
else
    nopotential(foldername, filename, [], xij, xik, xil, xi, xj, xk, xl, qi, qj, qk, ql, ti, tj, tk, tl, ai, aj, ak, al, mu, eta, kappa, potform);
    genpotentialgrad(foldername, filename * "Gradient", [], xij, xik, xil, xi, xj, xk, xl, qi, qj, qk, ql, ti, tj, tk, tl, ai, aj, ak, al, mu, eta, kappa, potform);
end

filename = "Tripletc";
mu = [SymPy.symbols("mu$i") for i=1:app.ncmu3c];
rho = [SymPy.symbols("rho$i") for i=1:1];
if isdefined(potential, Symbol(filename))
    u = potential.Tripletc(xij, xik, qi, qj, qk, ti, tj, tk, rho, mu, eta, kappa);
    genbo3(foldername, filename, u, xij, xik, qi, qj, qk, ti, tj, tk, ai, aj, ak, rho, mu, eta, kappa, app.natomtype);
    genbo3grad(foldername, filename * "Gradient", u, xij, xik, qi, qj, qk, ti, tj, tk, ai, aj, ak, rho, mu, eta, kappa, app.natomtype);
else
    ifile = 0;
    gen = 0;
    fp = "TripletcPair";
    fn = "TripletcDensity";
    stropu, strcpu, strgpu = gentriplet(foldername, filename, [], xij, xik, qi, qj, qk, ti, tj, tk, ai, aj, ak, mu, eta, kappa, gen, ifile);
    tmqopu, tmqcpu, tmqgpu = genpair(foldername, fp, [], xij, qi, qj, ti, tj, ai, aj, mu, eta, kappa, gen, ifile);
    tmpopu, tmpcpu, tmpgpu = gendensity(foldername, fn, [], mu, mu, eta, kappa, gen, ifile);
    stropu = stropu  * "\n" * tmqopu * "\n" * tmpopu * "\n";
    strcpu = strcpu  * "\n" * tmqcpu * "\n" * tmpcpu * "\n";
    strgpu = strgpu  * "\n" * tmqgpu * "\n" * tmpgpu * "\n";
    cppfiles(foldername, filename, stropu, strcpu, strgpu, 1);
    genbo3grad(foldername, filename * "Gradient", [], xij, xik, qi, qj, qk, ti, tj, tk, ai, aj, ak, rho, mu, eta, kappa, app.natomtype);
end

potform = 4;
filename = "Quadrupleta";
mu = [SymPy.symbols("mu$i") for i=1:app.ncmu4a];
if isdefined(potential, Symbol(filename))
    u = potential.Quadrupleta(xij, xik, xil, qi, qj, qk, ql, ti, tj, tk, tl, mu,  eta, kappa);
    genpotential(foldername, filename, u, xij, xik, xil, xi, xj, xk, xl, qi, qj, qk, ql, ti, tj, tk, tl, ai, aj, ak, al, mu, eta, kappa, potform);
    genpotentialgrad(foldername, filename * "Gradient", u, xij, xik, xil, xi, xj, xk, xl, qi, qj, qk, ql, ti, tj, tk, tl, ai, aj, ak, al, mu, eta, kappa, potform);
else
    nopotential(foldername, filename, [], xij, xik, xil, xi, xj, xk, xl, qi, qj, qk, ql, ti, tj, tk, tl, ai, aj, ak, al, mu, eta, kappa, potform);
    genpotentialgrad(foldername, filename * "Gradient", [], xij, xik, xil, xi, xj, xk, xl, qi, qj, qk, ql, ti, tj, tk, tl, ai, aj, ak, al, mu, eta, kappa, potform);
end

filename = "Quadrupletb";
mu = [SymPy.symbols("mu$i") for i=1:app.ncmu4b];
if isdefined(potential, Symbol(filename))
    u = potential.Quadrupletb(xij, xik, xil, qi, qj, qk, ql, ti, tj, tk, tl, mu,  eta, kappa);
    genpotential(foldername, filename, u, xij, xik, xil, xi, xj, xk, xl, qi, qj, qk, ql, ti, tj, tk, tl, ai, aj, ak, al, mu, eta, kappa, potform);
    genpotentialgrad(foldername, filename * "Gradient", u, xij, xik, xil, xi, xj, xk, xl, qi, qj, qk, ql, ti, tj, tk, tl, ai, aj, ak, al, mu, eta, kappa, potform);
else
    nopotential(foldername, filename, [], xij, xik, xil, xi, xj, xk, xl, qi, qj, qk, ql, ti, tj, tk, tl, ai, aj, ak, al, mu, eta, kappa, potform);
    genpotentialgrad(foldername, filename * "Gradient", [], xij, xik, xil, xi, xj, xk, xl, qi, qj, qk, ql, ti, tj, tk, tl, ai, aj, ak, al, mu, eta, kappa, potform);
end

end

end

max_stars_repo_name: cesmix-mit/MDP

max_issues_repo_name: cesmix-mit/MDP

max_forks_repo_name: cesmix-mit/MDP