#include "i3_mcmc.h"
#include "i3_load_data.h"
#include "i3_math.h"
#include "i3_model.h"
#include "i3_minimizer.h"
#include "i3_fisher.h"
#include "i3_mcmc.h"
#include "i3_image_fits.h"
#ifdef MPI
#include "mpi.h"
#endif


#include "i3_analyze.h"
#include "i3_model_tools.h"
#include "i3_data_set.h"
#include "i3_great.h"
#include "i3_psf.h"
#include "i3_image_fits.h"
#include "i3_utils.h"
#include "i3_options.h"
#include <time.h>



int main(){
	
	//Initialize FFTW
	i3_fftw_load_wisdom();
	atexit(i3_fftw_save_wisdom);

	int N=48;
	
	i3_flt e1 = 0.1;
	i3_flt e2 = -0.4;
	i3_flt x0 = N/2.;
	i3_flt y0 = N/2.;
        i3_flt SN_scale = 1000.;
	i3_flt bulge_amplitude = 1.0;
	i3_flt disc_amplitude = 0.5;
	i3_flt bulge_radius = 3.0;
	i3_flt disc_radius = 6.0;

	i3_flt noise_sigma = 1.;
	
	i3_image * image = i3_image_create(N,N);
	i3_image * weight = i3_image_create(N,N); 
        i3_image_fill(weight, 1. / noise_sigma / noise_sigma);

	// Build options
	i3_options * options = i3_options_default();
	options->stamp_size = N;
	options->minimizer_tolerance = 1.e-5;
	options->minimizer_loops = 1;
	options->milestone_radius_ratio_max = 10.;
	
	//Create basic information and space for the chosen model (e.g. milestone)
	i3_model * model = i3_model_create("milestone", options);


	//Set up the PSF
	i3_moffat_psf psf_form;
	psf_form.beta = 3.;
	psf_form.e1 = 0.;
	psf_form.e2 = 0.;
	psf_form.fwhm = 4.;
        
	i3_data_set * dataset = i3_build_dataset(options, image, weight, &psf_form, GREAT10_PSF_TYPE_MOFFAT);

	// Build the model parameters
	i3_milestone_parameter_set * params_gal = i3_model_option_starts(model, options);
	params_gal->radius = disc_radius;
	params_gal->radius_ratio = bulge_radius / disc_radius;
	params_gal->flux_ratio = bulge_amplitude / (bulge_amplitude + disc_amplitude);
	params_gal->bulge_index = 4.0;
	params_gal->disc_index = 1.0;
        params_gal->x0 = x0;
	params_gal->y0 = y0;
	params_gal->e1 = e1;
       	params_gal->e2 = e2;

	// Build the image
	i3_milestone_model_image(params_gal,dataset,image);
	i3_image_scale(image, SN_scale / i3_image_sum(image));

	// Add noise
	i3_init_rng();
	i3_image_add_white_noise(image, noise_sigma);
	i3_image_save_fits(image, "optest_tmp.fits");

	// Analyse the test data
	// i3_analyze_dataset_method(dataset, model, options, stdout, 1, i3_minimizer_method_praxis);
	// i3_analyze_dataset_method(dataset, model, options, stdout, 2, i3_minimizer_method_levmar);

	i3_image * image_dvc = i3_image_create(n_pix*n_sub,n_pix*n_sub);	i3_image_zero(image_dvc);
	i3_image * image_exp = i3_image_create(n_pix*n_sub,n_pix*n_sub);	i3_image_zero(image_exp);

	i3_sersics_components( params_gal,dataset,image_dvc,image_exp );  

	i3_flt * logL_xy0 = i3_sersics_likelihood_in_xy0(image_dvc,image_exp,params_gal,dataset );

	



}