; Copyright (c) 1998-2007 A.P. Hitchcock  All rights reserved
;+
;NAME:
;	AX_INTERP
;
;LAST CHANGED: ----------------------------------- 04-Aug-07
;
;PURPOSE:
;	This function changes meshes of 1-d (spectra) or
; 2-d (images) AXIS data buffers, based on user input
; used to force image or spectra to a particular spacing
;
;CATEGORY:
;	AXIS: image and spectral processing
;
;CALLING SEQUENCE:
;	RESULT = AX_INTERP( tmp, mesh=mesh )
;
;CALLED FROM AXIS:
;	(utilities~change mesh) TMP -  spectrum or image
;
;KEYWORDS:
;	MESH 	[x-pixels, y-pixels]
;
; OUTPUT
;	RESULT is the input data converted to modified mesh
;
;COMMON BLOCKS:
;	@AXIS_COM	standard set of common blocks
;
;MODIFICATION HISTORY:
; (01-nov-00 aph) first version
; (31-dec-01 aph) add MESH keyword for auto operation
; (02-feb-05 aph) unknown update
; (04-aug-07 aph) corrected for distortion (see 70804-utilities-mesh-bug)
;-

FUNCTION ax_interp,tmp, mesh=mesh
@axis_com
on_error,2

if widget_info(/active) then axis_on=1 else axis_on=0
; ------------------------- images -----------
if tmp.t EQ '2d' then begin
	nx = n_elements(tmp.x)
	ny = n_elements(tmp.y)
	dx = (tmp.x(n_elements(tmp.x)-1) - tmp.x(0))/nx
	dy = (tmp.y(n_elements(tmp.y)-1) - tmp.y(0))/ny

; --------- by-pass if mesh defined from keyword
	if NOT keyword_Set(mesh) then begin
		text = string(format='(i4,"x",i4," pixels. ",f6.3,"x",f6.3," um.")', $
		   nx,ny,dx,dy)
	   	if axis_on then axis_log, text else print, text
	   	text = string(nx, format='("# of X pixels. Old = ",i4)')
		IF N_ELEMENTS(mesh_x) EQ 0 THEN mesh_x = nx
		IF N_ELEMENTS(mesh_y) EQ 0 THEN mesh_y = ny
		if (axis_on) then $
		     nx = get_num(prompt=text, val=mesh_x, group=axis_ID) $
		else nx = get_num(prompt=text, val=mesh_x)
		dx = (tmp.x(n_elements(tmp.x)-1) - tmp.x(0))/nx
		new_ny = fix((tmp.y(n_elements(tmp.y)-1) - tmp.y(0))/dx)
		text = string(format= $
		    '("New values: dx = ",f6.3, ".   # y pts for square pixels = ",i4)', dx, new_ny)
	   	if axis_on then axis_log, text else print, text
	   	text = string(ny, format='("# of Y pixels. Old = ",i4)')
		if (axis_on) then ny = get_num(prompt=text,val=mesh_y,group=axis_ID) $
			else ny = get_num(prompt=text,val=mesh_y)
	endif else begin
		nx = mesh(0)
		ny = mesh(1)
	endelse

; interpolate tmp data to new mesh defined by (nx, ny)
	dx = (tmp.x(n_elements(tmp.x)-1) - tmp.x(0))/nx
	dy = (tmp.y(n_elements(tmp.y)-1) - tmp.y(0))/ny
	xn = findgen(nx)*dx + tmp.x(0)			; reversed order of dx versus xn calc
	yn = findgen(ny)*dy + tmp.y(0)

;	xn = congrid(tmp.x,nx,/interp) ; original code - caused rounding errors
;	yn = congrid(tmp.y,ny,/interp)
	dn = congrid(tmp.d,nx,ny, /interp, cubic=-0.5)
	tmp1 = {t:'2d', x:xn, y:yn, xl:tmp.xl, yl:tmp.yl, d:dn, e:tmp.e, dl:'M '+tmp.dl}
	return, tmp1	; returns modified data file
endif


; ----------------- re-mesh 1-d (spectral) data ----------
if tmp.t EQ 1 then begin
	nx = n_elements(tmp.x)
	if (axis_on) then $
		 nx = get_num(prompt='# of E points',val=nx,group=axis_ID) $
	else nx = get_num(prompt='# of E points',val=nx)
; interpolate tmp data to new mesh
	dx = (tmp.x(n_elements(tmp.x)-1) - tmp.x(0))/nx
	xn = findgen(nx)*dx + tmp.x(0)
;	xn = congrid(tmp.x,nx,/interp)       ; original
	dn = congrid(tmp.d,nx,/interp)
	tmp1 = {t:'1d', x:xn, d:dn, xl:tmp.dl, dn:dn, dl:'M '+tmp.dl}
	return, tmp1	; returns modified data file
endif
end

