text:
#
# April 2020 Benefit Rates: CAREFUL
#
sys.it.savings_rates   =  [0.0, 20.0, 40.0, 45.0]
sys.it.savings_thresholds   =  [5_000.0, 37_500.0, 150_000.0]
sys.it.savings_basic_rate = 2 # above this counts as higher rate

sys.it.non_savings_rates  =  [19.0,20.0,21.0,41.0,46.0]
sys.it.non_savings_thresholds  =  [2_085, 12_658, 31_931, 150_000.0]
sys.it.non_savings_basic_rate = 2 # above this counts as higher rate

sys.it.dividend_rates  =  [0.0, 7.5,32.5,38.1]
sys.it.dividend_thresholds  =  [2_000.0, 37_500.0, 150_000.0]
sys.it.dividend_basic_rate = 2 # above this counts as higher rate

sys.it.personal_allowance          = 12_500.00
sys.it.personal_allowance_income_limit = 100_000.00
sys.it.personal_allowance_withdrawal_rate = 50.0
sys.it.blind_persons_allowance     = 2_500.00

sys.it.married_couples_allowance   = 9_075.00
sys.it.mca_minimum                 = 3_510.00
sys.it.mca_income_maximum          = 29_600.00
sys.it.mca_credit_rate             = 10.0
sys.it.mca_withdrawal_rate        = 50.0

sys.it.marriage_allowance          = 1_260.00 # check 20/21
sys.it.personal_savings_allowance  = 1_000.00

# FIXME better to have it straight from
# the book with charges per CO2 range
# and the data being an estimate of CO2 per type
merge( sys.it.company_car_charge_by_CO2_emissions,
    Dict([
        Missing_Fuel_Type=>0.1,
        No_Fuel=>0.1,
        Other=>0.1,
        Dont_know=>0.1,
        Petrol=>0.25,
        Diesel=>0.37,
        Hybrid_use_a_combination_of_petrol_and_electricity=>0.16,
        Electric=>0.02,
        LPG=>0.02,
        Biofuel_eg_E85_fuel=>0.02 ]))

sys.it.fuel_imputation  = 24_500.00 # 20/21

#
# pensions
# CHECK THESE FOR UPRATING
sys.it.pension_contrib_basic_amount = 3_600.00
sys.it.pension_contrib_annual_allowance = 40_000.00
sys.it.pension_contrib_annual_minimum = 4_000.00
sys.it.pension_contrib_threshold_income = 240_000.00
sys.it.pension_contrib_withdrawal_rate = 50.0

# sys.it.non_savings_income = NON_SAVINGS_INCOME
# sys.it.all_taxable = ALL_TAXABLE_INCOME
# sys.it.savings_income = SAVINGS_INCOME
# sys.it.dividend_income = DIVIDEND_INCOME

# sys.it.mca_date = MCA_DATE

sys.ni.abolished = false
sys.ni.primary_class_1_rates  = [0.0, 0.0, 12.0, 2.0 ]
sys.ni.primary_class_1_bands  = [120.0, 183.0, 962.0, 9999999999999.9] # the '-1' here is because json can't write inf
sys.ni.secondary_class_1_rates  = [0.0, 13.8, 13.8 ] # keep 2 so
sys.ni.secondary_class_1_bands  = [170.0, 967.0, 99999999999999.9 ]
sys.ni.state_pension_age = 66; # fixme move
sys.ni.class_2_threshold = 6_475.0;
sys.ni.class_2_rate = 3.05;
sys.ni.class_4_rates  = [0.0, 9.0, 2.0 ]
sys.ni.class_4_bands  = [9_500.0, 50_270.0, 99999999999999.9 ]
# sys.ni.class_1_income = IncludedItems([WAGES],[PENSION_CONTRIBUTIONS_EMPLOYER])
# sys.ni.class_4_income = [SELF_EMPLOYMENT_INCOME]

# sys.uc.

sys.lmt.abolished = false

## FIXME we can't turn off pension credit individually here..

sys.lmt.premia.family = 17.60
sys.lmt.premia.family_lone_parent = 22.20
sys.lmt.premia.carer_single = 37.50
sys.lmt.premia.carer_couple = 2*37.50
sys.lmt.premia.disabled_child = 65.52
sys.lmt.premia.disability_single = 34.95
sys.lmt.premia.disability_couple = 49.80
sys.lmt.premia.enhanced_disability_child = 26.60
sys.lmt.premia.enhanced_disability_single = 17.10
sys.lmt.premia.enhanced_disability_couple = 24.50
sys.lmt.premia.severe_disability_single = 66.95
sys.lmt.premia.severe_disability_couple = 133.90
sys.lmt.premia.pensioner_is = 148.40

sys.lmt.allowances.age_18_24 = 58.90
sys.lmt.allowances.age_25_and_over = 74.35
sys.lmt.allowances.age_18_and_in_work_activity = 74.35
sys.lmt.allowances.over_pension_age = 187.75 #
sys.lmt.allowances.lone_parent = 74.70
sys.lmt.allowances.lone_parent_over_pension_age = 191.15
sys.lmt.allowances.couple_both_under_18 = 74.35
sys.lmt.allowances.couple_both_over_18 = 74.35
sys.lmt.allowances.couple_over_pension_age =  280.85
sys.lmt.allowances.couple_one_over_18_high =  116.80
sys.lmt.allowances.couple_one_over_18_med = 74.35
sys.lmt.allowances.pa_couple_one_over_18_low = 58.90
sys.lmt.allowances.child = 68.27
sys.lmt.allowances.pc_mig_single = 173.75
sys.lmt.allowances.pc_mig_couple = 265.20

# sys.lmt.income_rules.
sys.lmt.income_rules.permitted_work = 140.00
sys.lmt.income_rules.lone_parent_hb = 25.0
sys.lmt.income_rules.high = 20.0
sys.lmt.income_rules.low_couple = 10.0
sys.lmt.income_rules.low_single = 5.0
sys.lmt.income_rules.hb_additional = 17.10
sys.lmt.income_rules.childcare_max_1 = 175.0
sys.lmt.income_rules.childcare_max_2 = 300.0
sys.lmt.income_rules.incomes     = LEGACY_MT_INCOME
sys.lmt.income_rules.hb_incomes  = LEGACY_HB_INCOME
sys.lmt.income_rules.pc_incomes  = LEGACY_PC_INCOME
sys.lmt.income_rules.sc_incomes  = LEGACY_SAVINGS_CREDIT_INCOME
sys.lmt.income_rules.capital_min = 6_000.0
sys.lmt.income_rules.capital_max = 16_000.0
sys.lmt.income_rules.pc_capital_min = 10_000.0
sys.lmt.income_rules.pc_capital_max = 99999999999999.9
sys.lmt.income_rules.pensioner_capital_min = 10_000.0
sys.lmt.income_rules.pensioner_capital_max = 16_000.0

