package Vector;

use strict;
use vars qw($VERSION);

#@ISA = qw( AutoLoader);

$VERSION = "0.05";

# Autoload methods go after =cut, and are processed by the autosplit program.


# Preloaded methods go here.

# the constructor.  can be called as a class method or
# an object method.  
sub new
{
	my $proto = shift;
	my $class = ref($proto) || $proto  ;
	my $array_ref;
	my %args;
	if (ref($_[0]) eq "")
	{
		%args = @_;
		$array_ref = $args{element};
	}
	else
	{
		$array_ref = shift;
	}
	my $self = { element=> $array_ref ||  [] };

	bless $self, $class;
}

#get/set a particular component.  Will
#grow the array as required
sub component
{
	my $self = shift;
	my $index = shift;
	my $value = shift;

	if ( defined $value )
	{
		$$self{element}[$index] = $value;
	}
	else
	{
		return $$self{element}[$index];
	}

}

sub as_array
{
	my $self = shift;
	return @{$$self{element}} ;
}

# returns something like <0,0,3>
sub as_string
{
	my $self = shift;
	return '<'. join(',', $self->as_array ) .'>';
}

# how many elements?
sub length
{
	my $self = shift;
	return scalar( $self->as_array );
}

# "times" a scalar
# returns a new thingy of the same class as the calling object
sub times
{
	my $self = shift;
	my $scalar = shift;
	my @temp_el =  $self->as_array;
	my $element;
	foreach $element ( @temp_el )
	{
		$element *= $scalar; # remember, $element is by reference
	}
	return $self->new(\@temp_el);
}

# "plus" a vector
# returns a new thingy of the same class as the calling object
# which is the sum of the object and the argument
sub plus
{
	my $self = shift;
	my $other_guy = shift;
	my @temp1 =  $self->as_array;
	my @temp2 =  $other_guy->as_array;
	my $last = $#temp1 > $#temp2 ? $#temp1 : $#temp2;
	my $i;

	foreach $i ( 0..$last )
	{
		$temp1[$i] += $temp2[$i]; 
	}
	return $self->new(\@temp1);
}

# dot product is pretty simple--sum of the products of
# corresponding pieces.  
sub dot
{
	my $self = shift;
	my $other_guy = shift;
	my $result = 0;

	my @temp1 = $self->as_array;
	my @temp2 = $other_guy->as_array;

	my $last = $#temp1 ;#note:  it doesn't matter if the other one is
						#longer--if it is, this will just act as if we are
						# padding the shorter one with zeroes

	my $i;

	foreach $i (0..$last)
	{
		$result += $temp1[$i] * $temp2[$i];
	}

	return $result;

}

sub projected_onto
{
	my $self = shift;

	#my @basis;
	#my $possible_normal;
	#my $normal;
	if ( ref ($_[0]) ne '' )
	{
		my $other_guy = shift;
		if ( ref($other_guy) eq 'ARRAY' )
		{
			$other_guy = $self->new($other_guy);
		}
		return $other_guy->normalized->times($self->dot($other_guy->normalized));
	}
	else
	{
		my %args = @_;
		
		if ( exists($args{plane}) )
		{
			my @basis;
			my $normal;

			my $possible_normal = $args{plane};
			if (ref($possible_normal) eq 'ARRAY')
			{
				$normal = $self->new( $possible_normal );
			}
			else
			{
				$normal = $possible_normal;
			}

			$basis[0] = $self->cross($normal)->normalized;
			$basis[1] = $normal->cross($basis[0])->normalized;

			#@_=();
			return $self->projected_onto($basis[0])->plus($self->projected_onto($basis[1]));
		}
	}
}
# magnitude of a vector is the square root of it dotted with itself
sub magnitude
{
	my $self = shift;

	return sqrt( $self->dot($self) );
}

