File Coverage

blib/lib/Physics/Particles.pm

Criterion	Covered	Total	%
statement	99	113	87.6
branch	6	14	42.8
condition	3	9	33.3
subroutine	14	17	82.3
pod	8	8	100.0
total	130	161	80.7

# See the POD documentation at the end of this
# document for detailed copyright information.
# (c) 2002-2003 Steffen Mueller, all rights reserved.

package Physics::Particles;

use 5.006;
use strict;
use warnings;

use constant C_VACUUM         => 299792458;
use constant C_VACUUM_SQUARED => 299792458 * 299792458;

use Carp;

use Data::Dumper;

use vars qw/$VERSION/;
$VERSION = '1.02';


# constructor new
#
# Does not require any arguments. All arguments
# directly modify the object as key/value pairs.
# returns freshly created simulator object.

sub new {
	my $proto = shift;
	my $class = ref $proto || $proto;

	# Make a new object with default values.
	my $self = {
		forces => [], # All forces (callbacks) will be stored here.
		p      => [], # All particles (hashrefs) will be stored here.
		p_attr => {
			x => 0,
			y => 0,
			z => 0,
			vx => 0,
			vy => 0,
			vz => 0,
			m => 1,
			n => '',
		},
		@_
	};

	bless $self => $class;
}


# method set_particle_default
#
# Takes a hash reference whichis then used as is for
# particle default attributes such as position, velocity, mass,
# unique ID, or whatever properties you fancy.

sub set_particle_default {
	my $self = shift;
	my $hashref = shift;

	ref $hashref eq 'HASH'
	  or croak "Must pass hash reference to set_particle_default().";

	$self->{p_attr} = $hashref;
	return 1;
}

# method add_particle
#
# Takes key/value pairs of particle attributes.
# Attributes default to whatever has been set using
# set_particle_default.
# A new particle represented by the attributes is then
# created in the simulation and its particle number is
# returned. Attributes starting with an underscore are
# reserved to internal attributes (so don't use any of them).

sub add_particle {
	my $self = shift;

	my $particle = {
		%{$self->{p_attr}},
		@_,
		_fx => 0,
		_fy => 0,
		_fz => 0,
	};

	my $particle_no = $self->_make_particle();

	$self->{p}[$particle_no] = $particle;

	return $particle_no;
}


# private method _make_particle
#
# Returns a currently unused particle number or
# appends an empty particle to the particle list.

sub _make_particle {
	my $self = shift;
	my $count = 0;
	foreach (@{$self->{p}}) {
		return $count unless ref $_;
		$count++;
	}

	push @{ $self->{p} }, undef;
	return $count;
}


# method clear_particles
#
# Removes all particles from the simulation.

sub clear_particles {
	my $self = shift;
	$self->{p} = [];
	return 1;
}


# method remove_particle
#
# Removes a specific particle from the simulation.
# Takes the particle number as argument.
# Returns 1 on success, 0 otherwise.

sub remove_particle {
	my $self = shift;
	my $particle_no = shift;

	return 0 if $particle_no < 0 || $particle_no > $#{$self->{p}};

	$self->{p}[$particle_no] = undef;

	return 1;
}


# method add_force
#
# Adds a new force to the simulation.
# Takes a subroutine reference (the force) as argument,
# appends it to the list of forces and returns the force it
# just appended. Second argument is optional and defaults to false.
# It is a boolean indicating whether or not the force is symmetric.
# Symmetric means that the force particle1 excerts on particle2
# is exactly zero minus the force particle2 excerts on particle1.

sub add_force {
	my $self = shift;
	my $force = shift;
	my $symmetric = shift;
	$symmetric = 0 unless defined $symmetric;

	push @{$self->{forces}}, [$force, $symmetric];
	return $force;
}


# method iterate_step
#
# Applies all forces (excerted by every particle) to all particles.
# That means this is of complexity no_particles*no_particles*forces.
# Takes a list of additional parameters as argument that will be
# passed to every force subroutine.

sub iterate_step {
	my $self = shift;
	my @params = @_;

	my $forces    = [];
	foreach (@{ $self->{p} }) {
		push @$forces, [$_->{_fx}, $_->{_fy}, $_->{_fz}];
		$_->{_fx} = $_->{_fy} = $_->{_fz} = 0;
	}

	my $new_state = [];

	my $p_no = 0;
	foreach my $p (@{ $self->{p} }) {
		my $new_p = { %$p };
		push @$new_state, $new_p;
		my $f = $forces->[$p_no];
		foreach my $force_def (@{$self->{forces}}) {
			my $symm = $force_def->[1];
			my $force = $force_def->[0];
			my $exc_no = 0;
			foreach my $excerter (@{$self->{p}}) {
				last if $symm and $exc_no >= $p_no;
				$exc_no++, next if $exc_no == $p_no;
				my @this_f = $force->(
					$p,
					$excerter,
					\@params
				);
				$f->[0] += $this_f[0];
				$f->[1] += $this_f[1];
				$f->[2] += $this_f[2];
				if ($symm) {
					$forces->[$exc_no][0] -= $this_f[0];
					$forces->[$exc_no][1] -= $this_f[1];
					$forces->[$exc_no][2] -= $this_f[2];
				}
				$exc_no++;
			}
		}
		$p_no++;
	}

	$p_no = 0;
	foreach my $new_p (@$new_state) {
		my $f = $forces->[$p_no++];
		# accel_i = force_i / mass (i-th component)
		my $m = $new_p->{m} /
			sqrt(1 -
				($new_p->{vx}**2 +
				 $new_p->{vy}**2 +
				 $new_p->{vz}**2) /
			 	C_VACUUM_SQUARED
			);
		my @a = ($f->[0] / $m, $f->[1] / $m, $f->[2] / $m);

		@a = map $_ * $params[0], @a;

		$new_p->{x} += $new_p->{vx} * $params[0] +
				0.5 * $a[0] * $params[0];
		$new_p->{y} += $new_p->{vy} * $params[0] +
				0.5 * $a[1] * $params[0];
		$new_p->{z} += $new_p->{vz} * $params[0] +
				0.5 * $a[2] * $params[0];

		$new_p->{vx} += $a[0];
		$new_p->{vy} += $a[1];
		$new_p->{vz} += $a[2];
	}

	$self->{p} = $new_state;

	return 1;
}


# method dump_state
#
# Returns a Data::Dumper dump of the state of all particles.

sub dump_state {
	my $self = shift;
	return Dumper($self->{p});
}

1;


__END__