sys.lmt.income_rules.capital_tariff = 250
sys.lmt.income_rules.pensioner_tariff = 500
# FIXME why do we need a seperate copy of HoursLimits here?

sys.lmt.hours_limits.lower = 16
sys.lmt.hours_limits.med = 24
sys.lmt.hours_limits.higher = 30

sys.lmt.savings_credit.abolished = false
sys.lmt.savings_credit.withdrawal_rate = 60.0
sys.lmt.savings_credit.threshold_single = 150.47
sys.lmt.savings_credit.threshold_couple = 239.17
sys.lmt.savings_credit.max_single = 13.97
sys.lmt.savings_credit.max_couple = 15.62
sys.lmt.savings_credit.available_till = Date( 2016, 04, 06 )

sys.lmt.child_tax_credit.abolished = false
sys.lmt.child_tax_credit.family = 545.0
sys.lmt.child_tax_credit.child  = 2_830.0
sys.lmt.child_tax_credit.disability = 3_415
sys.lmt.child_tax_credit.severe_disability = 4825
sys.lmt.child_tax_credit.threshold = 16_385.0

sys.lmt.working_tax_credit.basic = 3_040 # with the 10pw uplift
sys.lmt.working_tax_credit.lone_parent = 2_045
sys.lmt.working_tax_credit.couple  = 2_045
sys.lmt.working_tax_credit.hours_ge_30 = 825
sys.lmt.working_tax_credit.disability = 3_220
sys.lmt.working_tax_credit.severe_disability = 1_390
sys.lmt.working_tax_credit.age_50_plus  = 1_365.00  # not used?
sys.lmt.working_tax_credit.age_50_plus_30_hrs = 2_030.00 # not used?
sys.lmt.working_tax_credit.childcare_max_2_plus_children = 300.0 # pw
sys.lmt.working_tax_credit.childcare_max_1_child  = 175.0
sys.lmt.working_tax_credit.childcare_proportion = 70.0
sys.lmt.working_tax_credit.taper = 41.0
sys.lmt.working_tax_credit.threshold = 6_530.0
sys.lmt.working_tax_credit.non_earnings_minima = 300.0 # FIXME check

sys.lmt.hb.taper = 65.0
sys.lmt.hb.ndd_deductions = [15.85,36.45,50.05,81.90,93.25,102.35]
sys.lmt.hb.ndd_incomes = [149.0,217.0,283.0,377.0,469.0,99999999999999.9]

sys.lmt.ctr.taper = 20.0
sys.lmt.ctr.ndd_deductions = []
sys.lmt.ctr.ndd_incomes = []

sys.uc.abolished = false
sys.uc.threshold = 2_500.0 ## NOT USED
sys.uc.age_18_24 = 342.72
sys.uc.age_25_and_over = 409.89

sys.uc.couple_both_under_25 = 488.59
sys.uc.couple_oldest_25_plus = 594.04

sys.uc.first_child  = 281.25
sys.uc.subsequent_child = 235.83
sys.uc.disabled_child_lower = 128.25
sys.uc.disabled_child_higher = 400.29
sys.uc.limited_capcacity_for_work_activity = 341.92
sys.uc.carer = 163.73

sys.uc.ndd = 75.15

sys.uc.childcare_max_2_plus_children  = 1_108.04 # pm
sys.uc.childcare_max_1_child  = 646.35
sys.uc.childcare_proportion = 85.0 # pct

sys.uc.minimum_income_floor_hours = 35*WEEKS_PER_MONTH

sys.uc.work_allowance_w_housing = 292.0
sys.uc.work_allowance_no_housing = 512.0
sys.uc.other_income = UC_OTHER_INCOME
# sys.uc.earned_income :: IncludedItems = UC_EARNED_INCOME
sys.uc.capital_min = 6_000.0
sys.uc.capital_max = 16_000.0
# £1 *per week* ≆ 4.35 pm FIXME make 4.35 WEEKS_PER_MONTH?
sys.uc.capital_tariff = 250.0/4.35
sys.uc.taper = 63.0
sys.uc.ctr_taper = 20.0

sys.age_limits.state_pension_ages = pension_ages()
sys.age_limits.savings_credit_to_new_state_pension = Date( 2016, 04, 06 )

sys.hours_limits.lower = 16
sys.hours_limits.med = 24
sys.hours_limits.higher = 30

sys.child_limits.max_children = 2
# TILL APRIL 22
sys.minwage.ages = [16,18,21,23]
sys.minwage.wage_per_hour = [4.62, 6.56, 8.36, 8.91]
sys.minwage.apprentice_rate = 4.30

sys.hr.maximum_rooms = 4
sys.hr.rooms_rent_reduction = [14.0,25.0]
sys.hr.single_room_age = 35
#
# These are unchanged in 3 years; see:
# https://www.gov.scot/publications/local-housing-allowance-rates-2021-2022/
#
sys.hr.brmas = loadBRMAs( 4, T, DEFAULT_BRMA_2021 )

sys.nmt_bens.attendance_allowance.abolished = false
sys.nmt_bens.attendance_allowance.higher = 89.60
sys.nmt_bens.attendance_allowance.lower = 60.00

sys.nmt_bens.child_benefit.abolished = false
sys.nmt_bens.child_benefit.first_child = 21.05
sys.nmt_bens.child_benefit.other_children = 13.95
sys.nmt_bens.child_benefit.high_income_thresh = 50_000.0
sys.nmt_bens.child_benefit.withdrawal = 1/100
sys.nmt_bens.child_benefit.guardians_allowance = 17.90

sys.nmt_bens.dla.abolished = false
sys.nmt_bens.dla.care_high = 89.15
sys.nmt_bens.dla.care_middle  = 59.70
sys.nmt_bens.dla.care_low  = 23.60
sys.nmt_bens.dla.mob_high  = 62.25
sys.nmt_bens.dla.mob_low  = 23.60

sys.nmt_bens.carers.abolished = false
sys.nmt_bens.carers.allowance =  67.25
sys.nmt_bens.carers.scottish_supplement = 231.40
sys.nmt_bens.carers.hours :: Int = 35
sys.nmt_bens.carers.gainful_employment_min = 128.00

sys.nmt_bens.pip.abolished = false
sys.nmt_bens.pip.dl_standard = 59.70
sys.nmt_bens.pip.dl_enhanced = 89.15
sys.nmt_bens.pip.mobility_standard = 23.60
sys.nmt_bens.pip.mobility_enhanced = 62.25

sys.nmt_bens.esa.abolished = false
sys.nmt_bens.esa.assessment_u25 = 58.90
sys.nmt_bens.esa.assessment_25p = 74.35
sys.nmt_bens.esa.main           = 74.35
sys.nmt_bens.esa.work           = 29.55
sys.nmt_bens.esa.support        = 39.20

sys.nmt_bens.jsa.abolished = false
sys.nmt_bens.jsa.u25 = 58.90
sys.nmt_bens.jsa.o24 = 74.35

