/*----------------------------------------------------------------------------
   
   File name 	:	encircl.c
   Author 		:	Christian Drouet d'Aubigny
   Created on	:	21 Oct 1996
   Description	:	compute the radius corresponding to a percentage of
   					encircled energy

 ---------------------------------------------------------------------------*/

/*
	$Id: encircl.c,v 1.24 2002/01/29 12:49:11 ndevilla Exp $
	$Author: ndevilla $
	$Date: 2002/01/29 12:49:11 $
	$Revision: 1.24 $
 */

/*---------------------------------------------------------------------------
  								Includes
 ---------------------------------------------------------------------------*/

#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include <math.h>

#include "eclipse.h"


/*---------------------------------------------------------------------------
   							Private declarations	
 ---------------------------------------------------------------------------*/

static void
usage(char *pname) ;

double
get_radius_on_image(
    image_t    *   image_in,
    int             x_expect,
    int             y_expect,
    int             half_size,
    double          Plate_Scale,
    double          TotalRadius,
    int             percent
);


static char prog_desc[] =
"arcsecond radius for given percentage encircled energy" ;

/*---------------------------------------------------------------------------
  							Main code	
 ---------------------------------------------------------------------------*/

int main(int argc, char *argv[])
{
	cube_t		*	cube_in ;
	char			inname[FILENAMESZ+1] ;
	int				c ;
  
	/* Parameters to fwhm computation   */
	int				np ;
	char	*		op_scal ;
	double			platescale ;
	int				x_expect,
					y_expect,
					half_size ;
	double			total_radius ;
	int				percent ;
	double			radius ;
	char			format[200] ;
	int				tot,
  					dec ;

	/* Initialize variables */
	platescale = -1 ;
	total_radius = 1.4 ;
	half_size = 20 ;
  
	if (argc<2)
		usage(argv[0]) ;
	  
	/* Command line parsing by getopt() */
	while ((c = getopt(argc, argv, "Lh:p:r:")) != EOF)
    switch(c) { 

		/* Standard option: display license (not documented in usage)   */
		case 'L':
		eclipse_display_license() ;
		return 0 ;
	
		/* Special options : plate scale */
		case 'p':
		platescale = (double) atof (optarg) ;
		break ; 
	
		/* Special option : 100 % radius size (arc-sec) */ 
		case 'r':
		total_radius = (double) atof (optarg) ;
		break ;

		/* Special option: modify half size defaulted to 20 pix */
		case 'h':
		half_size = (int) atoi (optarg) ;
		break ;

		default:
		usage(argv[0]) ;
		break ;
	}

	/* Initialize eclipse environment */
	eclipse_init();
 
	/* Get arguments    */
	if ((argc-optind)!=4) {
		e_error("invalid number of arguments") ;
		return -1 ;
	}
  
	strncpy(inname, argv[optind++], FILENAMESZ) ;
	x_expect = (int) atoi (argv[optind++]) ;
	y_expect = (int) atoi (argv[optind++]) ;
	percent  = (int) atoi (argv[optind]);

  
	/* load the cube and check for errors */
	if (platescale<0){
		op_scal = qfits_query_hdr(inname, "OP_SCAL");
		if (op_scal==NULL) {
			e_error("cannot find OP_SCAL in header");
			return -1 ;
		}
	}
  
	cube_in = cube_load(inname) ;
	if (cube_in == NULL){
		e_error("in loading cube [%s]: aborting", inname) ;
		return -1 ;
	}
  
	for (np=0 ; np<cube_in->np ; np++) { 
		/* working on user coordinate format 1,2,3...n and then calling */
		/*  routines in 0,1,2...n-1 eclipse format */    
		radius = get_radius_on_image(cube_in->plane[np],
									 x_expect,
									 y_expect,
									 half_size,
									 platescale,
									 total_radius,
									 percent) ;

		if (-1 == radius){
			e_error("cannot compute radius");
			return -1;
		}
		if (platescale < 1){
			dec = 2 - (int) (log10(platescale)) ;
			tot = 2 + dec ;
		} else if (log10(platescale) < 1) {
   			dec = 1  ;
   		} else {
   			dec = 0  ;
   		}
   		tot = 2 + (int) (log10(platescale)) ;

		sprintf(format,
				"radius for %d percent is: %%%d.%df arc-sec\n",
				percent,tot,dec);
		printf(format,radius);
	}
 
	cube_del(cube_in) ;
  
	if (debug_active())
		xmemory_status() ;
	return 0 ;
}

/*
 * This function only gives the usage for the program
 */

static void
usage(char *pname)
{
	hello_world(pname, prog_desc) ;
	printf(
"use : %s [options] incube x_expect y_expect percent \n"
"options are :\n"
"\t[-h halfsize] sets confidence window half size\n"
"\t[-r radius] to override the 1.4 arc-sec hundred percent energy radius \n"
"\t[-p plate_scale] in arc-sec per pixel (overrides FITS header)\n"
"\n\n", pname) ;
    exit(0) ;
}

 
/*---------------------------------------------------------------------------
   Function :   get_radius_on_image()
   In       :   1 image
                x,y position of expected peak
                half size of zone containing peak
                plate scale in arc second per pixel
                TotalRadius (100% Energy radius) in arc seconds
                percentage of energy for which the radius is wanted
   Out      :   double radius
   Job      :   computes the radius which corresponds to a percentage
                of the encircled energy
   Notice   :   returns -1 if something goes wrong
 ---------------------------------------------------------------------------*/
 
