#include <iostream>
#include <omp.h>
#include "Science.h"
#include "taos.h"

Science::Science() {
	data = NULL;
}

Science::Science(string name): SourceFits(name) {
	data = NULL;
}

Science::~Science() {
	if (data) {
		std::cout << "Liberando Science data ... ";
		free(data);
		data = NULL;
		std::cout << "Done" << std::endl;
	}
}

void Science::Init() {

	//Conversion de string a char[] para funcion de fits_open_file
	char* file = new char[filename.size() + 1];
	filename.copy(file, filename.size() + 1);
	file[filename.size()] = '\0';

	//Abrimos el archivo
	if (fits_open_file(&fptr, file, READONLY, &status)) {
		fits_report_error(stderr, this->status); //Imprime error de status
		exit(EXIT_FAILURE);
	}

	//Ancho y Largo de la imagen en pixeles
	if (fits_read_keys_lng(fptr, "NAXIS", 1, 3, naxes, &nfound, &status)) {
		fits_report_error(stderr, this->status); //Imprime error de status
		exit(EXIT_FAILURE);
	}

	Xs = naxes[0];
	Ys = naxes[1];
	Frames = naxes[2];

	FitsReadFile();

	//Cerramos el archivo
	if (fits_close_file(fptr, &status)) {
		fits_report_error(stderr, this->status); //Imprime error de status
		exit(EXIT_FAILURE);
	}

	delete []file;

}

void Science::FitsReadFile() {
	
	//Declaramos desde donde queremos comenzar a leer la imagen en fpixel:
	long fpixel[] = { 1, 1, 1 };

	//Debemos obtener el total de pixeles que contiene la imagen y la cantidad de frames del archivo
	long int nelements = Xs;	//Pixeles a leer por renglon

	//Tamaņo del arreglo donde se guardaran los datos:
	data = (double*)calloc(Xs * Ys * Frames, sizeof(double));
	long int* tempdata = (long int*)calloc(nelements, sizeof(long int));

	//nulval se de clara como short -1 y anynul se trabaja con un int 0; datatype para leer TULONG
	short nulval = -1;

	int frame, row;
	int XYs = Xs * Ys;

	if (Frames > 0) {

		//Recorremos cada frame
		for (int i = 0; i < naxes[2]; i++)
		{

			fpixel[2] = i + 1;
			frame = i;

			//Recorremos cada renglon
			for (int j = 0; j < nelements; j++)
			{
				fpixel[1] = j + 1;
				row = j;

				if (fits_read_pix(fptr, TLONG, fpixel, nelements, &nulval, tempdata, 0, &this->status)) {
					fits_report_error(stderr, this->status); //Imprime error de status
					exit(EXIT_FAILURE);
				}

				for (int k = 0; k < nelements; k++)
				{
					data[(frame * XYs) + (row * nelements) + k] = (double)tempdata[k];
				}
			}
		}
	}
	else {

		fpixel[2] = 1;

		//Recorremos cada renglon
		for (int j = 0; j < Ys; j++)
		{
			fpixel[1] = j + 1;
			row = j;

			if (fits_read_pix(this->fptr, TLONG, fpixel, nelements, &nulval, tempdata, 0, &this->status)) {
				fits_report_error(stderr, this->status); //Imprime error de status
				exit(EXIT_FAILURE);
			}

			for (int k = 0; k < nelements; k++)
			{
				data[(row * nelements) + k] = (double)tempdata[k];
			}
		}
	}

	free(tempdata);
}

void Science::CleanData(long* BiasMaster, long* DarksMaster, double* FlatsMaster) {
	
	int xs = Xs;
	int ys = Ys;
	int xys = xs * ys;
	int frames = Frames;
	int row;
	long indx = 0;
	long frame = 0;
	double * tempptr = data;

	std::cout << "Comenzando a limpiar pixeles: \n";

#pragma omp parallel for shared(xs, ys, frames, xys, tempptr) private(indx, frame, row) num_threads(THREADNUM)
	//Recorremos Frames
	for (int pf = 0; pf < frames; pf++)
	{

		frame = pf * xys;

		//Recorremos Y's
		for (int py = 0; py < ys; py++)
		{
			row = py * xs;

			//Recorremos X's
			for (int px = 0; px < xs; px++)
			{
				indx = px + row;

				tempptr[indx + frame] = (data[indx + frame] - BiasMaster[indx] - DarksMaster[indx]) / FlatsMaster[indx];
			}
		}
	}

}

double* Science::Data() {
	return data;
}

void Science::free_data() {
	if (data) {
		std::cout << "Liberando SCIENCE data ... ";
		free(data);
		data = NULL;
		std::cout << "Done \n";
	}
}

//Establecer ejes manualmente
void Science::setAxes(long xs, long ys, long frames) {
	Xs = xs;
	Ys = ys;
	Frames = frames;
}

//Cambio de archivo
void Science::ChangeFile(string newfile) {
	filename = newfile;

	if (data) {
		free(data);
		data = NULL;
	}

	Xs = 0;
	Ys = 0;
	Frames = 0;
}

void Science::setMemory(unsigned int pixels) {
	data = (double*)calloc(pixels, sizeof(double));
}

void Science::Print(int limx, int limy, int limframe) {

	std::cout << "\t" << filename << "\n";

	//Recorremos Frames
	for (int pf = 0; pf < limframe; pf++)
	{
		//Recorremos Y's
		for (int py = 0; py < limy; py++)
		{
			//Recorremos X's
			for (int px = 0; px < limx; px++)
			{
				std::cout << "[" << px << "]" << "[" << py << "]" << "[" << pf << "]: "
					<< data[px + (py * Ys) + (pf * Ys * Xs)] << "\n";
			}
		}
	}

	std::cout << "\n\n";

}