#pragma once
#include <stdio.h>
#include <math.h>
#include "taos.h"

using namespace std;

void Photometry(string filename, Star star, int lz);
void add_to_file(string filename, vector<double> curve_by_frame, vector<double> snr);
double* readFile(string filename, long int &total_frames);
double get_curve(double* data, long int frame, int outter_ring_diameter, int* circle, int* ring);
int* create_mask_circle(int diameter1, int diameter2);
int* create_mask_ring(int diameter1, int diameter2);
vector<long int> size_snr(long int total_frames, int lz);
vector<vector<double>> sign_noise(double* data, long int total_frames, int outter_ring_diameter, vector<long int> z, int lz, int* circle, int* ring);
double sign_curve(vector<double> data, long int total_frames);
void supr_outfl(vector<double>& curve_by_frame);

void Photometry(string filename, Star star, int lz) {

	//Leemos la información guardada en el archivo fits y obtenemos la cantidad de frames a trabajar
	long int total_frames;
	double* data = readFile(filename, total_frames);

	//Definimos el tamaño de los anillos interior y exterior
	int outter_ring_diameter = star.getSquareSide();
	int inner_ring_diameter = (outter_ring_diameter / sqrt(2)) + 0.5;

	//Obtenemos las matrices de mascara para obtener la curva
	int* circle = create_mask_circle(outter_ring_diameter, inner_ring_diameter);
	int* ring = create_mask_ring(outter_ring_diameter, inner_ring_diameter);

	//Creamos el vector donde guardaremos los valores de la curva de luz
	vector<double> curve_by_frame;

	for (long int frame = 0; frame < total_frames; frame++)
	{
		curve_by_frame.push_back(get_curve(data, frame, outter_ring_diameter, circle, ring));
	}

	supr_outfl(curve_by_frame);

	vector<double> snr;
	snr.push_back(sign_curve(curve_by_frame, total_frames));

	//Si existe LZ >= 1, entonces realizamos lo siguiente:
	if (lz >= 1) {

		//Obtenemos la señal a ruido
		vector<long int> z = size_snr(total_frames, lz);
		vector<vector<double>> C_SNR = sign_noise(data, total_frames, outter_ring_diameter, z, lz, circle, ring);

		for (int i = 0; i < lz; i++)
		{
			snr.push_back(sign_curve(C_SNR[i], (total_frames / z[i])));
		}

	}

	//Liberamos la memoria de las mascaras
	free(circle);
	free(ring);

	//Agregamos los datos al archivo hdf5 correspondiente
	add_to_file(filename, curve_by_frame, snr);

}

void add_to_file(string filename, vector<double> curve_by_frame, vector<double> snr) {

	herr_t status;

	//Abrimos archivo
	hid_t file_id;

	file_id = H5Fopen(filename.c_str(), H5F_ACC_RDWR, H5P_DEFAULT);

	//Abrimos los grupos que usaremos
	hid_t rgroup1_id = H5Gopen2(file_id, CURVE_CENTER, H5P_DEFAULT);
	hid_t rgroup3_id = H5Gopen2(file_id, DATA, H5P_DEFAULT);

	hdf5_save_vector2(rgroup1_id, CURVE, curve_by_frame.data(), curve_by_frame.size());
	hdf5_save_vector2(rgroup3_id, SNR, snr.data(), snr.size());

	//Cerramos grupos
	status = H5Gclose(rgroup1_id);
	status = H5Gclose(rgroup3_id);

	//Cerramos archivo
	status = H5Fclose(file_id);

}

double* readFile(string filename, long int &total_frames) {

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

	hid_t file_id, dataset_id;
	herr_t status;

	//Abrimos el archivo
	file_id = H5Fopen(file, H5F_ACC_RDONLY, H5P_DEFAULT);

	//Abrimos el grupo
	hid_t group_id = H5Gopen(file_id, STAR_CENTER, H5P_DEFAULT);

	//Abrimos el dataset que nos interesa
	dataset_id = H5Dopen2(group_id, STAR_DATA, H5P_DEFAULT);

	//Conseguimos dimensiones de la informacion
	hid_t dataspace_id;
	hsize_t dim[5];
	int ndims;

	dataspace_id = H5Dget_space(dataset_id);
	ndims = H5Sget_simple_extent_dims(dataspace_id, dim, NULL);

	total_frames = (long int)dim[0]; //Asignamos los frames totales a trabajar

	//Apartamos memoria para guardar la informacion
	double* data = (double*)calloc(dim[0] * dim[1] * dim[2], sizeof(double));
	//Guardamos la informacion en data
	status = H5Dread(dataset_id, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, data);

	//Cerramos los id's
	status = H5Sclose(dataspace_id);
	status = H5Dclose(dataset_id);
	status = H5Gclose(group_id);
	status = H5Fclose(file_id);

	return data;
}

double get_curve(double * data, long int frame, int outter_ring_diameter, int * circle, int * ring) {

	double A, B;
	double curve;

	A = 0.0;
	B = 0.0;

	//Sumamos todos los datos que solo pertenecen a la estrella
	int matrix_dims = outter_ring_diameter * outter_ring_diameter;
	
	for (int i = 0; i < matrix_dims; i++)
	{
		A += data[(frame * matrix_dims) + i] * circle[i];
		B += data[(frame * matrix_dims) + i] * ring[i];
	}

	curve = -2.5 * log10(A - B);

	return curve;
}

