#include<iostream>
#include<vector>
#include<fstream>
#include<algorithm>
#include<string>
#include<cmath>
using namespace std;

#define pi 3.14159265358979
#define eps 1e-5

const int ni = 100,nj = 80;
const int n_pt = 3;
const double L = 20,H = 5;

const double d_xi = 1.0/(ni-1),d_eta = 1.0/(nj-1);

// PT có dạng U E F

double data_t[n_pt+1][ni+1][nj+1]; // data tại thời điểm t hay vector U
double data_t_dt[n_pt+1][ni+1][nj+1]; // data tại thời điểm t + d_t hay vector U*
double G[n_pt+1][ni+1][nj+1]; // mang do nhot gia
double E[n_pt+1][ni+1][nj+1]; // vector E
double F[n_pt+1][ni+1][nj+1]; // vector F

// Các thông số vật lý cho các điểm lưới
double u[ni+1][nj+1], v[ni+1][nj+1], Density[ni+1][nj+1], Pressure[ni+1][nj+1], sound_speed[ni+1][nj+1];
double R_gas_constant = 333.33,Cx = 0.2,Cy = 0.2,Cp = 1007; // thông số nhiệt động

// Grid OXY
double X[ni+1][nj+1], Y[ni+1][nj+1];

//Grid O xi eta
double Jacobi_1[ni+1][nj+1]; //Jacobi trong biến đổi từ Oxy -> O xi eta
double Y_eta[ni+1][nj+1], X_eta[ni+1][nj+1], Y_xi[ni+1][nj+1], X_xi[ni+1][nj+1];// Cần đc tính toán
double xi_X[ni+1][nj+1], xi_Y[ni+1][nj+1], eta_X[ni+1][nj+1], eta_Y[ni+1][nj+1];// Grid CFD

//Thông số Inlet
double Temperature_0,Pressure_0,Density_0,M_0,u_0;

//Input
double DH = 1.0, d_t = 1e-4;
double NA,NAF,ITOUT,ITCLD,ITPRI;

// PT dang U E F

void InputData();
void MaccorMack();
void output();
void f_sound_speed();
void U_t(); // vector U
void U_t_dt(); // vector U*
void E_comma(); // ma trận E' sau biến đổi tọa độ
void F_comma(); // ma trận F' sau biến đổi tọa độ
void fake_viscosity(); // tính độ nhớt giả
void return_variables(); // trả lại ma trận U gồm density U V với dòng 2D , ...
void boundary_condition(); // khai báo điều kiện biên
void boundary_wall(); //
void transform_xi(double old_[ni+1][nj+1],double new_[ni+1][nj+1]); // đổi tọa độ theo d_xi
void transform_eta(double old_[ni+1][nj+1],double new_[ni+1][nj+1]); // đổi tọa độ theo e_ta

int main()
{

    return 0;
}

void InputData()
{
    ifstream openfile("D:/VSCode/HW/MCM/INPUT.dat");
    string tmp;
    openfile >> tmp >> NA;
    openfile >> tmp >> NAF;
    openfile >> tmp >> ITOUT;
    openfile >> tmp >> ITCLD;
    openfile >> tmp >> ITPRI;
    openfile >> tmp >> d_t;
    openfile >> tmp >> Pressure_0;
    openfile >> tmp >> Temperature_0;
    openfile >> tmp >> u_0;
    openfile.close();

    /**
     * @brief tạo file chia lưới ở đây: X và Y
     * 
     */

    for(int j=1;j<=nj;j++)
    {
        for(int i=2;i<=ni;i++)
        {
            if(i<=50)
            {
                Pressure[i][j] = 100000;
                Density[i][j] = 1;
                u[i][j] = 0;
            }
            else
            {
                Pressure[i][j] = 10000;
                Density[i][j] = 0.1;
                u[i][j] = 0;
            }
        }
    }

    // tinh cac gia tri bo tro
    transform_xi(X,X_xi);
    transform_xi(Y,Y_xi);
    transform_eta(X,X_eta);
    transform_eta(Y,Y_eta);

    // buoc chuyen doi he toa do
    for (int j=1;j<=nj;j++)
	{
		for (int i=1;i<=ni;i++)
		{
			Jacobi_1[i][j] = X_xi[i]*Y_eta[i][j] - Y_xi[i][j]*X_eta[i][j];
            xi_X[i][j] = Y_eta[i][j] / Jacobi_1[i][j];
            xi_Y[i][j] = -1.0*X_eta[i][j] / Jacobi_1[i][j];
            eta_X[i][j] = -1.0*Y_xi[i][j] / Jacobi_1[i][j];
            eta_Y[i][j] = X_xi[i][j] / Jacobi_1[i][j];
		}
	}

}