sys.nmt_bens.pensions.abolished = false
sys.nmt_bens.pensions.new_state_pension = 175.20
# pension_start_date = Date( 2016, 04, 06 )
sys.nmt_bens.pensions.cat_a     = 134.25
sys.nmt_bens.pensions.cat_b     = 134.25
sys.nmt_bens.pensions.cat_b_survivor = 80.45
sys.nmt_bens.pensions.cat_d     = 80.45

sys.nmt_bens.bereavement.abolished = false
# higher effectively just means 'with children';
sys.nmt_bens.bereavement.lump_sum_higher = 3_500
sys.nmt_bens.bereavement.lump_sum_lower  = 2_500
sys.nmt_bens.bereavement.higher = 350
sys.nmt_bens.bereavement.lower  = 100

sys.nmt_bens.widows_pension.abolished = false
sys.nmt_bens.widows_pension.industrial_higher = 134.25
sys.nmt_bens.widows_pension.industrial_lower = 40.28
sys.nmt_bens.widows_pension.standard_rate =  121.95
sys.nmt_bens.widows_pension.parent =  121.95
sys.nmt_bens.widows_pension.ages = collect(54:-1:45)
sys.nmt_bens.widows_pension.age_amounts = [113.41,104.88,96.34,87.80,79.27,70.73,62.19,53.66,45.12,36.59]

#
# young carer grant
sys.nmt_bens.maternity.abolished = false
sys.nmt_bens.maternity.rate = 151.20

sys.nmt_bens.smp = 151.97 ## 90% of earn cpag 21/2 812
# = XX

sys.bencap.abolished = false
sys.bencap.outside_london_single = 257.69
sys.bencap.outside_london_couple = 384.62
# not really needed, but anyway ..
sys.bencap.inside_london_single = 296.35
sys.bencap.inside_london_couple  = 442.31
sys.bencap.uc_incomes_limit  = 617

sys.scottish_child_payment.amount = 10.0
sys.scottish_child_payment.maximum_age = 5

sys.ubi.abolished = true
sys.ubi.adult_amount = 4_800.0
sys.ubi.child_amount= 3_000.0
sys.ubi.universal_pension = 8_780.0
sys.ubi.adult_age = 17
sys.ubi.retirement_age = 66

weeklyise!( sys )


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

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local iv = try Base.loaded_modules[Base.PkgId(Base.UUID("6e696c72-6542-2067-7265-42206c756150"), "AbstractPlutoDingetjes")].Bonds.initial_value catch; b -> missing; end
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : iv(el)
        el
    end
end

# ╔═╡ 2f0ce4ef-4c45-4175-a9a7-df70cb268629
begin
    CURRENT_DIR = @__DIR__
    ROOT_DIR = if basename(CURRENT_DIR) == "pluto"
        dirname(dirname(CURRENT_DIR))
    else
        CURRENT_DIR
    end

    using Pkg
    Pkg.activate(ROOT_DIR; io = devnull)

    using Logging
    disable_logging(Logging.Warn)

    println('\n', " "^4, "> Loading the packages...")

    using LaTeXStrings
    using Plots
    using PlutoUI

    # Use the GR backend for plots
    gr()

    # Change some of the default parameters for plots
    default(fontfamily = "Computer Modern", dpi = 300, legend = nothing)

    # Define the paths to output directories
    PLOTS_DIR = joinpath(ROOT_DIR, "plots")
    NOTEBOOKS_PLOTS_DIR = joinpath(PLOTS_DIR, "notebooks")
    M1_DIR = joinpath(NOTEBOOKS_PLOTS_DIR, "m1")

    # Make sure the needed directories exist
    mkpath(M1_DIR)

    # An alias for a floating point type
    F = Float64

    md"This is a Pluto notebook to test the parameters of the first model."
end

# ╔═╡ 52342993-f23a-4a50-b3e0-e9e4b032f1cf
md"""
Masses are in units of ``M_0 = 2.325 \times 10^7 M_\odot``, other parameters are in ``\text{kpc}``.

Initial values are taken from Bajkova, Bobylev (2020), ranges are taken from Granados et al. (2021).
"""

# ╔═╡ 1fa86bf7-e1dc-4cd1-9ff0-73c14225ef2c
begin
    # Bulge
    M_B_S = @bind M_B NumberField(
        43.0:0.001:2580.0,
        default = 443.0,
    )
    A_B_S = @bind A_B NumberField(0.0:0.001:0.05, default = 0.0)
    B_B_S = @bind B_B NumberField(0.0:0.001:2.0, default = 0.2672)

    # Disk
    M_D_S = @bind M_D NumberField(
        1500.0:0.001:12903.0,
        default = 2798.0,
    )
    A_D_S = @bind A_D NumberField(1.0:0.001:20.0, default = 4.40)
    B_D_S = @bind B_D NumberField(0.1:0.001:16.0, default = 0.3084)

    # Halo
    M_H_S = @bind M_H NumberField(
        43.0:0.001:12903.0,
        default = 12474.0,
    )
    A_H_S = @bind A_H NumberField(0.1:0.001:30.0, default = 7.7)

    md"""
    ``M_b``: $(M_B_S) ``a_b``: $(A_B_S) ``b_b``: $(B_B_S) \
    ``M_d``: $(M_D_S) ``a_d``: $(A_D_S) ``b_d``: $(B_D_S) \
    ``M_h``: $(M_H_S) ``a_h``: $(A_H_S)
    """
end

# ╔═╡ a5d896ff-570f-49a9-adf0-19f11cb899e0
begin
    "Calculate the value of the R derivative of the Miyamoto & Nagai potential
    ``[ 100 \\, \\text{km} \\, \\text{s}^{-2} ]``"
    function miyamoto_nagai_phi_dr(r, z, m, a, b)::F
        return m * r / (r^2 + (a + √(z^2 + b^2))^2)^(1.5)
    end

    "Calculate the value of the R derivative of the Navarro-Frenk-White potential
    ``[ 100 \\, \\text{km} \\, \\text{s}^{-2} ]``"
    function navarro_frenk_white_phi_dr(r, z, m, a)::F
        sq_sum = r^2 + z^2
        k_1 = √(sq_sum) / a + 1.0
        return m * r * log(k_1) / sq_sum^(1.5) - m * r / a / sq_sum / k_1
    end

    "Calculate the value of circular velocity according to
    ``v_c^2(r) = R \\left. \\frac{\\partial \\Phi(R, Z)}{\\partial R} \\right|_{Z=0}``
    ``[ \\text{km}^2 \\, \\text{s}^{-2} ]``"
    circ_vel(r) = √(r * phi_dr(r, 0) * 100)

    "Calculate the value of ``\\partial \\Phi(R, Z) / \\partial R``
    ``[ 100 \\, \\text{km} \\, \\text{s}^{-2} ]``"
    function phi_dr(r, z)::F
        return miyamoto_nagai_phi_dr(r, z, M_B, A_B, B_B) +
               miyamoto_nagai_phi_dr(r, z, M_D, A_D, B_D) +
               navarro_frenk_white_phi_dr(r, z, M_H, A_H)
    end

    println(" "^4, "> Plotting the rotation curve...")

    # Plot the rotation curve
    p1 = plot(
        circ_vel, 0.2, 20;
        title = "Rotation curve",
        xlabel = L"R \;\, [\mathrm{kpc}]",
        ylabel = L"v_c \; [\mathrm{km^2 \, s^{-2}}]",
        xminorticks = 8,
        yminorticks = 8
    )

    # Save the figure
    savefig(p1, joinpath(M1_DIR, "Rotation curve.pdf"))
    savefig(p1, joinpath(M1_DIR, "Rotation curve.png"))

    p1