# returns a three element vector.  Ignores everything
# after the first element.  Shorter vectors will be
# padded with zeroes (undefs, actually)
sub cross
{
	my $self = shift;
	my $other_guy = shift;

	my @temp1 = $self->as_array;
	my @temp2 = $other_guy->as_array;

# this is calculated via the standard determinant
# for the cross product:
#
#             | i  j  k |
# a x b = det | ax ay az|
#             | bx by bz|

	return $self->new ( 
		[$temp1[1]*$temp2[2] - $temp2[1]*$temp1[2],
		 $temp1[2]*$temp2[0] - $temp2[2]*$temp1[0],
		 $temp1[0]*$temp2[1] - $temp2[0]*$temp1[1]
		]
	)
}

sub normalized
{
	my $self = shift;
	my @temp_el = @{$$self{element}};
	my $mag = $self->magnitude;
	my $len = $self->length;
	my $i;
	# The zero vector gives the zero vector
	if ( $mag == 0 )
	{
		foreach $i (0..$len-1 )
		{
			$temp_el[$i] = 0;
		}
	}
	else
	# we know the magnitude isn't zero, so it's safe to
	# divide each element by it (unless it's really, really
	# close to zero--but we'll make an exception for that in the warranty
	{
		foreach $i (0..$len-1)
		{
			$temp_el[$i] /= $mag;
		}
	}
	return $self->new( \@temp_el );
}

# rotate the object counterclockwise through the given angle, 
# about the given vector (or one going through the given array)
sub rotate
{
	my $self = shift;
	my $alpha = shift;
	my $axis = shift;

	#stuff for the rotation matrix (see below)
	my $cos_alpha = cos($alpha);
	my $one_minus = 1 - $cos_alpha;
	my $sin_alpha = sin($alpha);

	# vectorify the array if they gave us one
	if ( ref($axis) eq 'ARRAY' )
	{
		$axis = $self->new( $axis );
	}

	# these are the famous "direction cosines" of the rotation axis
	my ($a, $b, $c) = $axis->normalized->as_array;

	# now build the rows of the rotation matrix
	#See http://wwwinfo.cern.ch/asd/lhc++/clhep/manual/UserGuide/Vector/vector.html
	#or http://freeabel.geom.umn.edu/docs/reference/CRC-formulas/node45.html
	#for more information.

	# we'll want to rearrange how things are done if 
	# this ever goes production, because it will be
	# wasteful if we rebuild this every time for 
	# ten thousand vectors
	my $row1 = $self->new( 
		[
		$a*$a*$one_minus + $cos_alpha,
		$b*$a*$one_minus - $c*$sin_alpha,
		$c*$a*$one_minus + $b*$sin_alpha 
		]);


	my $row2 = $self->new( 
		[
		$a*$b*$one_minus + $c*$sin_alpha,
		$b*$b*$one_minus + $cos_alpha,
		$c*$b*$one_minus - $a*$sin_alpha 
		]);

	my $row3 = $self->new( 
		[
		$a*$c*$one_minus - $b*$sin_alpha,
		$b*$c*$one_minus + $a*$sin_alpha,
		$c*$c*$one_minus + $cos_alpha 
		]);

	# and return a vector which is the product of that matris
	# and the calling object
	return $self->new(  [ 	$self->dot($row1) ,
									$self->dot($row2),
									$self->dot($row3) ] );

}

sub tilt
{
	my $self = shift;
	my $tilt_from = shift;

	if ( ref($tilt_from) eq 'ARRAY' )
	{
		$tilt_from = new Vector ($tilt_from);
	}

	my $tilt_to = shift;
	unless (defined ($tilt_to))
	{
		$tilt_to = new Vector ( [0,0,1] );
	}

	if ( ref($tilt_to) eq 'ARRAY' )
	{
		$tilt_to = new Vector ($tilt_to);
	}

	my $tilt_axis = $tilt_from->cross($tilt_to);
	my $from_mag = $tilt_from->magnitude;
	my $to_mag = $tilt_to->magnitude;
	my $product_thereof = $from_mag*$to_mag;
	my $cos_theta = ( $tilt_from->dot($tilt_to) ) / $product_thereof;
	my $sin_theta = ( $tilt_from->cross($tilt_to)->magnitude ) / $product_thereof;
	my $tilt_angle = atan2( $sin_theta, $cos_theta );

	#@_=();
	return $self->rotate($tilt_angle, $tilt_axis);
}

