import math

# Общие константы
DELTA_CONST = 1e-6
MIN_PART_DELTA_VAL = 1e-4
MIN_DUML_INGRL_VAL = 1e-4

# Единичный интеграл дюамеля
class DuhamelPart():

    def __init__(self, _beg_time, _amp, _t_per):
        self.beg_time = _beg_time
        self.amp = _amp
        self.t_per = _t_per
        self.is_const = False

    def part_step(self, _calc_time):
        if self.is_const:
            return self.amp
        exp_val = math.exp( (self.beg_time - _calc_time)/self.t_per )
        self.is_const = exp_val < DELTA_CONST
        return self.amp * (1. - exp_val)

    def const_part(self):
        return self.is_const


# совокупность интегралов дюамеля
class DuhamelIntegral():

    def __init__(self, _up_per = 1.0, _down_per = 20.0):
        # параметры интеграла
        # вроде как вначале одинаковые для всех
        self.curr_time = -1.0
        self.curr_src_val = 0.0

        # непосредственно влияющие на интеграл
        self.up_per = _up_per
        self.down_per = _down_per

        self.parts_set = []
        self.duml_val = 0.0

    def add_part(self, delta):
        if abs(delta) < MIN_PART_DELTA_VAL:
            return
        if delta > 0:
            self.parts_set.append(DuhamelPart(self.curr_time, delta, self.up_per))
        else:
            self.parts_set.append(DuhamelPart(self.curr_time, delta, self.down_per))

    def calc_step(self, calc_time):
        self.duml_val = 0.
        for part in self.parts_set:
            self.duml_val += part.part_step(calc_time)

        return self.duml_val

    def duml_step(self, time, value):
        self.add_part(value - self.curr_src_val)
        self.curr_time = time
        self.curr_src_val = value
        return self.calc_step(self.curr_time)

    def clear_parts_set(self):
        if abs(self.duml_val) >= MIN_DUML_INGRL_VAL:
            return
        for part in self.parts_set:
            if not part.const_part():
                return
        self.parts_set = []

# Полная модель интегралов дюамеля с учетом корректирующих коэффициентов
class DumlRiverModel:

    def __init__(self, _k_n = 1.0, _k_korr = 25.0, _res_add_val = 20.0, _dmi_up_per = 1.0, _dmi_down_per = 20.0):

        self.duml_intgrl = DuhamelIntegral(_down_per=_dmi_down_per, _up_per=_dmi_up_per)

        self.k_n = _k_n             # константный коэффициент коррекции результата
        self.save_k_n = self.k_n

        self.k_korr = _k_korr       # настраиваемый коэффициент коррекции результата
        self.save_k_korr = self.k_korr

        self.res_add_val = _res_add_val # аддитивная величина коррекции результата
        self.save_add_val = self.res_add_val

    def model_step(self, time, value):
        return self.k_n * self.k_korr * self.duml_intgrl.duml_step(time, value) + self.res_add_val

# Расчёт стока с водосборов
class CatchmentArea():
    def __init__(self, _capacity = 22.2, _in_filter = 1.0, _out_filter = 0.009):
        # параметры водосбора
        self.capacity = _capacity  # вместимость водосбора, мм
        self.in_filter = _in_filter  # инфильтрация в почву за ед. периода, мм (за 10 мин)
        self.out_filter = _out_filter  # фильтрация из почвы за ед. периода, мм (за 10 мин)
        self.level_delta = self.capacity

    def calc_surface_flow(self, _rainfall):
        мах_in = min(_rainfall, self.in_filter)
        if мах_in <= self.level_delta:
            self.level_delta -= мах_in
            self.level_delta += self.out_filter
            if self.level_delta > self.capacity:
                self.level_delta = self.capacity
        else:
            self.level_delta += self.out_filter
            if self.level_delta > self.capacity:
                self.level_delta = self.capacity
            мах_in = min(мах_in, self.level_delta)
            self.level_delta -= мах_in

        return _rainfall - мах_in

# реализация метода ньютона-рафсона
def newton_raphson(study_func, target_val, start_x, max_iter_count=10, max_err=1e-3, delta_x=1e-6):

    # начальное значение аргумента
    finded_x = start_x

    # Цикл нахождения решения
    iter_count = 0
    err_y = 1e20

    # цикл поиска решения
    while (iter_count <= max_iter_count) and (abs(err_y) > max_err):

        # Вычисляем производную dy/dx в точке x
        # dfdx в принципе нам может быть известен аналитически, и в таком случае
        # предпочтительно использовать аналитическое выражение.
        dfdx = (study_func(finded_x + 0.5 * delta_x) - study_func(finded_x - 0.5 * delta_x)) / delta_x

        # Находим приращение по y
        err_y = study_func(finded_x) - target_val
        # и по x
        err_x = err_y / dfdx
        # Получаем новое решение
        finded_x = finded_x - err_x

        # меняем к - во итераций
        iter_count += 1

    return finded_x

# для расчёта уровня воды по расходу в створе
class RiverCrossSection():

    def __init__(self, _left_slope = 5.0, _bottom = 18.0, _right_slope = 5.0, _max_h = 3.6, _zero_bsv = -1.06, _n_roughness = 0.035, _i_bias = 0.00244):
        # параметры трапеции
        self._zero_wlevel = _zero_bsv
        self._trap_a = _left_slope  # ширина левого склона трапеции русла
        self._trap_b = _bottom  # ширина дна трапеции русла
        self._trap_c = _right_slope # ширина правого склона трапеции русла
        self._trap_d = _max_h  # высота трапеции русла
        if self._trap_d >= 0:
            self._trap_leftCtg = self._trap_a / self._trap_d  # котангенс левого угла трапеции русла
            self._trap_rightCtg = self._trap_c / self._trap_d  # котангенс правого угла трапеции русла
        else:
            self._trap_leftCtg = 0  # котангенс левого угла трапеции русла
            self._trap_rightCtg = 0  # котангенс правого угла трапеции русла

        self._i_rb_part = _i_bias  # уклон русла
        self._n_rb_part = _n_roughness  # шероховатость русла

    # площадь сечения (м^2)
    def f_trap_area(self, h):
      return self._trap_b * h + (self._trap_leftCtg + self._trap_rightCtg) * h * h / 2

    # смоченный периметр (м)
    def f_trap_perimetr(self, h):
      return self._trap_b + \
             math.sqrt(self._trap_leftCtg ** 2 + 1) * h + \
             math.sqrt(self._trap_rightCtg ** 2 + 1) * h

    # гидрологический радиус
    def f_hydro_r(self, h):
        return self.f_trap_area(h) / self.f_trap_perimetr(h)

    # скорость потока (м/с)
    def f_flow_speed(self, h):
        return self.f_hydro_r(h)**(2./3) * math.sqrt(self._i_rb_part) / self._n_rb_part

    # расход воды (м^3/с)
    def f_flow_rate(self, h):
        if h <= 0.:
            return 0.
        return self.f_trap_area(h) * self.f_flow_speed(h)

    # расход воды (м^3/с)
    def f_flow_rate_by_bs_h(self, bs_h):
        return self.f_flow_rate(bs_h - self._zero_wlevel)

    # уровень воды (м), определяемый по расходу воды (м^3/с)
    def high_by_flow_rate(self, flow_rate):
        if flow_rate <= 0.:
            return 0.
        return newton_raphson(self.f_flow_rate, flow_rate, 1.)

    # уровень воды (м), определяемый по расходу воды (м^3/с) по балтийской системе высот
    def bs_water_level(self, flow_rate):
        return self.high_by_flow_rate(flow_rate) + self._zero_wlevel