end

# ╔═╡ d4ba276c-a578-4a64-9f54-c33c81225a0b
begin
    # Plot the rotation curve
    p2 = plot(
        circ_vel, 0.2, 200;
        title = "Rotation curve",
        xlabel = L"R \;\, [\mathrm{kpc}]",
        ylabel = L"v_c \; [\mathrm{km^2 \, s^{-2}}]",
        xticks = [1e-1, 1, 1e1, 1e2],
        xlims = (0.1, 390),
        xscale = :log10,
        xminorticks = 9,
        yminorticks = 8
    )

    # Save the figure
    savefig(p2, joinpath(M1_DIR, "Rotation curve (log).pdf"))
    savefig(p2, joinpath(M1_DIR, "Rotation curve (log).png"))

    println()

    p2
end

# ╔═╡ bb727eda-9d7c-48e2-bac6-23e90aa6aab1
md"""
From Bajkova, Bobylev (2020):

$(PlutoUI.LocalResource(joinpath(ROOT_DIR, "assets", "Rotation curve of the NFWBB model (Bajkova, Bobylev, 2020).png")))
"""

# ╔═╡ Cell order:
# ╟─2f0ce4ef-4c45-4175-a9a7-df70cb268629
# ╟─52342993-f23a-4a50-b3e0-e9e4b032f1cf
# ╟─1fa86bf7-e1dc-4cd1-9ff0-73c14225ef2c
# ╟─a5d896ff-570f-49a9-adf0-19f11cb899e0
# ╟─d4ba276c-a578-4a64-9f54-c33c81225a0b
# ╟─bb727eda-9d7c-48e2-bac6-23e90aa6aab1

ext: jl


text:
# Implementation of the show function for Model, Kpoints, PlaneWaveBasis
# to avoid clutter in their respective files. For "shorter" files
# (e.g. elements.jl, Energies.jl) the preference is still to keep
# show in the same file to avoid forgetting to add a print statement
# for an added field.

function Base.show(io::IO, model::Model)
    nD = model.n_dim == 3 ? "" : "$(model.n_dim)D, "
    print(io, "Model(", model.model_name, ", ", nD,
          "spin_polarization = ", model.spin_polarization, ")")
end

function Base.show(io::IO, ::MIME"text/plain", model::Model)
    println(io, "Model(", model.model_name, ", $(model.n_dim)D):")
    for i = 1:3
        header = i==1 ? "lattice (in Bohr)" : ""
        showfieldln(io, header, (@sprintf "[%-10.6g, %-10.6g, %-10.6g]" model.lattice[i, :]...))
    end
    showfieldln(io, "unit cell volume", @sprintf "%.5g Bohr³" model.unit_cell_volume)

    if !isempty(model.atoms)
        println(io)
        showfieldln(io, "atoms", chemical_formula(model))
        elements = first.(model.atoms)
        for (i, el) in enumerate(elements)
            header = i==1 ? "atom potentials" : ""
            showfieldln(io, header, el)
        end
    end

    println(io)
    showfieldln(io, "num. electrons",    model.n_electrons)
    showfieldln(io, "spin polarization", model.spin_polarization)
    showfieldln(io, "temperature",       @sprintf "%.5g Ha" model.temperature)
    if model.temperature > 0
        showfieldln(io, "smearing",      model.smearing)
    end

    println(io)
    for (i, term) in enumerate(model.term_types)
        header = i==1 ? "terms" : ""
        showfield(io, header, sprint(show, "text/plain", term))
        i < length(model.term_types) && println(io)
    end
end

# TODO show function for the Kpoint struct

function Base.show(io::IO, basis::PlaneWaveBasis)
    print(io, "PlaneWaveBasis(model = ", basis.model, ", Ecut = ", basis.Ecut, " Ha")
    if !isnothing(basis.kgrid)
        print(io, ", kgrid = ", basis.kgrid)
        if !isnothing(basis.kshift) && !iszero(basis.kshift)
            print(io, ", kshift = ", basis.kshift)
        end
    else
        print(io, ", num. irred. kpoints = ", length(basis.kcoords_global))
    end
    print(io, ")")
end

function Base.show(io::IO, ::MIME"text/plain", basis::PlaneWaveBasis)
    println(io, "PlaneWaveBasis discretization:")
    showfieldln(io, "Ecut",     basis.Ecut, " Ha")
    showfieldln(io, "fft_size", basis.fft_size)
    if !basis.variational
        showfieldln(io, "variational", "false")
    end
    showfieldln(io, "kgrid type", "Monkhorst-Pack")
    showfieldln(io, "kgrid",    basis.kgrid)
    if !isnothing(basis.kshift) && !iszero(basis.kshift)
        showfieldln(io, "kshift",   basis.kshift)
    end
    showfieldln(io, "num. irred. kpoints", length(basis.kcoords_global))

    println(io)
    modelstr = sprint(show, "text/plain", basis.model)
    indent = " " ^ SHOWINDENTION
    print(io, indent, "Discretized " * replace(modelstr, "\n" => "\n" * indent))
end


text:
# include("E:\\Programs\\blockpop\\TSSOS\\src\\TSSOS.jl")
using TSSOS
include("E:\\Programs\\blockpop\\TSSOS\\examples\\modelopf.jl") # include the file modelopf.jl

cd("E:\\Programs\\PolyOPF\\pglib") # You need to download the problem data from PGLiB first.

silence()

case = "793_goc"
AC = 260200
data = parse_file("pglib_opf_case" * case * ".m")

# the first order relaxation
model,_ = pop_opf_real(data, normal=true, AngleCons=true, LineLimit="relax")
n = model.n
m = model.m
numeq = model.numeq
supp = model.supp
coe = model.coe
mc = maximum(abs.(coe[1]))
coe[1]=coe[1]./mc

time = @elapsed begin
opt,sol,popd = cs_tssos_first(supp, coe, n, 1, numeq=numeq, CS=false, TS="MF", MomentOne=false)
end
opt *= mc
mb = maximum(maximum.(popd.sb)) # maximal block size
gap = (AC-opt)*100/AC # optimality gap
println("n = $n, m = $m")
println("opt = $opt, time = $time, mb = $mb, gap = $gap%")

