NAME:
	ibclrf
 PURPOSE:
	Takes a longword (32-bit) integer, and clears (sets to zero) a bit
 CATEGORY:
	Bits and bytes
 CALLING SEQUENCE:
	Result = ibclrf(bitset, bitpos)
 INPUT:
	bitset	longword    A longword integer, one of it
		    bits will be clear or zero
	bitpos	longword    The position of the bit that
		    will be cleared; 0 <= bitpos <= 31
 OUTPUT:
	Result	longword    Same as bitset, (bitwise) but
		    at bit position, bitpos, it will
		    definitely be a zero
 INCLUDE:
	@compile_opt.pro    ; On error, return to caller
 CALLED BY:
	ice_pack
 SEE ALSO:
	ibsetf, ibtestf
 MODIFICATION HISTORY:
	8/23/00 	Kevin Nguyen

 NAME:
	ibsetf
 PURPOSE:
	Sets a bit in a longword (32-bit) integer
 CATEGORY:
	Bits and bytes
 CALLING SEQUENCE:
	Result = ibsetf(bitset, bitpos)
 INPUT:
	bitset longword integer 32 bits longword integer.
		    Its bit at bit position bitpos
		    will be set at 1
	bitpos longword integer The bit position of bit to be
		    set as 1 (0 <= bitpos <= 31)
 OUTPUT:
	Result longword integer Same at bitset bitwise except
		    at bit position bitpos
 INCLUDE:
	@compile_opt.pro    ; On error, return to caller
 CALLED BY:
	ice_pack, ice_unpack, smei_frm_flag
 SEE ALSO:
	ibclrf, ibtestf
 MODIFICATION HISTORY:
	8/23/00     Kevin Nguyen

 NAME:
	ibtestf
 PURPOSE:
	Tests if the bit at bitpos of a longword (32-bit) integer is set
 CATEGORY:
	Bits and bytes
 CALLING SEQUENCE:
	Result = ibtestf(bitset, bitpos)
 INPUT:
	bitset	longword integer  will be used to test if
		      bit position, bitpos is zero or one
	bitpos	longword integer  the bit postion, bit pos, being tested
 OUTPUT:
	ibtestf     scalar; type:  byte
		1B if bit is set;  0B bit is clear
 INCLUDE:
	@compile_opt.pro	; On error, return to caller
 CALLED BY:
	ice_pack, ice_unpack, smei_frm_flag
 SEE ALSO:
	ibclrf, ibsetf
 MODIFICATION HISTORY:
	8/23/00     Kevin Nguyen

 NAME:
	ice_pack
 PURPOSE:
	Compresses an integer array using a modified Rice algorithm
 CATEGORY:
	Bits and bytes
 CALLING SEQUENCE:
	FUNCTION ice_pack, Orig, Pack, kmax=kmax, kshift=kshift, lenbit=lenbit, sign=sign, perm=perm
 INPUTS:
	Orig	array; type: long integer
		array to be compressed
 OPTIONAL INPUTS:
	kmax	scalar; type: integer; default: 8
		cutoff for applying compression (see PROCEDURE) (kmax > 0)
	kshift	scalar; type: integer; default: 0
		minimum bits copied for each difference
	lenbit	scalar; type: integer; default: 16
		# bits encoded when full value instead of a difference
		is used (see PROCEDURE)
	sign	scalar; type: integer; default: 1
		-1: all values in Orig < zero
		+1: all values in Orig >= zero
		 0: both positive and negative values are present
	perm	array[0:kmax-kshift]; type: integer
		permutation table (see PROCEDURE)
 OUTPUTS:
	ibit	scalar; type: long integer
		# bits in Pack used to hold the compressed data;
		on failure ibit = 0 is returned; this happens only if
		the compressed array Pack is bigger than the input
		array Orig.
	Pack	array; type: long integer
		compressed data;
		The # elements in Pack is 1+(ibit-1)/nbit.
		Bits 0..ibit-1 are filled with compressed data.
		For the last element only bits 0 .. ib-1 are
		used [ib = 1+ ((ibit-1) mod nbit)].
 INCLUDE:
	@compile_opt.pro	; On error, return to caller
 CALLS: ***
	InitVar, ibclrf, ibsetf, ibtestf
 SEE ALSO:
	ice_unpack
 PROCEDURE:
	The technique is based on the algorithm described by Michael W. Richmond and Nancy E.
	Ellman, in March 1996 paper which we have in preprint form.
	COMPRESSION ALGORITHM
	---------------------
	The Rice algorithm stores an integer in the following bit pattern:
	    0.......0 1 x.......x
	    |_______|	|_______|
	      ipad	itop + 1
	a string of ipad 0-bits, followed by a terminating 1 bit, followed by a string of itop + 1
	bits encoding the integer value.

	The actual number encoded is either the full value from the input array, or the difference
	with the previous element in the array (!!! for the first element in the array, the 'previous
	element' is assumed to be zero).

	The two cases are distinguished  using the kmax value: if the absolute difference is less than
	2^(kmax - 1) then the difference is encoded; if not then the full value is encoded

	CASE A: Non-zero absolute difference less than 2^(kmax - 1)
	-----------------------------------------------------------
	2^(kmax -1) has only bit kmax - 1 set. I.e. for absolute differences less than 2^(kmax - 1) this
	bit, and all higher bits, are NOT set, i.e. only the first kmax - 1 bits 0..kmax - 2 might be set.

	Let itop be the highest 1-bit of the absolute difference. For a non-zero difference
	0 <= itop <= kmax-2. Since bit itop is by definition set, the absolute difference
	is fully described by the first itop bits 0...itop-1. An additional bit is needed
	to encode the sign of the difference. So itop+1 bits are needed to store the difference:
	first ipad=itop+1 0-bits, followed by a terminator 1-bit, indicate how many bits are
	used. The terminator bit is then followed by the first itop bits of the absolute difference
	followed by a sign bit. The total number of bits used is 2*itop+3.

	    0..........0 1 x..........x
	    |__________|   |__________|
	    ipad=itop+1       itop+1

	CASE B: Zero difference
	-----------------------
	No bit at all is needed to store a zero, except for a terminator bit. Not that this
	is CASE A with ipad = 0, itop = -1.

	CASE C:  Absolute difference greater/equal 2^(kmax-1)
	-----------------------------------------------------
	In this case the full value instead of the difference is encoded.
	For case A the maximum number of leading 0-bits is kmax-1. A string of kmax 0-bits
	is used to indicate that a full value is encoded. After the terminator bit follow
	the first lenbit bits. The total number of bits used is kmax+1+lenbit

	    0..........0 1 x..........x
	    |__________|   |__________|
	    ipad=kmax	      lenbit

	MINIMUM NUMBER OF BITS
	----------------------

	The above describes the the compression algorithm for the case kshift=0.
	Setting kshift to non-zero value modifies the way differences are stored: for each
	difference the kshift lowest bits are always stored. The differences are now encoded
	as follows:

	Zero difference (ipad = 0)
	    1 0
	The 1-bit is the terminating bit (i.e. no leading 0-bits), followed by a 0-bit to
	indicate the presence of a zero difference.

	Difference needing 1..kshift bits (1 <= ipad <= kshift)
	    1 1 s x.....x
		  |_____|
		   kshift
	The leading 1-bit is the terminating bit (i.e. no leading 0-bits). It is followed by
	a 1-bit to be able to distinguish it from a zero difference. Then a bit follows
	storing the sign of the difference (s-bit). Then follow the lowest kshift bits.

	Differences needing kshift+1...kmax-1 bits (kmax-1 <= ipad <= kmax-1 above)
	    0.......0 1 x........x
		|_______|   |________|
	       ipad-kshift     ipad

	Full value (ipad=kmax)
	    0.......0 1 x........x
		|_______|   |________|
	       kmax-kshift    lenbit


	PERMUTATION TABLE
	-----------------

	The number of leading 0-bits varies from izero=0 to izero=kmax-kshift. These differences
	serve merely as identifiers for the type of information stored after the terminating
	1-bit. Rather than using them directly as described above, it is more efficient to
	use a permutation table based on the histogram of the number of bits needed for to store
	the differences. E.g. if for some sky image 3-bit difference occur most, followed by 2-bit,
	then 1-bit, then 0-bit, then 4-bit, 5-bit, etc. then the permutation table would be
	    perm = 3,2,1,0,4,5...kmax-kshift
 MODIFICATION HISTORY:
	AUG-2000, Paul Hick (UCSD; pphick@ucsd.edu)
	    Rewrite of Andy Buffington's ricearc.for.
	    Main modifications:
	    - large intermediate byte array for storing individual bits is removed
	    - # leading zeroes used to store full values was reduced from kmax+1 to kmax
	    - use of kshift allows for slightly better compression
	8/24/2000, Kevin Nguyen
	    - Rewrite the code from fortran to idl

 NAME:
	ice_read
 PURPOSE:
	Reads a compressed file written by ice_write into an integer array
 CATEGORY:
	Bits and bytes
 CALLING SEQUENCE:
	FUNCTION ice_read, cFile, Orig, info=info
 INPUTS:
	cFile	    scalar; type: string
		    file name of compressed file
 OUTPUTS:
	Result	    scalar; type: integer
		    0: failure (an error message is displayed)
		    1: success
	Orig	    array; type: long integer
 OPTIONAL OUTPUT PARAMETERS:
	info=info   scalar; type: string
		    trailer string extracted from end of file
 INCLUDE:
	@compile_opt.pro    ; On error, return to caller
 CALLS: ***
	IsType, ice_unpack
 SEE ALSO:
	ice_pack, ice_write
 MODIFICATION HISTORY:
	AUG-2000, Paul Hick (UCSD/CASS; pphick@ucsd.edu)
	AUG-2000, Kevin Nguyen translated from Fortran to IDL

 NAME:
	ice_unpack
 PURPOSE:
	Deompresses an integer array compressed using ice_pack
 CATEGORY:
	Bits and bytes
 CALLING SEQUENCE:
	FUNCTION ice_unpack, ibit, Pack, Orig, kmax=kmax, kshift=kshift, lenbit=lenbit, sign=sign, perm=perm
 INPUTS:
	ibit	    scalar; type: integer
		    # bits in Pack containing compressed data.

 OPTIONAL INPUTS:
	kmax	    scalar; type: integer; default: 8
		    cutoff for applying compression
	kshift	    scalar; type: integer; default: 0
		    minimum # bits copied for each difference
	lenbit	    scalar: type: integer; default: 16
		    # bits encoded  when the full value instead
		    of a difference is used
	sign	    scalar; type: integer
		    -1: all values in Orig are <  zero
		    +1: all values in Orig are >= zero
		     0: both positive and negative values are present
	perm	    array[0:kmax-kshift]; type: integer; default: lindgen(kmax-kshift+1)
		    permutation table
	Pack	    array; type:long integer
		    compressed data
 OUTPUTS:
	Result	    scalar; type: long integer
		    # bits in Pack used to hold the compressed data
		    as determined during the decompression. This
		    should match the input value ibit.
	Orig	    array; long integer
		    decompressed array
 INCLUDE:
	@compile_opt.pro    ; On error, return to caller
 CALLS: ***
	InitVar, ORIG, ibsetf, ibtestf
 CALLED BY:
	ice_read
 SEE ALSO:
	ice_pack
 PROCEDURE:
 >	The input values for ibit, kmax, lenbit, sign, perm must be the same as for
	the ice_pack call that created the compressed array Pack.
 >	Inverts the compression of ice_pack.
 MODIFICATION HISTORY:
	AUG-2000, Paul Hick (UCSD; pphick@ucsd.edu)

 NAME:
	InfiniteValue
 PURPOSE:
	Provide the infinity value for any data type
 CATEGORY:
	sat/idl/toolbox
 CALLING SEQUENCE:
	FUNCTION InfiniteValue, X
 INPUTS:
	X	    variable of any type
 OUTPUTS:
	Result	    scalar; type: as implied by 'size' vector
			NaN value
 INCLUDE:
	@compile_opt.pro    ; On error, return to caller
 CALLS: ***
	IsType
 CALLED BY:
	ThomsonLOSFar, ThomsonLOSRomb
 PROCEDURE:
	For integers the number zero is returned.
	For floating point the appropriate value from !values is returned.
 MODIFICATION HISTORY:
	MAR-2000, Paul Hick (UCSD/CASS; pphick@ucsd.edu)

 NAME:
	Inside_Wedge
 PURPOSE:
	Tests whether image locations are inside a wedge
 CATEGORY:
 CALLING SEQUENCE:
	result = Inside_Wedge(p_box,p)
 INPUTS:
	p_box	    array[2,2]; type: float
			two opposite corners of the wedge in polar coordinates
			in the form [ [angle1,radius1],[angle2,radius2] ]
	p	    array[2,n,m,...]; type: float
			image locations in polar coordinates
			(same form as p_box)
 OPTIONAL INPUT PARAMETERS:
	exclude_p_box=exclude_p_box
		    array[2,2]; type: float
			two opposite corners of another wedge in polar coordinates
			in the form [ [angle1,radius1],[angle2,radius2] ]
			This second wedge must lie entirely inside the first one.
			If specified the coordinate array p is tested for inside
			p_box and outside exclude_p_box.
 OUTPUTS:
	result	    array[n,m,...]; type: byte
			0 for locations p outside the wedge
			1 for locations inside the wedge
 CALLS: ***
	IsType, SubArray
 CALLED BY:
	TMO_tracksky, cvsmei, even_light, qGlitch_Run, qImage_cw_BoxCosine
	qImage_cw_SmeiMask, qImage_cw_ZEllipse, qView_PlotSeries, smei_ccd2sky
	smei_sky2ccd, wedge_content
 RESTRICTIONS:
	If p contains !values.f_nan the value 1 is returned
 PROCEDURE:
 >	Formerly called qImage_cw_WTest
 >	The array p can have more than two dimensions, e.g. [2,n,m].
 MODIFICATION HISTORY:
	FEB-2000, Paul Hick (UCSD/CASS; pphick@ucsd.edu)

 NAME:
	InSitu
 PURPOSE:
	Plots in-situ data from various spacecraft
 CATEGORY:
 CALLING SEQUENCE:
	PRO InSitu, T,	$
	    delt    = delt,	$
	    weight  = weight,	$

	    source  = source,	$
	    step    = step,	$

	    xtime   = xtime,	$
	    xlng    = xlng,	$
	    xlat    = xlat,	$

	    ylng    = ylng,	$
	    ylat    = ylat,	$
	    ydis    = ydis,	$
	    yden    = yden,	$
	    yvel    = yvel,	$
	    ybr     = ybr ,	$
	    ybt     = ybt ,	$
	    ybn     = ybn ,	$
	    ybb     = ybb ,	$
	    ykp     = ykp ,	$
	    ydst    = ydst,	$
	    ypdyn   = ypdyn,	$
	    yvlng   = yvlng,	$
	    yvlat   = yvlat,	$

	    traceback=traceback,$
	    xrange  = xrange,	$
	    charsize=charsize,	$

	    file    = file,	$
	    nocvsky = nocvsky,	$

	    _extra  = _extra
 INPUTS:
	T	    array[2], type: standard time structure
			    time range for time series
 OPTIONAL INPUT PARAMETERS:
	/xtime, /xlng,	/xlat
			select one to set the quantity plotted on the x-axis: time, heliographic longitude
			or heliographic latitude. Only one can be specified; if none is specified
			then /xtime is assumed
	source=source
		scalar; type: integer
			identifies in situ instrument (see function Instrument)
	/ylng,	/ylat,	/ydis,	/yden,	/yvel, /ybr, /ybt, /ybn, /ybb, /ykp, /ydst, /ypdyn, /yvlng, yvlat
			detemines quantity plotted on y-axis. Multiple keywords can be specified.
			If none is specified then velocity and density is plotted.
	traceback=traceback
			scalar; type: any; default:0 (no traceback)
			heliocentric traceback distance in AU; adjusts time and heliographic longitude to
			give a 'time and location of origin' at the traceback distance
 OUTPUTS:
	Output to screen
 OPTIONAL OUTPUT PARAMETERS:
	(none)
 INCLUDE:
	@compile_opt.pro		; On error, return to caller
 CALLS: ***
	Carrington, INSITUTIMESERIES0, InitVar, Instrument, IsType, PlotCurve, TimeGet
	TimeSet, TimeUnit, TimeXAxis
 PROCEDURE:
 MODIFICATION HISTORY:
	SEP-1999, Paul Hick (UCSD/CASS)
	DEC-2001, Paul Hick (UCSD/CASS)
	    added call to Carrington to allow input of time as Carrington variable
	SEP-2006, Paul Hick (UCSD/CASS)
	    Added /ypdyn keyword to process Mars Global Surveyor dynamic pressure data
	SEP-2006, Paul Hick (UCSD/CASS; pphick@ucsd.edu)
	    Added /yvlng and /yvlat keywords to plot velocity angles from the
	    SOMNI and EOMNI databases.

 NAME:
	InsituTimeSeries
 PURPOSE:
	Extract in-situ data in specified interval, and average if necessary
 CALLING SEQUENCE:
	FUNCTION InsituTimeseries,$
	    data_request	, $
	    t_request		, $
	    t_given		, $
	    data_given		, $
	    delt    = delt	, $
	    weight  = weight	, $
	    traceback=traceback , $
	    source  = source_	, $
	    nocvsky = nocvsky	, $
	    silent  = silent
 INPUTS:
	data_request	array, type: string
			    determines which data are extracted.
			    Identifiers are:
			    'lng','lat','dis','den','vel','br','bt','bn','bb','pdyn','vlng','vlat'
			    'lng'   = heliographic longitude (deg)
			    'lat'   = heliographic latitude (deg)
			    'dis'   = heliocentric distance (AU)
			    'den'   = density (cm^-3)
			    'vel'   = velocity (km/s)
			    'br'    = radial magnetic field (nT)
			    'bt'    = tangential magnetic field (nT)
			    'bn'    = normal magnetic field (nT)
			    'bb'    = field strength (nT)
			    'pdyn'  = dynamic pressure (nPa)
			    'vlng'  = azimuthal flow angle
			    'vlat'  = elevation flow angle
	t_request	array[2],array[N] type: standard time structure
			    If array[2]: data are extracted between t_request[0] and t_request[1]
				at times extracted from s/c insitu files
			    If array[N]: data are extracted at times t_request
			    Should be in chronological order
 OPTIONAL INPUT PARAMETERS:
	source=source	scalar; type: integer or string
			    identifies insitu instrument either by name or integer
			    as provided by function Instrument.
	delt=delt	scalar; type: time difference structure, or float; default: 0
			    if present, and non-zero, time window used for averaging.
			    in-situ data inside [t-delt/2,t+delt/2] are averaged.
			    If a scalar (float or integer) is specified, then it is
			    assumed to be the time window in days.
	/weight 	if set, standard deviations are used to weight the averages (if available
			    for selected instrument and data)
			    if not set, an unweighted average is calculated
	traceback=traceback
			scalar; type: any; default:0
			    trace back distance (AU)
			if set and nonzero then the in situ data are traced back to a distance of
			'traceback' AU, i.e. the time and heliographic longitude reflect the time of
			origin at the traceback distance.
			This option explicitly assumes that the distances in the in situ files are
			in AU, the velocities in km/s, and the heliographic longitudes in degrees.
			If this option is used date_request MUST contain 'lng','dis','vel'
 OUTPUTS:
	Result		scalar; type: integer
			0: if something went wrong
			1: time series succesfully set up
			(first check this status indicator before using the output arrays!!!)
	t_given 	array[*]; type: time structure
			    times for output time series
			    If delt is set to a positive value and t_request contains more than two
			    elements then t_given=t_request and the time series contains data
			    averaged over delt days.
			    If delt is zero, or t_request contains only 2 elements then
			    t_given contains all times for which spacecraft data were found
			    (averaged over delt days if delt > zero).
			    (if /traceback is used then the time range is shifted relative to
			    the input range)
	data_given	array[*,n_elements(data_request)], type: float
			    data_given[*,i] are the data for the quantity indicated by data_request[i]
			    missing data are marked by the value !values.f_nan
 OPTIONAL OUTPUT PARAMETERS:
	source=source	scalar; type: integer
			    if not set on input, the return value is set to whatever default
			    is set by the function Instrument (currently /somni).
 INCLUDE:
	@compile_opt.pro	; On error, return to caller
 CALLS: ***
	BadValue, CvSky, FILEPATH, InitVar, Instrument, IsTime, IsType, TimeFixYear, TimeGet
	TimeLimits, TimeOp, TimeSet, TimeUnit, UNIQ, big_body, big_eph, boost, destroyvar
	flt_read, hide_env, jpl_body, say
 CALLED BY:
	vu_insitu
 RESTRICTIONS:
	Needs system variables: !sun.spiral and !sun.au
	The input time array should be in chronological order
 PROCEDURE:
 >	Input files for density and velocity are hourly-averages. The filename is:
	    fil = filepath(root=getenv('SSW_SMEI_DAT'),sub='insitu','sw'+spacecraft+'_'+strcompress(year,/rem)+'.hravg')
	    e.g. file = 'imp8_2000.hravg'
 >	The files are assumed to be chronologically ordered.
 >	The file structure is described internally as column numbers for all quantities.
 > Heliographic longitudes are returned in monotonic decreasing order by adding a multiple of 360
	degrees when the timeseries crosses from one Carrington rotation to the next.
 MODIFICATION HISTORY:
	AUG-1998, Paul Hick (UCSD/CASS)
	JUL-1999, Paul Hick (UCSD/CASS); documentation added
	SEP-2001, Paul Hick (UCSD/CASS)
	    Added GSE & GSM to RTN conversion; currently this is only used
	    for the real-time ACE data.
	JUL-2002, Paul Hick (UCSD/CASS)
	    Fixed a bug with the delt keyword. When delt was input as a
	    difference time structure it was converted to hours, instead of days.
	DEC-2003, Paul Hick (UCSD/CASS)
	    Changed /omni to /somni (accesses the files with B in RTN)
	    Added /eomni (accesses the files with B in GSM)
	    The insitu data files can now also be located in a subdirectory
	    of $SSW_SMEI_DATA/insitu, and the filenames can be more general
	    than the default template described in PROCEDURE.
	SEP-2006, Paul Hick (UCSD/CASS)
	    Added Dana Crider Mars Orbiter data
	OCT-2006, Paul Hick (UCSD/CASS)
	    All ephemeris calculations now done by big_eph.
	    Added 'vlng','vlat' to pull the off-radial flow angles of the solar
	    wind speed out of the OMNI databases. Note that ONMI2 (EOMNI) and
	    modified OMNI (SOMNI) use different coordinates systems and different
	    definitions for the flow angles.
	FEB-2007, Paul Hick (UCSD/CASS)
	    Replaced "from" keyword by "source" keyword. "from" is still accepted
	    but is now considered obsolete.
	AUG-2007, Paul Hick (UCSD/CASS)
	    Input keyword can now also be a string (instead of just an integer)
	    In addition made the order of the various spacecraft independent
	    of the order in which function Instruments sets them up.
	    Removed keyword 'from'.
	SEP-2007, Paul Hick (UCSD/CASS)
	    Added instrument vhm to plot Ulysses magnetic field data.
	    Added instrument windb to plot Wind magnetic field data
	    Rename instrument wind to windsw (plots n,v data)
	FEB-2012, Paul Hick (UCSD/CASS)
	    Bugfix. For magnetic field data requiring a conversion from GSM
	    or GSE, the conversion would be silently skipped if the data request
	    did not include all three magnetic field components, br, bt, and bn.
	    There now is a check which results in an empty data array (and
	    return value 0) if only one or two magnetic field component are
	    specified.
	JAN-2013, Paul Hick (UCSD/CASS)
	    Added instrument 'rcelias' (real-time SOHO/Celias data)
	    Bugfix. Test for presence of doy instead of month,day of month
	    did not always work.
	JUN-2013, Paul Hick (UCSD/CASS; pphick@ucsd.edu)
	    Bumped 'atleast' value for rcelias from 5 to 6 after
	    minor change in header lines of insitu data file.

 NAME:
	InsituTimeSeries
 PURPOSE:
	Extract in-situ data in specified interval, and average if necessary
 CALLING SEQUENCE:
	FUNCTION InsituTimeseries0,$
	    data_request	, $
	    t_request		, $
	    t_given		, $
	    data_given		, $
	    delt    = delt	, $
	    weight  = weight	, $
	    traceback=traceback , $
	    source  = source_	, $
	    nocvsky = nocvsky	, $
	    silent  = silent
 INPUTS:
	data_request	array, type: string
			    determines which data are extracted.
			    Identifiers are:
			    'lng','lat','dis','den','vel','br','bt','bn','bb','pdyn','vlng','vlat'
			    'lng'   = heliographic longitude (deg)
			    'lat'   = heliographic latitude (deg)
			    'dis'   = heliocentric distance (AU)
			    'den'   = density (cm^-3)
			    'vel'   = velocity (km/s)
			    'br'    = radial magnetic field (nT)
			    'bt'    = tangential magnetic field (nT)
			    'bn'    = normal magnetic field (nT)
			    'bb'    = field strength (nT)
			    'pdyn'  = dynamic pressure (nPa)
			    'vlng'  = azimuthal flow angle
			    'vlat'  = elevation flow angle
	t_request	array[2],array[N] type: standard time structure
			    If array[2]: data are extracted between t_request[0] and t_request[1]
				at times extracted from s/c insitu files
			    If array[N]: data are extracted at times t_request
			    Should be in chronological order
 OPTIONAL INPUT PARAMETERS:
	source=source	scalar; type: integer or string
			    identifies insitu instrument either by name or integer
			    as provided by function Instrument.
	delt=delt	scalar; type: time difference structure, or float; default: 0
			    if present, and non-zero, time window used for averaging.
			    in-situ data inside [t-delt/2,t+delt/2] are averaged.
			    If a scalar (float or integer) is specified, then it is
			    assumed to be the time window in days.
	/weight 	if set, standard deviations are used to weight the averages (if available
			    for selected instrument and data)
			    if not set, an unweighted average is calculated
	traceback=traceback
			scalar; type: any; default:0
			    trace back distance (AU)
			if set and nonzero then the in situ data are traced back to a distance of
			'traceback' AU, i.e. the time and heliographic longitude reflect the time of
			origin at the traceback distance.
			This option explicitly assumes that the distances in the in situ files are
			in AU, the velocities in km/s, and the heliographic longitudes in degrees.
			If this option is used date_request MUST contain 'lng','dis','vel'
 OUTPUTS:
	Result		scalar; type: integer
			0: if something went wrong
			1: time series succesfully set up
			(first check this status indicator before using the output arrays!!!)
	t_given 	array[*]; type: time structure
			    times for output time series
			    If delt is set to a positive value and t_request contains more than two
			    elements then t_given=t_request and the time series contains data
			    averaged over delt days.
			    If delt is zero, or t_request contains only 2 elements then
			    t_given contains all times for which spacecraft data were found
			    (averaged over delt days if delt > zero).
			    (if /traceback is used then the time range is shifted relative to
			    the input range)
	data_given	array[*,n_elements(data_request)], type: float
			    data_given[*,i] are the data for the quantity indicated by data_request[i]
			    missing data are marked by the value !values.f_nan
 OPTIONAL OUTPUT PARAMETERS:
	source=source	scalar; type: integer
			    if not set on input, the return value is set to whatever default
			    is set by the function Instrument (currently /somni).
 INCLUDE:
	@compile_opt.pro	; On error, return to caller
 CALLS: ***
	BadValue, CvSky, FILEPATH, INSITUTIMESERIES0, InitVar, Instrument, IsTime, IsType
	TimeFixYear, TimeGet, TimeLimits, TimeOp, TimeSet, TimeUnit, UNIQ, big_body, big_eph
	boost, destroyvar, flt_read, hide_env, jpl_body, say
 CALLED BY:
	vu_insitu
 RESTRICTIONS:
	Needs system variables: !sun.spiral and !sun.au
	The input time array should be in chronological order
 PROCEDURE:
 >	Input files for density and velocity are hourly-averages. The filename is:
	    fil = filepath(root=getenv('SSW_SMEI_DAT'),sub='insitu','sw'+spacecraft+'_'+strcompress(year,/rem)+'.hravg')
	    e.g. file = 'imp8_2000.hravg'
 >	The files are assumed to be chronologically ordered.
 >	The file structure is described internally as column numbers for all quantities.
 > Heliographic longitudes are returned in monotonic decreasing order by adding a multiple of 360
	degrees when the timeseries crosses from one Carrington rotation to the next.
 MODIFICATION HISTORY:
	AUG-1998, Paul Hick (UCSD/CASS)
	JUL-1999, Paul Hick (UCSD/CASS); documentation added
	SEP-2001, Paul Hick (UCSD/CASS)
	    Added GSE & GSM to RTN conversion; currently this is only used
	    for the real-time ACE data.
	JUL-2002, Paul Hick (UCSD/CASS)
	    Fixed a bug with the delt keyword. When delt was input as a
	    difference time structure it was converted to hours, instead of days.
	DEC-2003, Paul Hick (UCSD/CASS)
	    Changed /omni to /somni (accesses the files with B in RTN)
	    Added /eomni (accesses the files with B in GSM)
	    The insitu data files can now also be located in a subdirectory
	    of $SSW_SMEI_DATA/insitu, and the filenames can be more general
	    than the default template described in PROCEDURE.
	SEP-2006, Paul Hick (UCSD/CASS)
	    Added Dana Crider Mars Orbiter data
	OCT-2006, Paul Hick (UCSD/CASS)
	    All ephemeris calculations now done by big_eph.
	    Added 'vlng','vlat' to pull the off-radial flow angles of the solar
	    wind speed out of the OMNI databases. Note that ONMI2 (EOMNI) and
	    modified OMNI (SOMNI) use different coordinates systems and different
	    definitions for the flow angles.
	FEB-2007, Paul Hick (UCSD/CASS)
	    Replaced "from" keyword by "source" keyword. "from" is still accepted
	    but is now considered obsolete.
	AUG-2007, Paul Hick (UCSD/CASS)
	    Input keyword can now also be a string (instead of just an integer)
	    In addition made the order of the various spacecraft independent
	    of the order in which function Instruments sets them up.
	    Removed keyword 'from'.
	SEP-2007, Paul Hick (UCSD/CASS)
	    Added instrument vhm to plot Ulysses magnetic field data.
	    Added instrument windb to plot Wind magnetic field data
	    Rename instrument wind to windsw (plots n,v data)
	FEB-2012, Paul Hick (UCSD/CASS)
	    Bugfix. For magnetic field data requiring a conversion from GSM
	    or GSE, the conversion would be silently skipped if the data request
	    did not include all three magnetic field components, br, bt, and bn.
	    There now is a check which results in an empty data array (and
	    return value 0) if only one or two magnetic field component are
	    specified.
	JAN-2013, Paul Hick (UCSD/CASS)
	    Added instrument 'rcelias' (real-time SOHO/Celias data)
	    Bugfix. Test for presence of doy instead of month,day of month
	    did not always work.
	JUN-2013, Paul Hick (UCSD/CASS; pphick@ucsd.edu)
	    Bumped 'atleast' value for rcelias from 5 to 6 after
	    minor change in header lines of insitu data file.

 NAME:
	Instrument
 PURPOSE:
	Provides a means for consistently dealing with all the different instrument for which
	in situ observations are available
 CATEGORY:
	Auxilliary
 CALLING SEQUENCE:
	FUNCTION Instrument, which, $
	    imp8    = imp8	, $
	    celias  = celias	, $
	    rcelias = rcelias	, $
	    windsw  = windsw	, $
	    windb   = windb	, $
	    acesw   = acesw	, $
	    swace   = swace	, $
	    aceb    = aceb	, $
	    somni   = somni	, $
	    eomni   = eomni	, $
	    helios1 = helios1	, $
	    helios2 = helios2	, $
	    stereoa = stereoa	, $
	    stereob = stereob	, $
	    swoops  = swoops	, $
	    vhm     = vhm	, $
	    mgs     = mgs	, $

	    name    = name	, $
	    label   = label	, $

	    nearearth = nearearth,$
	    L1	    = L1	, $
	    list    = list
	Result = Instrument( which [, /name, /label, /nearearth )
	Result = Instrument( /imp8 [, /name, /label, /nearearth )
 INPUTS:
	which	    scalar; string or integer
			integer ID of primary name of instrument
 OPTIONAL INPUT PARAMETERS:
	The following instrument keywords are available
	(instead of these keywords the argument 'which' can be set to 'imp8','celias', etc.)
	/imp8
	/celias
	/windsw
	/windb
	/acesw
	/aceb
	/somni
	/helios1
	/helios2
	/stereoa	Stereo Ahead
	/stereob	Stereo Behind
	/swoops 	Ulysses plasma data
	/vhm		Ulysses magnetic field data
	/mgs		Mars global surveyor

	/name		retrieves primary name
	/label		retrieves alternative name (usually same as primary name)
	/nearearth	identifies instrument as near-earth or deep-space
	/L1		identifies instrument as at L1.
 OUTPUTS:
	Result		scalar; integer, byte or string
			/name set: string containing the primary name
			/label set: string containing the alternative name
			/nearearth set: 1B if instrument is a nearearth instrument;
			    0B if it is a deep-space instrument
			no keyword set: integer identifier for the instrument
 INCLUDE:
	@compile_opt.pro	; On error, return to caller
 CALLED BY:
	InSitu, InsituTimeSeries [1], InsituTimeSeries [2], vu_insitu, vu_movie
 PROCEDURE:
	Each of the in situ instruments is associated with an integer identifier.
	Currently listed by name are
	    ['imp8','celias','windsw','somni','helios1','helios2','stereoa','stereob','swoops','acesw','aceb']
 MODIFICATION HISTORY:
	FEB-2000, Paul Hick (UCSD/CASS)
	FEB-2000, Paul Hick (UCSD/CASS)
	    Changed /omni to /somni; added /eomni
	FEB-2000, Paul Hick (UCSD/CASS)
	    Added Mars Global Surveyor
	AUG-2007, John Clover (UCSD/CASS; jclover@ucsd.edu)
	    Added STEREO A/B
	AUG-2007, Paul Hick (UCSD/CASS)
	    Added keyword list to pull out list of all instruments.
	SEP-2007, Paul Hick (UCSD/CASS)
	    Added instrument 'vhm'
	JAN-2013, Paul Hick (UCSD/CASS; pphick@ucsd.edu)
	    Added instrument 'rcelias' for real-time Celias data
	    Also marked both types of Celias data as 'nearL1'

 NAME:
	IntegrateLOS
 PURPOSE:
	Calculate a sky map by integrating through a 3D heliospheric matrix.
	Two forms of integrated expressions are implemented:
	(1) A weighted integration (used for g-level and IPS velocity):
	    F = Integral(F*W ds)/Integral(W ds)
	(2) Straight integration (used for Thomson scattering):
	    F = Integral(F ds)
	See PROCEDURE for more details.
 CATEGORY:
	sat/idl/util
 CALLING SEQUENCE:
	FUNCTION IntegrateLOS,	R, FF, WW,	$

	    ; Keywords passed to InterpolateHeliosphere

	    R3D     = R3D	, $
	    dR3D    = dR3D	, $
	    xcrange = xcrange	, $
	    cv2carrington=cv2carrington,$
	    xlrange = xlrange	, $
	    opengrid= opengrid	, $
	    xcgrid  = xcgrid	, $
	    xlgrid  = xlgrid	, $
	    rrgrid  = rrgrid	, $
	    fillbad = fillbad	, $

	    degrees = degrees	, $
	    minRR   = minRR	, $

	    f3dfnc  = f3dfnc	, $
	    f3darg  = f3darg	, $
	    w3dfnc  = w3dfnc	, $
	    w3darg  = w3darg	, $
	    rr_earth= rr_earth	, $
	    ut_earth= ut_earth	, $
	    pa_earth= pa_earth	, $
	    elo_earth=elo_earth , $
	    elo_sun = elo_sun	, $
	    silent  = silent
 INPUTS:
	R	array[3,*,n]; type: float
		    Locations in spherical coordinates for all points along all
		    lines of sight. '*' can be absent or represent more than one dimension.
		    E.g. for a typical skymap there will be two dimensions l,m
		    representing longitude and latitude. 'n' represents the
		    'depth' dimension along the lines of sight.

		    R[0,*,n] contains longitude angles in [0,360]
		    R[1,*,n] contains latitudes in [-90,90]
		    R[2,*,n] contains heliocentric distances in AU

	The last dimension is the 'depth' dimension along the line of sight over
	which the integration is performed.

	'Densities' and 'Weights' can be specified in two ways.

	If R3D and xcgrid are NOT specified then FF and WW directly specify densities
	and weights along all lines of sight:

	FF	    array[*,n]; type: float
			volume data at locations 'R'.
	WW	    array[*,n]; type: float
			weights at locations 'R'.

	If R3D or xcgrid are specified than calls to InterpolateHeliosphere
	are made with arrays FF and WW to determine the function values and
	weight along each line of sight.

	In this case the following keywords are passed to InterpolateHeliosphere
	    R
	    FF or WW
	    R3D     = R3D
	    dR3D    = dR3D
	    xcrange = xcrange
	    cv2carrington=cv2carrington
	    xlrange = xlrange
	    opengrid= opengrid
	    xcgrid  = xcgrid
	    xlgrid  = xlgrid
	    rrgrid  = rrgrid
	    fillbad = fillbad

	See InterpolateHeliosphere for more information.
 CALLED BY:
	RemoteView_Display2D, vu_coronagraph, vu_earthskymap, vu_elotime, vu_lineofsight
	NOTE: R[0,*] MUST BE longitude angles, even though InterpolateHeliosphere
	    in principle allows both R[0,*] and xcrange to be Carrington variables.
	    R[0,*] must be a longitude angle because R is also passed as
	    argument to f3dfnc and w3dfnc, where it is used in combination
	    with the longitude angles in rr_earth

 OPTIONAL INPUT PARAMETERS:
	/degrees    if set then all angles are in degrees (default: radians)

	The integrands are set up using a number of additional keywords:

	ut_earth=ut_earth
	rr_earth=rr_earth
		position Earth in same coordinates as R.
	pa_earth=pa_earth
		array[*,n]; type: float
		    line of sight position angle measured counterclockwise
		    from ecliptic north
	elo_earth=elo_earth
		array[*,n]; type: float
		    angle Sun-observer-segment on line of sight (i.e. the elongation)
		    for all line of sight segments (note that array[*,0] = array[*,1] etc.
	elo_sun=elo_sun
		array[*,n]; type: float
		    angle observer-Sun-segment on line of sight
 OUTPUTS:
	Result	array[*]; type: float
		    sky map with line-of-sight integrated quantities
 INCLUDE:
	@compile_opt.pro		; On error, return to caller
 CALLS: ***
	BadValue, InitVar, InterpolateHeliosphere, IsType
 PROCEDURE:
 > Two types of expressions are evaluated

	(1) <F> = Integral[ w3dfnc(s,W(s))*f3dfnc(s,F(s)) ds]/ Integral[ w3dfnc(s,W(s)) ds]
	and
	(2) <F> = Integral[ w3dfnc(s,W(s))*f3dfnc(s,F(s)) ds]

	The line of sight integrations is done for all lines of sight specified.
	The integration dimension always is the last dimension in the input arrays.
	The integration is performed simply by summing over the integration dimension,
	implicitly assuming a constant stepsize ds:

	(1) <F> = Sum[ w3dfnc(s,W(s))*f3dfnc(s,F(s)) ]/ Sum[ w3dfnc(s,W(s)) ]
	and
	(2) <F> = Sum[ w3dfnc(s,W(s))*f3dfnc(s,F(s)) ]

	Note that for expression (1) the stepsize ds will cancel. For expression (2)
	it is the users responsibility to absorb the stepsize ds into one of the
	input arrays F or W.

 >	W(s) and F(s) represent two physical quantities known on a heliographic grid
	(for example solar wind velocity and density, and the s-dependence reflects some
	geometrical aspect of the integration.

	In all applications I can think of the dependence of f3dfnc on F(s) will be linear:
	f3dfnc(s,t) = f3dfnc(s)*t. In that case the integration looks like:
	    <F> = Integral[ w3dfnc(s,W(s))*f3dfnc(s)*F(s) ds ]/ Int[ w3dfnc(s,W(s)) ds ]
	This form applies to the calculation of IPS velocities, where the weight factor
	depends on the density: w3dfnc is 'IPS_velocity_w', W(s) = n(s), and
	the perpendicular velocity is integrated: f3dfnc(s) = sin(angle(s)) (is
	'IPS_velocity') and F(s) = V(s).

	In its simplest form w3dfnc(s,t) = w3dfnc(s) and f3dfnc(s,t) = t. In that case:
	    <F> = Integral[ w3dfnc(s)*F(s) ds ]/ Int[ w3dfnc(s) ds ]
	i.e. <F> is a weighted average of F(s) along the line of sight. This form applies
	to the calculation of IPS g-level with a weight factor that depends only on the
	line of sight geometry: w3dfnc is 'IPS_hlevel_w', and the quantity F(s) is
	the 'gamma function' (introduced in the tomography manifesto).

	A special case occurs when there is no denominator:
	    <F> = Integral[ w3dfnc(s)*F(s) ds ]
	This applies to the Thomson scattering brightness where w3dfnc is the scattered
	intensity from a single electron ('ThomsonBrightness'), and F(s) is the
	solar wind density n(s). The stepsize ds is absorbed into ThomsonBrightness.

 >	Bad values in the input array F3D are detected with the IDL 'finite' function.
	This is done InterpolateHeliosphere. The treatment of these bad values depends on
	the setting of 'fillbad'.

 MODIFICATION HISTORY:
	AUG-1999, Paul Hick (UCSD/CASS)
	AUG-2004, Paul Hick (UCSD/CASS)
	    Bug fix: w3dfnc was not be used if WW did not exist. This is wrong
	    when w3dfnc only uses R and not WW (as for the g-level determination).
	    Bug fix: if neither w3dfnc nor WW is defined (as for Thomson
	    scattering the average rather than the total along the line of
	    sight was calculated.
	    Also rearranged some of the arrangement. Now the presence of R3D
	    triggers a call to InterpolateHeliosphere.
	JUN-2006, Paul Hick (UCSD/CASS)
	    Fixed bug in the summation expressions. It was explicitly set to
	    the third dimension. Now it is set to the last dimension present.
	    in F and W. Added grid keywords xcgrid, xlgrid and rrgrid to call to
	    InterpolateHeliosphere.
	FEB-2007, Paul Hick (UCSD/CASS; pphick@ucsd.edu)
	    Added test for existence of xcgrid to trigger call to
	    InterpolateHeliosphere.

 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: ***
	BadValue, Carrington, InitVar, IsType, MEAN, REVERSE, ToRadians, boost, gridfill
 CALLED BY:
	IntegrateLOS, qvu_draw, vu_atlocation, vu_timeseries
 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).

 NAME:
	IPS_hlevel_w
 PURPOSE:
	Calculates weights for integrating IPS g-level
 CATEGORY:
	sat/idl/toolbox
 CALLING SEQUENCE:
	FUNCTION IPS_hlevel_w, rr_earth, R, F, w3darg,	$
	    ut_earth	= ut_earth  , $
	    pa_earth	= pa_earth  , $
	    elo_earth	= elo_earth , $
	    elo_sun	= elo_sun   , $
	    degrees	= degrees
 INPUTS:
	rr_earth    array[3]		for 'plain' sky map
		    array[3,N]		for 'transit' sky map
					location of Earth for all lines of sight
					in same spherical coordinates as R.
	R	    array[3,N,L,M]	(N,L may be 1 or absent)
					locations of all segments for all lines of sight
					in same spherical coordinates as rr_earth
					(usually heliographic)
					This can be a grid of NxL lines of sight with M
					segments along each line of sight
					R[0,*,*,*] : longitude
					R[1,*,*,*] : latitude
					R[2,*,*,*] : heliocentric distance
	F				not used
	w3darg	    array[4]		w3darg[0] = observing freq (MHz)
					w3darg[1] = source size (arcsec)
					w3darg[2] = extra source size scale factor
					w3darg[3] = beta_r
 OPTIONAL INPUT PARAMETERS:
	/degrees			if set all input angles should be in degrees (default: radians)
	ut_earth			not used
	pa_earth			not used
	elo_earth			not used
	elo_sun 			not used
 OUTPUTS:
	R	    array[N,L,M]	weight factors
 INCLUDE:
	@compile_opt.pro	    ; On error, return to caller
 EXTERNAL BY:
	vu_coronagraph, vu_earthskymap, vu_lineofsight
 CALLS: ***
	CV_COORD, IPS_WeightFnc, SubArray, SuperArray
 CALLED BY:
	IPS_velocity_w
 RESTRICTIONS:
	For each line of sight segment the weight is the product of a factor
	which is a function heliocentric distance and a factor which is a
	function of the geocentric distance. The second factor is calculated
	for one line of sight only on the assumption that it's the same for
	all lines of sight.
 PROCEDURE:
 MODIFICATION HISTORY:
	SEP-1999, Paul Hick (UCSD/CASS)
	FEB-2007, Paul Hick (UCSD/CASS; pphick@ucsd.edu
	    Modified to work with list of los specified in any configuration
	    (not just as a 2D skymap)

 NAME:
	IPS_params
 PURPOSE:
	Sets parameters for calculating IPS velocities and g-levels
 CATEGORY:
	Physics
 CALLING SEQUENCE:
	R = IPS_params( /nagoya )
	R = IPS_params( freqmhz=FreqMHz, arcsec=ArcSec, sourcescale=SourceScale)
 OPTIONAL INPUT PARAMETERS:
	/nagoya 	if set, returns [327 , 0.1, 1.0]
	/cambridge	if set, returns [81.5, 0.3, 1.0]
	/sandiego	if set, returns [73.8, 0.3, 1.0]
 OUTPUTS:
	R	    array[3]; type: float; default: [327, 0.1, 1.0] (Nagoya values)
			R[0]: observing frequency in MHz
			R[1]: source size (arcsec)
			R[2]: source scale (multiplicative factor to be applied to source size)
 CALLED BY:
	IPS_WeightFnc, RemoteView_Display2D, vu_coronagraph, vu_earthskymap, vu_elotime
	vu_lineofsight
 MODIFICATION HISTORY:
	JUN-2000, Paul Hick (UCSD/CASS; pphick@ucsd.edu)

 NAME:
	IPS_velocity
 PURPOSE:
	Calculates component of solar wind velocity perpendicular to line of sight
 CATEGORY:
	sat/idl/toolbox
 CALLING SEQUENCE:
	FUNCTION IPS_velocity, rr_earth, R, F, w3darg,	$
	    ut_earth	= ut_earth  , $
	    pa_earth	= pa_earth  , $
	    elo_earth	= elo_earth , $
	    elo_sun	= elo_sun   , $
	    degrees	= degrees
 INPUTS:
	rr_earth			not used
	R				not used
	F	    array[N,L,M]	solar wind velocity matrix
	w3darg				not used
 OPTIONAL INPUT PARAMETERS:
	/degrees			if set all input angles are in degrees
	ut_earth			not used
	pa_earth			not used
	elo_earth			angle Sun-Earth-LOS segment
	elo_sun 			angle Sun-LOS segment-Earth
 OUTPUTS:
	R	    array[N,L,M]	weight factors
 INCLUDE:
	@compile_opt.pro	; On error, return to caller
 EXTERNAL BY:
	vu_coronagraph, vu_earthskymap, vu_lineofsight
 CALLS: ***
	ToRadians
 MODIFICATION HISTORY:
	SEP-1999, Paul Hick (UCSD/CASS; pphick@ucsd.edu)

 NAME:
	IPS_velocity_w
 PURPOSE:
	Calculates weights for integrating IPS g-level
 CATEGORY:
	Physics: IPS
 CALLING SEQUENCE:
	FUNCTION IPS_velocity_w, rr_earth, R, F, w3darg,    $
	    ut_earth	= ut_earth  , $
	    pa_earth	= pa_earth  , $
	    elo_earth	= elo_earth , $
	    elo_sun	= elo_sun   , $
	    degrees	= degrees
 INPUTS:
	rr_earth    array[3]		for 'plain' sky map
		    array[3,N]		for 'transit' sky map
					location of Earth for all lines of sight
					in same spherical coordinates as R.
	R	    array[3,N,L,M]	(N,L may be 1 or absent)
					locations of all segments for all lines of sight
					in same spherical coordinates as rr_earth
					(usually heliographic)
					This can be a grid of NxL lines of sight with M
					segments along each line of sight
					R[0,*,*,*] : longitude
					R[1,*,*,*] : latitude
					R[2,*,*,*] : heliocentric distance
	F	    array[N,L,M]	(N,L may be 1 or absent)
					g-values at all los segments
	w3darg	    array[4]		w3darg[0] = observing freq (MHz)
					w3darg[1] = source size (arcsec)
					w3darg[2] = extra source size scale factor
					w3darg[3] = betar
 OPTIONAL INPUT PARAMETERS:
	/degrees			if set all input angles should be in degrees (default: radians)
	ut_earth			not used
	pa_earth			not used
	elo_earth			not used
	elo_sun 			not used
 OUTPUTS:
	R	    array[N,L,M]	weight factors
 INCLUDE:
	@compile_opt.pro	; On error, return to caller
 EXTERNAL BY:
	vu_coronagraph, vu_earthskymap, vu_lineofsight
 CALLS: ***
	IPS_hlevel_w
 RESTRICTIONS:
	For each line of sight segment the weight is the product of a
	factor which is a function heliocentric distance and a factor
	which is a function of the geocentric distance. The second
	factor is calculated for one line of sight only on the assumption
	that it's the same for all lines of sight.
	These weights are calculated by ips_hlevel_w, and are here
	multiplied by the F array.
 MODIFICATION HISTORY:
	SEP-1999, Paul Hick (UCSD/CASS; pphick@ucsd.edu)

 NAME:
	IPS_WeightFnc
 PURPOSE:
	Calculates weight function along the line of sight for IPS
 CATEGORY:
	Physics
 CALLING SEQUENCE:
	W = IPS_WeightFnc(ZAu, params)
 INPUTS:
	ZAu	    array[N]; type float
			distance from Earth along line of sight (AU)
	params	    array[3]; type: float
			params[0]: observing freq in MHz
			params[1]: source size (arcsec)
			params[2]: extra scale factor for source size
 OUTPUTS:
	W	    array[N]; type: float
			weight factors
 CALLS: ***
	IPS_params, SPECTRALFNC, nrSimpson
 CALLED BY:
	IPS_hlevel_w
 PROCEDURE:
	CALCULATE INTEGRAL FOR WEIGHTING FUNCTION OF Z.
	SIMPSON METHOD DOUBLE INTEGRAL
	SpectralFnc = q*fresnel(q)*size(q)*spectr(q)
	fresnel= fresnel filter term
	size   = source size term
	spectr = spectrum fluctuation term
 MODIFICATION HISTORY:
	MAR-1999, Paul Hick (pphick@ucsd.edu)
	    Adapted from STELab Fortran code