<?php
// ============================================================
// HumidityUncertaintyEngine
//
// PHP port of the GSA Humidity Laboratory calibration Python script.
// Implements all three uncertainty sub-functions exactly as in
// the original:
//   - calc_uncertainty_rh()     → expanded U for %RH
//   - calc_uncertainty_temp()   → expanded U for temperature
//   - partial_derivative_rh_wrt_temp()  → Magnus formula
//
// All method signatures mirror the Python function signatures.
// Default parameter values are identical to the originals.
// ============================================================
class HumidityUncertaintyEngine {

    // ── Magnus formula constants ─────────────────────────────
    const MAGNUS_A = 6.112;
    const MAGNUS_B = 17.62;
    const MAGNUS_C = 243.12;

    // ============================================================
    // partial_derivative_rh_wrt_temp()
    //
    // Calculates ∂RH/∂T using the Magnus formula.
    // Identical to the Python function of the same name.
    //
    // @param float $temp_c      Air temperature (°C)
    // @param float $dew_point_c Dew point temperature (°C)
    // @return float  ∂RH/∂T in %RH/°C
    // ============================================================
    public static function partialDerivativeRhWrtTemp(float $temp_c, float $dew_point_c): float {
        $A = self::MAGNUS_A;
        $B = self::MAGNUS_B;
        $C = self::MAGNUS_C;

        // Saturation vapour pressures
        $alpha = ($B * $temp_c)      / ($C + $temp_c);
        $beta  = ($B * $dew_point_c) / ($C + $dew_point_c);

        $e_s = $A * exp($alpha);
        $e_d = $A * exp($beta);

        // Derivative of e_s with respect to T
        $d_es_dT = $e_s * $B * $C / (($C + $temp_c) ** 2);

        // ∂RH/∂T = -100 × (e_d / e_s²) × (de_s/dT)
        $partial = -100.0 * ($e_d / ($e_s ** 2)) * $d_es_dT;

        return round($partial, 4);
    }

    // ============================================================
    // calc_uncertainty_rh()
    //
    // Calculates expanded uncertainty for %RH measurement.
    // Mirrors the Python calc_uncertainty_rh() exactly.
    //
    // @param float $std_dev_ref  Standard deviation of reference readings
    // @param float $std_dev_eut  Standard deviation of EUT readings
    // @param float $res_eut      Resolution of EUT (%RH)
    // @param float $k            Coverage factor (default 2)
    // @return float Expanded uncertainty U (%RH)
    // ============================================================
    public static function calcUncertaintyRh(
        float $std_dev_ref,
        float $std_dev_eut,
        float $res_eut,
        float $k = 2.0
    ): float {

        // Sensitivity coefficients (all 1 for RH, per original)
        // Using temp=10°C, dew_point=-6.77°C as in original defaults
        $sens_coef_std   = self::partialDerivativeRhWrtTemp(10, -6.77);
        $sens_coef_eut   = self::partialDerivativeRhWrtTemp(10, -6.77); // same value, mirroring original
        $sens_coef_chamb = self::partialDerivativeRhWrtTemp(10, -6.77); // same value, mirroring original

        // ── Uncertainty of Standard ──────────────────────────
        $s_unc_std_sq = self::uncertaintyOfStandardRh(
            $std_dev_ref,
            $sens_coef_std,
            $sens_coef_eut
        );

        // ── Uncertainty of Chamber ───────────────────────────
        $s_unc_chamb_sq = self::uncertaintyOfChamberRh();

        // ── Uncertainty of EUT ───────────────────────────────
        $s_unc_eut_sq = self::uncertaintyOfEutRh(
            $std_dev_eut,
            $res_eut,
            $sens_coef_std
        );

        // Combined uncertainty (RSS)
        $com_unc = sqrt($s_unc_std_sq + $s_unc_eut_sq + $s_unc_chamb_sq);

        // Expanded uncertainty
        $exp_unc = $k * $com_unc;

        return round($exp_unc, 3);
    }

