//#import "i3_great.h"
#include "i3_great.h"
#include "i3_global.h"
#include "i3_model_tools.h"
#include "i3_psf.h"
#include "tools/i3_options.h"


#ifdef I3_USE_DOUBLE
#define CAT_ROW_FORMAT "%lf %lf %lf %lf %d %d"
#else
#define CAT_ROW_FORMAT "%f %f %f %f %d %d"
#endif

//#define GREAT08_STAR_SET_1_G1 -0.019
//#define GREAT08_STAR_SET_1_G2 -0.007
// I had to flip the sign to match the GREAT08 simulations code output
// I haven't checked more than just set 1 yet. SLB


i3_moffat_psf * i3_great10_psf_catalog_read(char * filename, int *n)
{
  /* 
   * read in catalogue file 
   */

  FILE *fl;
  char line[i3_options_max_line_length];
  int counter;

  *n=GREAT10_IMAGES_PER_FILE;
  i3_moffat_psf *output = malloc(sizeof(i3_moffat_psf)*(*n));

  fl = fopen(filename,"r");	     
  if (fl == NULL) {
     printf("\nim3shape : Failed to open file: %s", filename);
     exit(-1);
  }
   /* get each line */
  counter =0;
  /* loop forever */
  for(;;){
    /* if first element of a line = # */
     fgets(line, i3_options_max_line_length, fl);
     if (line[0] == '#'){
       counter += 1;
       // printf("%i\n",counter);
     } else {
      /* else break */
     break;
     }
  }
  int i=0;
	
  sscanf(line,CAT_ROW_FORMAT,&((output+i)->beta),&((output+i)->fwhm),&((output+i)->e1),&((output+i)->e2),&((output+i)->x),&((output+i)->y));
  for (i=1;i<*n;i++){
    fgets(line, i3_options_max_line_length, fl);
  	sscanf(line,CAT_ROW_FORMAT,&((output+i)->beta),&((output+i)->fwhm),&((output+i)->e1),&((output+i)->e2),&((output+i)->x),&((output+i)->y));
  }
  return output;
};

i3_fourier * i3_build_great10_kernel(i3_moffat_psf * psf,int npix, int upsampling, int padding, int psf_type){
	return i3_fourier_conv_kernel_moffat(npix, upsampling, padding, psf->beta, psf->fwhm, psf->e1, psf->e2, 1e30, psf_type);
};




i3_fourier * i3_build_great08_kernel(int set, int nx, int ny, int upsampling)
{
	int n1 = nx + 1;
	int N1 = n1 * upsampling;
	i3_flt x0 = nx/2.0;
	i3_flt y0 = ny/2.0;
	i3_flt X0 = x0 * upsampling;
	i3_flt Y0 = y0 * upsampling;
	int padding = N1;
	i3_image * psf = i3_great08_psf(set, N1, X0, Y0, upsampling);
	i3_image * padded_psf = i3_image_padded_corner(psf,N1);

	i3_fourier * psf_fourier = i3_image_fourier(padded_psf);
	i3_fourier * downsampler_fourier = i3_fourier_downsampling_psf(n1, n1, upsampling, n1);
	i3_fourier * kernel = i3_fourier_create(N1+padding,N1+padding);
	i3_multiply_fourier_fourier_into(psf_fourier, downsampler_fourier, kernel);
	i3_fourier_destroy(psf_fourier);
	i3_fourier_destroy(downsampler_fourier);
	return kernel;
	
}

i3_moffat_psf * i3_great08_psf_form(int start_set, int stamp_size)
{
	i3_moffat_psf * psf = malloc(sizeof(i3_moffat_psf));
	psf->beta = GREAT08_MOFFAT_BETA;
	psf->fwhm = GREAT08_PSF_FWHM_PIXELS;
	psf->x = stamp_size/2.0;
	psf->y = stamp_size/2.0;
	
	switch(start_set){
		case 1:
			psf->e1 = GREAT08_STAR_SET_1_G1;
			psf->e2 = GREAT08_STAR_SET_1_G2;
			break;
		case 2:
			psf->e1 = GREAT08_STAR_SET_2_G1;
			psf->e2 = GREAT08_STAR_SET_2_G2;	
			break;
		case 3:
			psf->e1 = GREAT08_STAR_SET_3_G1;
			psf->e2 = GREAT08_STAR_SET_3_G2;	
			break;
		default:
			I3_FATAL("Unknown Great08 star psf set",1);
			break;
	}
	return psf;
}


i3_image * i3_great08_psf(int start_set, int npix, i3_flt x0, i3_flt y0, int upsampling){
	
	
	if (start_set<1 || start_set>3) I3_FATAL("Requested Great08 PSF set <1 or >3.  Only three available.",1);
		
	i3_flt e1 = 0 , e2 = 0;
	if (start_set==1){
		e1 = GREAT08_STAR_SET_1_G1;
		e2 = GREAT08_STAR_SET_1_G2;
	}
	else if (start_set==2){
		e1 = GREAT08_STAR_SET_2_G1;
		e2 = GREAT08_STAR_SET_2_G2;	
	}
	else if (start_set==3){
		e1 = GREAT08_STAR_SET_3_G1;
		e2 = GREAT08_STAR_SET_3_G2;		
	}
	else{
		I3_FATAL("Unknown Great08 star psf set",1);
	}
	
	//i3_flt fwhm_pixels = GREAT08_PSF_FWHM_PIXELS * (1.0*npix)/GREAT08_ORIGINAL_STAMP_SIZE; //Rescale the fwhm to fit the number pf pixels
	i3_flt fwhm_pixels = GREAT08_PSF_FWHM_PIXELS *upsampling; 
	i3_flt moffat_beta = GREAT08_MOFFAT_BETA;
	i3_flt moffat_radius = fwhm_pixels*MOFFAT_RADIUS_FOR_UNIT_FWHM;
	//i3_flt truncation_radius = fwhm_pixels;
	// testing with truncation switched off for now! 20110617SLB
	i3_flt truncation_radius = fwhm_pixels*2.;
	
	i3_image * moffat = i3_image_create(npix,npix);
	i3_image_zero(moffat);
	
	i3_image_add_truncated_moffat(moffat,x0,y0, moffat_radius, e1, e2, moffat_beta, truncation_radius);

	i3_flt sum = i3_image_sum(moffat);
	i3_image_scale(moffat,1./sum);

	return moffat;
}