FUNCTION vectorproduct, r1, r2, unit=unit, cylin=cylin, sphere=sphere, rect=rect,   $
    degrees=degrees, scalar=scalar, angle=angle
    @compile_opt.pro            ; On error, return to caller

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