NAME:
	Rank
 PURPOSE:
	Construct rank table from given index table.
 CATEGORY:
	Math: sorting
 CALLING SEQUENCE:
	subroutine Rank(N1,N2,M1,M2,INDX,IRANK)
 INPUTS:
	N1,N2	integer		physical dimensions of arrays INDX and IRANK
	M1,M2	integer		logical dimensions of arrays INDX and IRANK
		(only the logical part of INDX is ranked)
	INDX(M1:M2)
		integer		index table for a one-dimensional array
				ARRIN(M1:M2) (obtained from subroutine Index??)
 OUTPUTS:
	IRANK(M1:M2)
		integer		rank table
 SIDE EFFECTS:
	The part of the rank array IRANK outside the logical range
	(i.e. [M1:N1-1] and [N2+1,M2]) remains untouched
 RESTRICTIONS:
	The index range [N1,N2] must be a subset of [M1,M2]
 SEE ALSO:
	IndexI4, Sort2I4, SortI4
 PROCEDURE:
	The index values [N1,N2] are arranged in IRANK(N1:N2) in such a way
	that IRANK(J) gives the rank of the element ARRIN(J), i.e. N1 if 
	ARRIN(J) is the smallest element in the range [N1,N2]; and N2 if it is
	the largest.
	See Press et al., "Numerical Recipes", Cambridge UP (1989), Chapter 8, 
	p. 234
 MODIFICATION HISTORY:
	1990, Paul Hick (UCSD)

 NAME:
	RClean
 PURPOSE:
	Cleans 'dirty' or 'NaN' REAL*4 zeros
 CATEGORY:
	gen/for/lib
 CALLING SEQUENCE:
	subroutine Rclean (idirty, iclean)
 INPUTS:
	X	real	value to be checked
 OUTPUTS:
	X	real	set to zero if it was a dirty zero
	I	integer		0: if the was not modified
				1: if the input was set to zero
 CALLS:
	(none)
 SEE ALSO:
	Clean0, iClean0, iCleanNaN
 RESTRICTIONS:
	Should be called with REAL*4 argument!
 PROCEDURE:
	A dirty zero is tested by the condition
		iand(idirty,Z'0000FF80') .eq. 0

	The NaN value is tested by
		iand(idirty,Z'0000FF80') .eq. 0
 MODIFICATION HISTORY:
	SEP-1999, Paul Hick (UCSD/CASS; pphick@ucsd.edu)

 NAME:
	RebinSphere
 PURPOSE:
	(NOT WORKING YET)
	Fills 'holes' and/or smooths a 3D array using a Gaussian weighting function.
	Two of the array dimensions represent a regular grid in spherical
	coordinates (longitude and latitude). The third represents a regular
	grid in radial distance.
 CALLING SEQUENCE:
	subroutine RebinSphere(dXC,nLng,nLat,nR,Z,WidDeg,R1,iBadZ,WThreshold,ClipLng, mLng,mLat,mR,ZZ)
 INPUTS:
	dXC		real		range of Carrington variable covered by nLng
					i.e. number of 360 deg rotations.
	nLng		integer		# longitudes
	nLat		integer		# latitudes
	nR		integer		# radial distances
	Z(nLng,nLat,nR)	real		3D array of function values
	WidDeg		real		Width of Gaussian used for angular
					smoothing (in degrees) (WidDeg<0 has
					special meaning; see PROCEDURE)
	R1		real		(3D version only)
				R1 > 0:Radial distance for inner boundary (in
					grid spacings), i.e. location of Z(*,*,1)
					R1=0 IS INVALID; MUST BE R1>0
				R1 < 0:(value does not matter)
					suppresses the radial smoothing, i.e.
					each subarray Z(*,*,i),i=1,nR is
					2D smoothed over longitude and latitude only.
	iBadZ		integer	     0: Invalid elements remain invalid;
					valid elements are replaced by a
					smoothed value
				     1: All elements (valid and invalid) are
					replaced by smoothed values.
				     2: Invalid elements are replaced by
					smoothed values; valid elements remain
					untouched.

					The following two options are useful if all
					bad bins have to be filled in with something.
					Use them at your own risk (see PROCEDURE).

				     3: Same as iBadZ=1, but if any invalid
					elements remain at the end, these are
					all filled in by a call to GridFill
				     4: Same as iBadZ=2, but if any invalid
					elements remain at the end, these are
					all filled in by a call to GridFill

				  >=10: Use open grid for longitude and latitude (?)
				 >=100: Use lookup table for Gaussian

	WThreshold	real		The replacement by a smoothed value is
					made only if total weight is 
					larger than/equal to the threshold value.
					(see PROCEDURE).
	ClipLng		real		Longitude difference (see PROCEDURE)
					(usually set to zero).
 OUTPUTS:
	Z(nLng,nLat,nR)	real		smoothed array of function values
 INCLUDE:
	include		'math.h'
 CALLS: ***
	ArrR4Copy, BadR4, GaussLookup, GridFill, Say, iArrR4ValuePresent
 RESTRICTIONS:
 >	GridSphere3D uses internal scratch arrays, which are
	dimensioned using 3 parameters, LNG,LAT and LNR. If the scratch
	arrays are too small, program execution stops.
 >	iBadZ=3 and iBadZ=4:
	GridFill uses a 2D algorithm to fill in holes. If nR>1 then GridFill
	is applied to each level k=1,nR separately.
 PROCEDURE:
 >	Invalid or 'Bad' array elements are indicated by the value BadR4()
 >	The 1st dimension represents a regular grid in longitude covering
	[0,360] degrees, i.e. Lng(i) = 360*(i-1)/(nLng-1), i=1,nLng.
	The 2nd dimension represents a regular grid in latitude covering
	[-90,+90] degrees, i.e. Lat(j) = -90*[-1+2*(j-1)/(nLat-1)], j=1,nLat.
	The 3rd dimension represents a regular grid in radial distance,
	i.e. R(k) = R1+k-1
 >	Obviously, 2D arrays can be smoothed by setting nR=1.
 >	GridSphere2D is just a special call to GridSphere3D with R1=-1 (the
	value does not matter as long as it's negative).
 >	The weighting is Gaussian with weight=1 at the origin. To avoid loosing
	valid elements for iBadZ=0 or 1, set the threshold weight to a
	value less than 1.
 >	The angular weighting factor involves the angular distance between
	two locations on a sphere. The points are connected by two arcs,
	one is the minimum angular distance A (0<A<180), the other is the
	complementary angle AC=360-A (180<AC<360).
	If ClipLng=0 the angle A is always used. As a result, the smoothed
	function value for small/large longitudes will depend on the input
	function values at large/small longitudes with fairly large weighting
	factors (since the angular distance is small).
	This effect can be suppressed by setting ClipLng to a non-zero value.
	If the longitude difference abs(Lng(i1)-Lng(i2)) is larger than
	ClipLng, then the corresponding element is excluded.
	BE CAREFUL WHEN USING THIS OPTION. In particular, setting
	ClipLng to values around 180 degrees affects the smoothing for
	points near the poles (where points with a wide range of longitude
	differences are close together). As a guideline 
	ClipLng > 180+3*WidDeg is probably not a bad place to start.
 >	If WidDeg > 0, the smoothing procedure is speeded up by looping
	over only part of the array when calculating the smoothed function 
	value. Excluded are all elements which are more than 3 times the
	Gaussian width away from the element for which the smoothed value is
	calculated (the weight would be less than exp(-9)). The exclusion
	is based on an approximate calculation of the distance between points.
	Near the poles there may be a problem (see the explanation in the code),
	but I think it's under control now.
	In case you don't trust the result set WidDeg to a negative value.
	If WidDeg <0, then the absolute value is used as Gaussian width, and
	the whole array is used in the calculation of the smoothed average.
	The procedure will be a LOT slower though.
 >	If ClipLng=0 then the averaging algorithm is mod(XCend-XCbeg), i.e. if
	in the input array column 1 (XC = XCend) and column nLng (XC = XCbeg)
	are the same this is also true for the output array.
	HOWEVER this symmetry is broken if iBadZ=3 or 4 because of the extra
	call to GridFill.
 MODIFICATION HISTORY:
	FEB-1997, Paul Hick (UCSD/CASS; pphick@ucsd.edu)

 NAME:
	RemoveTag
 PURPOSE:
	Remove a tag (substring) from a list of tags
 CATEGORY:
	Strings: tag manipulation
 CALLING SEQUENCE:
	subroutine RemoveTag(cTags,cSep,cTag)
 INPUTS:
	cTags		character*(*)	list of tags, separated by the cSep character
	cSep		character	character used as separater (usually comma)
	cTag		character*(*)	tag (substring) to be removed
 OUTPUTS:
	cTags		character*(*)	updated list of tags
 CALLS: ***
	LocateTag, itrim
 SEE ALSO:
	ReplaceTag, SwitchTag
 PROCEDURE:
 MODIFICATION HISTORY:
	NOV-1994, Paul Hick (UCSD/CASS; pphick@ucsd.edu)
	JUN-1995, Paul Hick (UCSD; pphick@ucsd.edu), added SwitchTag

 NAME:
	ReplaceTag
 PURPOSE:
	Replace a tag (substring) in a list of tags by another tag
 CATEGORY:
	Strings: tag manipulation
 CALLING SEQUENCE:
	subroutine ReplaceTag(cTags,cSep,cOldTag,cNewTag)
 INPUTS:
	cTags		character*(*)	list of tags, separated by the cSep character
	cSep		character	character used as separater (usually comma)
	cOldTag		character*(*)	tag (substring) to be replaced
	cNewTag		character*(*)	tag inserted at location of old tag
 OUTPUTS:
	cTags		character*(*)	updated list of tags (RemoveTag, SwitchTag only)
 CALLS: ***
	LocateTag, itrim
 CALLED BY:
	SwitchTag
 SEE ALSO:
	RemoveTag, SwitchTag
 PROCEDURE:
 MODIFICATION HISTORY:
	NOV-1994, Paul Hick (UCSD/CASS; pphick@ucsd.edu)
	JUN-1995, Paul Hick (UCSD; pphick@ucsd.edu), added SwitchTag

 NAME:
	rice
 PURPOSE:
	Compress/decompress *.nic files
 CATEGORY:
	I/O
 CALLING SEQUENCE:
	program rice

	rice file1 [file2 ...] [/directory[=directory]]
	rice /help	(displays help)
			(on Linux use a '-' instead of '/')
 INPUTS:
	file1 file2	file specifications for files to be processed. Only
			files with type .nic or .ice are accepted. The file
			name can contain wildcards. Directories are optional.

			If no files are specified then the user is prompted for
			a single file specification
	/directory=directory
			output directory for compressed/decompressed files
			Specifying /directory (without an '=directory') will
			put output files in the same directory as the input
			files.

			If no directory is specified then the user is prompted
			to specify a directory. The default directory presented
			will be the directory of the first input file.
 CALLS: ***
	AskChar, ForeignArg, Int2Str, LocFirst, LocFirstLen, Say, Str2Str, bCompressNic
	bUncompressNic, cInt2Str, iCheckDirectory, iGetFileSpec, iSearch [1], iSearch [2]
	iSearch [3], iSearch [4], iSetFileSpec, lowercase
 INCLUDE:
	include		'dirspec [1].h'
	include		'dirspec [2].h'
	include		'dirspec [3].h'
	include		'dirspec [4].h'
	include		'openfile.h'	
	include 	'filparts.h'
 RESTRICTIONS:
	This program is only intended to compress/decompress *.nic files from
	the SMEI camera. If a *.ice file is decompressed which is not a
	compressed *.nic file then the decompressed file will contain garbage.
 PROCEDURE:
 >	On NT this program will usually be executed by clicking on the program
	icon, i.e. without any command line arguments. When prompted for
		File(s) to be processed
	specify a file name and type (the directory is optional).
	The name may include a wild card; the type must be .nic or .ice

 >	In NT it is possible to select up to 100 files, and 'drag and drop' them
	on the program icon (the value 100 is defined as a parameter nVar=100)

 >	When prompted for
		Output directory
	specify a destination directory for the output files. The default
	(selected by hitting return) is the same directory as the input file(s).

 >	.nic files are compressed to .ice files
 >	.ice files are decompressed to .nic files.

 >	Files implied by the file1, file2, etc. arguments are accumulated in
	an array of file names. Currently a maximum of 2000 files are processed.
 MODIFICATION HISTORY:
	DEC-2000, Kevin Nguyen, Paul Hick (UCSD; pphick@ucsd.edu)

 NAME:
	rotate
 PURPOSE:
	Calculate polar angles in rotated coordinate system
 CATEGORY:
	Math: coordinate transformation
 CALLING SEQUENCE:
	subroutine rotate(ALFA,BETA,GAMMA,PHI,RLAT)
 INPUTS:
	ALFA	real	phase angle of the pole Z(new) in old coordinate system
	BETA	real	polar angle of the pole Z(new) in old coordinate system
	GAMMA	real	phase angle of X(new) in coordinate system (2)
	PHI	real	phase angle in old coordinate system
	RLAT	real	latitude in old coordinate system (-90<=RLAT<=90)
 OUTPUTS:
	PHI	real	phase angle in new coordinate system (0<=PHI<360)
	RLAT	real	latitude in new coordinate system
 CALLS: ***
	Say, dacosd, datan2d, dcosd, dsind
 CALLED BY:
	ECLIPTIC_EQUATOR, ECLIPTIC_HELIOGRAPHIC, KeplerOrbit, PRECESSION_ROT
 PROCEDURE:
	Given phase and polar angles of a point in some coordinate system
	[axes: X(old),Y(old),Z(old)], phase and polar angles are calculated in
	the new coordinate system [axes: X(new),Y(new),Z(new)] rotated with
	respect to the original one over the Euler angles ALFA, BETA, GAMMA.
	The rotation from old to new system is realized by 
	(1) rotation around Z(old) over ALFA  ---> X(1),Y(1),Z(1)=Z(old)
	(2) rotation around Y(1)   over BETA  ---> X(2),Y(2)=Y(1),Z(2)=Z(new) 
	(3) rotation around Z(new) over GAMMA ---> X(new),Y(new),Z(new). 
	All angles are in degrees.
 MODIFICATION HISTORY:
	1989, Paul Hick (MPAE,UCSD/CASS; pphick@ucsd.edu)

 NAME:
	rotate_dcm2quat
 PURPOSE:
	Convert direction cosine matrix to quaternion
 CATEGORY:
	gen/for/lib
 CALLING SEQUENCE:
	subroutine rotate_dcm2quat(dcm,qq)
 INPUTS:
	dcm(3,3)	double precision	direction cosine matrix
 OUTPUTS:
	q(4)		double precision	quaternion
 CALLED BY:
	Time2smei_quaternion
 PROCEDURE:
 MODIFICATION HISTORY:
	JUL-2006, Paul Hick (UCSD/CASS; pphick@ucsd.edu)

 NAME:
	rotate_euler2dcm
 PURPOSE:
	Convert Euler angles to direction cosine matrix
 CATEGORY:
	gen/for/lib
 CALLING SEQUENCE:
	subroutine rotate_euler2dcm(alfa,beta,gamma,dcm)
 INPUTS:
	alfa		double precision	rotation around z-axis
	beta		double precision	rotation around y-axis
	gamma		double precision	rotation arounc z-axis
				Euler angles
 OUTPUTS:
	dcm(3,3)	double precision	direction cosine matrix
 CALLS: ***
	dcosd, dsind
 CALLED BY:
	Time2smei_quaternion
 PROCEDURE:
 MODIFICATION HISTORY:
	JUL-2006, Paul Hick (UCSD/CASS; pphick@ucsd.edu)

 NAME:
	rotated
 PURPOSE:
	Calculate polar angles in rotated coordinate system
 CATEGORY:
	Math: coordinate transformation
 CALLING SEQUENCE:
	subroutine rotated(ALFA,BETA,GAMMA,PHI,RLAT)
 INPUTS:
	ALFA	double precision	phase angle of the pole Z(new) in old coordinate system
	BETA	double precision	polar angle of the pole Z(new) in old coordinate system
	GAMMA	double precision	phase angle of X(new) in coordinate system (2)
	PHI	double precision	phase angle in old coordinate system
	RLAT	double precision	latitude in old coordinate system (-90<=RLAT<=90)
 OUTPUTS:
	PHI	double precision	phase angle in new coordinate system (0<=PHI<360)
	RLAT	double precision	latitude in new coordinate system
 CALLS: ***
	Say, dacosd, datan2d, dcosd, dsind
 CALLED BY:
	Time2EclipticEquatorial, Time2EclipticHeliographic, Time2KeplerOrbit
	Time2Precession
 PROCEDURE:
	Given phase and polar angles of a point in some coordinate system
	[axes: X(old),Y(old),Z(old)], phase and polar angles are calculated in
	the new coordinate system [axes: X(new),Y(new),Z(new)] rotated with
	respect to the original one over the Euler angles ALFA, BETA, GAMMA.
	The rotation from old to new system is realized by 
	(1) rotation around Z(old) over ALFA  ---> X(1),Y(1),Z(1)=Z(old)
	(2) rotation around Y(1)   over BETA  ---> X(2),Y(2)=Y(1),Z(2)=Z(new) 
	(3) rotation around Z(new) over GAMMA ---> X(new),Y(new),Z(new). 
	All angles are in degrees.
 MODIFICATION HISTORY:
	JAN-2004, Paul Hick (UCSD/CASS; pphick@ucsd.edu)