const ASM3 = require('./asm3_class')
const { create, all } = require('mathjs')

const config = {
  matrix: 'Array' // Choose 'Matrix' (default) or 'Array'
}

const math = create(all, config)

class Reactor_CSTR {
    
    constructor(volume, n_inlets, kla, initial_state) {
        this.state = initial_state;
        this.asm = new ASM3();

        this.Vl = volume; // fluid volume reactor [m3]
        this.Fs = Array(n_inlets).fill(0.0); // fluid debits per inlet [m3 d-1]
        this.Cs_in = Array.from(Array(n_inlets), () => new Array(13).fill(0.0)); // composition influents
        this.OTR = 0.0; // oxygen transfer rate [g O2 d-1]
        
        this.kla = kla; // if NaN, use external OTR [d-1]

        this.currentTime = Date.now(); // milliseconds since epoch [ms]
        this.timeStep = 1/(24*60*15); // time step [d]
        this.speedUpFactor = 1;
    }

    set setInfluent(input) { // setter for C_in (WIP)
        let index_in = input.payload.inlet;
        this.Fs[index_in] = input.payload.F;
        this.Cs_in[index_in] = input.payload.C;
    }

    set setOTR(input) { // setter for OTR (WIP) [g O2 d-1]
        this.OTR = input.payload;
    }

    get getEffluent() { // getter for Effluent, defaults to inlet 0
        return {topic: "Fluent", payload: {inlet: 0, F: math.sum(this.Fs), C:this.state}, timestamp: this.currentTime};
    }

    calcOTR(S_O, T=20.0) { // caculate the OTR using basic correlation, default to temperature: 20 C
        let S_O_sat = 14.652 - 4.1022e-1*T + 7.9910e-3*T*T + 7.7774e-5*T*T*T;
        return this.kla * (S_O_sat - S_O);
    }

    // expect update with timestamp
    updateState(newTime) {

        const day2ms = 1000 * 60 * 60 * 24; 

        let n_iter = Math.floor(this.speedUpFactor*(newTime - this.currentTime) / (this.timeStep * day2ms));
        if (n_iter) {
            let n = 0;
            while (n < n_iter) {
                this.tick_fe(this.timeStep);
                n += 1;
            }
            this.currentTime += n_iter * this.timeStep * day2ms / this.speedUpFactor;
        }
    }

    tick_fe(time_step) { // tick reactor state using forward Euler method
        const r = this.asm.compute_dC(this.state);
        const dC_in = math.multiply(math.divide([this.Fs], this.Vl), this.Cs_in)[0];
        const dC_out = math.multiply(-1*math.sum(this.Fs)/this.Vl, this.state);
        const t_O = Array(13).fill(0.0);
        t_O[0] = isNaN(this.kla) ? this.OTR : this.calcOTR(this.state[0]); // calculate OTR if kla is not NaN, otherwise use externaly calculated OTR

        const dC_total = math.multiply(math.add(dC_in, dC_out, r, t_O), time_step);

        // clip value element-wise to each subarray to avoid negative concentrations
        this.state = math.add(this.state, dC_total).map(val => val < 0 ? 0 : val);
        return this.state;
    }
}

class Reactor_PFR {
    
    constructor(volume, length, resolution_L, n_inlets, kla, initial_state) {
        this.asm = new ASM3();

        this.Vl = volume; // fluid volume reactor [m3]
        this.length = length; // reactor length [m]
        this.n_x = resolution_L; // number of slices
        this.d_x = length / resolution_L;

        this.A = volume / length; // crosssectional area [m2]

        this.state = Array.from(Array(this.n_x), () => initial_state.slice())
        
        // console.log("Initial State: ")
        // console.log(this.state)

        this.Fs = Array(n_inlets).fill(0.0); // fluid debits per inlet [m3 d-1]
        this.Cs_in = Array.from(Array(n_inlets), () => new Array(13).fill(0.0)); // composition influents
        this.OTR = 0.0; // oxygen transfer rate [g O2 d-1]
        this.D = 0.0; // axial dispersion [m2 d-1]
        
        this.kla = kla; // if NaN, use external OTR [d-1]

        this.currentTime = Date.now(); // milliseconds since epoch [ms]
        this.timeStep = 1/(24*60*15); // time step [d]
        this.speedUpFactor = 60;

        this.D_op = this.makeDoperator(true);
        this.D2_op = this.makeD2operator();
    }

