import numpy as np
from input.user_input import *


def cd_Palumbo(fr, re, a_top):
    k_cd = 19.8                                          # +/- 1.2
    k1 = 1.50 * 10 ** 7                                  # +/- 0.78
    k2 = 0.86 * 10 ** (-7)                               # +/- 0.43
    a_top0 = np.pi / 4 * 1.383 ** 2 * V_design ** (2/3)  # a_top0 for zero-pressure shape

    value = 0.2 * k_cd/fr * (k1/re + k2 * re) * a_top / (a_top0 ** (3/2))

    if value > 1.6:
        value = 1.6
    elif value <= 0.1:
        value = 0.1

    return value


def cd_PalumboLow(fr, re, a_top):
    k_cd = 19.8 - 1.2
    k1 = (1.50 - 0.78) * 10 ** 7
    k2 = (0.86 - 0.43) * 10 ** (-7)
    a_top0 = np.pi / 4 * 1.383 ** 2 * V_design ** (2/3)  # a_top0 for zero-pressure shape

    value = 0.2 * k_cd/fr * (k1/re + k2 * re) * a_top / (a_top0 ** (3/2))

    if value > 1.6:
        value = 1.6
    elif value <= 0.1:
        value = 0.1

    return value


def cd_PalumboHigh(fr, re, a_top):
    k_cd = 19.8 + 1.2
    k1 = (1.50 + 0.78) * 10 ** 7
    k2 = (0.86 + 0.43) * 10 ** (-7)
    a_top0 = np.pi / 4 * 1.383 ** 2 * V_design ** (2/3)  # a_top0 for zero-pressure shape

    value = 0.2 * k_cd/fr * (k1/re + k2 * re) * a_top / (a_top0 ** (3/2))

    if value > 1.6:
        value = 1.6
    elif value <= 0.1:
        value = 0.4

    return value


def drag(c_d, rho_air, A_proj, v):
    return 0.5 * c_d * rho_air * A_proj * v ** 2


def cd_sphere(Re):
    "This function computes the drag coefficient of a sphere as a function of the Reynolds number Re."
    # Curve fitted after fig . A -56 in Evett and Liu: "Fluid Mechanics and Hydraulics"

    from numpy import log10, array, polyval

    if Re <= 0.0:
        CD = 0.0
    elif Re > 8.0e6:
        CD = 0.2
    elif Re > 0.0 and Re <= 0.5:
        CD = 24.0 / Re
    elif Re > 0.5 and Re <= 100.0:
        p = array([4.22, -14.05, 34.87, 0.658])
        CD = polyval(p, 1.0 / Re)
    elif Re > 100.0 and Re <= 1.0e4:
        p = array([-30.41, 43.72, -17.08, 2.41])
        CD = polyval(p, 1.0 / log10(Re))
    elif Re > 1.0e4 and Re <= 3.35e5:
        p = array([-0.1584, 2.031, -8.472, 11.932])
        CD = polyval(p, log10(Re))
    elif Re > 3.35e5 and Re <= 5.0e5:
        x1 = log10(Re / 4.5e5)
        CD = 91.08 * x1 ** 4 + 0.0764
    else:
        p = array([-0.06338, 1.1905, -7.332, 14.93])
        CD = polyval(p, log10(Re))
    return CD