//! Rain attenuation.
//!
//! The specific attenuation of rain is modelled as a power law in the rain
//! rate, `gamma = k R^alpha`, where `k` and `alpha` depend on the frequency,
//! on the polarisation, and on the elevation. The coefficients are given by
//! ITU-R Recommendation P.838 as closed form curve fits rather than as a
//! table, which is what the current revision of that recommendation uses.
//!
//! The attenuation along a slant path follows ITU-R Recommendation P.618.
//! The path length through the rain is limited by the height of the rain
//! layer above the station, and two reduction factors account for the fact
//! that the rain rate is not uniform along the path: a horizontal reduction
//! and a vertical adjustment. Both depend on the specific attenuation, so
//! the calculation is not a simple product of a rate and a length.
//!
//! Sources:
//!
//! - Specific attenuation coefficients: ITU-R Recommendation P.838-3,
//!   *Specific attenuation model for rain for use in prediction methods*,
//!   equations 2 and 3 and tables 1 and 2.
//! - Slant path method: ITU-R Recommendation P.618-13, *Propagation data and
//!   prediction methods required for the design of Earth-space
//!   telecommunication systems*, equations 1 through 10.

use apsis_rf::Decibel;

/// The effective radius of the Earth used by the path geometry, in
/// kilometers.
///
/// The value accounts for the refraction of the ray in the lower atmosphere,
/// which makes the Earth appear larger than it is.
const EFFECTIVE_EARTH_RADIUS_KM: f64 = 8_500.0;

/// A mid latitude rain height, in kilometers above sea level.
///
/// Recommendation P.618 takes the rain height from the zero degree isotherm
/// height of Recommendation P.839, which is tabulated as a function of
/// latitude and longitude. Prefer [`rain_attenuation_at`], which looks the
/// value up, over [`rain_attenuation`] with this constant.
pub const DEFAULT_RAIN_HEIGHT_KM: f64 = 3.0;

/// One term of a P.838-3 curve fit.
#[derive(Debug, Clone, Copy)]
struct CurveTerm {
    /// The amplitude coefficient.
    amplitude: f64,
    /// The centre of the term, on a base ten logarithmic frequency axis.
    centre: f64,
    /// The width of the term, on a base ten logarithmic frequency axis.
    width: f64,
}

/// The coefficients of the `k` curve fit for horizontal polarisation.
const K_HORIZONTAL_TERMS: &[CurveTerm] = &[
    CurveTerm {
        amplitude: -5.339_80,
        centre: -0.100_08,
        width: 1.130_98,
    },
    CurveTerm {
        amplitude: -0.353_51,
        centre: 1.269_70,
        width: 0.454_00,
    },
    CurveTerm {
        amplitude: -0.237_89,
        centre: 0.860_36,
        width: 0.153_54,
    },
    CurveTerm {
        amplitude: -0.941_58,
        centre: 0.645_52,
        width: 0.168_17,
    },
];

/// The coefficients of the `k` curve fit for vertical polarisation.
const K_VERTICAL_TERMS: &[CurveTerm] = &[
    CurveTerm {
        amplitude: -3.805_95,
        centre: 0.569_34,
        width: 0.810_61,
    },
    CurveTerm {
        amplitude: -3.449_65,
        centre: -0.229_11,
        width: 0.510_59,
    },
    CurveTerm {
        amplitude: -0.399_02,
        centre: 0.730_42,
        width: 0.118_99,
    },
    CurveTerm {
        amplitude: 0.501_67,
        centre: 1.073_19,
        width: 0.271_95,
    },
];

/// The coefficients of the `alpha` curve fit for horizontal polarisation.
const ALPHA_HORIZONTAL_TERMS: &[CurveTerm] = &[
    CurveTerm {
        amplitude: -0.143_18,
        centre: 1.824_42,
        width: -0.551_87,
    },
    CurveTerm {
        amplitude: 0.295_91,
        centre: 0.775_64,
        width: 0.198_22,
    },
    CurveTerm {
        amplitude: 0.321_77,
        centre: 0.637_73,
        width: 0.131_64,
    },
    CurveTerm {
        amplitude: -5.376_10,
        centre: -0.962_30,
        width: 1.478_28,
    },
    CurveTerm {
        amplitude: 16.172_1,
        centre: -3.299_80,
        width: 3.439_90,
    },
];

