"""Restricted Paardekooper2010v1 torque prescription, independently implemented.

G=Mstar=AU=1; time unit sqrt(AU^3/GMstar), not one year.
This module is public scientific material, not an interactive Lab Engine.
"""

import math

EARTH_SOLAR_MASS_RATIO = 3.0035e-6
SECONDS_PER_CODE_TIME = math.sqrt(149597870700.0**3 / 1.3271244004193938e20)
YEARS_PER_CODE_TIME = SECONDS_PER_CODE_TIME / (86400 * 365.25)

DEFAULT_DISK = {
    "surfaceDensityAtOneAu": 1e-4,
    "densitySlope": 0.5,
    "aspectRatioAtOneAu": 0.055,
    "temperatureSlope": 0.5,
    "viscosity": 1e-6,
    "prandtlNumber": 0.5,
    "adiabaticGamma": 1.4,
    "temperatureTransition": False,
    "transitionRadiusAu": 1.2,
    "transitionWidthLogRadius": 0.5,
}


def effective_gamma(thermal_parameter, adiabatic_gamma=1.4):
    """Equation46, rationalized to retain the Q->0 limit without cancellation."""
    if thermal_parameter <= 0:
        return adiabatic_gamma
    gamma = adiabatic_gamma
    squared = thermal_parameter**2
    radical = math.sqrt((gamma * gamma * squared + 1) ** 2 - 16 * squared * (gamma - 1))
    divided_inner = 2 * (
        (2 * gamma * gamma - 16 * (gamma - 1) + gamma**4 * squared) / (radical + 1)
        + gamma * gamma
    )
    return 2 * gamma / (gamma + 0.5 * math.sqrt(max(0, divided_inner)))


def saturation_function(parameter):
    return 1 / (1 + (parameter / 1.3) ** 2)


def diffusion_cutoff(parameter, turnover_parameter):
    threshold = math.sqrt(4 * turnover_parameter / (45 * math.pi))
    if parameter < threshold:
        return (
            16 / 25 * (45 * math.pi / (4 * turnover_parameter)) ** 0.75 * parameter**1.5
        )
    return 1 - 9 / 25 * (4 * turnover_parameter / (45 * math.pi)) ** (
        4 / 3
    ) * parameter ** (-8 / 3)


def local_disk(radius, configuration):
    density_slope = configuration["densitySlope"]
    slope = configuration["temperatureSlope"]
    log_temperature = -slope * math.log(radius)
    if configuration["temperatureTransition"]:
        # A smooth declared temperature profile, not an opacity/thermal-equilibrium solution.
        inner_slope, outer_slope = 1.4, 0.4
        midpoint = (inner_slope + outer_slope) / 2
        contrast = (inner_slope - outer_slope) / 2
        width = configuration["transitionWidthLogRadius"]
        centre = configuration["transitionRadiusAu"]
        coordinate = math.log(radius / centre) / width
        reference = math.log(1 / centre) / width
        slope = midpoint - contrast * math.tanh(coordinate)
        log_temperature = -midpoint * math.log(radius) + contrast * width * (
            math.log(math.cosh(coordinate)) - math.log(math.cosh(reference))
        )
    aspect = configuration["aspectRatioAtOneAu"] * math.sqrt(
        radius * math.exp(log_temperature)
    )
    density = configuration["surfaceDensityAtOneAu"] * radius ** (-density_slope)
    frequency = radius ** (-1.5)
    return density, aspect, frequency, slope


def torque_components(planet_mass, semimajor_axis, eccentricity, configuration):
    density, aspect, frequency, temperature_slope = local_disk(
        semimajor_axis, configuration
    )
    if density == 0:
        return {
            "wave": 0.0,
            "corotation": 0.0,
            "total": 0.0,
            "normalization": 0.0,
            "eccentricityTime": math.inf,
            "aspect": aspect,
            "thermalMassRatio": planet_mass / aspect**3,
            "toomreQ": math.inf,
            "effectiveGamma": configuration["adiabaticGamma"],
            "viscousParameter": 0.0,
            "thermalParameter": 0.0,
        }
    viscosity = configuration["viscosity"]
    diffusivity = viscosity / configuration["prandtlNumber"]
    gamma = configuration["adiabaticGamma"]
    gamma_effective = effective_gamma(
        2 * diffusivity / (3 * aspect**3 * semimajor_axis**2 * frequency), gamma
    )
    horseshoe_width = 1.1 / gamma_effective**0.25 * math.sqrt(planet_mass / aspect)
    viscous_parameter = (2 / 3) * math.sqrt(
        semimajor_axis**2 * frequency * horseshoe_width**3 / (2 * math.pi * viscosity)
    )
    thermal_parameter = math.sqrt(
        semimajor_axis**2 * frequency * horseshoe_width**3 / (2 * math.pi * diffusivity)
    )
    viscous_f = saturation_function(viscous_parameter)
    thermal_f = saturation_function(thermal_parameter)
    viscous_g = diffusion_cutoff(viscous_parameter, 2)
    thermal_g = diffusion_cutoff(thermal_parameter, 2)
    viscous_k = diffusion_cutoff(viscous_parameter, 7)
    thermal_k = diffusion_cutoff(thermal_parameter, 7)
    density_slope = configuration["densitySlope"]
    entropy_slope = temperature_slope - (gamma - 1) * density_slope
    wave = (-2.5 - 1.7 * temperature_slope + 0.1 * density_slope) / gamma_effective
    barotropic = (
        (1.1 * viscous_f * viscous_g + 0.7 * (1 - viscous_k))
        * (1.5 - density_slope)
        / gamma_effective
    )
    entropy = (
        7.9
        * entropy_slope
        / gamma_effective**2
        * viscous_f
        * thermal_f
        * math.sqrt(viscous_g * thermal_g)
    )
    entropy += (
        (2.2 - 1.4 / gamma_effective)
        * entropy_slope
        / gamma_effective
        * math.sqrt((1 - viscous_k) * (1 - thermal_k))
    )
    scaled_eccentricity = eccentricity / aspect
    wave *= (1 - (scaled_eccentricity / 2.02) ** 4) / (
        1 + (scaled_eccentricity / 2.25) ** 1.2 + (scaled_eccentricity / 2.84) ** 6
    )
    corotation = (barotropic + entropy) * math.exp(
        -eccentricity / (0.5 * aspect + 0.01)
    )
    normalization = (
        (planet_mass / aspect) ** 2 * density * semimajor_axis**4 * frequency**2
    )
    wave_time = aspect**4 / (planet_mass * density * semimajor_axis**2 * frequency)
    eccentricity_time = (
        wave_time
        / 0.780
        * (1 - 0.14 * scaled_eccentricity**2 + 0.06 * scaled_eccentricity**3)
    )
    return {
        "wave": wave * normalization,
        "corotation": corotation * normalization,
        "total": (wave + corotation) * normalization,
        "normalization": normalization,
        "eccentricityTime": eccentricity_time,
        "aspect": aspect,
        "thermalMassRatio": planet_mass / aspect**3,
        "toomreQ": aspect / (math.pi * density * semimajor_axis**2),
        "effectiveGamma": gamma_effective,
        "viscousParameter": viscous_parameter,
        "thermalParameter": thermal_parameter,
    }