//! Solution of the radiative transfer equation by Feautrier method
//! for two continuum points.
//!
//! Translated from SYNSPEC54.FOR subroutine RTECD (line 12952).
//!
//! Used when one employs RTEDFE (DFE method) for the inner frequency points.
//! The solver uses variable Eddington factors and a 3×3 matrix inversion.

use crate::synspec::math::inibla::{BN, HK};

// ============================================================================
// Constants
// ============================================================================

/// Gauss quadrature points (μ values) for 3-point scheme
const AMU3: [f64; 3] = [0.887_298_334_620_742, 0.5, 0.112_701_665_379_258];

/// Gauss quadrature weights
const WTMU3: [f64; 3] = [0.277_777_777_777_778, 0.444_444_444_444_444, 0.277_777_777_777_778];

/// Reference optical depth (τ = 2/3)
const TAUREF: f64 = 0.666_666_666_666_7;

/// Speed of light (Å/s) for wavelength conversion
const CL_ANGSTROM: f64 = 2.997_925e18;

// ============================================================================
// RTECD parameters
// ============================================================================

/// Input parameters for the RTECD solver.
pub struct RtecdParams<'a> {
    /// Number of depth points
    pub nd: usize,
    /// Mass depth array (g/cm²) [nd]
    pub dm: &'a [f64],
    /// Density array (g/cm³) [nd]
    pub dens: &'a [f64],
    /// Temperature array (K) [nd]
    pub temp: &'a [f64],
    /// Frequency array (Hz) [2] (two continuum points)
    pub freq: &'a [f64],
    /// Absorption coefficient [2 × nd] (continuum)
    pub ch: &'a [Vec<f64>],
    /// Emission coefficient [2 × nd] (continuum)
    pub et: &'a [Vec<f64>],
    /// Scattering coefficient [2 × nd] (continuum)
    pub sc: &'a [Vec<f64>],
    /// Number of μ points for specific intensity (NMU0)
    pub nmu0: usize,
    /// Angles for specific intensity output [nmu0]
    pub angl: &'a [f64],
    /// Weangles for specific intensity output [nmu0]
    pub wangl: &'a [f64],
    /// Iflux flag: 0 = no specific intensity, >=1 = compute
    pub iflux: i32,
    /// Iprin flag for diagnostic output
    pub iprin: i32,
}

/// Output from the RTECD solver.
pub struct RtecdResult {
    /// Flux at the two continuum frequencies [2]
    pub flux: [f64; 2],
    /// Scattering source function for continuum 1 [nd]
    pub scc1: Vec<f64>,
    /// Scattering source function for continuum 2 [nd]
    pub scc2: Vec<f64>,
    /// Specific intensities at the surface [nmu0] (if iflux >= 1)
    pub rint_surface: Vec<f64>,
    /// Emergent flux from specific intensities (if iflux >= 1)
    pub flx: f64,
}

// ============================================================================
// 3×3 matrix inversion helper
// ============================================================================

/// Inlined 3×3 matrix inversion (MINV3).
///
/// Replaces the Fortran code that inlines MATINV for a specific 3×3 case.
/// Modifies bb in-place.
fn minv3(bb: &mut [[f64; 3]; 3]) {
    // Forward elimination
    bb[1][0] /= bb[0][0];
    bb[1][1] -= bb[1][0] * bb[0][1];
    bb[1][2] -= bb[1][0] * bb[0][2];
    bb[2][0] /= bb[0][0];
    bb[2][1] = (bb[2][1] - bb[2][0] * bb[0][1]) / bb[1][1];
    bb[2][2] -= bb[2][0] * bb[0][2] - bb[2][1] * bb[1][2];

    // Back substitution
    bb[2][1] = -bb[2][1];
    bb[2][0] = -bb[2][0] - bb[2][1] * bb[1][0];
    bb[1][0] = -bb[1][0];

    bb[2][2] = 1.0 / bb[2][2];
    bb[1][2] = -bb[1][2] * bb[2][2] / bb[1][1];
    bb[1][1] = 1.0 / bb[1][1];
    bb[0][2] = -(bb[0][1] * bb[1][2] + bb[0][2] * bb[2][2]) / bb[0][0];
    bb[0][1] = -bb[0][1] * bb[1][1] / bb[0][0];
    bb[0][0] = 1.0 / bb[0][0];

    // Final transformation
    bb[0][0] = bb[0][0] + bb[0][1] * bb[1][0] + bb[0][2] * bb[2][0];
    bb[0][1] += bb[0][2] * bb[2][1];
    bb[1][0] = bb[1][1] * bb[1][0] + bb[1][2] * bb[2][0];
    bb[1][1] += bb[1][2] * bb[2][1];
    bb[2][0] *= bb[2][2];
    bb[2][1] *= bb[2][2];
}