double
get_radius_on_image(
    image_t    *   image_in,
    int             x_expect,
    int             y_expect,
    int             half_size,
    double          Plate_Scale,
    double          TotalRadius,
    int             percent
)
{
  double         Radius ;
  double         TotalEnergy ;
  double         TotalEnergyPercentValue ;
  int           TotalRadiusInPixel ;
  double        *EnergyArray;
  image_t      *sub_image ;
  int           x_min, x_max, y_min, y_max ;
  image_stats    *sub_ima_stats ;
  int           i ;
  int           x1, x2 ;
  double        y1, y2 ;
  int			PassedValue ;
  int           x, y;
 
  /* initialization */
 
  PassedValue = 0 ;
 
  /* Error handling: test entries */
 
  if (image_in == NULL) return -1.0 ;
 
  /* Check that the peak position estimate is in the frame */
  if ((x_expect < 1) ||
      (x_expect > image_in->lx) ||
      (y_expect < 1) ||
      (y_expect > image_in->ly)) {
    e_error("peak position estimate out of frame: [%d %d]",
			x_expect, y_expect);
	return -1.0 ;
  }
 
  /* defines the expectation window */
 
  x_min = x_expect - half_size ;
  y_min = y_expect - half_size ;
  x_max = x_expect + half_size ;
  y_max = y_expect + half_size ;
 
  if(x_min < 1) x_min = 1 ;
  if(y_min < 1) y_min = 1 ;
  if(x_max > image_in->lx) x_max = image_in->lx ;
  if(y_max > image_in->ly) y_max = image_in->ly ;
 
  /* extracts subframe */
  sub_image = image_getvig(image_in,x_min,y_min,x_max,y_max) ;
 
  /* do some stats on that subimage */
  sub_ima_stats = image_getstats(sub_image) ;
  image_del(sub_image);
 
  /* get the position of the peak and check for positivity */
 if( 0 > sub_ima_stats->min_pix){
    e_error("some pixels have negative value (minimum is %g)",
            sub_ima_stats->min_pix);
    return -1.0;
  }
 
  x = sub_ima_stats->max_x + x_min ;
  y = sub_ima_stats->max_y + y_min ;
  free(sub_ima_stats) ;
 
  /* computes the 100 % radius in pixels */
 
  if ( Plate_Scale < 1e-7){
    e_error("plate scale too small: %g arcsec/pix", Plate_Scale);
    return -1.0;
  }
 
  TotalRadiusInPixel = (int)(TotalRadius/Plate_Scale)  ;
 
  if ( 0 == TotalRadiusInPixel ){
    e_error("total radius in pixel is 0...");
    return -1.0;
  }
 
  /* extract a new re-centered re scaled sub image */
 
  x_min = x - TotalRadiusInPixel ;
  y_min = y - TotalRadiusInPixel ;
  x_max = x + TotalRadiusInPixel ;
  y_max = y + TotalRadiusInPixel ;
 
  sub_image = image_getvig(image_in,x_min,y_min,x_max,y_max) ;
  if (sub_image == NULL) {
    e_error("100%% energy radius window out of image\n");
    return -1.0;
  }
 
  sub_ima_stats = image_getstats(sub_image);
  x = sub_ima_stats->max_x ;
  y = sub_ima_stats->max_y ;
  free(sub_ima_stats);
 
  /* computes the 100 % energy  */
  TotalEnergy = image_get_radenergy(sub_image, x, y, TotalRadiusInPixel);
 
  /* computes the value of x % of the total encircled energy */
  TotalEnergyPercentValue = TotalEnergy * percent / 100 ;
 
  /* reserve a double array to hold the encircled  energy values */
  EnergyArray = malloc(TotalRadiusInPixel * sizeof(double));
 
  /* fills up the array with encicled energies */
 
  for(i=0 ; i < TotalRadiusInPixel ; i++){
    EnergyArray[i] =  image_get_radenergy(sub_image, x, y, i+1) ;
  }
  image_del(sub_image) ;
 
  /* finds where the encicled energy equals TotalEnergyPercentValue */
  if ( EnergyArray[0] >  TotalEnergyPercentValue ){
    e_warning("extrapolating. can be really imprecise");
    if (imstat_x_for_y_between_2_points(1, EnergyArray[0], 2, EnergyArray[1],
                                 TotalEnergyPercentValue, &Radius) != 0) {
        e_error("in extrapolation: aborting") ;
        return -1.0 ;
    }
    if( 0 >= Radius){
      e_error("extrapolation returned negative value, percentage too small");
      return -1.0 ;
    }
    return Radius * Plate_Scale;
  }
  else{
    x1 = 1 ;
    x2 = 2 ;
    y1 = EnergyArray[0] ;
    y2 = EnergyArray[1] ;
    for(i=3 ; i < TotalRadiusInPixel ; i++){
      if (!PassedValue){
        if ( y2 > TotalEnergyPercentValue ){
          PassedValue = 1 ;
        }
        else {
          x1++ ;
          x2++ ;
          y1 = y2 ;
          y2 = EnergyArray[i];
        }
      }
    }
  }
 
  free(EnergyArray);
 
  /* given two points around TotalEnergyPercentValue get subpixel precision */
  if (imstat_x_for_y_between_2_points(x1, y1, x2, y2,
                               TotalEnergyPercentValue, &Radius) != 0) {
        e_error("in interpolation: aborting") ;
        return -1.0 ;
    }
 
  return Radius * Plate_Scale ;
}