;+
;
; NAME:
;       PIX_SAMPLE_1D
;
;
; PURPOSE:
;   	Computes the 1D pixel sampled vector of a y(x) profile given on a regular or irregular grid.
;
;
; DESCRIPTION:
;		This function does conceptually the same as the IDL function INTERPOL, i.e. it computes
;		a regularly gridded representation of the given function y=f(x) where the x vector
;		can also be non uniform.
;
;		The difference is that INTERPOL just interpolates the function onto a 1D
;		regularly spaced grid, while PIX_SAMPLE_1D gives the pixel sampled array of the function, i.e.
;		each point in the final array contains the integral (normalized to the pixel size) of the
;		input function over the sampling interval.
;
;		This makes difference when y(x) contains peak features narrower than the choosed
;		sample interval.
;		In this case PIX_SAMPLE_1D will conserve the energy in each sampling interval whatever the
;		choosed sampling resolution is.
;
;		Very useful when Fourier transforms are required for irregularly spacing of x: the transform of
;		the 1-D vector could depend on the choosed sampling unless pixel sampling is used.
;
;		NOTE: if the y(x) input function represent the radial profile of a polar surface and this function is
;		called to compute the final 2D polar surface, use PIXEL_SAMPLED_POLAR_SURFACE instead of this function.
;
;
; CALLING SEQUENCE:
;
;		Result = PIX_SAMPLE_1D(x, y, SAMPLE=SAMPLE, _EXTRA=EXTRA)
;
; INPUTS:
;       x:		1d array of the independent variable; may be irregularly spaced.
;
;		y:		1d array of the dependent variable, of the same number of elements as r.
;
;
; KEYWORDS:
;
;		SAMPLE:		Sampling interval (i.e. 1D pixel size).
;					If the whole range of 'x' cannot be sampled with an integer number of
;					samples, the larger integer multiple of the 'sample' value will be taken as max(x).
;					The minimum value of x, min(x), will remain the same.
;
;		_EXTRA:		Any keyword of INTERPOL can be used in calling this function.
;
;
; OUTPUTS:
;
;       Result:	The function returns a 1D array with the 1D pixel sampled function.
;
;
; CALLS:
;
;		The function TSUM from W.B. Landsman for trapezoidal summation is called and is included
;		inside this function.
;
; MODIFICATION HISTORY:
;
;       September 2004, G. Li Causi, Rome Astronomical Observatory
;						licausi@mporzio.astro.it
;						http://www.mporzio.astro.it/~licausi/
;
; PLANNED IMPROVEMENTS:
;
;	- Add CENTRAL_POINT keyword to specify an origin to keep fixed.
;-


FUNCTION Pix_Sample_1D, x, y, SAMPLE=SAMPLE, _EXTRA=EXTRA


;***********
;Input Check
;***********
IF SIZE(x, /N_DIMENSIONS) NE 1 THEN MESSAGE, 'Independent variable x must be a 1-D array!'
IF SIZE(y, /N_DIMENSIONS) NE 1 THEN MESSAGE, 'Function value y must be a 1-D array!'
IF n_elements(y) NE n_elements(x) THEN MESSAGE, 'Function values y must have same number of elements as independent variable x!'



;**********************************************************
;Ricampionamento ed integrazione del vettore y(x) sui pixel
;**********************************************************

IF NOT KEYWORD_SET(SAMPLE) THEN BEGIN
	dx = x - shift(x, 1)
	dx = dx[1:*]
	sample = min(abs(dx))
ENDIF ELSE sample = ABS(sample)

num_x_samples = long( (max(x) - min(x)) / sample)			;number of pixels

xs = dindgen(num_x_samples + 1) * sample + min(x)			;x of pixel left borders

y_res_interp = INTERPOL(y, x, xs, _EXTRA=EXTRA)				;interpolation of data at pixel left borders

