//             by: Alexander Reeves

const ctx = document.getElementById("myChart");
const ctx2 = document.getElementById("myChart2");

var myChart = new Chart(ctx, {});
var myChart2 = new Chart(ctx2, {});

var labels = []; // labels for x axis

var t = 0; // time variable in s

// ----------- Assign Parameters -----------

// must be 'var' if parameters are changed during runtime ( see updateParameters() )
var T = 723.15; // Temperature in K -> isothermal case
var xN2_in = document.getElementById("input_x_N2").value; // inlet composition of component N2
var xH2_in = document.getElementById("input_x_H2").value; // inlet composition of component H2
var xNH3_in = document.getElementById("input_x_NH3").value; // inlet composition of component NH3
var V = 10; // reactor volume in m^3
var R = 8.314; // universal gas constant J/K⋅mol
var P = 100000 * 150; // pressure in Pa -> isobaric
var F_in = 100000; //(P*V)/(R*T);       // inlet molar flowrate in mol/s (approx defined as total reactor content)

function updateParameters() {
  // parameters can be re-assigned dynamically via if statements related to the decision variable
  // example:
  //if (t >= 100) { F_in = F_in*1.2 }  // 20% increase in inlet flowrate at time 100s
  //if (t >= 2) { P = (t-2)/(4-2)*100000*10 }  // Gradual pressure change from P_initial to 10 bar over time from 2s to 4s { (t-t_start)/(t_end-t_start) * target_value }
}

// ----------- Initialize Model variables -----------

// algebraic:
var var_x = 0; // example
var x_N2 = 1; // molar fraction of N2
var x_H2 = 0; // molar fraction of H2
var x_NH3 = 0; // molar fraction of NH3
var M_t = 0; // total mass in system
//var M_N2 = 0;       // molar holdup of N2
//var M_H2 = 0;       // molar holdup of H2
//var M_NH3 = 0;      // molar holdup of NH3
var r = 0; // reaction rate (Ammonia production: N_2 + 3H_2 => 2NH_3) in kmol/(m3*s)
var p_N2 = 0; // partial pressure of component i
var p_H2 = 0; // partial pressure of component i
var p_NH3 = 0; // partial pressure of component i
//var P = 100000;          // total pressure in system
var rho = 0; // molar density in system
var F_out = 0; //F_in/100;   // constant flowrate

// differential vbles arrays:
var M_N2 = [F_in * 1]; // molar holdup of N2
var M_H2 = [0]; // molar holdup of H2
var M_NH3 = [0]; // molar holdup of N2

// ----------- Initialize dataset arrays (for graphing) -----------

var graph_x_N2 = [];
var graph_x_H2 = [];
var graph_x_NH3 = [];

// ############### Algebraic Equations ###############

function solveAlgebraic(i) {
  M_t = M_N2[i] + M_H2[i] + M_NH3[i];
  x_N2 = M_N2[i] / M_t;
  x_H2 = M_H2[i] / M_t;
  x_NH3 = M_NH3[i] / M_t;
  p_N2 = x_N2 * P;
  p_H2 = x_H2 * P;
  p_NH3 = x_NH3 * P;
  r =
    (Math.exp(-91000 / (R * T)) * Math.pow(p_N2, 0.5) * Math.pow(p_H2, 1.5) -
      2.57 * Math.pow(10, 5) * Math.exp(-140000 / (R * T)) * p_NH3) /
    1000; // rate equation kmol/(m3*s) -> /1000 for mol

  rho = M_t / V;
  F_out = M_t; // outlet flowrate = reactor content -> no increase or loss

  console.log("r: " + r);
  console.log("x_H2: " + x_H2);
  console.log("p_H2: " + p_H2);

  //console.log('Algebraic variables: ');
  //console.log('Variable X: '+var_x);
}

// ############### ODE Equations ###############

function f_diff(t, y, functionSelect) {
  t == 0;
  y == 0;

  // return the value of y' for the selected function from the ODE function array
  if (functionSelect == 1) {
    return F_in * xN2_in - F_out * x_N2 + -0.5 * r * V;
  }
  if (functionSelect == 2) {
    return F_in * xH2_in - F_out * x_H2 + -1.5 * r * V;
  }
  if (functionSelect == 3) {
    return F_in * xNH3_in - F_out * x_NH3 + 1 * r * V;
  } else {
    return 0;
  }
}

// ############### ODE Solver ###############

function rk4iter(t, y, h, functionSelect) {
  // compute the four RK coefficients using the RHS of the diff eq.
  k1 = f_diff(t, y, functionSelect);
  k2 = f_diff(t + h / 2, y + (k1 * h) / 2, functionSelect);
  k3 = f_diff(t + h / 2, y + (k2 * h) / 2, functionSelect);
  k4 = f_diff(t + h, y + k3 * h, functionSelect);

  // compute the y value at the end of the iteration
  return y + ((k1 + 2 * k2 + 2 * k3 + k4) * h) / 6;
}

