use crate::prelude::*;
use crate::shared::ecs::{Part, Radiator, Temperature};

const STEFAN_BOLTZMANN: f64 = 5.670374419E-8;
const T_ENV: f64 = 4.0; // units: Kelvin
const NEWTON_ITERATIONS: u32 = 3;

pub fn heat_radiation_plugin(app: &mut App) {
    app.add_systems(Update, part_radiation);
}

fn part_radiation(
    time: Res<Time>,
    parts: Query<(&mut Temperature, &Radiator, &Part)>
) {
    for (mut temperature, radiator, part) in parts {
        // time for the implicit euler method
        // e = emissivity, o = stefan-boltzmann constant, m = mass,
        // c = specific heat, A = surface area
        // k = (eoA)/(mc)
        // dT/dt = -k(T^4 - T_env^4)
        // T_env = 200.0 K
        // T_n+1 = T_n + dt*f(T_n+1)
        // T_n+1 = T_n - k*dt*((T_n+1)^4 - T_env^4)
        // g(T_n+1) = T_n+1 + k*dt*(T_n+1)^4 - T_n - k*dt*(T_env)^4
        // g'(T_n+1) = 1 + 4*k*dt*(T_n+1)^3
        // newton's method:
        // x_n+1 = x_n - g(x_n)/g'(x_n)
        // where x_n is the previous guess for T_n+1
        let initial_temp = temperature.0;
        let k = (radiator.emissivity * STEFAN_BOLTZMANN * radiator.surface_area)
            / (part.config.physics.mass * part.config.part.specific_heat);
        let dt = time.delta_secs() as f64;

        // initial guess
        let mut next_temp = initial_temp;
        for _ in 0..NEWTON_ITERATIONS {
            let g = next_temp + k*dt*next_temp.powi(4) - initial_temp - k*dt*T_ENV;
            let g_prime = 1.0 + 4.0*k*dt*(next_temp).powi(3);

            next_temp = next_temp - g/g_prime;
        }
        temperature.0 = next_temp;
    }
}