/// The coefficients of the `alpha` curve fit for vertical polarisation.
const ALPHA_VERTICAL_TERMS: &[CurveTerm] = &[
    CurveTerm {
        amplitude: -0.077_71,
        centre: 2.338_40,
        width: -0.762_84,
    },
    CurveTerm {
        amplitude: 0.567_27,
        centre: 0.955_45,
        width: 0.540_39,
    },
    CurveTerm {
        amplitude: -0.202_38,
        centre: 1.145_20,
        width: 0.268_09,
    },
    CurveTerm {
        amplitude: -48.299_1,
        centre: 0.791_669,
        width: 0.116_226,
    },
    CurveTerm {
        amplitude: 48.583_3,
        centre: 0.791_459,
        width: 0.116_479,
    },
];

/// The linear coefficients of the `k` curve fits.
const K_SLOPE_HORIZONTAL: f64 = -0.189_61;
/// The constant of the horizontal `k` curve fit.
const K_OFFSET_HORIZONTAL: f64 = 0.711_47;
/// The linear coefficient of the vertical `k` curve fit.
const K_SLOPE_VERTICAL: f64 = -0.163_98;
/// The constant of the vertical `k` curve fit.
const K_OFFSET_VERTICAL: f64 = 0.632_97;
/// The linear coefficient of the horizontal `alpha` curve fit.
const ALPHA_SLOPE_HORIZONTAL: f64 = 0.678_49;
/// The constant of the horizontal `alpha` curve fit.
const ALPHA_OFFSET_HORIZONTAL: f64 = -1.955_37;
/// The linear coefficient of the vertical `alpha` curve fit.
const ALPHA_SLOPE_VERTICAL: f64 = -0.053_739;
/// The constant of the vertical `alpha` curve fit.
const ALPHA_OFFSET_VERTICAL: f64 = 0.834_33;

/// A rain attenuation result at a given elevation.
#[derive(Debug, Clone, Copy, PartialEq)]
pub struct RainAttenuation {
    /// The specific attenuation, in decibels per kilometer.
    pub specific_db_per_km: f64,
    /// The effective path length through the rain, in kilometers.
    pub path_length_km: f64,
    /// The total attenuation along the path, in decibels.
    pub attenuation_db: f64,
    /// The horizontal reduction factor, dimensionless and at most one.
    pub horizontal_reduction: f64,
    /// The vertical adjustment factor, dimensionless.
    ///
    /// The factor is the ratio of the effective path length to the longer
    /// of the two geometric path estimates. It is at most one when the path
    /// meets the rain layer at a steep angle and can exceed one at a shallow
    /// angle, which is the behaviour Recommendation P.618 gives.
    pub vertical_adjustment: f64,
}

impl RainAttenuation {
    /// Returns the total attenuation.
    pub fn attenuation(self) -> Decibel {
        Decibel::from_db(self.attenuation_db)
    }

    /// Returns `true` when every value is finite and not negative.
    pub fn is_valid(self) -> bool {
        self.specific_db_per_km.is_finite()
            && self.specific_db_per_km >= 0.0
            && self.path_length_km.is_finite()
            && self.path_length_km >= 0.0
            && self.attenuation_db.is_finite()
            && self.attenuation_db >= 0.0
    }
}

/// Evaluates one P.838-3 curve term.
///
/// The term is a Gaussian on a base ten logarithmic frequency axis, so the
/// frequency must be in gigahertz and positive.
fn curve_term(term: CurveTerm, frequency_ghz: f64) -> f64 {
    let logarithm = frequency_ghz.log10();
    term.amplitude * (-((logarithm - term.centre) / term.width).powi(2)).exp()
}

/// Returns the `k` for a polarisation, before the elevation correction.
fn base_k(terms: &[CurveTerm], slope: f64, offset: f64, frequency_ghz: f64) -> f64 {
    let gaussian: f64 = terms
        .iter()
        .map(|term| curve_term(*term, frequency_ghz))
        .sum();
    10f64.powf(gaussian + slope * frequency_ghz.log10() + offset)
}

/// Returns the `alpha` for a polarisation, before the elevation correction.
fn base_alpha(terms: &[CurveTerm], slope: f64, offset: f64, frequency_ghz: f64) -> f64 {
    let gaussian: f64 = terms
        .iter()
        .map(|term| curve_term(*term, frequency_ghz))
        .sum();
    gaussian + slope * frequency_ghz.log10() + offset
}

