FUNCTION AlignHeliosphere, ff, dOff, aligneps=aligneps, nowrap=nowrap

;+
; NAME:
;	AlignHeliosphere
; PURPOSE:
;	For a given set of 3D heliospheric matrices covering different ranges of
;	heliographic longitudes (or, equivalently, Carrington variables) 'align'
;	all matrices so they start and end at the same heliographic longitude,
;	interpolating linearly where necessary.
; CATEGORY:
;	sat/idl/util
; CALLING SEQUENCE:
;	result = AlignHeliosphere(ff, dOff)
; INPUTS:
;	ff			array[n,l,m,k] or array[n,l,m,d,k]; type: float
;					3D heliospheric matrix specified on a heliographic grid
;					n: # heliographic longitudes (can cover more than one rotation)
;					l: # latitudes
;					m: # heliocentric distances
;					d: (optional) # data types (e.g. k=2 if velocity and density
;						are stored in one array)
;					k: # matrices to be aligned (usually this represents time)
;
;	dOff		array[k]; type: float
;					offset between k matrices in grid spacings
; OPTIONAL INPUT PARAMETERS:
;	aligneps=aligneps
;				scalar; type: float
;					if dOff is less than eps from an integer number of grid spacings
;					than it is rounded to the nearest integer. This reduces the
;					alignment to a cyclic shift (bypassing the call to the IDL
;					interpolate function). This should be faster.
;	/nowrap 	by default the alignment is done by applying a mod 360 operation.
;					If /nowrap is set this is suppressed and the missing part of a
;					map after alignment is set to BadValue(0.0)
; OUTPUTS:
;	result		same as input 'ff', but now aligned
; INCLUDE:
	@compile_opt.pro			; On error, return to caller
; CALLS:
;	InitVar, BadValue, destroyvar
; PROCEDURE:
;	This procedure was written to deal with (a sequence of) 3D matrices (longitude, latitude,
;	radial distance). For a matrix with only longitude, or only longitude and latitude, this
;	procedure can be used after inserting a dummy dimension of 1 for the missing coordinate(s).
; MODIFICATION HISTORY:
;	APR-2001, Paul Hick (UCSD/CASS)
;	JUL-2002, Paul Hick (UCSD/CASS)
;		Allowed for special case of 3D array ff[n,l,m] (i.e. k=1 but without degenerate dimension)
;	NOV-2002, Paul Hick (UCSD/CASS; pphick@ucsd.edu)
;		Added /nowrap keyword
;-

InitVar, aligneps, 0.0
InitVar, nowrap  , /key

sz = size(ff)

nLng = sz[1]
nLat = sz[2]
nRad = sz[3]

CASE sz[0] OF
3   : nTim = 1					; Intercept 3D array; set nTim = 1
ELSE: nTim = sz[sz[0]]			; Last dimension is time dimension
ENDCASE

bDat = sz[0] EQ 5				; 4th dimension is data index (e.g. velocity and density)
IF bDat THEN nDat = sz[4]

Lng = indgen(nLng)
Lat = indgen(nLat)
Rad = indgen(nRad)

Bad = BadValue(0.0)

CASE bDat OF					; Create output array
0: ff_align = make_array(dim=[nLng,nLat,nRad	 ,nTim], value=Bad)
1: ff_align = make_array(dim=[nLng,nLat,nRad,nDat,nTim], value=Bad)
ENDCASE

FOR i=0,nTim-1 DO BEGIN

	L = Lng-dOff[i]

	CASE nowrap OF
	0: BEGIN

		Lbeg = 0
		Lend = nLng-1

		s = where(L LT 0)
		IF s[0] NE -1 THEN L[s] = L[s]+nLng-1

		s = where(L GT nLng-1)
		IF s[0] NE -1 THEN L[s] = L[s]-nLng+1

	END
	1: BEGIN
	
		s = where(0 LE L AND L LE nLng-1)
		IF s[0] NE -1 THEN BEGIN
			L = L[s]
			Lbeg = s[0]
			Lend = s[n_elements(s)-1]
		ENDIF ELSE BEGIN
			destroyvar, L				; Redundant, but leave it in anyway
			Lbeg =  0
			Lend = -1
		ENDELSE
	
	END
	ENDCASE

	IF Lend GT Lbeg THEN BEGIN

		; If dOff[i] is an integer then interpolation is not necessary.
		; A cyclic shift probably is much faster.

		IF abs(dOff[i]-round(dOff[i])) LT aligneps THEN BEGIN

			L = round(L)

			CASE bDat OF
			0: ff_align[Lbeg:Lend,*,*,  i] = ff[L,*,*,  i]
			1: ff_align[Lbeg:Lend,*,*,*,i] = ff[L,*,*,*,i]
			ENDCASE

		ENDIF ELSE BEGIN

			CASE bDat OF
			0: ff_align[Lbeg:Lend,*,*,i] = interpolate(ff[*,*,*,i], L, Lat, Rad, /grid, missing=Bad)
			1: FOR j=0,nDat-1 DO ff_align[Lbeg:Lend,*,*,j,i] = interpolate(ff[*,*,*,j,i], L, Lat, Rad, /grid, missing=Bad)
			ENDCASE

		ENDELSE

	ENDIF

ENDFOR

RETURN, ff_align  &  END