    // ── uncertainty_of_standard_rh() ────────────────────────
    // Direct port of the Python inner function.
    // All default values match the Python originals exactly.
    private static function uncertaintyOfStandardRh(
        float $std_dev_std,
        float $sens_coef_std,
        float $sens_coef_eut,
        float $unc_std1          = 0.7,
        float $unc_std2          = 0.7,
        float $res_std1          = 0.1,
        float $res_std2          = 0.1,
        float $interp_cal_points = 0.5,
        float $drift_std1        = 0.1,
        float $drift_std2        = 0.1,
        float $hys_std1          = 0.31,
        float $hys_std2          = 0.21,
        float $diff_std          = 0.0,
        float $temp_dep_std      = 0.6,
        float $heat_diss_std     = 0.0,
        float $self_heat_std     = 0.25
    ): float {
        $sqrt3 = sqrt(3);

        $s_std_dev_std       = ($std_dev_std / 1)              * 1;                  // normal (k=1) × sens=1
        $s_unc_std1          = ($unc_std1 / 2)                 * 1;                  // normal (k=2) × sens=1
        $s_unc_std2          = ($unc_std2 / 2)                 * 1;                  // normal (k=2) × sens=1
        $s_res_std1          = ($res_std1 / (2 * $sqrt3))      * 1;                  // rectangular
        $s_res_std2          = ($res_std2 / (2 * $sqrt3))      * 1;                  // rectangular
        $s_interp_cal_points = ($interp_cal_points / $sqrt3)   * 1;                  // rectangular
        $s_drift_std1        = ($drift_std1 / $sqrt3)          * 1;                  // rectangular
        $s_drift_std2        = ($drift_std2 / $sqrt3)          * 1;                  // rectangular
        $s_hys_std1          = ($hys_std1 / $sqrt3)            * 1;                  // rectangular
        $s_hys_std2          = ($hys_std2 / $sqrt3)            * 1;                  // rectangular
        $s_diff_std          = ($diff_std / $sqrt3)            * 1;                  // rectangular
        $s_temp_dep_std      = ($temp_dep_std / $sqrt3)        * 1;                  // rectangular
        $s_heat_diss_std     = ($heat_diss_std / $sqrt3)       * $sens_coef_std;     // rectangular × sens_coef
        $s_self_heat_std     = ($self_heat_std / $sqrt3)       * $sens_coef_std;     // rectangular × sens_coef

        $sq = ($s_std_dev_std**2) + ($s_unc_std1**2) + ($s_unc_std2**2)
            + ($s_res_std1**2) + ($s_res_std2**2) + ($s_interp_cal_points**2)
            + ($s_drift_std1**2) + ($s_drift_std2**2) + ($s_hys_std1**2)
            + ($s_hys_std2**2) + ($s_diff_std**2) + ($s_temp_dep_std**2)
            + ($s_heat_diss_std**2) + ($s_self_heat_std**2);

        return round($sq, 4);
    }

    // ── uncertainty_of_chamber_rh() ─────────────────────────
    private static function uncertaintyOfChamberRh(
        float $temp_inst_chamb  = 0.01,
        float $spat_inhom_chamb = 0.29,
        float $rad_infl_chamb   = 0.09
    ): float {
        $sqrt3 = sqrt(3);

        $s_temp_inst  = ($temp_inst_chamb  / $sqrt3) * 1;   // rectangular × sens=1
        $s_spat_inhom = ($spat_inhom_chamb / $sqrt3) * 1;   // rectangular × sens=1
        $s_rad_infl   = ($rad_infl_chamb   / $sqrt3) * 1;   // rectangular × sens=1

        $sq = ($s_temp_inst**2) + ($s_spat_inhom**2) + ($s_rad_infl**2);

        return round($sq, 3);
    }

    // ── uncertainty_of_eut_rh() ─────────────────────────────
    private static function uncertaintyOfEutRh(
        float $std_dev_eut,
        float $res_eut,
        float $sens_coef_std,
        float $heat_diss_eut = 0.0,
        float $self_heat_eut = 0.1,
        float $hys_eut       = 0.0
    ): float {
        $sqrt3 = sqrt(3);

        $s_std_dev_eut  = ($std_dev_eut   / 1)           * 1;                // normal (k=1) × sens=1
        $s_res_eut      = ($res_eut       / (2 * $sqrt3))* 1;                // rectangular
        $s_heat_diss    = ($heat_diss_eut / $sqrt3)      * 1;                // rectangular × sens=1
        $s_self_heat    = ($self_heat_eut / $sqrt3)      * $sens_coef_std;   // rectangular × sens_coef
        $s_hys_eut      = ($hys_eut       / $sqrt3)      * 1;                // rectangular

        $sq = ($s_std_dev_eut**2) + ($s_res_eut**2) + ($s_heat_diss**2)
            + ($s_self_heat**2) + ($s_hys_eut**2);

        return round($sq, 3);
    }

