FUNCTION vectorproduct, r1, r2, unit=unit, cylin=cylin, sphere=sphere, rect=rect,	$
	degrees=degrees, scalar=scalar, angle=angle
;+
; NAME:
;	vectorproduct
; PURPOSE:
;	Calculate vector or scalar product for two vectors
; CATEGORY:
;	gen/idl/toolbox/math
; CALLING SEQUENCE:
;	result = vectorproduct(r1, r2 [, /unit, /cylin, /sphere, /rect, /degrees, /scalar, /angle])
; INPUTS:
;		r1 		arrays[3,*]; type: int or float
;		r1	 	arrays[3,*]; type: int or float
;					vectors for which the product is to be taken.
;					the first dimension identifies the three coordinates: Cartesian (default),
;					cylindrical or spherical.
; OPTIONAL INPUT PARAMETERS:
;		/scalar		by default the vector (cross-product) is calculated. If /scalar is set
;					then the scalar product is calculated instead.
;		/angle		if set then information about the angle between the two vectors is returned:
;						for vector products sin(angle) = |r1 x r2|/|r1||r2|
;						for scalar products cos(angle) = |r1.r2|/|r1||r2|
;		/unit		if set, then a unit vector along the cross-product vector is returned
;						(keyword not used if /scalar is set)
;		/cylin, /sphere, /rect
;					indicates whether input vectors are in cylindrical, spherical or rectangular
;						coordinates. If none is set then rectangular is assumed.
;						For cross-products the output vectors are in the same units as the input vectors.
;		/degrees	if set, and spherical or cylindrical coordinates are used, then the angles
;						are assumed to be in degrees (default: radians)
; OUTPUTS:
;		result	array[3,*] or array[*]	same type as input
;					cross-products are returned as array[3,*] (same as the input arrays)
;					scalar products or angle information are returned as array[n]
; OPTIONAL OUTPUT PARAMETERS:
; INCLUDE:
	@compile_opt.pro			; On error, return to caller
; CALLS:
;	SyncArgs, SyncDims, InitVar
; PROCEDURE:
; >	The arrays r1 and r2 do not necessarily have the same dimensions, i.e.
;	if r1 is an array[3,n] and r2 is an array[3] then r2 is interpreted as an array[3,n]
;	with all n vectors the same (SyncArgs is used to synchronize the array dimensions).
; >	Arrays r1 and r2 can have more than two dimensions, i.e. if r1 and r2 are both
;	arrays[3,n,m] then the output cross-vectors also will be an array[3,n,m] while scalar
;	products will be returned as array[n,m]
; MODIFICATION HISTORY:
;	AUG-1999, Paul Hick (UCSD/CASS; pphick@ucsd.edu)
;-

InitVar, unit	, /key
InitVar, degrees, /key
InitVar, scalar	, /key
InitVar, angle	, /key

InitVar, sphere , /key

InitVar, cylin	, /key
cylin = cylin AND 1-sphere

InitVar, rect	, /key
rect = rect OR (1-sphere AND 1-cylin)

rr1 = r1
rr2 = r2
SyncArgs, rr1, rr2

sz = size(rr1)

rr1 = reform(rr1,sz[1],sz[sz[0]+2]/sz[1], /overwrite)
rr2 = reform(rr2,sz[1],sz[sz[0]+2]/sz[1], /overwrite)

CASE 1 OF
sphere: BEGIN
	rr1 = cv_coord(from_sphere=rr1, /to_rect, degrees=degrees)
	rr2 = cv_coord(from_sphere=rr2, /to_rect, degrees=degrees)
END
cylin: BEGIN
	rr1 = cv_coord(from_cylin =rr1, /to_rect, degrees=degrees)
	rr2 = cv_coord(from_cylin =rr2, /to_rect, degrees=degrees)
END
ELSE:	; Nothing to do
ENDCASE

CASE scalar OF
0B: begin				; Vector product
	pp = rr1[[1,2,0],*]*rr2[[2,0,1],*]-rr1[[2,0,1],*]*rr2[[1,2,0],*]

	CASE angle OF
	0: BEGIN
		IF unit THEN pp = pp/([1.,1.,1.]#[sqrt(total(pp*pp,1))])
		IF NOT rect THEN pp = cv_coord(from_rect=pp, to_sphere=sphere, to_cylin=cylin, degrees=degrees)
	END
	1: BEGIN
		pp = sqrt(total(pp*pp,1)/(total(rr1*rr1,1)*total(rr2*rr2,1)))
		sz = [sz[0]-1 , sz[2:sz[0]+1] , sz[sz[0]+2]/sz[1]]
	END
	ENDCASE
END
1B: BEGIN				; Scalar product
	pp = total(rr1*rr2,1)
	if angle then pp = pp/sqrt(total(rr1*rr1,1)*total(rr2*rr2,1))
	sz = [sz[0]-1 , sz[2:sz[0]+1] , sz[sz[0]+2]/sz[1]]
END
ENDCASE

SyncDims, pp, sizeinfo=sz

RETURN, pp  &  END