# the minimum order relaxation
model,_ = pop_opf_real(data, normal=true, AngleCons=true, LineLimit=true)
n = model.n
m = model.m
numeq = model.numeq
supp = model.supp
coe = model.coe
mc = maximum(abs.(coe[1]))
coe[1] = coe[1]./mc

time = @elapsed begin
opt,sol,popd = cs_tssos_first(supp, coe, n, "min", numeq=numeq, CS="MF", TS="block", MomentOne=true)
end
opt *= mc
maxc = maximum(popd.cliquesize) # maximal clique size
mb = maximum(maximum.(popd.sb)) # maximal block size
gap = 100*(AC-opt)/AC # optimality gap
println("n = $n, m = $m")
println("mc = $maxc, opt = $opt, time = $time, mb = $mb, gap = $gap%")


text:
# Julia wrapper for header: /usr/local/ffmpeg/2.4/include/libavutil/display.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

export
    av_display_rotation_get,
    av_display_rotation_set,
    av_display_matrix_flip

function av_display_rotation_get(matrix)
    ccall((:av_display_rotation_get,libavutil),Cdouble,(Ptr{Int32},),matrix)
end

function av_display_rotation_set(matrix,angle::Cdouble)
    ccall((:av_display_rotation_set,libavutil),Void,(Ptr{Int32},Cdouble),matrix,angle)
end

function av_display_matrix_flip(matrix,hflip::Integer,vflip::Integer)
    ccall((:av_display_matrix_flip,libavutil),Void,(Ptr{Int32},Cint,Cint),matrix,hflip,vflip)
end


text:
#! /usr/bin/julia

# Rosetta Code, Sort using a custom comparator

st = """You will rejoice to hear that no disaster has accompanied the
commencement of an enterprise which you have regarded with such evil
forebodings."""

wl = filter(x->length(x)>0, split(st, r"\W+"))

println("Original List:")
for w in wl
    println("   ", w)
end

wl = sort(wl, by=x->(-length(x), lowercase(x)))

println("\nSorted List:")
for w in wl
    println("   ", w)
end

max_stars_repo_name: MichaeLeroy/rosetta-code

max_issues_repo_name: MichaeLeroy/rosetta-code

max_forks_repo_name: MichaeLeroy/rosetta-code


text:
import Base: intersect, reverse

function getindex(r::StaticUnitRange, s::Union{StepRange{<:Integer},StaticStepRange{<:Integer}})
    Base.@_inline_meta
    @boundscheck checkbounds(r, s)
    st = ofeltype(first(r), first(r) + first(s)-SOne)
    range(st, step=step(s), length=length(s))
end

function getindex(r::StaticStepRange, s::AbstractRange{<:Integer})
    Base.@_inline_meta
    @boundscheck checkbounds(r, s)
    st = oftype(first(r), first(r) + (first(s)-SOne)*step(r))
    range(st, step=step(r)*step(s), length=length(s))
end

function getindex(r::StaticStepRangeLen{T},
                  s::Union{OrdinalRange{<:Integer},StaticOrdinalRange{<:Integer}}) where {T}
    Base.@_inline_meta
    @boundscheck checkbounds(r, s)
    # Find closest approach to offset by s
    ind = LinearIndices(s)
    f = max(min(SOne + round(Int, (offset(r) - first(s))/step(s)), last(ind)), first(ind))
    ref = _getindex_hiprec(r, first(s) + (f-SOne)*step(s))
    return StaticStepRangeLen{T}(ref, step(r)*step(s), length(s), f)
end

function getindex(
    r::StaticLinRange,
    s::Union{OrdinalRange{<:Integer},StaticOrdinalRange{<:Integer}})
    Base.@_inline_meta
    @boundscheck checkbounds(r, s)
    return StaticLinRange(unsafe_getindex(r, first(s)), vlast,
                          unsafe_getindex(r, last(s)))
end

function getindex(r::Union{StaticStepRangeLen,StaticLinRange}, i::Integer)
    Base.@_inline_meta
    @boundscheck checkbounds(r, i)
    unsafe_getindex(r, i)
end

# TODO ensure this computes the right thing on LinMRange
unsafe_getindex(r::StaticLinRange{T}, i::I) where {T,I<:Integer} =
    lerpi(i-one(i), lendiv(r), first(r), last(r))

function lerpi(j::Integer, d::Integer, a::A, b::B) where {A,B}
    Base.@_inline_meta
    t = j/d
    ofeltype(a, (one(t) - t) * a + t * b)
end

@inline function unsafe_getindex(
    r::StaticStepRangeLen{T,B}, i::Integer) where {T,B<:Union{<:TPVal,TwicePrecision}}
    Base.@_inline_meta
    u = i - offset(r)
    shift_hi, shift_lo = u*getstephi(r), u*getsteplo(r)
    x_hi, x_lo = add12(getrefhi(r), shift_hi)
    ofeltype(T, x_hi + (x_lo + (shift_lo + getreflo(r))))
end

@inline function unsafe_getindex(r::StaticStepRangeLen{T,B}, i::Integer) where {T,B}
    Base.@_inline_meta
    ofeltype(T, reference(r) + (i - offset(r)) * step(r))
end

# TODO: just for ordinal range in base, make sure it doesn't do anything crazy
Base.iterate(r::StaticRange) = isempty(r) ? nothing : (first(r), first(r))

function Base.iterate(r::StaticRange{T}, i::I) where {T,I}
    Base.@_inline_meta
    i == last(r) && return nothing
    next = ofeltype(T, i + step(r))
    (next, next)
end

reverse(r::StaticRange) = (:)(last(r), -step(r), first(r))
function reverse(r::StaticStepRangeLen)
    # If `r` is empty, `length(r) - r.offset + 1 will be nonpositive hence
    # invalid. As `reverse(r)` is also empty, any offset would work so we keep
    # `r.offset`
    offset = isempty(r) ? offset(r) : length(r)-offset(r)+One
    StaticStepRangeLen(first(r), -step(r), length(r), offset)
end

reverse(r::StaticLinRange) = StaticLinRange(last(r), first(r), length(r))

function in(x::Real, r::StaticRange{T}) where {T<:Integer}
    isa(x, Integer) && !isempty(r) && x >= minimum(r) && x <= maximum(r) &&
        (mod(ofeltype(T,x), step(r)) - mod(first(r), step(r)) == 0)
end

function in(x::AbstractChar, r::StaticRange{<:AbstractChar})
    !isempty(r) && x >= minimum(r) && x <= maximum(r) &&
        (mod(Int(x) - Int(first(r)), step(r)) == 0)
end

in(x::Integer, r::StaticUnitRange{<:Integer}) = (first(r) <= x) & (x <= last(r))

function in(x, r::StaticRange)
    anymissing = false
    for y in itr
        v = (y == x)
        if ismissing(v)
            anymissing = true
        elseif v
            return true
        end
    end
    return anymissing ? missing : false
end