    set setInfluent(input) { // setter for C_in (WIP)
        let index_in = input.payload.inlet;
        this.Fs[index_in] = input.payload.F;
        this.Cs_in[index_in] = input.payload.C;
        console.log("Pe total " + this.length*math.sum(this.Fs)/(this.D*this.A));
        console.log("Pe local " + this.d_x*math.sum(this.Fs)/(this.D*this.A));
        console.log("Co ad " + math.sum(this.Fs)*this.timeStep/(this.A*this.d_x));
        console.log("Co D " + this.D*this.timeStep/(this.d_x*this.d_x));
    }

    set setOTR(input) { // setter for OTR (WIP) [g O2 d-1]
        this.OTR = input.payload;
    }

    set setDispersion(input) { // setter for Axial dispersion [m2 d-1]
        this.D = input.payload;
    }

    get getEffluent() { // getter for Effluent, defaults to inlet 0
        return {topic: "Fluent", payload: {inlet: 0, F: math.sum(this.Fs), C:this.state.at(-1)}, timestamp: this.currentTime};
    }

    calcOTR(S_O, T=20.0) { // caculate the OTR using basic correlation, default to temperature: 20 C
        let S_O_sat = 14.652 - 4.1022e-1*T + 7.9910e-3*T*T + 7.7774e-5*T*T*T;
        return this.kla * (S_O_sat - S_O);
    }

    // expect update with timestamp
    updateState(newTime) {
        const day2ms = 1000 * 60 * 60 * 24; 

        let n_iter = Math.floor(this.speedUpFactor*(newTime - this.currentTime) / (this.timeStep * day2ms));
        if (n_iter) {
            let n = 0;
            while (n < n_iter) {
                this.tick_fe(this.timeStep);
                n += 1;
            }
            this.currentTime += n_iter * this.timeStep * day2ms / this.speedUpFactor;
        }
    }

    tick_fe(time_step) { // tick reactor state using forward Euler method
        const dispersion = math.multiply(this.D / (this.d_x*this.d_x), this.D2_op, this.state);
        const advection = math.multiply(-1*math.sum(this.Fs)/(this.A*this.d_x), this.D_op, this.state);
        const reaction = this.state.map((state_slice) => this.asm.compute_dC(state_slice));
        const transfer = Array.from(Array(this.n_x), () => new Array(13).fill(0.0));

        if (dispersion.some(row => row.some(Number.isNaN))) {
            throw new Error("NaN detected in dispersion!");
        }
        if (advection.some(row => row.some(Number.isNaN))) {
            throw new Error("NaN detected in advection!");
        }
        if (reaction.some(row => row.some(Number.isNaN))) {
            throw new Error("NaN detected in reaction!");
        }
        
        if (isNaN(this.kla)) { // calculate OTR if kla is not NaN, otherwise use externally calculated OTR
            transfer.forEach((x) => { x[0] = this.OTR; });
        } else {
            transfer.forEach((x, i) => { x[0] = this.calcOTR(this.state[i][0]); });
        }

        const dC_total = math.multiply(math.add(dispersion, advection, reaction, transfer), time_step);
        const new_state = math.add(this.state, dC_total);
        if (new_state.some(row => row.some(Number.isNaN))) {
            throw new Error("NaN detected in new_state after dC_total update!");
        }

        // apply boundary conditions
        if (math.sum(this.Fs) > 0) { // Danckwerts BC
            const BC_C_in = math.multiply(1/math.sum(this.Fs), [this.Fs], this.Cs_in)[0];
            const BC_gradient = Array(this.n_x).fill(0.0);
            BC_gradient[0] = -1;
            BC_gradient[1] = 1;
            const BC_dispersion = math.multiply(this.D * this.A / (math.sum(this.Fs)*this.d_x), [BC_gradient], new_state)[0];
            console.log(math.add(BC_C_in, BC_dispersion));
            new_state[0] = math.add(BC_C_in, BC_dispersion).map(val => val < 0 ? 0 : val);
            
        } else { // Neumann BC (no flux)
            new_state[0] = new_state[1];
        }
        // Neumann BC (no flux)
        new_state[this.n_x-1] = new_state[this.n_x-2]

        if (new_state.some(row => row.some(Number.isNaN))) {
            throw new Error("NaN detected in new_state after enforcing boundary conditions!");
        }

        this.state = new_state.map(row => row.map(val => val < 0 ? 0 : val)); // apply the new state
        return new_state;
    }

