; Copyright (c) 1998-2011 A.P. Hitchcock  All rights reserved
;+
;NAME:
;	AX_MESH
;
;LAST CHANGED: ----------------------------------- 09-Feb-11
;
;PURPOSE:
;	This function meshes two 1-d (spectra) or 2-d (images) AXIS data buffers
; to a common (x) or (x,y) scale
;
;CATEGORY:
;	AXIS: image and spectral processing
;
;CALLING SEQUENCE:
;	Result = AX_MESH( tmp1, tmp2, nd, first=first )
;
;CALLED FROM AXIS:
;	(ax_image_add, image_multiply, spectral_add)

;	TMP1 - first spectrum or image
;	TMP2 - second spectrum or image
;	ND   - 1 = spectra; 2 = images
;   FIRST - if set, use mesh of the first file (spectra only)
;
;KEYWORDS: none
;
; OUTPUT
;	Result = 1 signifies successful completion
;	TMP1 and TMP2 are returned converted to common scales
;
;COMMON BLOCKS:
;	@AXIS_COM	standard set of common blocks
;
;MODIFICATION HISTORY:
; ( 3-Dec-98 aph) first generated
; (11-Jan-99 aph) removed print lines
; (04-jul-99 aph) correct for images
; (28-dec-99 aph) correct for images (after extracting ax_image_add from axis_c)
; (31-dec-99 aph) spectral meshing ; AXIS standard documentation
; (09-jun-04 aph) include energy in output structures
; (21-apr-05 aph) add option to force mech to the first file (spectra only)
; (03-feb-11 aph) output axis format image files have a defined energy (default, the original one)
; (09-feb-11 aph) stop E-define applied to spectra
;-

FUNCTION ax_mesh,tmp1,tmp2,nd, first=first
@axis_com
on_error,2

if nd EQ 2 then begin
; interpolate each 2d file to same (x,y) range and mesh to finest of the pixel sizes
	xl=fltarr(2) & xh = xl  & yl = xl  & yh = xl
	xl(0) = min(tmp1.x,max=t) & xh(0)=t    &  xl(1) = min(tmp2.x,max=t) & xh(1)=t
	yl(0) = min(tmp1.y,max=t) & yh(0)=t    &  yl(1) = min(tmp2.y,max=t) & yh(1)=t
	xmin = max(xl) & xmax = min(xh)
	ymin = max(yl) & ymax = min(yh)
	if xmax LE xmin OR ymax LE ymin then goto, bailout
	xl(0) = min(tmp1.x(1:n_elements(tmp1.x)-1) - tmp1.x(0:n_elements(tmp1.x)-2))
	xl(1) = min(tmp2.x(1:n_elements(tmp2.x)-1) - tmp2.x(0:n_elements(tmp2.x)-2))
	yl(0) = min(xl)            ; smallest x-step
	xl(0) = min(tmp1.y(1:n_elements(tmp1.y)-1) - tmp1.y(0:n_elements(tmp1.y)-2))
	xl(1) = min(tmp2.y(1:n_elements(tmp2.y)-1) - tmp2.y(0:n_elements(tmp2.y)-2))
	yl(1) = min(xl)             ; smallest y-step
	step = abs(min(yl))         ; force to square pixels NB ABS accounts for reverse order !
;		print, 'Ax_Mesh: Common (x,y) range:'
;		print,' X ', xmin, xmax
;		print, 'Y ', ymin, ymax
;		print, 'MINIMUM x-step, y-step, step', yl(0), yl(1), step
	nx = (1 + (xmax - xmin) / step) & x = findgen(nx)  & x = xmin + step*x
	ny = (1 + (ymax - ymin) / step) & y = findgen(ny)  & y = ymin + step*y
	if n_elements(tmp1.e) GT 0 then e1 = tmp1.e else e1 = 0.  ; force output to have E defined
	if n_elements(tmp2.e) GT 0 then e2 = tmp2.e else e2 = 0.

; interpolate tmp1.d data to common mesh
	Ilo = DIndex(xmin,ymin,tmp1)
	Ihi = DIndex(xmax,ymax,tmp1)
	dn = tmp1.d(Ilo(0):Ihi(0),Ilo(1):Ihi(1))
	dn = congrid(dn,nx,ny,/interp)
	tmp1 = create_struct('t','2d','x',x,'y',y,'d',dn, 'e', tmp1.e, $
			'xl',tmp1.xl,'yl',tmp1.yl,'dl',tmp1.dl)
; interpolate tmp2.d data to common mesh
	Ilo = DIndex(xmin,ymin,tmp2)
	Ihi = DIndex(xmax,ymax,tmp2)
	dd = tmp2.d(Ilo(0):Ihi(0),Ilo(1):Ihi(1))
	dd = congrid(dd,nx,ny,/interp)
	tmp2 = create_struct('t','2d','x',x,'y',y,'d',dd,'e', tmp2.e,  $
			'xl',tmp2.xl,'yl',tmp2.yl,'dl',tmp2.dl)
	return, 1	; returns parameter = 1 to signify successful interpolation
endif

if nd EQ 1 then begin
	xmin = fltarr(2) & xmax = fltarr(2)               ; identify common range
	xmin(0) = min(tmp1.x) & xmin(1) = min(tmp2.x)
	xmax(0) = max(tmp1.x) & xmax(1) = max(tmp2.x)
	xmin = max(xmin) & xmax = min(xmax)
;	print, 'xmin,xmax ', xmin, xmax
	xval1 = where(tmp1.x GE xmin AND tmp1.x LE xmax, nt1)   ; check there is in-range data
	xval2 = where(tmp2.x GE xmin AND tmp2.x LE xmax, nt2)   ; in both buffers
;	n=where(tmp1.x,cnt1)  & n=where(tmp2.x,cnt2)
;	print, 'size-1, size-2, in-range sizes:' , cnt1, cnt2, nt1, nt2
;	print, ' steps', (xmax-xmin)/nt1, (xmax-xmin)/nt2
	if nt1 NE 0 AND nt2 NE 0 then BEGIN
		x1 = tmp1.x(where(tmp1.x GE xmin AND tmp1.x LE xmax))   ; trim arrays to common range
		d1 = tmp1.d(where(tmp1.x GE xmin AND tmp1.x LE xmax))   ; NB must use intermediate arrays!!
		x2 = tmp2.x(where(tmp2.x GE xmin AND tmp2.x LE xmax))
		d2 = tmp2.d(where(tmp2.x GE xmin AND tmp2.x LE xmax))
;		n=where(x1,nt1)  & n=where(x2,nt2)
		if keyword_set(first) then begin
				y1 = d1  & x = x1
				y2 = Interpol(d2, x2, x1)        ; mesh to first data set
		endif else begin
			if nt1 GE nt2 THEN BEGIN
				y1 = d1  & x = x1
				y2 = Interpol(d2, x2, x1)        ; mesh to densest data set
			endif else BEGIN
				y2 = d2 & x = x2
			    y1 = Interpol(d1, x1, x2)
			endelse
		endelse
	endif else goto, bailout
	tmp1 = {t:'1d', x: x, d: y1, xl: tmp1.xl,  dn: y1, dl: tmp1.dl}
	tmp2 = {t:'1d', x: x, d: y2, xl: tmp2.xl, dn: y2, dl: tmp2.dl}

	return, 1
endif
bailout:
print, ' ax_mesh: Incompatible arrays for interpolation'
return, 0
end