intersect(r::StaticUnitRange{<:Integer}, s::StaticUnitRange{<:Integer}) =
    max(first(r),first(s)):min(last(r),last(s))

intersect(i::Integer, r::StaticUnitRange{<:Integer}) =
    i < first(r) ? (first(r):i) :
    i > last(r) ? (i:last(r)) : (i:i)

intersect(r::StaticUnitRange{<:Integer}, i::Integer) = intersect(i, r)

function intersect(r::StaticUnitRange{<:Integer},
                   s::Union{StaticStepRange{<:Integer},StepRange{<:Integer}})
    if isempty(s)
        range(first(r), length=0)
    elseif step(s) == 0
        intersect(first(s), r)
    elseif step(s) < 0
        intersect(r, reverse(s))
    else
        sta = first(s)
        ste = step(s)
        sto = last(s)
        lo = first(r)
        hi = last(r)
        i0 = max(sta, lo + mod(sta - lo, ste))
        i1 = min(sto, hi - mod(hi - sta, ste))
        i0:ste:i1
    end
end

function intersect(r::StaticStepRange{<:Integer}, s::Union{AbstractUnitRange{<:Integer},<:StaticUnitRange{<:Integer}})
    if step(r) < 0
        reverse(intersect(s, reverse(r)))
    else
        intersect(s, r)
    end
end

function intersect(r::StaticStepRange, s::Union{<:StepRange,<:StaticStepRange})
    if isempty(r) || isempty(s)
        return range(first(r), step=step(r), length=0)
    elseif step(s) < 0
        return intersect(r, reverse(s))
    elseif step(r) < 0
        return reverse(intersect(reverse(r), s))
    end

    start1 = first(r)
    step1 = step(r)
    stop1 = last(r)
    start2 = first(s)
    step2 = step(s)
    stop2 = last(s)
    a = lcm(step1, step2)

    # if a == 0
    #     # One or both ranges have step 0.
    #     if step1 == 0 && step2 == 0
    #         return start1 == start2 ? r : AbstractRange(start1, 0, 0)
    #     elseif step1 == 0
    #         return start2 <= start1 <= stop2 && rem(start1 - start2, step2) == 0 ? r : AbstractRange(start1, 0, 0)
    #     else
    #         return start1 <= start2 <= stop1 && rem(start2 - start1, step1) == 0 ? (start2:step1:start2) : AbstractRange(start1, step1, 0)
    #     end
    # end

    g, x, y = gcdx(step1, step2)

    if rem(start1 - start2, g) != 0
        # Unaligned, no overlap possible.
        return range(start1, step=a, length=SZero)
    end

    z = div(start1 - start2, g)
    b = start1 - x * z * step1
    # Possible points of the intersection of r and s are
    # ..., b-2a, b-a, b, b+a, b+2a, ...
    # Determine where in the sequence to start and stop.
    m = max(start1 + mod(b - start1, a), start2 + mod(b - start2, a))
    n = min(stop1 - mod(stop1 - b, a), stop2 - mod(stop2 - b, a))
    m:a:n
end

Base.map(::Type{T}, r::StaticStepRange) where T =
    ofeltype(T, first(r)):ofeltype(T, setp(r)):ofeltype(T, last(r))
Base.map(::Type{T}, r::StaticUnitRange) where T =
    ofeltype(T, first(r)):ofeltype(T, last(r))
Base.map(::Type{T}, r::StaticStepRangeLen) where T =
    StaticStepRangeLen(ofeltype(T, first(r)), ofeltype(T, step(r)), length(r), offset(r))
Base.map(::Type{T}, r::StaticLinRange) where T =
    StaticLinRange(ofeltype(T, first(r)), ofeltype(T, last(r)), length(r))

Base.minimum(r::StaticRange) = min(first(r), last(r))
Base.maximum(r::StaticRange) = max(first(r), last(r))

Base.extrema(r::StaticRange) = (minimum(r), maximum(r))

function Base.minimum(r::StaticUnitRange)
    Base.@_inline_meta
    isempty(r) ? throw(ArgumentError("range must be non-empty")) : first(r)
end

function maximum(r::StaticUnitRange)
    Base.@_inline_meta
    isempty(r) ? throw(ArgumentError("range must be non-empty")) : last(r)
end

minimum(r::AbstractRange)  = isempty(r) ? throw(ArgumentError("range must be non-empty")) : min(first(r), last(r))
maximum(r::AbstractRange)  = isempty(r) ? throw(ArgumentError("range must be non-empty")) : max(first(r), last(r))

-(r1::StaticStepRangeLen, r2::StaticStepRangeLen) = +(r1, -r2)
-(r::StaticStepRange) = range(-first(r), step=-step(r), length=length(r))

-(r::StaticStepRangeLen) = StaticStepRangeLen(-first(r), -step(r), length(r), offset(r))
-(r::StaticLinRange) = StaticLinRange(-first(r), -last(r), length(r))

+(r1::StaticStepRangeLen{T,B1,S1,E1,L1}, r2::StaticStepRangeLen{T,B2,S2,E2,L2}) where {T,B1,S1,E1,L1<:SInteger,B2,S2,E2,L2<:SInteger} =
    throw(DimensionMismatch("argument dimensions must match"))
function +(r1::StaticStepRangeLen{T,B1,S1,E1,L}, r2::StaticStepRangeLen{T,B2,S2,E2,L}) where {T,B1,S1,E1,B2,S2,E2,L<:SInteger}
    StaticStepRangeLen(first(r1)+first(r2), step(r1)+step(r2), L())
end

function +(r1::StaticStepRangeLen{T}, r2::StaticStepRangeLen{T}) where {T}
    len = length(r1)
    (len == length(r2) || throw(DimensionMismatch("argument dimensions must match")))
    StaticStepRangeLen(first(r1)+first(r2), step(r1)+step(r2), len)
end

for f in (:+, :-)
    @eval begin
        function $f(r1::StaticStepRange, r2::StaticStepRange)
            r1l = length(r1)
            (r1l == length(r2) ||
             throw(DimensionMismatch("argument dimensions must match")))
            range($f(first(r1), first(r2)), step=$f(step(r1), step(r2)), length=r1l)
        end

        function $f(r1::StaticLinRange{T}, r2::StaticLinRange{T}) where T
            len = length(r1)
            (len == length(r2) ||
             throw(DimensionMismatch("argument dimensions must match")))
            StaticLinRange{T}(convert(T, $f(first(r1), first(r2))),
                              convert(T, $f(last(r1), last(r2))), len)
        end

        $f(r1::Union{StaticStepRangeLen, StaticOrdinalRange, StaticLinRange},
           r2::Union{StaticStepRangeLen, StaticOrdinalRange, StaticLinRange}) =
               $f(promote(r1, r2)...)
    end
end

function +(
    r1::StaticRange{T,L1},
    r2::StaticRange{T,L2}) where {T,L1,L2}
    throw(DimensionMismatch("argument dimensions must match"))
end

