#!/usr/bin/perl -w

# vector stuff

use strict;
#use varnings;
use POSIX;

package Vec;
#use Math::GSL::Matrix;
use Data::Dumper;

use Exporter 5.57 'import';
our @EXPORT_OK =
    qw(
	  PI RT2
	  Sin SinA SinR Cos CosA CosR Tan TanA TanR acos acosA acosR asin asinA asinR tan tanA tanR
	  Rad2Deg Rad2DegA Rad2DegR Deg2Rad Deg2RadA Deg2RadR
	  Acos AcosA AcosR Asin AsinA AsinR
	  RadNormPM RadNormPMA RadNormPMR RadNormPL RadNormPLA RadNormPLR DegNormPM DegNormPMA
	  DegNormPMR DegNormPL DegNormPLA DegNormPLR
	  Vec2Polar2d Polar2vec2d Normal2d nNormal2d AngleX2d diff_angle2d OuterProduct3d
	  MinMax Length AngleAxCos AngleAx Angle2Cos Angle2 Neg Mirror AddS Add Sub Distance
	  MulS Unit InnerProduct IsUndef NormPoint str2num str2vec str2veclist PrintXY vec2str ShowVeclist
  );

# vetor things, point or vector = [ $x, $y ]
# or of general length

###############
# missing math and rad/deg conversion

sub PI()  { 2 * atan2(1,0); }
sub RT2() { sqrt(2); }

sub Sin($) {
    my $x = shift;
    my $r;
    if (defined($x)) { $r = sin($x); }
    $r;
}
sub SinA(@) {
    my @r;
    for my $x (@_) { push @r, Sin($x); }
    @r;
}
sub SinR($) { [ SinA(@{$_[0]}) ]; }

sub Cos($) {
    my $x = shift;
    my $r;
    if (defined($x)) { $r = cos($x); }
    $r;
}
sub CosA(@) {
    my @r;
    for my $x (@_) { push @r, Cos($x); }
    @r;
}
sub CosR($) { [ CosA(@{$_[0]}) ]; }

sub Tan($) {
    my $x = shift;
    my $r;
    if (defined($x)) { $r = tan($x); }
    $r;
}
sub TanA(@) {
    my @r;
    for my $x (@_) { push @r, Tan($x); }
    @r;
}
sub TanR($) { [ TanA(@{$_[0]}) ]; }

sub acos($) {
    # not in perl, have to add it
    my $x = shift;
    my $r;
    if (defined($x) && abs($x) <= 1) { $r = atan2( sqrt(1 - $x * $x), $x ); }
    $r;
}
sub acosA(@) {
    my @r;
    for my $x (@_) { push @r, acos($x); }
    @r;
}
sub acosR($) {
    [ acosA(@{$_[0]}) ];
}

sub asin($) {
    # not in perl, have to add it
    my $x = shift;
    my $r;
    if (defined($x) && abs($x) <= 1) { $r = atan2(  $x, sqrt(1 - $x * $x) ); }
    $r;
}
sub asinA(@) {
    my @r;
    for my $x (@_) { push @r, asin($x); }
    @r;
}
sub asinR($) { [ asinA(@{$_[0]}) ]; }

sub tan($) {
    # not in perl, have to add it
    my $x = shift;
    my $r;
    if (defined($x)) {
	my $y = cos($x);
	if ($y != 0) { $r = sin($x) / $y; }
    }
    $r;
}
sub tanA(@) {
    my @r;
    for my $x (@_) { push @r, tan($x); }
    @r;
}
sub tanR($) { [ tanA(@{$_[0]}) ]; }

sub Rad2Deg($) {
    # conversion from radians to degrees (360° = one full turn)
    my $x = shift;
    my $r;
    if (defined($x)) { $r = 180 * $x / PI; }
    $r;
}
sub Rad2DegA(@) {
    my @r;
    for my $x (@_) { push @r, Rad2Deg($x); }
    @r;
}
sub Rad2DegR($) { [ Rad2DegA(@{$_[0]}) ]; }