    // ============================================================
    // calc_uncertainty_temp()
    //
    // Calculates expanded uncertainty for temperature measurement.
    // Mirrors the Python calc_uncertainty_temp() exactly.
    //
    // @param float $std_dev_ref  Standard deviation of reference readings
    // @param float $std_dev_eut  Standard deviation of EUT readings
    // @param float $res_eut      Resolution of EUT (°C)
    // @param float $k            Coverage factor (default 2)
    // @return float Expanded uncertainty U (K)
    // ============================================================
    public static function calcUncertaintyTemp(
        float $std_dev_ref,
        float $std_dev_eut,
        float $res_eut,
        float $k = 2.0
    ): float {

        // sens_coef = 1 for temperature (all three, per original)
        $s_unc_std_sq   = self::uncertaintyOfStandardTemp($std_dev_ref);
        $s_unc_chamb_sq = self::uncertaintyOfChamberTemp();
        $s_unc_eut_sq   = self::uncertaintyOfEutTemp($std_dev_eut, $res_eut);

        $com_unc = sqrt($s_unc_std_sq + $s_unc_chamb_sq + $s_unc_eut_sq);
        $exp_unc = $k * $com_unc;

        return round($exp_unc, 3);
    }

    // ── uncertainty_of_standard_temp() ──────────────────────
    private static function uncertaintyOfStandardTemp(
        float $std_dev_temp_std,
        float $res_std            = 0.01,
        float $interp_cal_points  = 0.2,
        float $drift_std          = 0.23,
        float $con_tech_std       = 0.001,
        float $par_therm_volt_std = 0.001,
        float $heat_diss_std      = 0.0,
        float $self_heat_std      = 0.03,
        float $hys_std            = 0.01,
        float $unc_std            = 0.2
    ): float {
        $sqrt3 = sqrt(3);

        $s_res_std            = ($res_std            / (2 * $sqrt3)) * 1;  // rectangular
        $s_interp_cal_points  = ($interp_cal_points  / $sqrt3)       * 1;  // rectangular
        $s_drift_std          = ($drift_std          / $sqrt3)       * 1;  // rectangular
        $s_con_tech_std       = ($con_tech_std       / $sqrt3)       * 1;  // rectangular
        $s_par_therm_volt_std = ($par_therm_volt_std / $sqrt3)       * 1;  // rectangular
        $s_heat_diss_std      = ($heat_diss_std      / $sqrt3)       * 1;  // rectangular
        $s_self_heat_std      = ($self_heat_std      / $sqrt3)       * 1;  // rectangular
        $s_hys_std            = ($hys_std            / $sqrt3)       * 1;  // rectangular
        $s_unc_std            = ($unc_std            / $sqrt3)       * 1;  // rectangular (note: original uses sqrt(3))
        $s_std_dev_temp_std   = ($std_dev_temp_std   / $sqrt3)       * 1;  // rectangular

        $sq = ($s_res_std**2) + ($s_interp_cal_points**2) + ($s_drift_std**2)
            + ($s_con_tech_std**2) + ($s_par_therm_volt_std**2) + ($s_heat_diss_std**2)
            + ($s_self_heat_std**2) + ($s_hys_std**2) + ($s_unc_std**2)
            + ($s_std_dev_temp_std**2);

        return round($sq, 3);
    }

    // ── uncertainty_of_chamber_temp() ───────────────────────
    private static function uncertaintyOfChamberTemp(
        float $spat_inhom_chamb = 0.3,
        float $rad_eff_chamb    = 0.14,
        float $temp_inst_chamb  = 0.01
    ): float {
        $sqrt3 = sqrt(3);

        $s_spat_inhom = ($spat_inhom_chamb / $sqrt3) * 1;  // rectangular
        $s_rad_eff    = ($rad_eff_chamb    / $sqrt3) * 1;  // rectangular
        $s_temp_inst  = ($temp_inst_chamb  / $sqrt3) * 1;  // rectangular

        $sq = ($s_spat_inhom**2) + ($s_rad_eff**2) + ($s_temp_inst**2);

        return round($sq, 3);
    }

    // ── uncertainty_of_eut_temp() ────────────────────────────
    private static function uncertaintyOfEutTemp(
        float $std_dev_temp_eut,
        float $res_temp_eut,
        float $heat_diss_eut = 0.0,
        float $self_heat_eut = 0.0,
        float $hys_eut       = 0.0
    ): float {
        $sqrt3 = sqrt(3);

        $s_std_dev_eut = ($std_dev_temp_eut / 1)           * 1;   // normal (k=1) × sens=1
        $s_res_eut     = ($res_temp_eut     / (2 * $sqrt3))* 1;   // rectangular
        $s_heat_diss   = ($heat_diss_eut    / $sqrt3)      * 1;   // rectangular
        $s_self_heat   = ($self_heat_eut    / $sqrt3)      * 1;   // rectangular (sens_coef_std = 1 for temp)
        $s_hys_eut     = ($hys_eut          / $sqrt3)      * 1;   // rectangular

        $sq = ($s_std_dev_eut**2) + ($s_res_eut**2) + ($s_heat_diss**2)
            + ($s_self_heat**2) + ($s_hys_eut**2);

        return round($sq, 4);
    }

