#pragma once

#include <span>
#include <fstream>

namespace export_format {

template <std::floating_point T>
void vtk(Cuboid c, std::span<T> rho, std::span<T> u, std::string path) {
	std::ofstream fout;
	fout.open(path.c_str());

	fout << "# vtk DataFile Version 3.0\n";
	fout << "lbm_output\n";
	fout << "ASCII\n";
	fout << "DATASET RECTILINEAR_GRID\n";
	fout << "DIMENSIONS " << c[0] << " " << c[1] << " " << c[2] << "\n";

	fout << "X_COORDINATES " << c[0] << " float\n";
	for (std::size_t iX=0; iX < c[0]; ++iX) {
		fout << iX << " ";
	}

	fout << "\nY_COORDINATES " << c[1] << " float\n";
	for (std::size_t iY=0; iY < c[1]; ++iY) {
		fout << iY << " ";
	}

	fout << "\nZ_COORDINATES " << c[2] << " float\n";
	for (std::size_t iZ=0; iZ < c[2]; ++iZ) {
		fout << iZ << " ";
	}
	fout << "\nPOINT_DATA " << c.volume() << "\n";

	fout << "VECTORS velocity float\n";
	for (std::size_t iZ=0; iZ < c[2]; ++iZ) {
		for (std::size_t iY=0; iY < c[1]; ++iY) {
			for (std::size_t iX=0; iX < c[0]; ++iX) {
				const std::size_t iCell = c.index(iX,iY,iZ);
				fout << u[iCell] << " " << u[c.volume() + iCell] << " " << u[2*c.volume() + iCell] << "\n";
			}
		}
	}

	fout << "SCALARS density float 1\n";
	fout << "LOOKUP_TABLE default\n";
	for (std::size_t iZ=0; iZ < c[2]; ++iZ) {
		for (std::size_t iY=0; iY < c[1]; ++iY) {
			for (std::size_t iX=0; iX < c[0]; ++iX) {
				const std::size_t iCell = c.index(iX,iY,iZ);
				fout << rho[iCell] << "\n";
			}
		}
	}

	fout.close();
}

}