Base.issorted(r::StaticUnitRange) = true

Base.sort(r::StaticUnitRange) = r
Base.sort!(r::StaticUnitRange) = r

Base.sort(r::StaticRange) = issorted(r) ? r : reverse(r)

Base.sortperm(r::StaticUnitRange) = SOne:length(r)
Base.sortperm(r::StaticRange) = issorted(r) ? (SOne:SOne:length(r)) : (length(r):-SOne:SOne)

# TODO these can easily be optimized with static params
Base.isempty(r::StaticStepRange) = (first(r) != last(r)) & ((step(r) > zero(step(r))) != (last(r) > first(r)))
Base.isempty(r::StaticUnitRange) = first(r) > last(r)
Base.isempty(r::StaticStepRangeLen) = length(r) == 0
Base.isempty(r::StaticLinRange) = length(r) == 0

Base.map(::Type{T}, r::StaticStepRange) where {T<:Real} = ofeltype(T, first(r)):ofeltype(T, step(r)):ofeltype(T, last(r))
Base.map(::Type{T}, r::StaticUnitRange) where {T<:Real} = ofeltype(T, first(r)):ofeltype(T, last(r))
Base.map(::Type{T}, r::StaticStepRangeLen) where {T<:AbstractFloat} =
    StaticStepRangeLen{T}(ofeltype(T, first(r)), ofeltype(T, step(r)), length(r), offset(r))
Base.map(::Type{T}, r::StaticLinRange) where T<:AbstractFloat =
    StaticLinRange{T}(ofeltype(T, first(r)), ofeltype(T, last(r)), length(r))

==(r::Union{StaticStepRangeLen{T},StaticLinRange{T}}, s::Union{StaticStepRangeLen{T},StaticLinRange{T}}) where T =
    (first(r) == first(s)) & (length(r) == length(s)) & (last(r) == last(s))
==(r::Union{StepRangeLen{T},LinRange{T}}, s::Union{StaticStepRangeLen{T},StaticLinRange{T}}) where T =
    (first(r) == first(s)) & (length(r) == length(s)) & (last(r) == last(s))
==(r::Union{StaticStepRangeLen{T},StaticLinRange{T}}, s::Union{StepRangeLen{T},LinRange{T}}) where T =
    (first(r) == first(s)) & (length(r) == length(s)) & (last(r) == last(s))


text:
import ObjFileBase: getSectionLoadAddress

r_symbolnum(r::relocation_info) = r.target & 0xFFFFFF
r_pcrel(r::relocation_info) = ((r.target & (1<<24)) >> 24) != 0
r_length(r::relocation_info) = (r.target & (UInt32(0b11) << 25)) >> 25
r_extern(r::relocation_info) = ((r.target & (1 << 27)) >> 27) != 0
r_type(r::relocation_info) = (r.target & (1 << 28)) >> 28
for f in (:r_symbolnum, :r_pcrel, :r_length, :r_extern, :r_type)
    @eval ($f)(r::RelocationRef) = ($f)(deref(r))
end

# The relocation to apply
immutable RelocToApply
    value::UInt64
    size::UInt8
end

function compute_value(h, symbols, sections, LOI, reloc)
    # Here one could to external symbol lookups, etc, but I'm not interested
    reloc = deref(reloc)
    if r_extern(reloc)
        symbols[r_symbolnum(reloc)]
    else
        @assert LOI !== nothing
        sections[r_symbolnum(reloc)]
    end
end

# Apply relocations in `buffer`. `h` should be the buffer being relocated
function relocate!(buffer, h; LOI = nothing, debug_only = true)
    scs = filter(LoadCmds(h)) do lc
        isa(lc.cmd, segment_commands)
    end
    symbols = Symbols(first(filter(LoadCmds(h)) do lc
        isa(lc.cmd, symtab_command)
    end))
    header = readheader(h)
    for lc in scs
        sects = Sections(lc)
        for sec in sects
            if debug_only && deref(sec).segname != "__DWARF"
                continue
            end
            for reloc in Relocations(sec)
              Value = compute_value(h, symbols, sects, LOI, reloc)
              Value = getSectionLoadAddress(LOI, Value) - deref(Value).addr
              rta = compute_relocation(header, reloc, Value)
              seek(buffer, sectionoffset(sec) + deref(reloc).address)
              write(buffer, rta.size == 8 ? rta.value :
                  rta.size == 4 ? convert(UInt32,rta.value) :
                  rta.size == 2 ? convert(UInt16,rta.value) :
                  rta.size == 1 ? convert(UInt8,rta.value) :
                  error("Unsupported Relocation Size"))
            end
        end
    end
end

function compute_relocation(header, reloc, Value)
    if header.cputype == CPU_TYPE_X86_64
        compute_X86_64(reloc, Value)
    end
end

function compute_X86_64(reloc, Value)
    if r_type(deref(reloc)) == X86_64_RELOC_UNSIGNED
        RelocToApply(Value, 1 << r_length(reloc))
    else
        error("Unknown Relocation $(X86_64_RELOC[r_type(deref(reloc))])")
    end
end

function relocateByLOI!(buffer, h, LOI)
    scs = filter(LoadCmds(h)) do lc
        isa(lc.cmd, segment_commands)
    end
    for lc in scs
      sects = Sections(lc)
      for sec in sects
          addr = try
              getSectionLoadAddress(LOI, sec)
          catch
              continue
          end
          ssize = sectionsize(sec)
          seek(buffer, sectionoffset(sec))
          write(buffer, pointer_to_array(
            reinterpret(Ptr{UInt8}, addr), (ssize,), false))
      end
    end
end


