function Sky_VolumeColors, rgbo,			$
	rgb_range		= rgb_range,			$

	full_range		= full_range,			$
	full_color		= full_color,			$
	full_opacity	= full_opacity,			$

	mid_range		= mid_range,			$
	mid_color		= mid_color,			$
	mid_opacity		= mid_opacity,			$

	opacity_power	= opacity_power,		$
	update			= update
;+
; NAME:
;	Sky_VolumeColors
; PURPOSE:
;	Sets up the RGBO vector for volume rendering
; CATEGORY:
;	Volume rendering: SkyImage
; CALLING SEQUENCE:
;	rgbo = Sky_VolumeColors(rgbo,			$
;			rgb_range		= rgb_range,			$
;
;			full_range		= full_range,			$
;			full_color		= full_color,			$
;			full_opacity	= full_opacity,			$
;
;			mid_range		= mid_range,			$
;			mid_color		= mid_color,			$
;			mid_opacity		= mid_opacity,			$
;
;			opacity_power	= opacity_power,		$
;			update			= update
; INPUTS:
;	full_range
;		scalar or array[2]; type: any; no default (MUST be specified)
;			A scalar is interpreted as [0,full_range]
;			The full range of values of rendered data in data units
;			(e.g. cm^-3 for a rendering of densities)
;			This range is converted to the byte range 'rgb_range' (using bytscl).
; OPTIONAL INPUT PARAMETERS:
;	rgbo
;		array[256,4]; type: byte
;			rgbo vector to be updated
;
;	/update		if set then the input rgbo vector is updated
;				(ignored if rgbo is not defined)
;
;	rgb_range	= rgb_range
;		scalar or array[2]; type: byte; default: [0B,255B]
;			Range of elements in rgbo assigned to the data set.
;
;	mid_range	= mid_range
;		scalar or array[2]; type: any; default: full_range
;			A scalar is interpreted as [0,full_range]
;			Internally this is restricted to be a subrange of 'full_range'
;
;			Range of data values defining 'intermediate data values'
;
;	full_color	= full_color
;		scalar or array[2]; default: [0,235]
;			A scalar is interpreted as [0,full_color].
;
;			Range of colors to be used in the rendering.
;			Internally this is restricted to the available colors
;			(currently 236).
;			color_range[0] is matched to data value full_range[0]
;			color_range[1] is matched to data value full_range[1]
;
;	mid_color	= mid_color
;		scalar or array[2]; default: full_color
;			Internally this is restricted to be a subrange of 'full_color'
;
;			Range of colors used for rendering the intermediate data values
;
;	full_opacity	= full_opacity
;		scalar or array[2]; type: float; default: [0,1]
;			range of opacity values
;
;	mid_opacity	= mid_opacity
;		scalar or array[2]; type: float; default: full_opacity
;			Internally this is restricted to be a subrange of 'full_opacity'
;
;			Opacity applied to data values to mid-range data values
;
;	opacity_power = opacity_power
;		scalar; type: float; default: 1
;			power applied to opacity for mid-range data values
;
; OUTPUTS:
;	rgbo		array[256,4]; type: bytarr
;					rgbo color vector needed for volume rendering
;					The first index [0,255] represent the data values.
;					Only the range 'rgb_range' is filled. The rest is zero
;					(if rgbo is created), or unmodified (if rgbo is updated)
; CALLS:
;	flt_read, filepath, SuperArray, gridgen, IsType, SetRange
; INCLUDE:
;	@skycolors.pro
; SIDE EFFECTS:
;	Needs files $SYS:skyinfo.txt and $SYS:skyrgb.txt, storing the color table information
; RESTRICTIONS:
; PROCEDURE:
;	The original 3D matrix will cover a certain range of data values. Volume rendering
;	can only handle a maximum of 256 different 'density values'.
;
;	We first need to decide how many of these we want to assign to the 3D density matrix.
;	If only one 3D volume data set is rendered then the full range [0,255] should be used,
;	but when several 3D volume data sets are merged each set is assigned its own subrange.
;	This is done with keyword 'rgb_range'
;
;	Next we need to set the range of data values in the 3D matrix to be rendered.
;	'full_range' sets the data range to be scaled to 'rgb_range'.
;	The 3D matrix will be converted to byte by the IDL bytscl command (in href=Sky_ColorView=)
;
;		rgb_range[0]+bytscl( Matrix, top=rgb_range[1]-rgb_range[0], min=full_range[0], max=full_range[1])
;
;	This scales the data values to the proper range 'rgb_range'
;	Note that data values below full_range[0] are treated the same as full_range[0], and data
;	values above full_range[1] the same as full_range[1]).
;
;	'mid_range' selects a subrange of 'full_range'. These two values separate the full range
;	in three sections:
;		low data values:		[full_range[0], mid_range[0]]
;		intermediate values:	[mid_range [0], mid_range[1]]
;		high data values:		[midrange  [1], full_range[0]]
;	These three ranges are treated separately.
;
;	To set up the rgbo vector a 8-bit color table and an opacity needs to be defined.
;	The color table information is stored in two files (see RESTRICTIONS).
;	Keywords 'full_color' and 'mid_color' select color ranges from this table.
;	Keywords 'full_opacity' and 'mid_opacity' set the opacities.
;
;	Data values, color indices and opacities are matched as follows:
;		full_range[0]	<----->		full_color[0]	<------>	full_opacity[0]
;		mid_range [0]	<----->		mid_color [0]	<------>	mid_opacity [0]
;		mid_range [1]	<----->		mid_color [1]	<------>	mid_opacity [1]
;		full_range[1]	<----->		full_color[1]	<------>	full_opacity[1]
;
;	For the low data values (between full_range[0] and mid_range[1]) both color index
;	and opacity are obtained by a linear scaling with the data value; the same is
;	the case for the high data values (between mid_range[1] and full_range[1]).
;
;	For the intermediate data values the color index scales linearly with the
;	data values. The opacity scales as a power law of data value with the power
;	specified in keyword opacity_power.
; MODIFICATION HISTORY:
;	JUN-2000, Paul Hick (UCSD/CASS; pphick@ucsd.edu)
;-
on_error, 2							; On error, return to caller
@skycolors.pro				; Dummy comment for Linux

