fixing the example

    output_data_path.mkdir(parents=True, exist_ok=True)

    # define the poro-depth and the poro-perm relationships:

    # porosity = a * exp(-b * depth) + c

    #  -with depth in km

    # porosity = a * exp(-b * depth), with depth in km

    poro_a = 28.97

    poro_b = 0.21

    poro_c = 0.0

    poro_std = 1

    poro_std = 5  # 1% uncertainty

    # ln(permeability) = a + b * poro

    # ln(permeability) = a * poro + b

    perm_a = 0.334

    perm_b = -1.655

    perm_std = 1.556

    # define range of depth interest:

    # number of depths and number of samples per depth

    n_samples = 500

    n_depths = 10

    thickness = 100

    ntg = 1.0

    porosity = poro_a * np.exp(-poro_b * (depths * 1e-3)) + poro_c

    porosity = poro_a * np.exp(-poro_b * (depths * 1e-3))

    # expand porosity to have n_samples per depth, with uncertainty in porosity

    porosity = np.random.normal(porosity[:, None], poro_std, size=(n_depths, n_samples))

    porosity = np.clip(porosity, 0.0, 100.0)

    permeability = np.exp(perm_a * porosity + perm_b)

    permeability = np.random.normal(permeability, perm_std)

    permeability = np.clip(permeability, 1.0, 1e4)

    reservoir_properties = xr.Dataset(

        {

        },

    )

    # For every sample, run a doublet simulation store the output values

    # For every sample, run a doublet simulation

    simulations = calculate_doublet_performance(

        reservoir_properties,

        chunk_size=100,

    )

    simulations_stoch = simulations.quantile([0.1,0.5,0.9], dim="samples")

    print(simulations)

    print(simulations_stoch)

    # make a plot showing there is a sweet spot between Doublet power, Cost of energy

    # and Transmissivity

    p50 = simulations_stoch.sel(quantile=0.5)

    p90 = simulations_stoch.sel(quantile=0.1)

    # UTC

    ax2.fill_betweenx(p50.depth, p10.utc, p90.utc, color="lightgrey", alpha=0.5)

    ax2.plot(p50.utc, p50.depth, label="UTC [€/kWh]", color="lightgrey")

    # power

    ax2.fill_betweenx(p50.depth, p10.power, p90.power, color="black", alpha=0.5)

    ax2.plot(p50.power, p50.power.depth, label="Power [MW]", color="black")

    ax.fill_betweenx(p50.depth, p10.transmissivity_with_ntg, p90.transmissivity_with_ntg, color="tab:green", alpha=0.5)

    ax.plot(p50.transmissivity_with_ntg, p50.transmissivity_with_ntg.depth, label="Transmissivity [Dm]", color="tab:green")

    # transmissivity Dm

    ax.fill_betweenx(

        p50.depth,

        p10.transmissivity_with_ntg,

        p90.transmissivity_with_ntg,

        color="tab:green",

        alpha=0.5,

    )

    ax.plot(

        p50.transmissivity_with_ntg,

        p50.transmissivity_with_ntg.depth,

        label="Transmissivity [Dm]",

        color="tab:green",

    )

    ax.set_ylabel("Depth [m]")

    ax2.set_xlabel("Power [MW] & UTC [€/kWh]")

    ax.set_xlabel("Transmissivity [Dm]", color="tab:green")

    ax.set_ylim([depths[-1], 0.0])

    plt.show()

    ax.set_ylim([depths[-1], depths[0]])

    # Collect handles and labels from both axes

    handles1, labels1 = ax.get_legend_handles_labels()

    handles2, labels2 = ax2.get_legend_handles_labels()

    handles = handles1 + handles2

    labels = labels1 + labels2

    ax.legend(handles, labels)

    plt.savefig(output_data_path / "power_utc_transmissivity.png")