y_res = dblarr(num_x_samples)
FOR jj = 0L, num_x_samples-2 DO BEGIN
	aindex = min(where(x GT xs[jj], count))
	IF count EQ 0 THEN aindex = 0
	xa = xs[jj]
	ya = y_res_interp[jj]
	bindex = max(where(x LT xs[jj+1], count))
	IF count EQ 0 THEN bindex = 0
	xb = xs[jj+1]
	yb = y_res_interp[jj+1]
	IF bindex LT aindex THEN BEGIN
		xvector = [xa, xb]
		yvector = [ya, yb]
	ENDIF ELSE BEGIN
		xvector = [xa, x[aindex:bindex], xb]
		yvector = [ya, y[aindex:bindex], yb]
	ENDELSE
	y_res[jj] = TSUM(xvector, yvector) / (xb - xa)
ENDFOR
xs = xs[0:num_x_samples-2]
y_res = y_res[0:num_x_samples-2]


RETURN, transpose([[xs], [y_res]])

END










FUNCTION TSUM,X,Y,IMIN,IMAX               ;Trapezoidal summation
;+
; NAME:
;       TSUM
; PURPOSE:
;       Trapezoidal summation of the area under a curve.
; EXPLANATION:
;       Adapted from the procedure INTEG in the IUE procedure library.
;
; CALLING SEQUENCE:
;       Result = TSUM(y)
;              or
;       Result = TSUM( x, y, [ imin, imax ] )
; INPUTS:
;       x = array containing monotonic independent variable.  If omitted, then
;               x is assumed to contain the index of the y variable.
;               x = lindgen( N_elements(y) ).
;       y = array containing dependent variable y = f(x)
;
; OPTIONAL INPUTS:
;       imin = scalar index of x array at which to begin the integration
;               If omitted, then summation starts at x[0].
;       imax = scalar index of x value at which to end the integration
;               If omitted then the integration ends at x[npts-1].
;
; OUTPUTS:
;       result = area under the curve y=f(x) between x[imin] and x[imax].
;
; EXAMPLE:
;       IDL> x = [0.0,0.1,0.14,0.3]
;       IDL> y = sin(x)
;       IDL> print,tsum(x,y)    ===>  0.0445843
;
;       In this example, the exact curve can be computed analytically as
;       1.0 - cos(0.3) = 0.0446635
; PROCEDURE:
;       The area is determined of individual trapezoids defined by x[i],
;       x[i+1], y[i] and y[i+1].
;
;       If the data is known to be at all smooth, then a more accurate
;       integration can be found by interpolation prior to the trapezoidal
;       sums, for example, by the standard IDL User Library int_tabulated.pro.
; MODIFICATION HISTORY:
;       Written, W.B. Landsman, STI Corp. May 1986
;       Modified so X is not altered in a one parameter call Jan 1990
;       Converted to IDL V5.0   W. Landsman   September 1997
;       Allow non-integer values of imin and imax  W. Landsman April 2001
;-
; Set default parameters
 On_error,2
 npar = N_params()

 if npar EQ 1 then begin
    npts = N_elements(x)
    yy = x
    xx = lindgen(npts)
    ilo = 0
    ihi = npts-1
 endif else begin

   if ( npar LT 3 ) then imin = 0
   npts = min( [N_elements(x), N_elements(y)] )
   if ( npar LT 4 ) then imax = npts-1
   ilo = long(imin)
   ihi = long(imax)
   xx = x[ilo:ihi]
   yy = y[ilo:ihi]
   npts = ihi - ilo + 1
 endelse
;
; Compute areas of trapezoids and sum result
;
  xdif = xx[1:*] - xx
  yavg =  ( yy[0:npts-2] + yy[1:npts-1] ) / 2.
  sum = total( xdif*yavg )

; Now account for edge effects if IMIN or IMAX parameter are not integers

  hi = imax - ihi
  lo = imin - ilo
  if (ihi LT imax) then sum = sum + (x[ihi+1]-x[ihi])*hi* $
              (y[ihi] + (hi/2.) *(y[ihi+1] - y[ihi]) )
  if (ilo LT imin) then sum = sum - (x[ilo+1]-x[ilo])*lo* $
              (y[ilo] + (lo/2.) *(y[ilo+1] - y[ilo]) )
  return, sum

  END