refresh = 1-keyword_set(update) or IsType(rgbo, /undefined)

if n_elements(opacity_power) eq 0 then	$
	opacity_power = 2.					$
else									$
	opacity_power = opacity_power[0]

; Index range in rgbo used for this data set

rgb_index = SetRange( [0B,255B], rgb_range)

; Full data value range

case n_elements(full_range) of
0   : message, 'no density range specified as first argument'
1   : data_index = [0,full_range[0]]
else: data_index = full_range[0:1]
endcase

; Set up range for intermediate data values. Check limits agains full_range

mid_data_index = SetRange(data_index, mid_range)

; Set opacity ranges

opacity_index	  = SetRange([0,1], full_opacity)
mid_opacity_index = SetRange(opacity_index, mid_opacity)

; Read color tables from file

if not flt_read(filepath(root=getenv('SYS'),'skyinfo.txt'),tmp,crumbs=Crumbs,	$
	/silent,error=error) then message, 'No color information available: '+error

tmp = fix( tmp(*, where(Crumbs eq cRgbo  )) )	; Indices into rgb_table

if not flt_read(filepath(root=getenv('SYS'),'skyrgb.txt'),rgb_table,	$
	/silent,error=errror) then message, 'No color information available: '+error

rgb_table = transpose( byte(rgb_table) )

; Color table available for the volume rendering
; Add a 4th row to the rows for R, G and B values filled with 255B

ncols = tmp[1]-tmp[0]+1				; Max # available colors
rgb_table = [ [rgb_table(tmp[0]:tmp[1],*)],	[replicate(255B, ncols)] ]

; Set limits on color index ranges

color_index		= SetRange([0,ncols-1], full_color)
mid_color_index = SetRange(color_index, mid_color )

; Set up rgbo vector

nrgb = rgb_index[1]-rgb_index[0]+1		; # entries in rgb to be set

color	= bytarr(nrgb)					; nrgb colors
opacity = fltarr(nrgb)

; Find index range into the opacity array occupied by the intermediate densities
; (in range [0, rgb_index[1]-rgb_index[0]])

mid_data_index = float(rgb_index[1]-rgb_index[0])/(data_index[1]-data_index[0])*(mid_data_index-data_index[0])
mid_data_index = round(mid_data_index)

if mid_data_index[0] gt 0 then begin
	color  [0:mid_data_index[0]] = round( gridgen(mid_data_index[0]+1, range=[color_index  [0], mid_color_index  [0]]) )
	opacity[0:mid_data_index[0]] =	      gridgen(mid_data_index[0]+1, range=[opacity_index[0], mid_opacity_index[0]])
endif

if mid_data_index[1] lt nrgb-1 then begin
	color  [mid_data_index[1]:*] = round( gridgen(nrgb-mid_data_index[1], range=[mid_color_index  [1], color_index  [1]]) )
	opacity[mid_data_index[1]:*] =	      gridgen(nrgb-mid_data_index[1], range=[mid_opacity_index[1], opacity_index[1]])
endif

nmid = mid_data_index[1]-mid_data_index[0]+1
color  [mid_data_index[0]:mid_data_index[1]] = round( mid_color_index  [0]+(mid_color_index  [1]-mid_color_index  [0])*gridgen(nmid, /one) )
opacity[mid_data_index[0]:mid_data_index[1]] =		  mid_opacity_index[0]+(mid_opacity_index[1]-mid_opacity_index[0])*gridgen(nmid, /one)^opacity_power


if refresh then rgbo = bytarr(256,4)

rgbo[rgb_index[0]:rgb_index[1],*] = round( SuperArray(opacity,4,/trail)*rgb_table[color,*])

return, rgbo  &  end