Resolve "cleanup the workflow/utc code"

from pythermogis.workflow.utc.flow import calculate_volumetric_flow

from pythermogis.workflow.utc.doublet_utils import calc_lifetime

from pythermogis.workflow.utc.flow import calculate_volumetric_flow

def calculate_cooling_temperature(

    props,

    input,

    drawdown_pressure: float,

    well_distance: float,

    injection_temp: float

    props, input, drawdown_pressure: float, well_distance: float, injection_temp: float

) -> float:

    results = calculate_volumetric_flow(

        props=props,

        input_data=input,

    if (not props.is_ates) and cooling_frac > 0.3:

        raise RuntimeError(

            f"Large cooling factor ({cooling_frac}) could result in less reliable calculation"

            f"Large cooling factor ({cooling_frac}) "

            f"could result in less reliable calculation"

        )

    return cooling_frac * (input.temperature - injection_temp)

EUR_PER_CT_PER_KWH = 0.36

NPV_SCALE = 1e-6

@dataclass

class DoubletInput:

    unknown_input_value: float

        DARCY_SI = 1.0e-12 / 1.01325

        return self.transmissivity / self.thickness * 1e-3 * DARCY_SI

@dataclass

class DoubletOutput:

    power: float

    production_temp: float

    injection_temp: float

def calculate_doublet_performance(props: UTCConfiguration, input: DoubletInput, verbose: bool = False) -> DoubletOutput | None:

def calculate_doublet_performance(

    props: UTCConfiguration, input: DoubletInput, verbose: bool = False

) -> DoubletOutput | None:

    timer = timeit.default_timer()

    # determine initial well distance

    well_distance = (

            props.hp_minimum_injection_temperature,

        )

    else:

        injection_temperature = max(input.temperature - props.max_cooling_temp_range,

            props.dh_return_temp)

        injection_temperature = max(

            input.temperature - props.max_cooling_temp_range, props.dh_return_temp

        )

    timer = print_time(timer, "injection temperature: ", verbose=verbose)

    # calculate maximum pressure

    # cooling temperature and well distance optimization

    if props.optim_well_dist:

        end_temperature_p = (

                props.optim_dist_cooling_fraction *

                (input.temperature - injection_temperature)

        end_temperature_p = props.optim_dist_cooling_fraction * (

            input.temperature - injection_temperature

        )

    else:

        end_temperature_p = calculate_cooling_temperature(

    # stimulation capex

    stimulation_capex = (

        0.0

        if (not props.use_stimulation or input.transmissivity_with_ntg > props.stim_kh_max)

        if (

            not props.use_stimulation

            or input.transmissivity_with_ntg > props.stim_kh_max

        )

        else props.stimulation_capex

    )

        props.salinity_surface,

        props.salinity_gradient,

        props.optim_dist_cp_rock,

        props.optim_dist_rho_rock

        props.optim_dist_rho_rock,

    )

    utc_cutoff = (

        props.utc_cutoff_deep if input.depth > props.utc_deep_depth else props.utc_cutoff

        props.utc_cutoff_deep

        if input.depth > props.utc_deep_depth

        else props.utc_cutoff

    )

    npv = (

        breakthrough=breakthrough_years,

        cooling=end_temperature_p,

        production_temp=production_temp,

        injection_temp=injection_temperature

        injection_temp=injection_temperature,

    )

import math

from numba import njit

from pythermogis.workflow.utc.water import density, heat_capacity, get_salinity, get_hydrostatic_pressure

from pythermogis.workflow.utc.water import (

    density,

    get_hydrostatic_pressure,

    get_salinity,

    heat_capacity,

)

@njit

def calculate_injection_temp_with_heat_pump(

    return max(

        needed_injection_temp,

        max(reservoir_temp - max_cooling_temp_range, hp_minimum_injection_temperature)

        max(reservoir_temp - max_cooling_temp_range, hp_minimum_injection_temperature),

    )

@njit

def get_orc_efficiency(Tx: float, Ts: float, etarel: float) -> float:

    return etarel * (Tx - Ts) / (Tx + Ts + 2 * 273.1)

@njit

def get_cop_carnot(eta: float, Tout: float, Tin: float) -> float:

    if Tout < 0.0 or Tin < 0.0:

    return eta * (Tcond + TKELVIN) / (Tcond - Tevap)

@njit

def calc_lifetime(

    well_distance: float,

    cp_rock: float,

    rho_rock: float,

) -> float:

    # --- Water properties at depth ---

    salinity = get_salinity(salinity_surface, salinity_gradient, depth)

    pressure = get_hydrostatic_pressure(depth)

    Aseg = Aseg1 - Aseg2

    Vseg = Aseg * thickness

    Eseg = 1e-9 * Vseg * (

        (1 - porosity) * cp_rock * rho_rock +

        porosity * cp_water * rho_water

    Eseg = (

        1e-9

        * Vseg

        * ((1 - porosity) * cp_rock * rho_rock + porosity * cp_water * rho_water)

    )

    flowseg = flowrate * 365 * 24 * dangle / math.pi

    curve_scaling_factor: float,

    max_true_vertical_depth_stepout_factor: float,

) -> float:

    if curve_scaling_factor == 0:

        return true_vertical_depth

    stepout -= horizontal_distance

    oblique_distance = math.sqrt(stepout**2 + true_vertical_depth**2) * curve_scaling_factor

    oblique_distance = (

        math.sqrt(stepout**2 + true_vertical_depth**2) * curve_scaling_factor

    )

    return oblique_distance + horizontal_distance

INCH_SI = 0.0254

@dataclass

class Doublet1DResults:

    geothermal_powers: float

    production_temp: float

    heat_power_produced: list[float]

def doubletcalc(

    props: UTCConfiguration,

    input,

        porosity=input.porosity,

        net_to_gross=input.ntg,

        thickness=input.thickness,

        salinity=get_salinity(props.salinity_surface, props.salinity_gradient, input.depth),

        salinity=get_salinity(

            props.salinity_surface, props.salinity_gradient, input.depth

        ),

        surface_temperature=props.surface_temperature,

        geothermal_gradient=get_geothermal_gradient(input.depth, input.thickness, input.temperature, props.surface_temperature),

        geothermal_gradient=get_geothermal_gradient(

            input.depth, input.thickness, input.temperature, props.surface_temperature

        ),

        rock_heat_capacity=855,

        rock_density=2720,

        thermal_conductivity=3,

    producer = Well(

        aquifer=aquifer,

        well_type="producer",

        pipe_segments=[WellPipeSegment(

            ah_depth=get_along_hole_length(input.depth, well_distance, props.well_curv_scaling, props.max_tvd_stepout_factor),

        pipe_segments=[

            WellPipeSegment(

                ah_depth=get_along_hole_length(

                    input.depth,

                    well_distance,

                    props.well_curv_scaling,

                    props.max_tvd_stepout_factor,

                ),

                tv_depth=input.depth,

                inner_diameter=props.inner_diameter * INCH_SI,

            roughness=props.roughness * 1e-3 * INCH_SI

        )],

                roughness=props.roughness * 1e-3 * INCH_SI,

            )

        ],

        aquifer_top_depth=input.depth,

        pipe_scaling=0,

        target_segment_length=props.segment_length,

        flowrate=doublet.flowrate,

        pump_power_required=doublet.pump_power,

        production_temp=doublet.production_temp,

        heat_power_produced=[doublet.geothermal_power,] * props.econ_lifetime_years,

        heat_power_produced=[

            doublet.geothermal_power,

        ]

        * props.econ_lifetime_years,

    )

def get_total_skin_injection(props, input):

    if (not props.use_stimulation) or (input.transmissivity_with_ntg > props.stim_kh_max):

    if (not props.use_stimulation) or (

        input.transmissivity_with_ntg > props.stim_kh_max

    ):

        stim_add_skin_inj = 0.0

    else:

        stim_add_skin_inj = props.stim_add_skin_inj()

    return props.skin_injector + stim_add_skin_inj

def get_total_skin_production(props, input):

    if (not props.use_stimulation) or (input.transmissivity_with_ntg > props.stim_kh_max):

    if (not props.use_stimulation) or (

        input.transmissivity_with_ntg > props.stim_kh_max

    ):

        stim_add_skin_prod = 0.0

    else:

        stim_add_skin_prod = props.stim_add_skin_prod

    return props.skin_producer + stim_add_skin_prod

def get_pump_production_depth(props, depth: float) -> float:

    return min(props.pump_depth, depth / 2)

@njit

def get_along_hole_length(

    true_vertical_depth: float,

from dataclasses import dataclass

import numpy as np

from dataclasses import dataclass

from pythermogis.workflow.utc.doublet_utils import get_along_hole_length

SECONDS_IN_YEAR = 365 * 24 * 3600

KWH_TO_MJ = 0.36 * 1e9

HOURS_IN_YEAR = 8760

@dataclass

class CapexCalculatorResults:

    sum_capex: float

    total_opex_ts: list[float]

    heat_power_per_year: list[float]

@dataclass

class UTCCalculatorResults:

    discounted_heat_produced: float

    utc: float

@dataclass

class EconomicsResults:

    capex: CapexCalculatorResults

    utc: UTCCalculatorResults

def calculate_economics(

    props,

    input,

    hp_added_power: float,

    pump_power_required: float,

) -> EconomicsResults:

    ah_length_array = [

        get_along_hole_length(

            true_vertical_depth=input.depth,

        utc=utc_results,

    )

def calculate_capex(

    props,

    heat_power_produced: list[float],

    ah_length_array: list[float],

    pump_power_required: float,

) -> CapexCalculatorResults:

    n_years = props.economic_lifetime

    heat_power_per_year = np.zeros(n_years)

    capex_well = calculate_capex_for_wells(props, ah_length_array, stimulation_capex)

    allocate_capex_for_wells(props, total_capex_ts, capex_well)

    hp_added_power = get_max_hp_added_power(initial_hp_added_power, hp_added_power_years)

    hp_after_hp = calculate_hp_added_power_after_heat_pump(props, hp_added_power, hp_cop)

    hp_added_power = get_max_hp_added_power(

        initial_hp_added_power, hp_added_power_years

    )

    hp_after_hp = calculate_hp_added_power_after_heat_pump(

        props, hp_added_power, hp_cop

    )

    installation_mw = calculate_installation_mw(shifted_heat_power, hp_added_power)

    total_capex_1year = calculate_total_capex_1year(

        for i in range(shift):

            series[i] = 0

def calculate_capex_for_wells(props, ah_length_array, stimulation_capex):

    inj_depth = ah_length_array[0]

    prod_depth = ah_length_array[0]

    )

    return capex * props.well_cost_scaling

def calculate_well_cost(props, depth):

    return (

        props.well_cost_z2 * depth**2

        + props.well_cost_z * depth

        props.well_cost_z2 * depth**2 + props.well_cost_z * depth

    ) * 1e-6 + props.well_cost_const

def allocate_capex_for_wells(props, total_capex_ts, capex_well):

    drilling_years = props.drilling_time

    yearly_cost = capex_well / max(drilling_years, 1)

    for y in range(drilling_years):

        total_capex_ts[y] += yearly_cost

def get_max_hp_added_power(initial, hp_power_years):

    if hp_power_years is not None:

        return max(hp_power_years)

    return initial

def calculate_hp_added_power_after_heat_pump(props, hp_added_power, hp_cop):

    alt_price = props.hp_alternative_heating_price

    return hp_added_power if alt_price < 0 else hp_added_power * (1 + 1 / (hp_cop - 1))

def calculate_installation_mw(heat_power_produced, hp_added_power):

    last_val = heat_power_produced[-1]

    return max(last_val - hp_added_power, 0)

def calculate_total_capex_1year(props, total_capex_ts, hp_after_hp, installation_mw):

    last_year = max(props.drilling_time - 1, 0)

        + (props.hp_capex * hp_after_hp + props.capex_variable * installation_mw) * 1e-3

    )

    total_capex_ts[last_year] *= (1 + props.capex_contingency / 100)

    total_capex_ts[last_year] *= 1 + props.capex_contingency / 100

    return total_capex_ts[last_year]

def process_annual_data(

    props,

    heat_power_produced,

    season_factor_years,

    pump_power_required,

):

    sum_capex = 0.0

    for year in range(props.economic_lifetime):

        # Season factor

        season_factor = (

            season_factor_years[year]

        # Annual heat (GJ/yr)

        heat_power_per_year[year] = (

            heat_power_produced[year]

            * season_factor

            * SECONDS_IN_YEAR

            * MJ_TO_GJ

            heat_power_produced[year] * season_factor * SECONDS_IN_YEAR * MJ_TO_GJ

        )

        sum_capex += total_capex_ts[year]

    return sum_capex

def calculate_variable_opex(

    props,

    season_factor,

):

    pump_cost = (

        -(pump_power * 1e3)

        * SECONDS_IN_YEAR * season_factor * elec_price

        / KWH_TO_MJ * 1e-6

        * SECONDS_IN_YEAR

        * season_factor

        * elec_price

        / KWH_TO_MJ

        * 1e-6

    )

    parasitic_cost = (

        -(heat_power_produced * heat_exchanger_parasitic * 1e6)

        * SECONDS_IN_YEAR * season_factor * elec_price

        / KWH_TO_MJ * 1e-6

        * SECONDS_IN_YEAR

        * season_factor

        * elec_price

        / KWH_TO_MJ

        * 1e-6

    )

    opex_per_energy = (

        -(1 / 0.36)

        * props.opex_per_energy * 1e-2

        * heat_power_gj_per_year * 1e4 / 36 * 1e-6

        * props.opex_per_energy

        * 1e-2

        * heat_power_gj_per_year

        * 1e4

        / 36

        * 1e-6

    )

    heat_pump_cost = 0.0

    if hp_cop != 0:

        heat_pump_cost = (

            -(hp_after_hp * 1e6 / hp_cop)

            * SECONDS_IN_YEAR * season_factor * elec_price

            / KWH_TO_MJ * 1e-6

            * SECONDS_IN_YEAR

            * season_factor

            * elec_price

            / KWH_TO_MJ

            * 1e-6

        )

    return pump_cost + parasitic_cost + heat_pump_cost + opex_per_energy

def calculate_fixed_opex(props, hp_after_hp, installation_mw, total_capex):

    return (

        -props.opex_base / 1e6

        - total_capex * props.opex_per_capex / 100

    )

def get_heat_exchanger_season_factor(props):

    return props.load_hours / HOURS_IN_YEAR

    total_opex_ts,

    total_capex,

) -> UTCCalculatorResults:

    net_cash_worth_eia = calculate_net_cash_worth_eia(props, total_capex)

    present_value = total_capex * props.debt_equity - net_cash_worth_eia

    yearly_payment = calculate_payment_value(props, present_value)

    discounted_heat_produced = 0.0

    for year in range(1, props.economic_lifetime):

        inflated_opex = total_opex_ts[year] * ((1 + props.inflation) ** (year - 1))

        future_value = calculate_future_value(

        utc=utc,

    )

def calculate_net_cash_worth_eia(props, total_capex):

    tax_rate = props.tax_rate

    project_interest_rate = calculate_project_interest_rate(props)

    return (tax_rate * (props.ecn_eia * total_capex)) / (1 + project_interest_rate)

def calculate_payment_value(props, present_value):

    if props.interest_loan == 0:

        return 0.0

    factor = (1 + props.interest_loan) ** props.economic_lifetime

    return (props.interest_loan / (factor - 1)) * (-(present_value * factor))

def calculate_depreciation_cost(props, total_capex):

    return -(total_capex / props.economic_lifetime)

def calculate_future_value(present_value, payment, interest, year):

    future_value = present_value * ((1 + interest) ** (year - 1))