    // ============================================================
    // processSetpoint()
    //
    // Processes one calibration setpoint exactly as the Python
    // main() loop does for each setpoint. Returns a result array
    // matching the Python `result` dict structure.
    //
    // @param array  $ref_temps        Reference temperature readings
    // @param float  $correction_temp  Correction for reference temperature
    // @param array  $ref_rh1          Ref RH Sensor 1 readings
    // @param float  $correction_rh1   Correction for Ref RH Sensor 1
    // @param array  $ref_rh2          Ref RH Sensor 2 readings
    // @param float  $correction_rh2   Correction for Ref RH Sensor 2
    // @param array  $eut_temps        EUT temperature readings
    // @param array  $eut_rhs          EUT RH readings
    // @param float  $res_rh_eut       EUT RH resolution (default 0.1)
    // @param float  $res_temp_eut     EUT Temp resolution (default 0.1)
    // @param int    $sp_number        Setpoint number (for labelling)
    // @return array Result dict matching Python output
    // ============================================================
    public static function processSetpoint(
        array $ref_temps,
        float $correction_temp,
        array $ref_rh1,
        float $correction_rh1,
        array $ref_rh2,
        float $correction_rh2,
        array $eut_temps,
        array $eut_rhs,
        float $res_rh_eut   = 0.1,
        float $res_temp_eut = 0.1,
        int   $sp_number    = 1
    ): array {

        // ── Reference calculations ───────────────────────────
        $avg_ref_temp       = array_sum($ref_temps) / count($ref_temps);
        $corrected_ref_temp = $avg_ref_temp + $correction_temp;

        $avg_ref_rh1        = array_sum($ref_rh1) / count($ref_rh1);
        $corrected_ref_rh1  = $avg_ref_rh1 + $correction_rh1;

        $avg_ref_rh2        = array_sum($ref_rh2) / count($ref_rh2);
        $corrected_ref_rh2  = $avg_ref_rh2 + $correction_rh2;

        $corrected_ref_rh   = ($corrected_ref_rh1 + $corrected_ref_rh2) / 2.0;

        // ── EUT statistics ───────────────────────────────────
        $n          = count($eut_temps);
        $avg_temp   = array_sum($eut_temps) / $n;
        $avg_rh     = array_sum($eut_rhs)   / $n;

        // Standard deviation (sample, ddof=1 — matches Python statistics.stdev)
        $std_temp = $n > 1 ? self::stdev($eut_temps) : 0.0;
        $std_rh   = $n > 1 ? self::stdev($eut_rhs)   : 0.0;

        // ── Deviations ───────────────────────────────────────
        $temp_dev = $avg_temp - $corrected_ref_temp;
        $rh_dev   = $avg_rh   - $corrected_ref_rh;

        // ── Uncertainties ────────────────────────────────────
        $temp_unc = self::calcUncertaintyTemp($std_temp, $std_temp, $res_temp_eut);
        $rh_unc   = self::calcUncertaintyRh($std_rh,  $std_rh,  $res_rh_eut);

        // ── Determine if temperature results should be included
        // (matches Python: if 45 < corrected_ref_rh < 55)
        $include_temp = ($corrected_ref_rh > 45 && $corrected_ref_rh < 55);

        return [
            'setpoint'            => $sp_number,
            'avg_ref_temp'        => round($avg_ref_temp, 4),
            'corrected_ref_temp'  => round($corrected_ref_temp, 4),
            'corrected_ref_rh'    => round($corrected_ref_rh, 4),
            'avg_eut_temp'        => round($avg_temp, 4),
            'std_eut_temp'        => round($std_temp, 6),
            'avg_eut_rh'          => round($avg_rh, 4),
            'std_eut_rh'          => round($std_rh, 6),
            'temp_deviation'      => round($temp_dev, 4),
            'rh_deviation'        => round($rh_dev, 4),
            'temp_uncertainty'    => $temp_unc,
            'rh_uncertainty'      => $rh_unc,
            'include_temp'        => $include_temp,
        ];
    }

    // ── Sample standard deviation (matches Python statistics.stdev) ──
    public static function stdev(array $values): float {
        $n    = count($values);
        if ($n < 2) return 0.0;
        $mean = array_sum($values) / $n;
        $sq   = array_reduce($values, fn($carry, $v) => $carry + ($v - $mean) ** 2, 0.0);
        return sqrt($sq / ($n - 1));
    }
}