/// Returns the P.838-3 coefficients of the specific attenuation power law.
///
/// The polarisation is an angle in radians between the electric field and
/// the horizontal, so zero is horizontal and a right angle is vertical. The
/// elevation is in radians.
///
/// The result is a pair of the coefficient `k` and the exponent `alpha` of
/// the relation `gamma = k R^alpha`, with the rain rate in millimeters per
/// hour and the result in decibels per kilometer.
///
/// Returns `None` when the frequency, polarisation, or elevation is not
/// usable. The recommendation is stated for frequencies from 1 GHz to
/// 1000 GHz.
///
/// Source: ITU-R Recommendation P.838-3, equations 2 and 3.
pub fn rain_coefficients(frequency: f64, polarisation: f64, elevation: f64) -> Option<(f64, f64)> {
    let frequency_ghz = frequency / 1.0e9;
    if !frequency_ghz.is_finite() || !(1.0..=1_000.0).contains(&frequency_ghz) {
        return None;
    }
    if !polarisation.is_finite() || !elevation.is_finite() {
        return None;
    }

    let k_horizontal = base_k(
        K_HORIZONTAL_TERMS,
        K_SLOPE_HORIZONTAL,
        K_OFFSET_HORIZONTAL,
        frequency_ghz,
    );
    let k_vertical = base_k(
        K_VERTICAL_TERMS,
        K_SLOPE_VERTICAL,
        K_OFFSET_VERTICAL,
        frequency_ghz,
    );
    let alpha_horizontal = base_alpha(
        ALPHA_HORIZONTAL_TERMS,
        ALPHA_SLOPE_HORIZONTAL,
        ALPHA_OFFSET_HORIZONTAL,
        frequency_ghz,
    );
    let alpha_vertical = base_alpha(
        ALPHA_VERTICAL_TERMS,
        ALPHA_SLOPE_VERTICAL,
        ALPHA_OFFSET_VERTICAL,
        frequency_ghz,
    );

    // The elevation and polarisation enter through the square of the cosine
    // of the elevation angle and the cosine of twice the polarisation tilt.
    // The elevation term vanishes at the zenith, where the distinction
    // between horizontal and vertical polarisation disappears.
    let geometry_factor = elevation.cos().powi(2) * (2.0 * polarisation).cos();
    let k = 0.5 * (k_horizontal + k_vertical + (k_horizontal - k_vertical) * geometry_factor);
    let alpha = (k_horizontal * alpha_horizontal
        + k_vertical * alpha_vertical
        + (k_horizontal * alpha_horizontal - k_vertical * alpha_vertical) * geometry_factor)
        / (2.0 * k);
    Some((k, alpha))
}

/// Returns the specific attenuation in rain, in decibels per kilometer.
///
/// The rain rate is in millimeters per hour.
///
/// Returns `None` when the rain rate is negative, or when the coefficients
/// are not available.
pub fn specific_rain_attenuation(
    rain_rate: f64,
    frequency: f64,
    polarisation: f64,
    elevation: f64,
) -> Option<f64> {
    if !rain_rate.is_finite() || rain_rate < 0.0 {
        return None;
    }
    let (k, alpha) = rain_coefficients(frequency, polarisation, elevation)?;
    Some(k * rain_rate.powf(alpha))
}

