#include <iostream>
#include <vector>
#include <math.h>
#include <iomanip>

using namespace std;

int main() 
{
    int l = 1925;               // Edge length
    int dx = 5;                 // Spatial step
    int kmax = l / dx;          // Number of spatial steps

    /**
     * Coefficients alpha and gamma are different for different 
     * edges. It depends on the edge's lenght, cross-sectional area etc.
    */
    double alpha = 0.744;
    double beta = 0.004;
    double gamma = 5829;

    /**
     * dt   - the number of approximation in Euler method
     * h    - discretization step in Euler method
     * 
     * Those two values work well with all existing steps and give us
     * result that satisfies the initial condition pretty close.
     * 
     * Changing of any of them might result in NaN error.
    */
    int dt = 1000;
    double h = 0.007;


    std::cout << "--------------------------------------------------" << std::endl;
    std::cout << "Simulation parameters" << std::endl;
    std::cout << "--------------------------------------------------" << std::endl;
    std::cout << "l:    " << l << std::endl;
    std::cout << "dx:   " << dx << std::endl;
    std::cout << "kmax: " << kmax << std::endl;
    std::cout << "dt:   " << dt << std::endl;
    std::cout << "h:    " << h << std::endl;
    std::cout << "--------------------------------------------------" << std::endl;

    // Stores the approximated values for each step 
    double Q[dt][kmax];
    double P[dt][kmax];

    // Prepare initial condition
    for (int i = 0; i < dt; i++) 
    {
        // Loop over spatial steps
        for (int k = 0; k < kmax; k++) 
        {
            // Boundary condition for the last spatial step
            if (k == kmax - 1) 
            {
                Q[i][k] = 0;
                P[i][k] = -202.0;   // P = -202 when Q = 10
            } else 
            {
                Q[i][k] = 0;
                P[i][k] = 0;
            }
        }
    }

    /**
     * In fact, we have two 'steps'. The first step represents the approximation 
     * criteria in Euler's method (dt) and the second one represents the actual
     * spatial step along the edge (dx).
    */
    for (int i = 2; i < dt; i++) 
    {
        // Loop over spatial steps for flow
        for (int k = 1; k < kmax; k++) {
            Q[i + 1][k] = Q[i][k] + h * (alpha * (P[i][k - 1] - P[i][k]) - beta * Q[i][k] * abs(Q[i][k]));
        }
        // Loop over steps for pressure, DO NOT calculate pressure for the last element as it is given as a boundary condition
        for (int k = 1; k < kmax - 1; k++)  {
            P[i + 1][k] = P[k][i - 2] + h * (gamma * (Q[i][k] - Q[i][k + 1]));
        }
    }
    
    // Print out the first 10 values of Q for each spatial step
    std::cout << "FIRST 10 VALUES\n";
    for(int i = 0; i < kmax; i++)
    {
        std::cout.width(5);
        std::cout.flags(std::ios::left);
        std::cout << i << ": ";
   
        for(int k = 0; k < 10; k++)
            std::cout << std::setw(10) << Q[k][i];
        std::cout << std::endl;
    }

    // Print out the last 10 values of Q for each spatial step
    std::cout << "LAST 10 VALUES:\n";
    for(int i = 1; i < kmax; i++)
    {
        std::cout.width(5);
        std::cout.flags(std::ios::left);
        std::cout << i << ": ";

        for(int k = 0; k < 10; k++)
            std::cout << std::setw(10) << Q[dt - 10 + k][i];
        std::cout << std::endl;
    }

    return 0;
}