sub Deg2Rad($) {
    # conversino from degrees to radians
    my $x = shift;
    my $r;
    if (defined($x)) { $r = PI * $x / 180; }
    $r;
}
sub Deg2RadA(@) {
    my @r;
    for my $x (@_) { push @r, Deg2Rad($x); }
    @r;
}
sub Deg2RadR($) { [ Deg2RadA(@{$_[0]}) ]; }

sub Acos($) {
    my $x = shift;
    Rad2Deg(acos($x));
}
sub AcosA(@) {
    my @r;
    for my $x (@_) { push @r, Acos($x); }
    @r;
}
sub AcosR($) { [ AcosA(@{$_[0]}) ]; }

sub Asin($) {
    my $x = shift;
    Rad2Deg(asin($x));
}
sub AsinA(@) {
    my @r;
    for my $x (@_) { push @r, Asin($x); }
    @r;
}
sub AsinR($) { [ AsinA(@{$_[0]}) ]; }

sub RadNormPM($) {
    # normalize angle, -PI < angle <= PI
    my $da = shift;
    if (defined($da)) {
	while ($da >   PI) { $da -= 2*PI; }
	while ($da <= -PI) { $da += 2*PI; }
    }
    $da;
}
sub RadNormPMA(@) {
    my @r;
    for my $x (@_) { push @r, RadNormPM($x); }
    @r;
}
sub RadNormPMR($) { [ RadNormPMA(@{$_[0]}) ]; }

sub RadNormPL($) {
    # normalize angle, 0 <= angle < 2*PI
    my $da = shift;
    if (defined($da)) {
	while ($da >= 2*PI) { $da -= 2*PI; }
	while ($da <  0    ) { $da += 2*PI; }
    }
    $da;
}
sub RadNormPLA(@) {
    my @r;
    for my $x (@_) { push @r, RadNormPL($x); }
    @r;
}
sub RadNormPLR($) { [ RadNormPLA(@{$_[0]}) ]; }

sub DegNormPM($) {
    # normalize angle, -PI < angle <= PI
    my $da = shift;
    if (defined($da)) {
	while ($da >   180) { $da -= 360; }
	while ($da <= -180) { $da += 360; }
    }
    $da;
}
sub DegNormPMA(@) {
    my @r;
    for my $x (@_) { push @r, DegNormPM($x); }
    @r;
}
sub DegNormPMR($) { [ DegNormPMA(@{$_[0]}) ]; }

sub DegNormPL($) {
    # normalize angle, 0 <= angle < 2*PI
    my $da = shift;
    if (defined($da)) {
	while ($da >= 360) { $da -= 360; }
	while ($da <  0  ) { $da += 360; }
    }
    $da;
}
sub DegNormPLA(@) {
    my @r;
    for my $x (@_) { push @r, DegNormPL($x); }
    @r;
}
sub DegNormPLR($) { [ DegNormPLA(@{$_[0]}) ]; }

sub Floor($) {
    # normalize angle, 0 <= angle < 2*PI
    my $x = shift;
    if (defined($x)) {
	return POSIX::floor($x);
    } else {
	return undef;
    }
}
sub FloorA(@) {
    my @r;
    for my $x (@_) { push @r, Floor($x); }
    @r;
}
sub FloorR($) { [ FloorA(@{$_[0]}) ]; }

sub Ceil($) {
    # normalize angle, 0 <= angle < 2*PI
    my $x = shift;
    if (defined($x)) {
	return POSIX::ceil($x);
    } else {
	return undef;
    }
}
sub CeilA(@) {
    my @r;
    for my $x (@_) { push @r, Ceil($x); }
    @r;
}
sub CeilR($) { [ CeilA(@{$_[0]}) ]; }

###############
# 2d stuff, arguments for thees are arrays: [ x, y ]