    makeDoperator(central=false, higher_order=false) { // create gradient operator
        if (higher_order) {
            if (central) {
                const I = math.resize(math.diag(Array(this.n_x).fill(1/12), -2), [this.n_x, this.n_x]);
                const A = math.resize(math.diag(Array(this.n_x).fill(-2/3), -1), [this.n_x, this.n_x]);
                const B = math.resize(math.diag(Array(this.n_x).fill(2/3), 1), [this.n_x, this.n_x]);
                const C = math.resize(math.diag(Array(this.n_x).fill(-1/12), 2), [this.n_x, this.n_x]);
                const D = math.add(I, A);
                const NearBoundary = Array(this.n_x).fill(0.0);
                NearBoundary[1] = -25/12;
                NearBoundary[2] = 4;
                NearBoundary[3] = -3;
                NearBoundary[4] = 4/3;
                NearBoundary[5] = -1/4;
                D[1] = NearBoundary;
                NearBoundary.reverse();
                D[this.n_x-2] = math.multiply(-1, NearBoundary)
            } else {
                throw new Error("Upwind higher order method not implemented! Use central scheme instead.");
            }
        } else {
            const I = math.resize(math.diag(Array(this.n_x).fill(1/(1+central)), central), [this.n_x, this.n_x]);
            const A = math.resize(math.diag(Array(this.n_x).fill(-1/(1+central)), -1), [this.n_x, this.n_x]);
            const D = math.add(I, A);
        }

        D[0] = Array(this.n_x).fill(0); // set by BCs elsewhere
        D[this.n_x-1] = Array(this.n_x).fill(0);
        return D;
    }

    makeD2operator() { // create the central second derivative operator
        const I = math.diag(Array(this.n_x).fill(-2), 0);
        const A = math.resize(math.diag(Array(this.n_x).fill(1), 1), [this.n_x, this.n_x]);
        const B = math.resize(math.diag(Array(this.n_x).fill(1), -1), [this.n_x, this.n_x]);
        const D2 = math.add(I, A, B);
        D2[0] = Array(this.n_x).fill(0); // set by BCs elsewhere
        D2[this.n_x-1] = Array(this.n_x).fill(0);
        return D2;
    }
}


// testing stuff
// state: S_O, S_I, S_S, S_NH, S_N2, S_NO, S_HCO, X_I, X_S, X_H, X_STO, X_A, X_TS
// let initial_state = [0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1];
// const Reactor = new Reactor_PFR(200, 10, 10, 1, 100, initial_state);
// Reactor.Cs_in[0] = [0.0, 30., 100., 16., 0., 0., 5., 25., 75., 30., 0., 0., 125.];
// Reactor.Fs[0] = 10;
// Reactor.D = 0.01;
// let N = 0;
// while (N < 5000) {
//     console.log(Reactor.tick_fe(0.001));
//     N += 1;
// }

module.exports = { Reactor_CSTR, Reactor_PFR };