//Recibimos de entrada el diametro de la caja y el diametro de la estrella
int* create_mask_circle(int diameter1, int diameter2) {

	float r1 = (diameter1 + 0.2) / 2;
	float r2 = (diameter2 + 0.2) / 2;
	int* mask = (int*)calloc(diameter1 * diameter1, sizeof(int));

	float a, b;
	float y1, y2;

	a = diameter1 / 2;
	b = a;

	for (int i = 0; i < diameter1; i++)
	{
		y2 = (b + sqrt((r2 * r2) - ((i - a) * (i - a))));
		y1 = (b - sqrt((r2 * r2) - ((i - a) * (i - a))));

		for (int j = 0; j < diameter1; j++)
		{
			if (j <= y2 && j >= y1) {
				mask[(i * diameter1) + j] = 1;
			}
			else {
				mask[(i * diameter1) + j] = 0;
			}
		}
	}

	return mask;
}

//Recibimos de entrada el diametro de la caja y el diametro de la estrella
int* create_mask_ring(int diameter1, int diameter2) {

	float r1 = (diameter1 + 0.2) / 2;
	float r2 = (diameter2 + 0.2) / 2;
	int* mask = (int*)calloc(diameter1 * diameter1, sizeof(int));

	float a, b;
	float y1, y2;
	float y1_b, y2_b;

	a = diameter1 / 2;
	b = a;

	for (int i = 0; i < diameter1; i++)
	{
		y2 = (b + sqrt((r1 * r1) - ((i - a) * (i - a))));
		y1 = (b - sqrt((r1 * r1) - ((i - a) * (i - a))));

		y2_b = (b + sqrt((r2 * r2) - ((i - a) * (i - a))));
		y1_b = (b - sqrt((r2 * r2) - ((i - a) * (i - a))));

		for (int j = 0; j < diameter1; j++)
		{
			if (j <= y2 && j >= y1) {
				mask[(i * diameter1) + j] = 1;
			}
			else {
				mask[(i * diameter1) + j] = 0;
			}

			if (j <= y2_b && j >= y1_b) {
				mask[(i * diameter1) + j] = 0;
			}
		}
	}

	return mask;
}

vector<long int> size_snr(long int total_frames, int lz) {

	vector<long int> snr_size;

	for (int i = 0; i < lz; i++)
	{
		snr_size.push_back((i + 1) * 10);
	}

	return snr_size;
}

vector<vector<double>> sign_noise(double* data,long int total_frames, int outter_ring_diameter, vector<long int> z, int lz, int* circle, int* ring) {

	vector<vector<double>> curves;
	vector<double> curve_by_z;

	double* aux_data;
	long int p;

	for (long int i = 0; i < lz; i++)
	{
		p = (long int)(total_frames / z[i]);
		int dims = outter_ring_diameter * outter_ring_diameter;
		aux_data = (double*)calloc(dims * p, sizeof(double));

		for (long int j = 0; j < p; j++)
		{
			for (long int k = 0; k < z[i]; k++)
			{
				for (long int l = 0; l < dims; l++)
				{
					aux_data[(j * dims) + l] += data[(j * z[i] * dims) + ((k * dims) + l)];
				}
			}

			curve_by_z.push_back(get_curve(aux_data, j, outter_ring_diameter, circle, ring));

		}

		supr_outfl(curve_by_z);
		curves.push_back(curve_by_z);
		curve_by_z.clear();

	}

	return curves;
}

double sign_curve(vector<double> data, long int total_frames) {

	double signal;
	double prom = 0;
	double std = 0;

	for (int i = 0; i < total_frames; i++)
	{
		prom += data[i];
	}

	prom = prom / total_frames;

	for (int i = 0; i < total_frames; i++)
	{
		std += pow(data[i] - prom, 2);
	}

	std = pow((std / (total_frames - 1)), 0.5);
	signal = prom / (2.5 * std);

	return signal;
}

void supr_outfl(vector<double> &curve_by_frame) {

	double p = 0.0;
	double s = 0.0;

	float a, b;

	long int frames = curve_by_frame.size();

	for (int i = 0; i < frames; i++)
	{
		if (isnan(curve_by_frame[i])) {
			continue;
		}

		p += curve_by_frame[i];
	}

	p = p / frames;

	for (int i = 0; i < frames; i++)
	{
		if (isnan(curve_by_frame[i])) {
			s += pow(p, 2);
		}
		else {
			s += pow((curve_by_frame[i] - p), 2);
		}
	}

	s = sqrt(s / (frames - 1));
	a = p - (3.5 * s);
	b = p + (3.5 * s);

	for (int i = 0; i < frames; i++)
	{
		if (isnan(curve_by_frame[i]) || curve_by_frame[i] == 0) {
			curve_by_frame[i] = p;
		}
	}

	for (int i = 0; i < frames; i++)
	{
		if (curve_by_frame[i] < a || curve_by_frame[i] > b) {

			curve_by_frame[i] = 0.0;

			if (i > 4) {

				for (int j = 1; j < 6; j++)
				{
					curve_by_frame[i] += curve_by_frame[i - j];
				}

			}
			else {

				for (int j = 1; j < 6; j++)
				{
					curve_by_frame[i] += curve_by_frame[i + j];
				}

			}

			curve_by_frame[i] = curve_by_frame[i] / 5;

		}
	}

}