sub Vec2Polar2d($) {
    my $p = shift;
    my $r = Length($p);
    my $a = AngleX2d($p);
    [ $r, $a ];
}

sub Polar2vec2d($) {
    my $p = shift;
    my $r = $$p[0];
    my $a = $$p[1];

    if (defined($r) && defined($a)) {
	return [ $r * cos($a), $r * sin($a) ];
    } else {
	return [ undef, undef ];
    }
}

sub Normal2d($) {
    # returns a rotated +90°
    my $a = shift;
    if (defined($$a[1])) {
	return [ -$$a[1], $$a[0] ];
    } else {
	return [ undef, $$a[0] ];
    }
}

sub nNormal2d($) {
    # returns a rotated -90°
    my $a = shift;
    if (defined($$a[0])) {
	return [ $$a[1], -$$a[0] ];
    } else {
	return [ $$a[1], undef ];
    }
}

sub AngleX2d($) {
    # gives a vectors angle to the x-axis
    # returns values in range [ -90, 90 ], but in radians
    my $v = shift;
    if (defined($$v[1]) && defined($$v[0])) {
	return atan2($$v[1], $$v[0]);
    } else {
	return undef;
    }
}

sub diff_angle2d($$) {
    # = angle between a and b
    my $a = shift;
    my $b = shift;

    my $s = AngleX2d($a);
    my $t = AngleX2d($b);

    my $d = $t - $s;
    RadNormPM($d);
}

###############
# 3d stuff, arguments for thees are arrays: [ x, y, z ]

sub OuterProduct3d($$) {
    my $a = shift;
    my $b = shift;

    return undef if (!defined($a) || !defined($b));
    for (my $ix = 0; $ix < 3; $ix++) {
	return undef if (!defined($$a[$ix]) || !defined($$b[$ix]));
    }
    my @c = (
	$$a[1] * $$b[2] - $$b[1] * $$a[2],
	$$a[2] * $$b[0] - $$b[2] * $$a[0],
	$$a[0] * $$b[1] - $$b[0] * $$a[1]
    );
    [ @c ];
}

###############
# generic length

sub MinMax(@) {
    my @veclist = @_;
    my @min;
    my @max;

    my $len = 0;
    # get length of longest vector in @veclist
    for my $vec (@veclist) {
	if (defined($vec) && $len < @$vec) {
	    $len = @$vec;
	}
    }

    for (my $ix = 0; $ix < $len; $ix++) {
	for (my $vx = 0; $vx < @veclist; $vx++) {
	    my $vec = $veclist[$vx];
	    next unless (defined($vec));
	    my $v = $$vec[$ix];
	    if (defined($v)) {
		if (!defined($min[$ix]) || $min[$ix] > $v) { $min[$ix] = $v; }
		if (!defined($max[$ix]) || $max[$ix] < $v) { $max[$ix] = $v; }
	    }
	}
    }

    (\@min, \@max);
}

sub Length($) {
    # returns Length of a vector
    # note, length is != dimension (number of elements in array)
    my $v = shift;

    my $sum = 0;
    return undef if (!defined($v));
    for (my $ix = 0; $ix < @$v; $ix++) {
	return undef if ( !defined($$v[$ix]) );
	$sum += $$v[$ix] * $$v[$ix];
    }
    sqrt($sum);
}

sub AngleAxCos($) {
    # gives the cosinus of a vectors angle to the x-axis
    my $v = shift;
    my $len = Length($v);

    my @res = ();
    if ($len) {
	for (my $ix = 0; $ix < @$v; $ix++) {
	    if ( defined($$v[$ix]) ) {
		$res[$ix] = $$v[$ix] / $len;
	    }
	}
    }

    [ @res ];
}
sub AngleAx($) {
    # get the angle to x,y,z,... axes from the cosines
    acosR(AngleAxCos($_[0]));
}