1;

__END__

# here's the docs

=head1 NAME

Vector.pm 0.05 - Simple vector stuff

=head1 SYNOPSIS

     use Vector;
     $a = new Vector ( [3,4] );
     $b = new Vector ( [1,0] );

     $ax = $a->dot($b);
     
     $c = $a->cross($b);
     
     print $c->cross($b)->times(3)->as_string;

=head1 DOWNLOADING

You can obtain Vector.pm from

=for text

http://fulcrum.org/projects/polyplus/Vector.pm



=head1 DESCRIPTION

A small class for doing things with vectors.  Mostly experimental--CPAN has 
much higher powered stuff in Math::Pari and PDL.

=over 5

=item C<new>

C<new( [1,0,0] ) > I<class/object method>

To create a new vector, send the "new" message to the Vector class, or call
it against an existing vector.  Initialize its elements with a reference to 
an array.

=item C<component>

C<component( $index [, $value])> I<object method>

If a $value is supplied, sets the $index-th component of the vector to
that value.  Otherwise, returns the existing value.  Will grow the
vector to accommodate if you assign off the end.

	$a = new Vector( [1,0,2] );
	$a_sub_x = $a->component(0);
	$a->component(1,5);
	

=item C<length>

C<length> I<object method>

Returns the length of the vector.

=item C<magnitude>

C<magnitude> I<object method>

Returns the magnitude of the vector (the square root of the sum of the squares
of the components).

=item C<times>

C<times( )> I<object method>

Returns the vector resulting from scaling the current vector by the 
indicated amount.

     $b = $a->times(3);

=item C<plus>

C<plus( $another_vector )> I<object method>

Returns the (vector) sum of the current vector and the argument:

     $sum = $a->plus($b);


=item C<dot>

C<dot( $another_vector )> I<object method>

Returns the dot product of the current vector and the argument:

     $square_mag = $a->dot($a);

=item C<projected_onto>

C<projected_onto( $another_vector )> I<object method>

Returns the vector "shadow" of the current vector projected onto the 
argument:

     $x_part = $a->projected_onto([1,0,0]);

The argument, as in this example, can be an array reference rather than
a vector.

=item C<cross>

C<cross( $another_vector ) > I<object method>

Returns the (vector) cross product of the object and the argument.

     $k = $i->cross($j);

=item C<normalized>

C<normalized> I<object method>

Returns a vector of magnitude one pointing in the same direction as the object:

     $n = $b->normalized;

=item C<as_array>

C<as_array> I<object method>

Returns the elements of the vector as an array:

     @dir_cos = $n->as_array;

=item C<as_string>

C<as_string> I<object method>

Returns a string of the form "<1,2>".

     print $n->as_string;
     print $a->cross($b)->normalized->times(3)->as_string;


=item C<rotate>

C<rotate ($angle, $axis | [$x,$y,$z] )> I<object method>

Returns the vector obtained by rotating the current vector through
$angle radians about the vector $axis (or, if an array ref is
given, a vector which goes through those points).  The $axis
need not be normalized (likewise $x, $y, $z can be any point).

	$pi = 4 * atan2(1,1);
	$a = new Vector ( [1,0,0] ) ;
	print $a->rotate( 2*$pi/3, [1,1,1] )->as_string ;

See http://wwwinfo.cern.ch/asd/lhc++/clhep/manual/UserGuide/Vector/vector.html
or http://freeabel.geom.umn.edu/docs/reference/CRC-formulas/node45.html
for more information.

=item C<tilt>

C<tilt($from_here, $to_here)> I<object method>

Take the given vector, apply the transformation to it which would rotate
the C<$from_here> vector so that it would be aligned with C<$to_here>,
and return the result.  (For example, if your camera is at C<$from_here>,
you can just tilt all your vectors with C<($from_here, [0,0,1])> and
then hand off the x- and y- components of those vectors to a renderer
to get the picture from that angle.)


=back
=cut