text:
"""
    ib_rh(set, ps, vs, gas)
    ib_rh(Ma, gam)

Initialize Rankine-Hugoniot relation
"""
function ib_rh(
    set::AbstractSetup,
    ps::AbstractPhysicalSpace,
    vs::Union{AbstractVelocitySpace,Nothing},
    gas::AbstractProperty,
)

    gam = gas.γ
    MaL = gas.Ma
    MaR = sqrt((MaL^2 * (gam - 1.0) + 2.0) / (2.0 * gam * MaL^2 - (gam - 1.0)))
    ratioT =
        (1.0 + (gam - 1.0) / 2.0 * MaL^2) * (2.0 * gam / (gam - 1.0) * MaL^2 - 1.0) /
        (MaL^2 * (2.0 * gam / (gam - 1.0) + (gam - 1.0) / 2.0))

    if set.nSpecies == 1

        primL = [1.0, MaL * sqrt(gam / 2.0), 1.0]
        primR = [
            primL[1] * (gam + 1.0) * MaL^2 / ((gam - 1.0) * MaL^2 + 2.0),
            MaR * sqrt(gam / 2.0) * sqrt(ratioT),
            primL[3] / ratioT,
        ]

        primL = ifelse(set.space[5:6] == "3v", [primL[1:2]; zeros(2); primL[end]], primL)
        primR = ifelse(set.space[5:6] == "3v", [primR[1:2]; zeros(2); primR[end]], primR)

        wL = prim_conserve(primL, gam)
        wR = prim_conserve(primR, gam)

        fw = function (x)
            if x <= (ps.x0 + ps.x1) / 2
                return wL
            else
                return wR
            end
        end

        bc = function (x)
            if x <= (ps.x0 + ps.x1) / 2
                return primL
            else
                return primR
            end
        end

        if set.space[1:4] == "1d0f"
            return fw, bc
        elseif set.space == "1d1f1v"
            ff = function (x)
                w = fw(x)
                prim = conserve_prim(w, gas.γ)
                h = maxwellian(vs.u, prim)
                return h
            end

            return fw, ff, bc
        elseif set.space == "1d2f1v"
            ff = function (x)
                w = fw(x)
                prim = conserve_prim(w, gas.γ)
                h = maxwellian(vs.u, prim)
                b = h * gas.K / 2 / prim[end]
                return h, b
            end

            return fw, ff, bc
        elseif set.space == "1d1f3v"
            ff = function (x)
                w = fw(x)
                prim = conserve_prim(w, gas.γ)
                h = maxwellian(vs.u, vs.v, vs.w, prim)
                return h
            end

            return fw, ff, bc
        end

    elseif set.nSpecies == 2

        primL = zeros(3, 2)
        primL[:, 1] .= [mi, MaL * sqrt(gam / 2.0), mi / 1.0]
        primL[:, 2] .= [me, MaL * sqrt(gam / 2.0), me / 1.0]

        primR = zeros(3, 2)
        for j = 1:2
            primR[1, j] = primL[1, j] * (gam + 1.0) * MaL^2 / ((gam - 1.0) * MaL^2 + 2.0)
            primR[2, j] = MaR * sqrt(gam / 2.0 / (mi * ni + me * ne)) * sqrt(ratioT)
            primR[3, j] = primL[3, j] / ratioT
        end

        wL = mixture_prim_conserve(primL, gam)
        wR = mixture_prim_conserve(primR, gam)

        fw = function (x)
            if x <= (ps.x0 + ps.x1) / 2
                return wL
            else
                return wR
            end
        end

        bc = function (x)
            if x <= (ps.x0 + ps.x1) / 2
                return primL
            else
                return primR
            end
        end

        if set.space == "1d2f1v"
            hL = mixture_maxwellian(vs.u, primL)
            hR = mixture_maxwellian(vs.u, primR)
            bL = similar(hL)
            bR = similar(hR)
            for j = 1:2
                bL[:, j] .= hL[:, j] .* gas.K ./ (2.0 .* primL[end, j])
                bR[:, j] .= hR[:, j] .* gas.K ./ (2.0 .* primR[end, j])
            end

            ff = function (x)
                if x <= (ps.x0 + ps.x1) / 2
                    return hL, bL
                else
                    return hR, bR
                end
            end

            return fw, ff, bc
        end

    end

    return nothing

end

function ib_rh(Ma::Real, gam::Real)
    MaL = Ma
    MaR = sqrt((MaL^2 * (gam - 1.0) + 2.0) / (2.0 * gam * MaL^2 - (gam - 1.0)))
    ratioT =
        (1.0 + (gam - 1.0) / 2.0 * MaL^2) * (2.0 * gam / (gam - 1.0) * MaL^2 - 1.0) /
        (MaL^2 * (2.0 * gam / (gam - 1.0) + (gam - 1.0) / 2.0))

    primL = [1.0, MaL * sqrt(gam / 2.0), 1.0]
    primR = [
        primL[1] * (gam + 1.0) * MaL^2 / ((gam - 1.0) * MaL^2 + 2.0),
        MaR * sqrt(gam / 2.0) * sqrt(ratioT),
        primL[3] / ratioT,
    ]

    return primL, primR
end


text:
"""
    hists_to_bars(hist1ds)

Given a vector of Hist1D, return `edges` (xs), `heights` (ys),
and `grps` (for grouping) that is useful for plotting stacked
histogram.
"""
function hists_to_bars(hist1ds)
    L = length(hist1ds)
    oneedge = binedges(hist1ds[1])[1:end-1]
    edges = repeat(oneedge, L)
    heights = mapreduce(bincounts, vcat, hist1ds)
    grps = repeat(1:L; inner=length(oneedge))

    edges, heights, grps
end

@inline function pearson_err(n::Real)
    s = sqrt(n+0.25)
    s+0.5, s-0.5
end
@inline function sqrt_err(n::Real)
    s = sqrt(n)
    s,s
end

_sturges(x) = StatsBase.sturges(length(x))

@inline function _is_uniform_bins(A::AbstractVector{T}) where T<:Real
    diffs = diff(A)
    diff1 = first(diffs)
    all(isapprox.(diff1, diffs; atol = 1e-9)) #hack
end
function _is_uniform_bins(A::AbstractRange{T}) where T<:Real
    true
end


text:
using BHAtp

ProjDir = @__DIR__

isdir(joinpath(ProjDir, "tp0")) && rm(joinpath(ProjDir, "tp0"), recursive=true)
isdir(joinpath(ProjDir, "tp")) && rm(joinpath(ProjDir, "tp"), recursive=true)
isdir(joinpath(ProjDir, "plots")) && rm(joinpath(ProjDir, "plots"), recursive=true)

ProjName = split(ProjDir, "/")[end]
bhaj = BHAJ(ProjName, ProjDir)
bhaj.ratio = 1.7

segs = [
# Element type,  Material,    Length,     OD,         ID,      fc
  (:bit,          :steel,     0.00,   2.75,   9.0,  0.0),
  (:collar,       :steel,    45.00,   2.75,   7.0,  0.0),
  (:stabilizer,   :steel,     0.00,   2.75,   9.0,  0.0),
  (:pipe,         :steel,    30.00,   2.75,   7.0,  0.0),
  (:stabilizer,   :steel,     0.00,   2.75,   9.0,  0.0),
  (:pipe,         :steel,    90.00,   2.75,   7.0,  0.0),
  (:stabilizer,   :steel,     0.00,   2.75,   9.0,  0.0),
  (:pipe,         :steel,   100.00,   2.75,   7.0,  0.0)
]

traj = [
#   Heading,      Diameter
  ( 60.0,      9.0)
]

wobs = 20:10:40
incls = 20:10:40             # Or e.g. incls = [5 10 20 30 40 45 50]

@time bhaj(segs, traj, wobs, incls)
println()

display("Fetch tp=0, wob=40, incl@bit=20 solution:")
df,df_tp = show_solution(ProjDir, 40, 20, show=false, tp=false)
df[:,[2; 5:6; 9:12]] |> display
println()

println("Fetch wob=40, incl@bit=20 solution:")
df,df_tp = show_solution(ProjDir, 40, 20, show=false);
df[:,[2; 5:6; 9:12]] |> display
println()

show_tp(ProjDir, wobs, incls) |> display