sub Angle2Cos($$) {
    # get the angle between to lines/vectors
    my $v1 = shift;
    my $v2 = shift;

    my $len1 = Length($v1);
    my $len2 = Length($v2);
    my $angleCos;
    if ($len1 && $len2) {
	my $sp   = InnerProduct($v1, $v2);
	$angleCos = $sp / ($len1 * $len2);
    } else {
	$angleCos = undef;
    }
    $angleCos;
}
sub Angle2($$) { Acos(Angle2Cos($_[0],$_[1])); }

sub Neg($) { [ map { defined() ? - $_ : undef } @{$_[0]} ]; }

sub Mirror($;$) {
    # Mirror($a,0) mirror on X-axel, Mirror($a,1) mirror on y-axis, etc.
    my $a = shift;
    my $ix = shift // 0;

    if (defined($$a[$ix])) {
	$$a[$ix] = - $$a[$ix];
    }
    $a;
}

sub AddS( $ $ ) {
    # Add a scalar, result[:] = a + b[:]
    my $a = shift;
    my $b = shift;

    my @res = ();

    if (!defined($a) || !defined($b)) { return []; }
    for (my $ix = 0; $ix < @$b; $ix++) {
	if (defined($$b[$ix])) {
	    $res[$ix] = $a + $$b[$ix];
	}
    }

    [ @res ];
}

sub Add(@) {
    # add vecors together, result = a + b + ...
    my @res = ();
    my $dim = 0;

    for my $a (@_) {
	if ($dim < @$a) {
	    $dim = @$a;
	}

	for (my $ix = 0; $ix < $dim; $ix++) {
	    if (defined($$a[$ix])) {
		if (defined($res[$ix])) {
		    $res[$ix] += $$a[$ix];
		} else {
		    $res[$ix] = $$a[$ix];
		}
	    }
	}
    }

    [ @res ];
}

sub Sub($$) {
    # result = a - b
    my $a = shift;
    my $b = shift;

    my @res = ();
    if (!defined($a) || !defined($b)) { return []; }
    my $dim = @$a;
    if (@$a > @$b) {
	$dim = @$b;
    }

    for (my $ix = 0; $ix < $dim; $ix++) {
	if (defined($$a[$ix]) && defined($$b[$ix])) {
	    $res[$ix] = $$a[$ix] - $$b[$ix];
	}
    }
    [ @res ];
}

sub Distance( $ $ ) {
    # = distance between point a and b
    my $a = shift;
    my $b = shift;

    Length(Sub($b, $a));
}

sub MulS($$) {
    # muliply vetctor with a scalar, returns $k * $b, where k is scalar and b vector
    my $k = shift;
    my $b = shift;

    my @res = ();
    for (my $ix = 0; $ix < @$b; $ix++) {
	if (defined($$b[$ix])) {
	    $res[$ix] = $k * $$b[$ix];
	}
    }
    [ @res ];
}

sub Unit($) {
    # returns a's unit vector
    my $a = shift;

    my $l = Length($a);
    my $e = [];
    if ($l) {
	$e = MulS( 1/$l, $a);
    }
    $e;
}

sub InnerProduct($ $) {
    # scalar product
    # result (a scalar) = a dot b
    my $a = shift;
    my $b = shift;

    my $sum = 0;
    my $len = @$a;
    if ($len < @$b) {
	$len = @$b;
    }

    for (my $ix = 0; $ix < $len; $ix++) {
	if ( defined($$a[$ix]) && defined($$b[$ix] ) ) {
	    $sum += $$a[$ix] * $$b[$ix];
	}
    }
    $sum;
}

sub IsUndef($) {
    my $v = shift;
    my $r = 0;
    if (defined($v)) {
	for my $e (@$v) {
	    if (!defined($e)) {
		$r = 1;
		last;
	    }
	}
    } else {
	$r = 1;
    }
    $r;
}

