;+
; NAME:
;	InterpolateHeliosphere
; PURPOSE:
;	Get interpolated function values on scalar heliospheric matrix at specified
;	locations in the heliosphere.
; CATEGORY:
;	sat/idl/util
; CALLING SEQUENCE:
	FUNCTION InterpolateHeliosphere, R, F3D, R3D, dR3D,		$
		xcrange 	= xcrange	, $
		xlrange 	= xlrange	, $
		is_carrington=is_carrington, $
		cv2carrington=cv2carrington, $
		opengrid	= opengrid	, $
		periodic	= periodic	, $
		xcgrid		= xcgrid	, $
		xlgrid		= xlgrid	, $
		rrgrid		= rrgrid	, $
		fillbad 	= fillbad	, $
		degrees 	= degrees

;	(1) Matrix specified on implicit regular grid in spherical coordinates:
;		(usually heliographic coordinates)
;
;	F = InterpolateHeliosphere(R, F3D, R3D, dR3D,	$
;		[xcrange=xcrange, xlrange=xlrange, /cv2carrington, /opengrid, /periodic, /fillbad, /degrees])
;
;	(2) Matrix on explicitly specified grid in spherical coordinates:
;		(usually heliographic coordinates)
;
;	F = InterpolateHeliosphere(R, F3D, xcgrid=xcgrid, xlgrid=xlgrid, rrgrid=rrgrid, 	$
;		[/fillbad, /degrees])
; INPUTS:
;	All distances (R[2,*], R3D, dR3D) must have the same units (usually AU).
;	The coordinate system used for R must be consistent with the grid on which
;	F3D is defined. For the tomography heliographic coordinates are used, but
;	other coordinate systems (e.g. ecliptic) are allowed.
;
;	(1) Matrix specified on regular grid in spherical coordinates:
;
;	R			array[3,*]; type: float
;					heliospheric locations where interpolated values are needed
;					R[0,*] longitudes (angle in radians or degrees, or Carrington variables)
;					R[1,*] latitudes (angle in radians or degrees)
;					R[2,*] heliocentric distance
;	F3D			array[L,M,N,K]; type: float
;					3D matrix of function values of scalar heliospheric quantity
;					1 dim: longitude angles or Carrington variables
;						The range is set by keyword xcrange. xcrange covers 360 degrees.
;						For longitude angles: l=0,L-1 covers [xcrange[0],xcrange[1]]
;							where xcrange[1]-xcrange[0] is 360 degrees.
;						For Carrington variables: l=0,L-1 covers [xcrange[1],xcrange[0]]
;							where xcrange[1]-xcrange[0] is 1
;					2 dim: latitude angles
;						The range is set by keyword xlrange.
;						m=0,M-1 covers [xlrange[0],xlrange[1] degrees
;
;					For both angles the grid is closed (with the outer grid points on
;					the edges of the ranges. If /opengrid is set the angular
;					grids are open (with the outer grid points half a grid spacing
;					away from the edges of the ranges).
;
;					3 dim: heliocentric distance: n=0,N-1 covers R3D+i*dR3D
;					4 dim: absent for scalar field; 3 for the 3 components for a vector field
;	R3D			scalar; type: float
;					heliocentric distance of inner boundary of F3D
;	dR3D		scalar; type: float
;					heliocentric distance resolution of F3D
;
;	See keywords xcrange, xlrange, opengrid, periodic, is_carrington and
;	cv_carrington for more control over the implicit grid for F3D.
;
;	(2) Matrix on explicitly specified grid in heliographic coordinates:
;
;	R			array[3,*]; type: float
;					heliospheric locations where interpolated values are needed
;					R[0,*] longitude angles (radians or degrees) or Carrington variables
;					R[1,*] latitude angles (radians or degrees)
;					R[2,*] heliocentric distance
;	F3D			array[L,M,N,K]; type: float
;					3D matrix of function values of scalar or vector heliospheric quantity
;
;	xcgrid=xcgrid
;				array[L]; type: any
;					grid of longitude angles or Carrington variables, matching
;					R[0,*] (covering 1st dimension of F3D)
;	xlgrid=xlgrid
;				array[M]; type: any
;					grid of latitude angles matching R[1,*]
;					(covering 2nd dimension of F3D)
;	rrgrid=rrgrid
;				array[N]; type: any
;					grid of heliocentric distances matching R[2,*]
;					(covering 3rd dimension of F3D)
;
; OPTIONAL INPUT PARAMETERS:
;	/degrees	if set then the longitude and latitude are assumed to be in
;					degrees (default: radians)
;	/fillbad	if set then GridFill is called first to fill in all 'bad' values
;
;	(1) Matrix specified on implicit regular grid in spherical coordinates:
;
;	/is_carrington
;				if SET then xcrange and R[0,*] are assumed to be
;				Carrington variables
;	/cv2carrington
;				if SET then xcrange is assumed to be specified as
;					as Carrington variable and R[0,*] is assumed to
;					be heliographic longitude.
;					The longitudes are converted to Carrington variables using
;						Carrington( mean(xcrange), R[0,*], /get_variable )
;					These values will be within 0.5 of mean(xcrange) if the longitudes
;					are in the range [0,360]
;					(the input R is NOT MODIFIED!!)
;
;	If neither /is_carrington nor /cv2carrington is set then xcrange and R[0,*]
;	are assumed to be longitude angles
;
;	xcrange=xcrange
;				scalar, array[2]; type: float; default: [0,360]
;					(or [0,1] if /is_carrington or /cv2carrington is SET). 
;					If xcrange is a scalar then xcrange+[0,360] is assumed.
;					(or xcrange+[0,1] if /is_carrington or /cv2carrington is SET).
;					Range of longitude angles or Carrington variables covered by
;					array F3D.
;					!!! the default [0,1] is useful only if /cv2carrington is
;						set and F3D corresponds to an integer rotation (longitude
;						range [0,360])
;	xlrange=xlrange
;				scalar; array[2]; type: float; default: [-90,90]
;					If xlrange is a scalar then xlrange*[-1,1] is assumed.
;					Range of heliographic latitudes covered by array F3D
;	/opengrid	if set then the angular grids (longitude and latitude) are assumed
;					to be open. By default the grids are assumed closed.
;	/periodic	maps all R[0,*] values back into the range defined by xcrange
;					Only useful if F3D represents exactly one rotation.
; OUTPUTS:
;	F			array[*]; type: float
;					interpolated function values
; INCLUDE:
	@compile_opt.pro			; On error, return to caller
; CALLS:
;	InitVar, ToRadians, GridFill, Carrington, BadValue, IsType, destroyvar
; PROCEDURE:
; >	Bad values are identified with the IDL 'finite' function
; >	The location R are translated into index values. The IDL function
;	interpolate is used to obtain an interpolated value. 'Missing' values
;	are set to BadValue(F3D).
;
;	(1) Matrix specified on implicit regular grid in spherical coordinates:
;
; >	If /is_carrington is SET then R[0,*] and xcrange MUST contain Carrington variables
; >	If /cv2carrington is SET then R[0,*] MUST contain heliographic longitudes and
;		xcrange MUST contain Carrington variables
; >	If neither are SET then R[0,*] and xcrange MUST contain longitude angles in
;		a spherical coordinate system.
;
; MODIFICATION HISTORY:
;	OCT-1999, Paul Hick (UCSD/CASS)
;	SEP-2002, Paul Hick (UCSD/CASS)
;		Added check for bad entries in positions R. If any one of the three
;		coordinates is bad the matching entry in F is set to BadValue(F3D)
;	NOV-2003, Paul Hick (UCSD/CASS)
;		Fixed problem with replicate statement (didn't work in IDL 5.3)
;	JUN-2006, Paul Hick (UCSD/CASS)
;		Changed order of arguments (F3D is now 2nd instead of 4th argument)
;		Added keywords xcgrid, xlgrid and rrgrid to handle matrices defined
;		on irregular grids.
;	OCT-2006, Paul Hick (UCSD/CASS)
;		Modified to allow interpolation on vector fields.
;	FEB-2007, Paul Hick (UCSD/CASS; pphick@ucsd.edu)
;		Generalized to work with spherical coordinate systems other the
;		just heliographic. Introduces keywords /is_carrington and /cv2carrington
;		to support the special case of heliographic coordinates
;		(used by the heliospheric tomography programs).
;-

InitVar, fillbad	, /key

rpm  = ToRadians(degrees=degrees)
BadV = BadValue(F3D)

szR  = size(R)						; szR[1] = 3
npnt = szR[szR[0]+2]/szR[1] 		; Number of points
IF size(F3D,/n_dim) EQ 4 THEN nvec = (size(F3D))[4] ELSE nvec = 1
									; Nr of components in vector (usually 3)
IF npnt GT 1 THEN dimFF = szR[2:szR[0]]
IF nvec GT 1 THEN boost, dimFF, nvec

FF = BadV							; Set up output array with proper dims
IF npnt*nvec GT 1 THEN FF = replicate(FF, npnt, nvec)

RR = reform(R, szR[1], npnt)		; Reform locations to 2D array; szR[1] = 3

GoodR = where( finite(RR[0,*]) AND finite(RR[1,*]) AND finite(RR[2,*]), n )

IF n NE 0 THEN BEGIN

	IF n NE npnt THEN RR = RR[*,GoodR]; Only keep valid locations

	CASE IsType(R3D,/defined) OF

	0: BEGIN						; Explicit grid specified for F3D

		XCI = dindgen(n_elements(xcgrid))
		XLI = dindgen(n_elements(xlgrid))
		RRI = dindgen(n_elements(rrgrid))

		Y2XC = spl_init(xcgrid, XCI, /double) 
		Y2XL = spl_init(xlgrid, XLI, /double) 
		Y2RR = spl_init(rrgrid, RRI, /double) 

		; Convert heliographic locations to array indices

		RR[0,*] = spl_interp(xcgrid, XCI, Y2XC, RR[0,*])
		RR[1,*] = spl_interp(xlgrid, XLI, Y2XL, RR[1,*])
		RR[2,*] = spl_interp(rrgrid, RRI, Y2RR, RR[2,*])

	END

	1: BEGIN						; Implicit grid specified for F3D

		InitVar, is_carrington, /key
		InitVar, cv2carrington, /key
		InitVar, periodic	  , /key
		InitVar, opengrid	  , /key

		; If is_carrington or cv2carrington is set then
		; xcrange is a Carrington variable.

		CASE is_carrington OR cv2carrington OF
		0: full_circle = 2*!dpi/rpm ; 360 degrees
		1: full_circle = 1			; One rotation
		ENDCASE

		CASE n_elements(xcrange) OF
		0   : xc = [0,full_circle]
		1   : xc = xcrange[0]+[0,full_circle]
		ELSE: xc = xcrange[0:1]
		ENDCASE

		; Latitude range

		CASE n_elements(xlrange) OF
		0   : xl = !dpi/2*[-1,1]/rpm; -90 to +90 degrees
		1   : xl = xlrange[0]*[-1,1]
		ELSE: xl = xlrange[0:1]
		ENDCASE

		; If /cv2carrington is set, convert the longitudes in RR[0,*]
		; from heliographic longitudes to Carrington variables.

		IF cv2carrington THEN	$
			RR[0,*] = Carrington( mean(xc), near_longitude=RR[0,*], degrees=degrees, /get_variable)

		; If keyword periodic is set, map RR[0,*] to range xc.

		IF periodic THEN BEGIN
			dx = xc[1]-xc[0]
			RR[0,*] = ( ( ( (RR[0,*]-xc[0]) mod dx ) + dx ) mod dx ) + xc[0]
		ENDIF

		szF3D = size(F3D)

		; Adjust angle ranges xc and xl for open grids.

		IF opengrid THEN BEGIN
			xc = xc+0.5d0*(xc[1]-xc[0])/szF3D[1]*[1,-1]
			xl = xl+0.5d0*(xl[1]-xl[0])/szF3D[2]*[1,-1]
		ENDIF

		; Carrington variables decrease with increasing
		; array index. Reverse xc to take this into account.

		IF is_carrington OR cv2carrington THEN xc = reverse(xc)

		; Convert heliographic locations to array indices

		RR[0,*] = (RR[0,*]-xc[0])/(xc[1]-xc[0])*(szF3D[1]-1)
		RR[1,*] = (RR[1,*]-xl[0])/(xl[1]-xl[0])*(szF3D[2]-1)	; [-90,90] -> [0,nDim[1]-1]
		RR[2,*] = (RR[2,*]-R3D  )/dR3D							; [R3D,.]  -> [0,.]

	END

	ENDCASE

	IF fillbad THEN BEGIN		; Remember bad values (needed to restore F3D to input value)
		BadF3D = where(1B-finite(F3D))
		IF BadF3D[0] NE -1 THEN 	$
			FOR n=0,nvec-1 DO		$
				FOR i=0,szF3D[3]-1 DO F3D[*,*,i,n] = GridFill(F3D[*,*,i,n])
	ENDIF

	; Interpolate on F3D array to get function values on all lines of sight

	FOR n=0,nvec-1 DO	$
		FF[GoodR,n] = interpolate(F3D[*,*,*,n], RR[0,*], RR[1,*], RR[2,*], missing=BadV)

	IF npnt*nvec GT 1 THEN FF = reform(FF, dimFF, /overwrite)

	IF fillbad THEN IF BadF3D[0] NE -1 THEN F3D[BadF3D] = BadV	; Restore input array to original state

ENDIF

RETURN, FF  &  END