/// Returns the total rain attenuation along a slant path.
///
/// The rain rate is the value exceeded for the percentage of time of
/// interest, in millimeters per hour. The station height above sea level and
/// the rain layer height are in kilometers; pass [`DEFAULT_RAIN_HEIGHT_KM`]
/// for the latter when no local value is available.
///
/// The method is the one of Recommendation P.618. The slant path length
/// follows from the rain height, the station height, and the elevation, with
/// a correction for the curvature of the Earth. Two reduction factors then
/// account for the non-uniformity of the rain: a horizontal factor that
/// depends on the specific attenuation, the frequency, and the horizontal
/// projection of the path, and a vertical factor that depends on the
/// latitude through the angle at which the path crosses the rain layer.
///
/// Returns `None` when the rain rate, frequency, or elevation is not usable,
/// when the station is above the rain layer, or when the elevation is at or
/// below the horizon.
///
/// Source: ITU-R Recommendation P.618-13, equations 1 through 10.
pub fn rain_attenuation(
    rain_rate: f64,
    frequency: f64,
    polarisation: f64,
    elevation: f64,
    station_height_km: f64,
    rain_height_km: f64,
    latitude: f64,
) -> Option<RainAttenuation> {
    if !rain_rate.is_finite() || rain_rate < 0.0 {
        return None;
    }
    if !elevation.is_finite() || elevation.sin() <= 0.0 {
        return None;
    }
    if !station_height_km.is_finite() || !rain_height_km.is_finite() {
        return None;
    }
    if !latitude.is_finite() {
        return None;
    }
    let specific = specific_rain_attenuation(rain_rate, frequency, polarisation, elevation)?;

    // Equation 2: the slant path length, with the Earth curvature term.
    let height_difference = rain_height_km - station_height_km;
    if height_difference <= 0.0 {
        // The station is at or above the rain layer, so the path is empty.
        return Some(RainAttenuation {
            specific_db_per_km: specific,
            path_length_km: 0.0,
            attenuation_db: 0.0,
            horizontal_reduction: 1.0,
            vertical_adjustment: 1.0,
        });
    }
    let sine = elevation.sin();
    let slant_length = 2.0 * height_difference
        / ((sine * sine + 2.0 * height_difference / EFFECTIVE_EARTH_RADIUS_KM).sqrt() + sine);

    // Equation 3: the horizontal projection of the slant path.
    let horizontal_length = (slant_length * elevation.cos()).abs();

    // Equation 6: the horizontal reduction factor.
    let frequency_ghz = frequency / 1.0e9;
    let horizontal_reduction = 1.0
        / (1.0 + 0.78 * (horizontal_length * specific / frequency_ghz).sqrt()
            - 0.38 * (1.0 - (-2.0 * horizontal_length).exp()));

    // Equation 7: the angle at which the path crosses the rain layer.
    let eta = (height_difference / (horizontal_length * horizontal_reduction)).atan();

    // Equation 8: the longer of the two path estimates.
    let path_length = if eta > elevation {
        horizontal_length * horizontal_reduction / elevation.cos()
    } else {
        height_difference / sine
    };

    // Equation 9: the vertical adjustment factor.
    //
    // The exponential takes the elevation in degrees, not radians, because
    // the coefficients of 31 and 0.45 were fitted with degrees. The sine
    // takes the elevation in radians, as usual.
    let latitude_magnitude = latitude.abs();
    let xi = if latitude_magnitude < 36.0 {
        36.0 - latitude_magnitude
    } else {
        0.0
    };
    let elevation_degrees = elevation.to_degrees();
    let vertical_adjustment = 1.0
        / (1.0
            + sine.sqrt()
                * (31.0
                    * (1.0 - (-(elevation_degrees / (1.0 + xi))).exp())
                    * (path_length * specific).sqrt()
                    / frequency_ghz.powi(2)
                    - 0.45));

    // Equation 10: the effective path length and the attenuation.
    let effective_length = path_length * vertical_adjustment;
    Some(RainAttenuation {
        specific_db_per_km: specific,
        path_length_km: effective_length,
        attenuation_db: specific * effective_length,
        horizontal_reduction,
        vertical_adjustment,
    })
}

/// Returns the angle in radians for a polarisation described by a tilt from
/// horizontal.
///
/// This is a convenience for the two cases that are quoted by name.
pub const fn polarisation_tilt(degrees_from_horizontal: f64) -> f64 {
    degrees_from_horizontal * core::f64::consts::PI / 180.0
}

/// Returns the total rain attenuation along a slant path, taking the rain
/// height from the Recommendation P.839 model.
///
/// This is the same calculation as [`rain_attenuation`], with the rain layer
/// height looked up from the tabulated zero degree isotherm height at the
/// given latitude and longitude rather than supplied by the caller. It is
/// the form to use unless a locally measured isotherm height is available.
///
/// The station height must be given; pass zero for a station at sea level.
///
/// Returns `None` when the coordinates are not usable, or when any of the
/// other arguments is rejected by [`rain_attenuation`].
///
/// Source: the rain height is from ITU-R Recommendation P.839-4 and the path
/// method from ITU-R Recommendation P.618-13.
pub fn rain_attenuation_at(
    rain_rate: f64,
    frequency: f64,
    polarisation: f64,
    elevation: f64,
    station_height_km: f64,
    latitude: f64,
    longitude: f64,
) -> Option<RainAttenuation> {
    let rain_height_km = crate::p839::rain_height(latitude, longitude)?;
    rain_attenuation(
        rain_rate,
        frequency,
        polarisation,
        elevation,
        station_height_km,
        rain_height_km,
        latitude,
    )
}