sub NormPoint($$$) {
    # returns point on line trough a and b where the lines normal goea through c
    my $a = shift;
    my $b = shift;
    my $c = shift;

    if (IsUndef($a) || IsUndef($b) || IsUndef($c)) { return []; }

    # translate so that oiunt a to origo
    my $B = Sub($b, $a);
    my $C = Sub($c, $a);
    my $D = Sub($c, $b);

    my $k = Length($B);
    my $l = Length($C);
    my $m = Length($D);
    my $r;

    return $a if ($l == 0 || $k == 0);
    return $b if ($m == 0);

    my $BU = Unit($B);
    my $CU = Unit($C);

    my $ip = InnerProduct($BU,$CU);
    print "ip: $ip\n";
    $ip *= $l;
    print "ip: $ip\n";

    $r = MulS($ip, $BU);
    $r = Add($a, $r);

    $r;
}

##################
# matrix * vector

##################
# supporting tools

sub bbInside($$;$);
sub bbInside($$;$) {
    my $bb = shift;
    my $pos = shift;
    my $ebb = shift // [];

    my $lenl = @{$$bb[0]};
    my $lenu = @{$$bb[1]};
    my $lenp = @$pos;

    return undef unless ($lenl == $lenu && $lenl == $lenp);
    return undef if (IsUndef($$bb[0]));
    return undef if (IsUndef($$bb[1]));
    return undef if (IsUndef($pos));

    my $res = 1;
    for (my $ix = 0; $ix < $lenl; $ix++) {
	if ($$bb[0][$ix] <= $$pos[$ix]  && $$pos[$ix] <= $$bb[1][$ix] ) {
	} else {
	    $res = 0;
	    last;
	}
    }
    if ($res == 1) {
	for my $BB (@$ebb) {
	    if (bbInside($BB, $pos)) {
		$res = 0;
		last;
	    }
	}
    }

    $res;
}

