#include "i3_image.h"
#include "i3_load_data.h"
#include "stdio.h"
#include "stdlib.h"
#include "i3_model.h"
#include "i3_model_tools.h"
#include "i3_psf.h"
#include <time.h>





int main(int argc, char *argv[]){

// init

	i3_fftw_load_wisdom();
	atexit(i3_fftw_save_wisdom);


// stamp parameters

	int N_pix = 39;
	int N_obs = 39; 
	int N_sub = 7;
	int N = N_pix*N_sub;
	int N_cut = N_pix - N_obs;
	

// psf parameters

	i3_flt psf_sersic_n = 0.5f; 
	i3_flt psf_A=1.0f;
	i3_flt psf_x0=N/2; 
	i3_flt psf_y0=N/2; 
	i3_flt psf_fwhm = 2.85f*N_sub;
	i3_flt psf_ab = i3_fwhm_to_ab(psf_fwhm,psf_sersic_n);
	i3_flt psf_e1 = 0.0f;
	i3_flt psf_e2 = 0.0f;
	i3_image * image_psf = i3_image_create(N,N);	i3_image_zero(image_psf);
	i3_add_real_space_sersic(image_psf, psf_ab,psf_e1,psf_e2, psf_A, psf_x0, psf_y0, psf_sersic_n );
	i3_image_save_text( image_psf,"image_psf.txt" );
	
	i3_image * image_psf_c = i3_image_centered( image_psf );
	i3_image_save_text( image_psf_c,"image_psf_c.txt" ); 

// get galaxy

	i3_flt gal_sersic_n = 1.f;
 	i3_flt gal_A=1.0f;
 	i3_flt gal_x0= N/2; 
	i3_flt gal_y0= N/2; 
 	i3_flt gal_fwhm = 2.0f*N_sub;
 	i3_flt gal_ab = i3_fwhm_to_ab(gal_fwhm,gal_sersic_n);
 	i3_flt gal_e1 = 0.5f;
 	i3_flt gal_e2 = 0.0f;
	i3_image * image_gal = i3_image_create(N,N);	i3_image_zero(image_gal);
	i3_add_real_space_sersic(image_gal , gal_ab, gal_e1 , gal_e2 , gal_A , gal_x0 , gal_y0 , gal_sersic_n );
	i3_image_save_text( image_gal,"image_gal.txt" );

	i3_image * image_gal_c = i3_image_centered( image_gal );
	i3_image_save_text( image_gal_c,"image_gal_c.txt" );	

// get sinc

	i3_image * image_pix = i3_image_get_top_hat( N_pix,N_pix,N_sub );
	i3_image_save_text( image_pix,"image_pix.txt" );

	i3_image * image_pix_c = i3_image_centered( image_pix );
	i3_image_save_text( image_pix_c,"image_pix_c.txt" );

// select your config here !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
// suffix _c means that it has been 180 rotated using i3_image_centered
	
	//i3_image * image_psf_used = image_psf;
	i3_image * image_psf_used = image_psf_c;

	//i3_image * image_psf_used = image_psf;
	i3_image * image_gal_used = image_gal;

	//i3_image * image_pix_used = image_psf;
	i3_image * image_pix_used = image_pix_c;


// get the FTs
	
		
	i3_fourier * fourier_psf = i3_image_fourier( image_psf_used );
	i3_fourier * fourier_pix = i3_image_fourier( image_pix_used );

// combine kernels
	
	i3_fourier * fourier_ker = i3_fourier_create(N,N);
	i3_multiply_fourier_fourier_into(fourier_pix,fourier_psf,fourier_ker);

// convolve
	
	i3_image * image_big = image_gal_used;
	
	i3_image_convolve_fourier( image_big, fourier_ker );

	i3_image * image_big_c = i3_image_centered( image_big ); 

	i3_image_save_text( image_big,"image_big.txt" );
	i3_image_save_text( image_big_c,"image_big_c.txt" );
		
// cut out the middle

	i3_image * image_obs = i3_image_copy_part( image_big , (N_cut/2)*N_sub, (N_cut/2)*N_sub, N_obs*N_sub, N_obs*N_sub );
	i3_image_save_text( image_obs,"image_obs.txt" );
		

	return 0;


/*
Results of checking symmetry of the images. Assumptions:
- psf is gaussian with zero ellipticity
- galaxy is sersic 1 with e1 > 0 and e2 = 0, symmetric on the x axis
- the image is subpixel integrated using a top hat
- both galaxy and psf are centered in the middle subpixel
This setup should return a convolved object that is symmetric around x axis.   
image_ - normal real space image
image_c - image with 180* phase shift (used function i3_centered)

gal * pix = phase shifted
gal * pix_c = OK
gal_c * pix = OK
gal * psf * pix = the center is shifted one pixel in both dimensions
gal_c * psf * pix = phase shifted 180
gal * psf_c * pix = phase shifted 180
gal * psf * pix_c = phase shifted 180
gal_c * psf_c * pix = OK
all other combinations containing 2x _c = OK
gal_c * psf_c * pix_c = phase shifted 180

see http://www.varioustopics.com/image-processing/425948-problems-with-phase-in-fftw.html
for reasonable explanation

*/

}