// ============================================================================
// RTECD main function
// ============================================================================

/// Solve the radiative transfer equation by Feautrier method for two
/// continuum frequency points.
///
/// # Fortran original
///
/// ```fortran
/// SUBROUTINE RTECD
///   Solution of the radiative transfer equation by Feautrier method
///   for two continuum points
///   used when one employs RTEDFE, ie. the DFE method for the
///   transfer equation for the inner frequency points
/// ```
pub fn rtecd(params: &RtecdParams) -> RtecdResult {
    let RtecdParams {
        nd,
        dm,
        dens,
        temp,
        freq,
        ch,
        et,
        sc,
        nmu0,
        angl,
        wangl,
        iflux,
        iprin,
    } = *params;

    let nmu: usize = 3;
    let nd1 = nd.saturating_sub(1);

    // Guard: return zeros if no frequency data
    if freq.is_empty() || nd < 2 {
        return RtecdResult {
            flux: [0.0; 2],
            scc1: vec![0.0; nd],
            scc2: vec![0.0; nd],
            rint_surface: vec![0.0; nmu0],
            flx: 0.0,
        };
    }

    // Output arrays
    let mut flux = [0.0f64; 2];
    let mut scc1 = vec![0.0f64; nd];
    let mut scc2 = vec![0.0f64; nd];
    let mut rint_surface = vec![0.0f64; nmu0];
    let mut flx_out = 0.0f64;

    // Working arrays
    let mut d = vec![[[0.0f64; 3]; 3]; nd]; // D[i][j][id]
    let mut anu = vec![[0.0f64; 3]; nd];     // ANU[i][id]
    let mut aanu = vec![0.0f64; nd];
    let mut ddd = vec![0.0f64; nd];
    let mut tau = vec![0.0f64; nd];
    let mut st0 = vec![0.0f64; nd];
    let mut ss0 = vec![0.0f64; nd];
    let mut fkk = vec![0.0f64; nd];
    let mut rdd = vec![0.0f64; nd];
    let mut rint = vec![vec![0.0f64; nmu0]; nd];
    let mut dt = vec![0.0f64; nd];

    // ========================================================================
    // Loop over two continuum frequencies
    // ========================================================================
    for ij in 0..2 {
        let taumin = ch[ij][0] / dens[0] * dm[0] * 0.5;
        tau[0] = taumin;

        for i in 0..nd1 {
            dt[i] = (dm[i + 1] - dm[i])
                * (ch[ij][i + 1] / dens[i + 1] + ch[ij][i] / dens[i])
                * 0.5;
            st0[i] = et[ij][i] / ch[ij][i];
            ss0[i] = -sc[ij][i] / ch[ij][i];
            tau[i + 1] = tau[i] + dt[i];
        }
        st0[nd - 1] = et[ij][nd - 1] / ch[ij][nd - 1];
        ss0[nd - 1] = -sc[ij][nd - 1] / ch[ij][nd - 1];

        let fr = freq[ij];
        let bnu = BN * (fr * 1.0e-15).powi(3);
        let pland = bnu / ((HK * fr / temp[nd - 1]).exp() - 1.0);
        let dplan = bnu / ((HK * fr / temp[nd - 2]).exp() - 1.0);
        let dplan = (pland - dplan) / dt[nd1 - 1];

        // Find reference depth (τ = 2/3)
        let mut iref: usize = 0;
        for i in 0..nd1 {
            if tau[i] <= TAUREF && tau[i + 1] > TAUREF {
                iref = i;
            }
        }

        // ================================================================
        // FIRST PART - Variable Eddington factors
        // ================================================================

        let alb1 = 0.0f64;

        // Upper boundary condition
        let mut dtp1 = dt[0];
        let mut q0 = 0.0f64;
        let p0: f64;

        // Allowance for non-zero optical depth at first depth point
        let tamm = taumin / AMU3[0];
        if tamm > 0.01 {
            p0 = 1.0 - (-tamm).exp();
        } else {
            p0 = tamm * (1.0 - 0.5 * tamm * (1.0 - tamm / 3.0 * (1.0 - 0.25 * tamm)));
        }
        let _ex = 1.0 - p0;
        q0 += p0 * AMU3[0] * WTMU3[0];

        let div = dtp1 / AMU3[0] / 3.0;
        let vl0 = div * (st0[0] + 0.5 * st0[1]) + st0[0] * p0;

        // Build BB matrix for upper boundary
        let mut bb = [[0.0f64; 3]; 3];
        let mut cc = [[0.0f64; 3]; 3];
        let mut vl = [0.0f64; 3];

        for i in 0..nmu {
            vl[i] = vl0;
            for j in 0..nmu {
                bb[i][j] = ss0[0] * WTMU3[j] * (div + p0) - alb1 * WTMU3[j];
                cc[i][j] = -0.5 * div * ss0[1] * WTMU3[j];
            }
            bb[i][i] += AMU3[i] / dtp1 + 1.0 + div;
            cc[i][i] += AMU3[i] / dtp1 - 0.5 * div;
            anu[i][0] = 0.0;
        }

        // 3×3 matrix inversion
        minv3(&mut bb);

        // Compute D and ANU at first depth
        for i in 0..nmu {
            for j in 0..nmu {
                let mut s = 0.0;
                for k in 0..nmu {
                    s += bb[i][k] * cc[k][j];
                }
                d[i][j][0] = s;
                anu[i][0] += bb[i][j] * vl[j];
            }
        }

        // Normal depth points
        for id in 1..nd1 {
            let dtm1 = dtp1;
            dtp1 = dt[id];
            let dt0 = 0.5 * (dtm1 + dtp1);
            let al = 1.0 / dtm1 / dt0;
            let ga = 1.0 / dtp1 / dt0;
            let be = al + ga;
            let a = (1.0 - 0.5 * al * dtp1 * dtp1) / 6.0;
            let c = (1.0 - 0.5 * ga * dtm1 * dtm1) / 6.0;
            let b = 1.0 - a - c;
            let vl0 = a * st0[id - 1] + b * st0[id] + c * st0[id + 1];

            let mut aa = [[0.0f64; 3]; 3];
            let mut bb = [[0.0f64; 3]; 3];
            let mut cc = [[0.0f64; 3]; 3];
            let mut vl = [0.0f64; 3];

            for i in 0..nmu {
                vl[i] = vl0;
                for j in 0..nmu {
                    aa[i][j] = -a * ss0[id - 1] * WTMU3[j];
                    cc[i][j] = -c * ss0[id + 1] * WTMU3[j];
                    bb[i][j] = b * ss0[id] * WTMU3[j];
                }
            }
            for i in 0..nmu {
                let div = AMU3[i] * AMU3[i];
                aa[i][i] += div * al - a;
                cc[i][i] += div * ga - c;
                bb[i][i] += div * be + b;
            }

            // Eliminate previous depth
            for i in 0..nmu {
                let mut s1 = 0.0;
                for j in 0..nmu {
                    let mut s = 0.0;
                    s1 += aa[i][j] * anu[j][id - 1];
                    for k in 0..nmu {
                        s += aa[i][k] * d[k][j][id - 1];
                    }
                    bb[i][j] -= s;
                }
                vl[i] += s1;
            }

            // 3×3 matrix inversion
            minv3(&mut bb);

            // Compute D and ANU at this depth
            for i in 0..nmu {
                anu[i][id] = 0.0;
                for j in 0..nmu {
                    let mut s = 0.0;
                    for k in 0..nmu {
                        s += bb[i][k] * cc[k][j];
                    }
                    d[i][j][id] = s;
                    anu[i][id] += bb[i][j] * vl[j];
                }
            }
        }

        // Lower boundary condition
        let id = nd - 1;
        let mut aa = [[0.0f64; 3]; 3];
        let mut bb = [[0.0f64; 3]; 3];
        let mut vl = [0.0f64; 3];

        for i in 0..nmu {
            aa[i][i] = AMU3[i] / dtp1;
            vl[i] = pland + AMU3[i] * dplan + aa[i][i] * anu[i][id - 1];
            for j in 0..nmu {
                bb[i][j] = -aa[i][i] * d[i][j][id - 1];
            }
            bb[i][i] += aa[i][i] + 1.0;
        }

        // 3×3 matrix inversion
        minv3(&mut bb);

        for i in 0..nmu {
            anu[i][id] = 0.0;
            for j in 0..nmu {
                d[i][j][id] = 0.0;
                anu[i][id] += bb[i][j] * vl[j];
            }
        }

        // Backsolution
        for iid in 0..nd1 {
            let id = nd1 - 1 - iid;
            for i in 0..nmu {
                for j in 0..nmu {
                    anu[i][id] += d[i][j][id] * anu[j][id + 1];
                }
            }
        }

        // Compute Eddington factors
        let mut aj = 0.0;
        let mut ak = 0.0;
        for i in 0..nmu {
            let div = WTMU3[i] * anu[i][0];
            aj += div;
            ak += div * AMU3[i] * AMU3[i];
        }
        fkk[0] = ak / aj;

        for id in 1..nd1 {
            aj = 0.0;
            ak = 0.0;
            for i in 0..nmu {
                let div = WTMU3[i] * anu[i][id];
                aj += div;
                ak += div * AMU3[i] * AMU3[i];
            }
            fkk[id] = ak / aj;
        }

        // Surface Eddington factor
        let mut ah = 0.0;
        for i in 0..nmu {
            ah += WTMU3[i] * AMU3[i] * anu[i][0];
        }
        let fh = ah / aj - 0.5 * alb1;

        fkk[nd - 1] = 1.0 / 3.0;

        // ================================================================
        // SECOND PART - Determination of mean intensities
        // ================================================================

        dtp1 = dt[0];
        let div = dtp1 / 3.0;
        let mut bbb = fkk[0] / dtp1 + fh + div + ss0[0] * (div + q0);
        let mut ccc = fkk[1] / dtp1 - 0.5 * div * (1.0 + ss0[1]);
        let vll = div * (st0[0] + 0.5 * st0[1]) + st0[0] * q0;
        aanu[0] = vll / bbb;
        ddd[0] = ccc / bbb;

        for id in 1..nd1 {
            let dtm1 = dtp1;
            dtp1 = dt[id];
            let dt0 = 0.5 * (dtp1 + dtm1);
            let al = 1.0 / dtm1 / dt0;
            let ga = 1.0 / dtp1 / dt0;
            let a = (1.0 - 0.5 * dtp1 * dtp1 * al) / 6.0;
            let c = (1.0 - 0.5 * dtm1 * dtm1 * ga) / 6.0;

            let aaa = al * fkk[id - 1] - a * (1.0 + ss0[id - 1]);
            ccc = ga * fkk[id + 1] - c * (1.0 + ss0[id + 1]);
            bbb = (al + ga) * fkk[id] + (1.0 - a - c) * (1.0 + ss0[id]);
            let vll = a * st0[id - 1] + c * st0[id + 1] + (1.0 - a - c) * st0[id];
            bbb -= aaa * ddd[id - 1];
            ddd[id] = ccc / bbb;
            aanu[id] = (vll + aaa * aanu[id - 1]) / bbb;
        }

        bbb = fkk[nd - 1] / dtp1 + 0.5;
        let aaa = fkk[nd1 - 1] / dtp1;
        bbb -= aaa * ddd[nd1 - 1];
        let vll = 0.5 * pland + dplan / 3.0;
        rdd[nd - 1] = (vll + aaa * aanu[nd1 - 1]) / bbb;

        for iid in 0..nd1 {
            let id = nd1 - 1 - iid;
            rdd[id] = aanu[id] + ddd[id] * rdd[id + 1];
        }

        flux[ij] = fh * rdd[0];

        // Store scattering source functions
        if ij == 0 {
            for id in 0..nd {
                scc1[id] = -rdd[id] * ss0[id] * ch[0][id];
            }
        } else {
            for id in 0..nd {
                scc2[id] = -rdd[id] * ss0[id] * ch[1][id];
            }
        }

        // Diagnostic output
        if iprin >= 3 {
            let t0 = (tau[iref + 1] / tau[iref]).ln();
            let x0 = (tau[iref + 1] / TAUREF).ln() / t0;
            let x1 = (TAUREF / tau[iref]).ln() / t0;
            let dmref = (dm[iref].ln() * x0 + dm[iref + 1].ln() * x1).exp();
            let tref = (temp[iref].ln() * x0 + temp[iref + 1].ln() * x1).exp();
            let stref = (st0[iref].ln() * x0 + st0[iref + 1].ln() * x1).exp();
            let scref = ((-ss0[iref]).ln() * x0 + (-ss0[iref + 1]).ln() * x1).exp();
            let ssref = ((-ss0[iref] * rdd[iref]).ln() * x0
                + (-ss0[iref + 1] * rdd[iref + 1]).ln() * x1)
                .exp();
            let sref = stref + ssref;
            let alm = CL_ANGSTROM / freq[ij];
            eprintln!(
                "RTECD: IJ={} ALM={:.1} IREF={} DMREF={:.3e} TREF={:.0} SCREF={:.3e} STREF={:.3e} SSREF={:.3e} SREF={:.3e}",
                ij + 1, alm, iref + 1, dmref, tref, scref, stref, ssref, sref
            );
        }

        // ================================================================
        // THIRD PART - Specific intensities
        // ================================================================

        if iflux == 0 {
            continue;
        }

        for imu in 0..nmu0 {
            let anx = angl[imu];
            dtp1 = dt[0];
            let div = dtp1 / 3.0 / anx;

            let tamm = taumin / anx;
            let p0: f64 = if tamm < 0.01 {
                tamm * (1.0 - 0.5 * tamm * (1.0 - tamm / 3.0 * (1.0 - 0.25 * tamm)))
            } else {
                1.0 - (-tamm).exp()
            };

            bbb = anx / dtp1 + 1.0 + div;
            ccc = anx / dtp1 - 0.5 * div;
            let vll = (div + p0) * (st0[0] - ss0[0] * rdd[0])
                + 0.5 * div * (st0[1] - ss0[1] * rdd[1]);
            aanu[0] = vll / bbb;
            ddd[0] = ccc / bbb;

            let div = anx * anx;
            for id in 1..nd1 {
                let dtm1 = dt[id - 1];
                dtp1 = dt[id];
                let dt0 = 0.5 * (dtp1 + dtm1);
                let al = 1.0 / dtm1 / dt0;
                let ga = 1.0 / dtp1 / dt0;
                let a = (1.0 - 0.5 * dtp1 * dtp1 * al) / 6.0;
                let c = (1.0 - 0.5 * dtm1 * dtm1 * ga) / 6.0;

                let aaa = div * al - a;
                ccc = div * ga - c;
                bbb = div * (al + ga) + 1.0 - a - c;
                let vll = a * (st0[id - 1] - ss0[id - 1] * rdd[id - 1])
                    + c * (st0[id + 1] - ss0[id + 1] * rdd[id + 1])
                    + (1.0 - a - c) * (st0[id] - ss0[id] * rdd[id]);
                bbb -= aaa * ddd[id - 1];
                ddd[id] = ccc / bbb;
                aanu[id] = (vll + aaa * aanu[id - 1]) / bbb;
            }

            // Lower boundary condition
            let aaa = anx / dtp1;
            bbb = aaa + 1.0;
            let vll = pland + anx * dplan;

            rint[nd - 1][imu] = (vll + aaa * aanu[nd1 - 1]) / (bbb - aaa * ddd[nd1 - 1]);
            for iid in 0..nd1 {
                let id = nd1 - 1 - iid;
                rint[id][imu] = aanu[id] + ddd[id] * rint[id + 1][imu];
            }
        }

        // Compute emergent flux from specific intensities
        let mut flx = 0.0;
        for imu in 0..nmu0 {
            rint[0][imu] /= 0.5;
            flx += angl[imu] * wangl[imu] * rint[0][imu];
        }
        flx *= 0.5;

        if iflux >= 1 {
            flx_out = flx;
            rint_surface.copy_from_slice(&rint[0][..nmu0]);
        }
    }

    RtecdResult {
        flux,
        scc1,
        scc2,
        rint_surface,
        flx: flx_out,
    }
}

#[cfg(test)]
mod tests {
    use super::*;

    #[test]
    fn test_minv3_identity() {
        let mut m = [[1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0]];
        minv3(&mut m);
        // Should be identity
        assert!((m[0][0] - 1.0).abs() < 1e-10);
        assert!((m[1][1] - 1.0).abs() < 1e-10);
        assert!((m[2][2] - 1.0).abs() < 1e-10);
        assert!(m[0][1].abs() < 1e-10);
        assert!(m[0][2].abs() < 1e-10);
        assert!(m[1][0].abs() < 1e-10);
    }

    #[test]
    fn test_minv3_diagonal() {
        let mut m = [[2.0, 0.0, 0.0], [0.0, 4.0, 0.0], [0.0, 0.0, 5.0]];
        minv3(&mut m);
        assert!((m[0][0] - 0.5).abs() < 1e-10);
        assert!((m[1][1] - 0.25).abs() < 1e-10);
        assert!((m[2][2] - 0.2).abs() < 1e-10);
    }
}