void transform_xi(double old_[ni+1][nj+1],double new_[ni+1][nj+1]) // 
{
    for(int j=1;j<=nj;j++)
    {
        for(int i=2;i<=ni-1;i++)
        {
            new_[i][j] = (old_[i+1][j] - old_[i-1][j])/(2.0*d_xi);
        }
    }

    for(int j=1;j<=nj;j++)
    {
        new_[1][j] = (old_[2][j] - old_[1][j])/(2.0*d_xi);
        new_[ni][j] = (old_[ni][j] - old_[ni-1][nj])/(2.0*d_xi);
    }
}

void transform_eta(double old_[ni+1][nj+1],double new_[ni+1][nj+1])
{
    for(int j=2;j<=nj-1;j++)
    {
        for(int i=1;i<=ni;i++)
        {
            new_[i][j] = (old_[i][j+1] - old_[i][j-1])/(2.0*d_eta);
        }
    }

    for(int i=1;i<=ni;i++)
    {
        new_[i][1] = (old_[i][2] - old_[i][1])/(2.0*d_eta);
        new_[i][nj] = (old_[i][nj] - old_[i][nj-1])/(2.0*d_eta);
    }
}

void f_sound_speed()
{
    for(int j=1;j<=nj;j++)
    {
        for(int i=1;i<=ni;i++)
        {
            sound_speed[i][j] = sqrtf(1.4*R_gas_constant*Temperature_0);
        }
    }
}

void U_t()
{
	for (int j = 1; j <= Ny; j++)
	{
		for (int i = 1; i <= Nx; i++)
		{
			data_t[1][i][j] = Density[i][j]*Jacobi_1[i][j];
			data_t[2][i][j] = data_t[1][i][j]*u[i][j]*Jacobi_1[i][j];
			data_t[3][i][j] = data_t[1][i][j]*v[i][j]*Jacobi_1[i][j];
		}
	}
}

void U_t_dt()
{
    for (int j=1;j<=nj;j++)
	{
		for (int i=1;i<=ni;i++)
		{
			data_t_dt[1][i][j] = Density[i][j]*Jacobi_1[i][j];
			data_t_dt[2][i][j] = data_t_dt[1][i][j]*u[i][j]*Jacobi_1[i][j];
			data_t_dt[3][i][j] = data_t_dt[1][i][j]*v[i][j]*Jacobi_1[i][j];
		}
	}
}

void E_comma() // hàm invisX hay vector E'
{
    
}

void F_comma() // hàm invisY hay vector F'
{

}

/**
 * @brief 
 * ? Giả sử dòng subsonic
 * ? 
 * 
 */
void boundary_condition() // khai báo điều kiện biên
{
    //Inlet
    for(int j=1;j<=nj;j++)
    {
        u[1][j] = u[2][j];
        v[1][j] = 0;
        Pressure[1][j] = Pressure[2][j];
        Density[1][j] = Pressure[1][j]/(R_gas_constant*Temperature_0);
    }

    //Outlet
    for (int j=1;j<=nj;j++)
	{
		u[ni][j] = u[ni - 1][j];
		V[ni][j] = v[ni - 1][j];
		Pressure[ni][j] = Pressure[ni - 1][j];
		Density[ni][j] = Pressure[ni][j]/(R_gas_constant*Temperature_0);
	}
}

/**
*? chất lỏng có nhớt thì tại tường u = 0, v = 0, w = 0 và đhr P / đhr (biến) = 0;
*? chất lỏng không nhớt thì u v w khác 0, cần đổi u v w sang hệ tọa độ tính toán
*? u = u_xi và v = v_eta
 * 
 */
void boundary_wall()
{

}