sub str2num($) {
    my $str = shift;
    $! = 0;
    my ($num, $n_unparsed) = POSIX::strtod($str);

    if (($str eq '') || ($n_unparsed != 0) || $!) {
	$num = undef;
    }
    $num;
}
sub str2vec($) {
    my $str = shift;

    my @list = split(/,/, $str);
    my @pos = ();
    for (my $ix = 0; $ix < @list; $ix++) {
	push @pos, str2num($list[$ix]);
    }
    @pos;
}
sub str2veclist($) {
    my $str = shift;
    my @list = split(/\//, $str);
    my @veclist = ();
    for (my $ix = 0; $ix < @list; $ix++) {
	push @veclist, [ str2vec($list[$ix]) ];
    }

    @veclist;
}

sub PrintXY(@) {
    my @point = @_;
    my $str = "";
    for my $p (@point) {
	if (!defined($p)) {
	    $str .= "\n";
	} elsif (ref($p) eq "ARRAY") {
	    $str .= sprintf "%7.2f %7.2f\n", @$p;
	} elsif (ref($p) eq "HASH" && defined($$p{xy})) {
	} else {
	}
    }
    $str;
}

sub vec2str($;$$) { # print out a vectors x and y value (mostly for debugging)
    my $v = shift;
    my $u = shift // " undef,";
    my $format = shift;

    my $str;
    if (defined($v)) {
	$str = "\[";
	if ($format) {
	    for (my $ix = 0; $ix < @$v; $ix++) {
		if (defined($$v[$ix])) {
		    $str .= sprintf($format, $$v[$ix]);
		} else {
		    $str .= $u;
		}
	    }
	} else {
	    for (my $ix = 0; $ix < @$v; $ix++) {
		if (defined($$v[$ix])) {
		    $str .= " $$v[$ix],";
		} else {
		    $str .= $u;
		}
	    }
	}
	#if (length($str) > 1) { substr($str, -1, 1, ""); } # drop the final ","
	$str .= " \]";
    } else {
	$str .= $u;
    }
    $str;
}

sub ShowVeclist(@) {
    my @argv = @_;

    for (my $ax = 0; $ax < @argv; $ax++) {
	my $e = $argv[$ax];
	if (@$e == 1) {
	    # just a single vecor
	    my $v = $$e[0];
	    if (@$v == 1) {
		my $s = $$v[0];
		print "scalar: $s\n";
	    } else {
		print "vector: ", vec2str($v), "\n";
	    }
	} else {
	    # list of vectors
	    print "list of vectors: ";
	    for my $v (@$e) {
		print vec2str($v), " ";
	    }
	    print "\n";
	}
    }
}

sub test($){
    my $func = shift;

    my @argv = ();
    my $u = undef;

    if ($func eq "math") {
	printf("pi : %10f\n", PI);
	printf("rt2: %10f\n", RT2);
	printf("Acos(%f): %10f\n", $_[0], Acos($_[0]));
	printf("Asin(%f): %10f\n", $_[0], Asin($_[0]));
	printf("Tan(%f): %10f\n", $_[0], Tan($_[0]));
	printf("rad_normalize1(%f): %10f\n", $_[1], rad_normalize1($_[1]));
	printf("rad_normalize2(%f): %10f\n", $_[1], rad_normalize2($_[1]));
    } else {
	#for (@_) { push @argv, [ str2veclist($_) ]; }
	for (@_) { push @argv, str2veclist($_); }
	#print Dumper(\@argv);

	if ($func eq "2d") {
	    my $a = $argv[0];
	    my $b = $argv[1];
	    my $c;
	    my $d;

	    print "a = ", vec2str($a), "\n";
	    print "b = ", vec2str($b), "\n";
	    $c = Vec2Polar2d($a);
	    print "Vec2Polar2d(a) = ", vec2str($c), " / ", Rad2Deg($$c[1]), "\n";
	    print "  and back: ", vec2str(Polar2vec2d($c)), "\n";

	    $c = Normal2d($a);
	    print "Normal2d(a): ", vec2str($c), "\n";
	    $c = nNormal2d($a);
	    print "nNormal2d(a): ", vec2str($c), "\n";

	    $c = AngleX2d($a);
	    my $dd = Rad2Deg($c);
	    print "Angle [$$a[0],$$a[1]] to x-axis:  ", $c, " rad, ", $dd, "° = ", DegNormPL($dd), "\n";

	    if (defined($b)) {
		$c = diff_angle2d($a, $b);
		print "diff_angle2d, how much to rot. a to line up with b: ", $c, " ", Rad2Deg($c), "\n";
	    }

	} elsif ($func eq "3d") {
	    my $a = $argv[0];
	    my $b = $argv[1];
	    if (!defined($a) && @$a < 3 && !defined($b) && @$b < 3) {
		print "need two veotors of length >= 3\n"
	    } else {
		my $c = OuterProduct3d($a, $b);
		print "a x b = ", vec2str($c), "\n";
	    }

	} elsif ($func eq "vmath") {
	    for (my $ix = 0; $ix < @argv; $ix++) {
		print "a[$ix] = ", vec2str($argv[$ix]), "\n";
	    }
	    print "\n";

	    for (my $ix = 0; $ix < @argv; $ix++) {
		my $s = AngleAxCos($argv[$ix]);
		my $t = acosR($s);
		my $u = Rad2DegR($t);
		print "cos (a[$ix]) = ", vec2str($s), "\n";
		print "rad (a[$ix]) = ", vec2str($t), "\n";
		print "deg (a[$ix]) = ", vec2str($u), "\n";
		print "\n";
	    }
	    print "\n";

	    my $td = [ 0, 30, 45, 60, 90, 120, 135, 150, 180, 210, 225, 240, 270, 300, 315, 330, 360 ];
	    my $tr = Deg2RadR($td);
	    my $cs = CosR($tr);
	    print " td       =\n", vec2str( $td      , $u, " %7.3f,"), "\n";
	    print " tr (rad) =\n", vec2str( $tr      , $u, " %7.3f,"), "\n";
	    print " sin(td ) =\n", vec2str( SinR($tr), $u, " %7.3f,"), "\n";
	    print "cs cos(td)=\n", vec2str( $cs      , $u, " %7.3f,"), "\n";
	    print " tan(td ) =\n", vec2str( TanR($tr), $u, " %7.3g,"), "\n";

	    print "acos(cs ) =\n", vec2str(Rad2DegR(acosR($cs)), $u, " %7.3g,"), "\n";
	    print "asin(cs ) =\n", vec2str(Rad2DegR(asinR($cs)), $u, " %7.3g,"), "\n";

	    print "Acos(cs ) =\n", vec2str(AcosR($cs), $u, " %7.3g,"), "\n";
	    print "Asin(cs ) =\n", vec2str(AsinR($cs), $u, " %7.3g,"), "\n";

	    my $d = RadNormPMR($tr);
	    print "norm pm   =\n", vec2str($d, $u, " %7.3g,"), "\n";
	    print "norm pl   =\n", vec2str(RadNormPLR($tr), $u, " %7.3g,"), "\n";
	    $d = DegNormPMR($td);
	    print "deg  pm   =\n", vec2str($d, $u, " %7.3g,"), "\n";
	    print "deg  pl   =\n", vec2str(DegNormPLR($td), $u, " %7.3g,"), "\n";
	} else {
	    my $a = $argv[0];
	    my $b = $argv[1];
	    my ($c, $d, $e, $f);
	    for (my $ix = 0; $ix < @argv; $ix++) {
		my $v = $argv[$ix];
		my $c = AngleAxCos($v);
		my $d = AngleAx($v);
		my $e = Rad2DegR($d);
		my $l = Length($v);
		print "a[$ix] = ", vec2str($v, "   undef,", " %7.3f,"), "\tLen: ", defined($l)? $l: "undef", "\n";
		print " neg = ", vec2str(Neg($v), "   undef,", " %7.3f,"), "\n";
		print "cos  = ", vec2str($c, "   undef,", " %7.3f,"), "\n";
		print " rad = ", vec2str($d, "   undef,", " %7.3f,"), "\n";
		print " deg = ", vec2str($e, "   undef,", " %7.3f,"), "\n";
		print "unit = ", vec2str(Unit($v), "   undef,", " %7.3f,"), "\n";
		print "\n";
	    }
	    print "\n";

	    ($c, $d) = MinMax(@argv);
	    print "minmax:\n ", vec2str($c), "\n ", vec2str($d), "\n\n";

	    $c = Angle2Cos($a, $b);
	    if (!defined($c)) {
		$c = $d = $e = $f = " undef";
	    } else {
		$d = acos($c);
		$e = Rad2Deg($d);
		$f = Angle2($a, $b);
	    }
	    print "Angle2Cos(a,b) = $c $d $e $f\n";

	    $d = $$b[0];
	    $c = AddS($d, $a);
	    print "a + $d = ", vec2str($c, "   undef,", " %7.3f,"), "\n";
	    $c = Add($a, $b);
	    print "a + b = ", vec2str($c, "   undef,", " %7.3f,"), "\n";
	    $c = Sub($a, $b);
	    print "a - b = ", vec2str($c, "   undef,", " %7.3f,"), "\n";
	    $c = Distance($a, $b);
	    print "dist a b = ", defined($c)? $c: "undef","\n";
	    $c = MulS($d, $a);
	    print "a * $d = ", vec2str($c, "   undef,", " %7.3f,"), "\n";

	    $c = InnerProduct($b, $argv[2]);
	    print "ip     = $c\n";
	    $c = NormPoint($a, $b, $argv[2]);
	    print "normp  = ", vec2str($c, "   undef,", " %7.3f,"), "\n";

	}
    }

}

#test(@ARGV);

1;


__END__