;+
; NAME:
;	vu_atlocation
; PURPOSE:
;	Extract a function values from heliospheric matrix at specified
;	locations (heliographic longitude, latitude, distance and time).
; CATEGORY:
;	WWW
; CALLING SEQUENCE:
	FUNCTION vu_atlocation, XC3D, R3D, dR3D, F3D, sequence=sequence,	$
		location=location, times=times, power=power, degrees=degrees, _extra=_extra
; INPUTS:
;	The 3D heliospheric matrix is described by the following input:
;
;	XC3D			array[2] or array[2,ntime]; type: float
;						Carrington range of reconstruction
;	R3D				scalar or array[ntime]; type: float
;						heliocentric distance at inner boundary of matrix
;	dR3D			scalar or array[ntime]; type: float
;						heliocentric distance resolution of the heliospheric matrix
;	F3D				array[n,m,l,ntype,ntime]; float
;						3D heliocentric matrix;
;						1st dim: Carrington variable, covering range XC3D
;						2nd dim: Heliographic latitude, covering [-90,90] degrees
;						3rd dim: Heliocentric distance, with inner boundary at R3D and resolution dR3D
;						4th dim: different data types (usually 2: velocity/density, Br/Bt
;						5th dim: # members in a (time-dependent) sequence
; OPTIONAL INPUT PARAMETERS:
;	sequence=sequence
;					array[ntime]; type: time structure of float; default: none
;						if specified, these are the times corresponding to the ntime
;						data arrays specified in argument F3D.
;	power=power		array[ntype]
;						power of heliocentric distance to be removed from the time series
;						(2 for normalized density; 2 for radial component of the magnetic field;
;						 1 for tangential component)
;	location=location
;					array[3,N]; type: float
;						heliographic locations: longitude or Carrington variable, latitude and
;						heliocentric distance.
;	times=times		array[N]; type: float or time structure
;						times associated with N locations.
;
;	/degrees		if set, then all angles are in degrees (default: radians)
;
;	/is_carrington	passed to href=InterpolateHeliosphere=
;	/cv2carrington	passed to href=InterpolateHeliosphere=
;	/periodic		passed to href=InterpolateHeliosphere=
; OUTPUTS:
;	F				array[N]; float
;						the function values for all N locations
; INCLUDE:
	@compile_opt.pro							; On error, return to caller
; CALLS:
;	Carrington, TimeUnit, SuperArray, SubArray, InterpolateHeliosphere, BadValue, SyncDims
;	TimeOp
; PROCEDURE:
; >	The time associated with a given Carrington variable is the time at which the
;	corresponding heliographic longitude crossed the center of the solar disk.
; MODIFICATION HISTORY:
;	AUG-2002, Paul Hick (UCSD/CASS; pphick@ucsd.edu)
;-

szf3d = size(F3D)
uday  = TimeUnit(/day)

CASE szf3d[0] OF
3: begin
	ntype = 1					; One matrix: one data type, one time
	ntime = 1
END
4: BEGIN
	IF szf3d[4] eq 2 THEN BEGIN
		ntype = szf3d[4]
		ntime = 1
	ENDIF ELSE BEGIN
		ntype = 1
		ntime = szf3d[4]
	ENDELSE
END
5: BEGIN
	ntype = szf3d[4]			; Multiple data types; multiple times
	ntime = szf3d[5]
END
ENDCASE

F3D = reform(F3D, [szf3d[1:3], ntype, ntime], /overwrite)


xc = XC3D
IF n_elements( xc) EQ 2 THEN xc = SuperArray(xc, ntime, /trail)

rr = R3D
IF n_elements( rr) EQ 1 THEN rr = replicate(rr[0], ntime)

drr = dR3D
IF n_elements(drr) EQ 1 THEN drr = replicate(drr[0], ntime)

npnt = (n_elements(location)/3) > n_elements(times)

; First make time series at the sequence times for all 'npnt' positions in 'location'

G = make_array( type=IsType(F3D), dim=[npnt, ntype, ntime] )

FOR n=0,ntime-1 DO			$			; Interpolate all npnt points at all matrices
 	FOR i=0,ntype-1 DO		$
		G[*,i,n] = InterpolateHeliosphere(location, F3D[*,*,*,i,n], rr[n], drr[n],	$
				xcrange=xc[*,n], degrees=degrees, _extra=_extra)

; Here we convert from 'normalized' quantities (read from file) to 'real' quantities:
;	nr^2 --> n;		r^2Br --> Br;	rBt --> Bt
; Note that 'power' needs to be set correctly.

IF IsType(power, /defined) THEN	$
	FOR i=0,ntype-1 DO			$
		IF power[i] NE 0 THEN G[*,i,*] = G[*,i,*]/SuperArray(SubArray(location,element=2)^power[i],ntime,/trail)


; Now interpolate at appropriate times for each point

F = replicate(BadValue(F3D), npnt, ntype)

CASE ntime OF
1   : F = G
ELSE: BEGIN
	ts = Carrington(sequence, /get_time)
	tc = Carrington(times   , /get_time)

	good = where(TimeOp(/subtract,tc,ts[0],uday) GE 0 AND TimeOp(/subtract,tc,ts[ntime-1],uday) LE 0, ngood)

	FOR n=0L,ngood-1 DO		$
		FOR i=0,ntype-1 DO	$
			F[good[n],i] = TimeInterpol(reform(G[good[n],i,*]), ts, tc[good[n]])
END
ENDCASE

SyncDims, F3D, size=szf3d

sz = size(location)

CASE npnt OF
1   : F = reform(F, ntype, /overwrite)
ELSE: F = reform(F, [sz[2:sz[0]],ntype], /overwrite)
ENDCASE

RETURN, F  &  END