;+
; NAME:
;	vu_syncgrid
; PURPOSE:
;	Change resolution of tomography arrays
; CATEGORY:
;	smei/sat/idl/util/vupack
; CALLING SEQUENCE:
	FUNCTION vu_syncgrid, old_hdr, old_ff, $
		ref_hdr			= ref_hdr_		, $
		ref_ff			= ref_ff_		, $
		new_hdr			= new_hdr_		, $
		new_ff			= new_ff_		, $
		silent			= silent
; INPUTS:
;	old_hdr		array[k]; type: structure
;					headers of tomography files
;	old_ff		array[n,l,m,2,k]; type: float
;					3d heliospheric arrays that need to be resampled
; OPTIONAL INPUT PARAMETERS:
;	ref_hdr = ref_hdr
;				array[K]; type: structure
;					reference headers
;	ref_ff  = ref_ff
;				array[N,L,M,2,k]; type: float
;					reference 3d heliospheric arrays used
; OUTPUTS:
;	result		scalar; type: integer
;				0: error
;				1: resampling successfull
; OPTIONAL OUTPUT PARAMETERS:
;	new_hdr = new_hdr
;				array[k]; type: structure
;					update headers matching the output array 'result'
;					These will be the same as the input 'old_hdr' array
;					with the grid parameters updated with values
;					from 'ref_hdr'
;	new_ff = new_ff
;				array[N,L,M,2,k]
;					resampled array
; INCLUDE:
	@compile_opt.pro		; On error, return to caller
; CALLS:
;	InitVar
; PROCEDURE:
;	Note that the last dimension MUST be the same
;	as for input array 'ff'.
; MODIFICATION HISTORY:
;	APR-2013, Paul Hick (UCSD/CASS; pphick@ucsd.edu)
;-

InitVar, silent  , 0

nhdr = n_elements(old_hdr)
nref = n_elements(ref_hdr_)

IF nhdr GT 1 AND nref EQ 1 THEN BEGIN

	ref_hdr = replicate(ref_hdr_,nhdr)
	ref_ff  = SuperArray(ref_ff_,nhdr,/trail)

ENDIF ELSE IF nhdr NE nref THEN BEGIN

	destroyvar, new_hdr_, new_ff_
	say, threshold=2, silent=silent, ['input has different nr of 3D arrays then reference set','... ABORTING']
	RETURN, 0

ENDIF ELSE BEGIN

	ref_hdr = ref_hdr_
	ref_ff	= ref_ff_

ENDELSE

old_dims = size(old_ff,/dim)	; Input dimensions
ref_dims = size(ref_ff,/dim)	; Reference dimensions

; Check that dimensions are really different

n = where(old_dims[0:2] NE ref_dims[0:2])

IF n[0] EQ -1 THEN BEGIN	; Same spatial dimensions: return input unmodified
	new_hdr_ = old_hdr
	new_ff_  = old_ff
	RETURN, 1
ENDIF

; Set up dimension of new array by taking the old
; dimensions and replacing the spatial dims that change

new_dims    = old_dims
new_dims[n] = ref_dims[n]

new_ff = make_array(dim=new_dims,type=IsType(old_ff))

; The new data are rebinned to match the spatial grid of the
; reference grid, so the new headers inherit all the spatial
; parameters from the reference headers.

new_hdr = old_hdr

new_hdr = vu_set(new_hdr, nlng				= vu_get(ref_hdr, /nlng				))
new_hdr = vu_set(new_hdr, array_range		= vu_get(ref_hdr, /array_range		))
new_hdr = vu_set(new_hdr, roi_range			= vu_get(ref_hdr, /roi_range		))

new_hdr = vu_set(new_hdr, nlat				= vu_get(ref_hdr, /nlat				))
new_hdr = vu_set(new_hdr, latitude_range	= vu_get(ref_hdr, /latitude_range	))

new_hdr = vu_set(new_hdr, nrad				= vu_get(ref_hdr, /nrad				))
new_hdr = vu_set(new_hdr, start_distance	= vu_get(ref_hdr, /start_distance	))
new_hdr = vu_set(new_hdr, distance_step		= vu_get(ref_hdr, /distance_step	))

FOR i=0,nhdr-1 DO BEGIN

	; Longitude (Carrington variable) grid

	new_range = reverse( vu_get(ref_hdr[i], /array_range) )
	old_range = reverse( vu_get(old_hdr[i], /array_range) )

print, new_range
print, old_range
	new_delta = (new_range[1]-new_range[0])/(vu_get(ref_hdr[i],/nlng)-1)
	old_delta = (old_range[1]-old_range[0])/(vu_get(old_hdr[i],/nlng)-1)
print, new_delta
print, old_delta
	lng_grid = (new_range[0]-old_range[0]+new_delta*findgen(new_dims[0]))/old_delta
	lng_grid = gridgen(new_dims[0],range=(new_range-old_range[0])/old_delta)

	; Latitude grid

	new_range = vu_get(ref_hdr[i], /latitude_range)
	old_range = vu_get(old_hdr[i], /latitude_range)

	new_delta = (new_range[1]-new_range[0])/(vu_get(ref_hdr[i],/nlat)-1)
	old_delta = (old_range[1]-old_range[0])/(vu_get(old_hdr[i],/nlat)-1)

	lat_grid = (new_range[0]-old_range[0]+new_delta*findgen(new_dims[1]))/old_delta
	lat_grid = gridgen(new_dims[1],range=(new_range-old_range[0])/old_delta)

	; Radial grid

	new_range = [vu_get(ref_hdr[i], /start_distance), vu_get(ref_hdr[i], /stop_distance)]
	old_range = [vu_get(old_hdr[i], /start_distance), vu_get(old_hdr[i], /stop_distance)]

	new_delta = vu_get(ref_hdr[i], /distance_step)
	old_delta = vu_get(old_hdr[i], /distance_step)

	rad_grid = (new_range[0]-old_range[0]+new_delta*findgen(new_dims[2]))/old_delta
	rad_grid = gridgen(new_dims[2],range=(new_range-old_range[0])/old_delta)

	FOR j=0,new_dims[3]-1 DO new_ff[*,*,*,j,i] = interpolate(old_ff[*,*,*,j,i],lng_grid,lat_grid,rad_grid,/grid)

ENDFOR

RETURN, 1  &  END