function rk4(n, t_end) {
  h = t_end / n; // compute the step size along the x-axis

  // add initial conditions to array
  t = 0; // initial value of x is 0

  labels.push(0);

  // loop until x reaches its target value
  for (var i = 1; i <= n; i++) {
    // check parameter states
    updateParameters();

    // re-solve algebraic values for each iteration
    solveAlgebraic(i - 1);

    // Perform iterations of the RK4 algorithm
    M_N2[i] = rk4iter(t, M_N2[i - 1], h, 1);
    M_H2[i] = rk4iter(t, M_H2[i - 1], h, 2);
    M_NH3[i] = rk4iter(t, M_NH3[i - 1], h, 3);

    // check for exceptions
    if (M_N2[i] < 0) {
      M_N2[i] = 0.0001;
      //r = 0;
    }
    if (M_H2[i] < 0) {
      M_H2[i] = 0.0001;
      //r = 0;
    }
    if (M_NH3[i] < 0) {
      M_NH3[i] = 0.0001;
      //r = 0;
    }

    t += h; // time step

    // record the x and y values for this iteration
    console.log(t + "s, " + M_N2[i] + ", " + M_H2[i] + ", " + M_NH3[i]);

    // Add data to the plot as well.
    labels.push(t.toFixed(2));

    // Add additional graphing variables
    graph_x_N2.push(x_N2);
    graph_x_H2.push(x_H2);
    graph_x_NH3.push(x_NH3);
  }
}

function run() {
  // check inputs, give warning
  if (
    Number(document.getElementById("input_x_N2").value) +
      Number(document.getElementById("input_x_H2").value) +
      Number(document.getElementById("input_x_NH3").value) !=
    1
  ) {
    for (
      i = 0;
      i <= document.getElementsByClassName("inputBox").length - 1;
      i++
    ) {
      document.getElementsByClassName("inputBox")[i].style.backgroundColor =
        "red";
    }
    setTimeout(function () {
      for (
        i = 0;
        i <= document.getElementsByClassName("inputBox").length - 1;
        i++
      ) {
        document.getElementsByClassName("inputBox")[i].style.backgroundColor =
          "initial";
      }
    }, 100);
    setTimeout(function () {
      for (
        i = 0;
        i <= document.getElementsByClassName("inputBox").length - 1;
        i++
      ) {
        document.getElementsByClassName("inputBox")[i].style.backgroundColor =
          "red";
      }
    }, 200);
    setTimeout(function () {
      for (
        i = 0;
        i <= document.getElementsByClassName("inputBox").length - 1;
        i++
      ) {
        document.getElementsByClassName("inputBox")[i].style.backgroundColor =
          "initial";
      }
      //location.reload(); //disabled by codepen
    }, 300);
  }

  // empty data
  clearAll();

  // runs 4th order runge kutta method (steps,y_initial,x_end)
  rk4(document.getElementById("myRange").value, 10);

  // re-loads chart.js graph
  loadGraph();
}

// ############### VISUALISATION ###############
//           graphing using Chart.js

function loadGraph() {
  // destroy first chart
  myChart.destroy();
  myChart2.destroy();

  // run Chart.js graph
  graph(labels);

  // empty data
  clearAll();
}

function clearAll() {
  labels = [];

  xN2_in = Number(document.getElementById("input_x_N2").value); // inlet composition of component N2
  xH2_in = Number(document.getElementById("input_x_H2").value); // inlet composition of component H2
  xNH3_in = Number(document.getElementById("input_x_NH3").value); // inlet composition of component NH3

  if (xN2_in + xH2_in + xNH3_in != 1) {
    // in case of incorrect inputs the values are scaled
    xN2_in = 1 - (xH2_in + xNH3_in);
  }

  M_N2 = [F_in]; // molar holdup of N2 (initial condition: 100% nitrogen)
  M_H2 = [0]; // molar holdup of H2
  M_NH3 = [0]; // molar holdup of NH3
  x_N2 = 1; // molar fraction of N2 (initial condition: 100% nitrogen)
  x_H2 = 0; // molar fraction of H2
  x_NH3 = 0; // molar fraction of NH3
  M_t = 0; // total mass in system
  r = 0; // reaction rate (Ammonia production: 0.5N_2 + 1.5H_2 => NH_3)
  p_N2 = 0; // partial pressure of component i
  p_H2 = 0; // partial pressure of component i
  p_NH3 = 0; // partial pressure of component i
  rho = 0; // molar density in system
  F_out = F_in; // constant flowrate

  // clear graphing data
  graph_x_N2 = [];
  graph_x_H2 = [];
  graph_x_NH3 = [];
}

function graph(inputLabels) {
  myChart = new Chart(ctx, {
    type: "line",
    data: {
      labels: inputLabels,
      datasets: [
        {
          label: "x_N2",
          data: graph_x_N2,
          borderWidth: 2,
          pointStyle: false
        },
        {
          label: "x_H2",
          data: graph_x_H2,
          borderWidth: 2,
          pointStyle: false
        },
        {
          label: "x_NH3",
          data: graph_x_NH3,
          borderWidth: 2,
          pointStyle: false
        }
      ]
    },
    options: {
      indexAxis: "x",
      scales: {
        y: {
          beginAtZero: true,
        title: {
                display: true,
                text: 'outlet composition'
              }
        },
        x: {
            title: {
                display: true,
                text: 'time (s)'
              }
        }
      }
    }
  });

  myChart2 = new Chart(ctx2, {
    type: "line",
    data: {
      labels: inputLabels,
      datasets: [
        {
          label: "M_N2",
          data: M_N2,
          borderWidth: 2,
          pointStyle: false
        },
        {
          label: "M_H2",
          data: M_H2,
          borderWidth: 2,
          pointStyle: false
        },
        {
          label: "M_NH3",
          data: M_NH3,
          borderWidth: 2,
          pointStyle: false
        }
      ]
    },
    options: {
      indexAxis: "x",
      scales: {
        y: {
          beginAtZero: true,
          title: {
                display: true,
                text: 'outlet flowrate (mol/s)'
              }
        },
        x: {
            title: {
                display: true,
                text: 'time (s)'
              }
        }
      }
    }
  });
}