#include <cstddef>
#include <tug/Boundary.hpp>
#include <iostream>

using namespace std;

auto calc_alpha_intercell(double alpha1, double alpha2, bool useHarmonic = false) {
    if (useHarmonic) {
        return 2 / ((1/alpha1) + (1/alpha2));
    } else {
        return 0.5 * (alpha1 + alpha2);
    }
}

MatrixXd FTCS_constant(Grid grid, Boundary bc, double timestep) {
    int rowMax = grid.getRow();
    int colMax = grid.getCol();
    double deltaRow = grid.getDeltaRow();
    double deltaCol = grid.getDeltaCol();

    // Matrix with concentrations at time t+1
    // TODO profiler / only use 2 matrices
    MatrixXd concentrations_t1 = MatrixXd::Constant(rowMax, colMax, 1);

    // inner cells
    cout << "Concentration 5,5: " << grid.getConcentrations()(5,5) << endl;
    cout << "Alpha Y 5,5: " << grid.getAlphaY()(5,5) << endl;
    cout << "calc alpha Y 5,5; 5,6: " << calc_alpha_intercell(grid.getAlphaY()(5,5), grid.getAlphaY()(5,6)) << endl;
    cout << "t1 Concentrations 5,5: " << concentrations_t1(5,5) << endl;
    for (int row = 1; row < rowMax-1; row++) {
        for (int col = 1; col < colMax-1; col++) {
            concentrations_t1(row, col) = grid.getConcentrations()(row, col) 
                + timestep / (deltaRow*deltaRow) * (
                    calc_alpha_intercell(grid.getAlphaY()(row+1,col), grid.getAlphaY()(row,col)) 
                    * grid.getConcentrations()(row+1,col)
                    - (calc_alpha_intercell(grid.getAlphaY()(row+1,col), grid.getAlphaY()(row,col))
                    + calc_alpha_intercell(grid.getAlphaY()(row-1,col), grid.getAlphaY()(row,col)))
                    * grid.getConcentrations()(row,col)
                    + calc_alpha_intercell(grid.getAlphaY()(row-1,col), grid.getAlphaY()(row,col))
                    * grid.getConcentrations()(row-1,col)
                )
                + timestep / (deltaCol*deltaCol) * (
                    calc_alpha_intercell(grid.getAlphaX()(row,col+1), grid.getAlphaX()(row,col)) 
                    * grid.getConcentrations()(row,col+1)
                    - (calc_alpha_intercell(grid.getAlphaX()(row,col+1), grid.getAlphaX()(row,col))
                    + calc_alpha_intercell(grid.getAlphaX()(row,col-1), grid.getAlphaX()(row,col)))
                    * grid.getConcentrations()(row,col)
                    + calc_alpha_intercell(grid.getAlphaX()(row,col-1), grid.getAlphaX()(row,col))
                    * grid.getConcentrations()(row,col-1)
                );
        }
    }

    // boundary conditions
    // left without corners / looping over rows
    int col = 0;
    // for (int row = 1; row < rowMax-1; row++) {
    //     concentrations_t1(row, col) = grid.getConcentrations()(row,col)
    //         + timestep / (deltaCol*deltaCol) 
    //             * (calc_alpha_intercell(grid.getAlphaX()(row,col+1), grid.getAlphaX()(row,col)) 
    //             * grid.getConcentrations()(row,col+1)
    //             - (calc_alpha_intercell(grid.getAlphaX()(row,col+1), grid.getAlphaX()(row,col)) 
    //             + 2 * grid.getAlphaX()(row,col)) * grid.getConcentrations()(row,col) 
    //             + 2 * grid.getAlphaX()(row,col) * bc.getBoundaryConditionValue(BC_SIDE_LEFT)(row))
    //         + timestep / (deltaRow*deltaRow)
    //             * (calc_alpha_intercell(grid.getAlphaY()(row+1,col), grid.getAlphaY()(row,col)) 
    //             * grid.getConcentrations()(row+1,col)
    //             - (calc_alpha_intercell(grid.getAlphaY()(row+1,col), grid.getAlphaY()(row,col))
    //             + calc_alpha_intercell(grid.getAlphaY()(row-1,col), grid.getAlphaY()(row,col)))
    //             * grid.getConcentrations()(row,col)
    //             + calc_alpha_intercell(grid.getAlphaY()(row-1,col), grid.getConcentrations()(row,col))
    //             * grid.getConcentrations()(row-1,col));
    // }

    // right without corners / looping over columns
    col = colMax-1;
    // for (int row = 1; row < rowMax-1; row++) {
    //     concentrations_t1(row,col) = grid.getConcentrations()(row,col)
    //         + timestep / (deltaCol*deltaCol) 
    //             * (2 * grid.getAlphaX()(row,col) * bc.getBoundaryConditionValue(BC_SIDE_RIGHT)(row)
    //             - (calc_alpha_intercell(grid.getAlphaX()(row,col-1), grid.getAlphaX()(row,col))
    //             + 2 * grid.getAlphaX()(row,col)) + 2 * grid.getAlphaX()(row,col) 
    //             * grid.getConcentrations()(row,col)
    //             + calc_alpha_intercell(grid.getAlphaX()(row,col-1), grid.getAlphaX()(row,col))
    //             * grid.getConcentrations()(row,col-1))
    //         + timestep / (deltaRow*deltaRow)
    //             * (calc_alpha_intercell(grid.getAlphaY()(row+1,col), grid.getAlphaY()(row,col))
    //             * grid.getConcentrations()(row+1,col)
    //             - (calc_alpha_intercell(grid.getAlphaY()(row+1,col), grid.getAlphaY()(row,col)) 
    //             + calc_alpha_intercell(grid.getAlphaY()(row-1,col), grid.getAlphaY()(row,col)))
    //             * grid.getConcentrations()(row,col)
    //             + calc_alpha_intercell(grid.getAlphaY()(row-1,col), grid.getAlphaY()(row,col))
    //             * grid.getConcentrations()(row-1,col));
    // }


    // top without corners / looping over cols
    // for(int col=1; col<colMax-1;col++){
    //     int row = 0;
    //     concentrations_t1(row, col) = grid.getConcentrations()(row, col)
    //         + timestep/(grid.getDeltaRow()*grid.getDeltaRow()) * (calc_alpha_intercell(grid.getAlphaY()(1, col), grid.getAlphaY()(0, col)) * grid.getConcentrations()(1,col)
    //             - (calc_alpha_intercell(grid.getAlphaY()(1, col), grid.getAlphaY()(0, col)) + 2 * grid.getAlphaY()(0, col)) * grid.getConcentrations()(0, col)
    //             + 2 * grid.getAlphaY()(0, col) * bc.getBoundaryConditionValue(BC_SIDE_TOP)(col))
    //         + timestep/(grid.getDeltaCol()*grid.getDeltaCol()) * (calc_alpha_intercell(grid.getAlphaX()(0, col+1), grid.getAlphaX()(0, col)) * grid.getConcentrations()(0, col+1)
    //             - (calc_alpha_intercell(grid.getAlphaX()(0, col+1), grid.getAlphaX()(0, col)) + calc_alpha_intercell(grid.getAlphaX()(0, col-1), grid.getAlphaX()(0, col))) * grid.getConcentrations()(0, col)
    //             + calc_alpha_intercell(grid.getAlphaX()(0, col-1), grid.getAlphaX()(0, col)) * grid.getConcentrations()(0, col-1));
    // }

    // bottom without corners / looping over cols
    int row = rowMax-1;
    // for(int col=1; row<colMax-1;col++){
    //     concentrations_t1(row, col) = grid.getConcentrations()(row, col)
    //         + timestep/(grid.getDeltaRow()*grid.getDeltaRow()) * (2 * grid.getAlphaY()(row, col) * bc.getBoundaryConditionValue(BC_SIDE_BOTTOM)(col)
    //             - (calc_alpha_intercell(grid.getAlphaY()(row, col), grid.getAlphaY()(row-1, col)) + 2 * grid.getAlphaY()(row, col)) * grid.getConcentrations()(row, col)
    //             + calc_alpha_intercell(grid.getAlphaY()(row, col), grid.getAlphaY()(row-1, col)))
    //         + timestep/(grid.getDeltaCol()*grid.getDeltaCol()) * (calc_alpha_intercell(grid.getAlphaX()(row, col+1), grid.getAlphaX()(row, col)) * grid.getConcentrations()(row, col+1)
    //             - (calc_alpha_intercell(grid.getAlphaX()(row, col+1), grid.getAlphaX()(row, col)) + calc_alpha_intercell(grid.getAlphaX()(row, col-1), grid.getAlphaX()(row, col))) * grid.getConcentrations()(row, col)
    //             + calc_alpha_intercell(grid.getAlphaX()(row, col-1), grid.getAlphaX()(row, col-1)) * grid.getConcentrations()(row, col-1));
    // }


    concentrations_t1(0,0) = 0;
    concentrations_t1(rowMax-1,0) = 0;
    concentrations_t1(0,colMax-1) = 0;
    concentrations_t1(rowMax-1,colMax-1) = 0;

    return concentrations_t1;
}

void FTCS_closed(Grid grid, Boundary bc, double timestep) {
    return;
}

MatrixXd FTCS(Grid grid, Boundary bc, double timestep) {
    if (bc.getBoundaryConditionType() == BC_TYPE_CONSTANT) {
        return FTCS_constant(grid, bc, timestep);
    } else if (bc.getBoundaryConditionType() == BC_TYPE_CLOSED) {
        FTCS_closed(grid, bc, timestep);
    }
}