File Coverage

blib/lib/Physics/Springs/Friction.pm

Criterion	Covered	Total	%
statement	21	54	38.8
branch	0	8	0.0
condition			n/a
subroutine	7	8	87.5
pod	1	1	100.0
total	29	71	40.8

package Physics::Springs::Friction;

use 5.006;
use strict;
use warnings;

use Carp;
use base 'Physics::Springs';

use constant FFORCES => '_PhSpringsFriction_friction_forces';
use Sub::Assert;

our $VERSION = '1.01';

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

	my $self = $class->SUPER::new();
	$self->{FFORCES} = [];
	return $self;
}

assert
	context => 'novoid',
	pre => '@PARAM == 1',
	post => '@RETURN == 1',
	sub => 'new';


sub add_friction {
	my $self   = shift;

	my $impl = shift;
	my $magn = shift;

	my $friction = {
		implementation => 'stokes',
		magnitude      => 0.1,
	};
	if (defined $impl) {
		unless (ref($impl) eq 'CODE' or
			$impl eq 'stokes' or
			$impl eq 'newton'
		       ) {

		    	$impl = <
sub {
	my \$P = shift;
	my \$M = shift;
	my \$F = [0, 0, 0];
	$impl
	return \@\$F;
}
ENDIMPL
			local $@;
			my $closure = eval $impl;
			croak "Error while evaluating friction " .
				"implementation.\n ($@)" if $@;
			$impl = $closure;
		}
		$friction->{implementation} = $impl;
	}

	$friction->{magnitude} = $magn if defined $magn;

	push @{$self->{FFORCES}}, $friction;
	return $#{$self->{FFORCES}};
}

assert
	pre =>	[
			'@PARAM <= 3 && @PARAM > 0',
			'@PARAM == 1 ? 1 : defined($PARAM[1])',
			'@PARAM != 3 ? 1 : defined($PARAM[2])',
		],
	post =>	'$VOID || (@RETURN == 1 and $RETURN >= 0)',
	sub => 'add_friction';



sub iterate_step {
	my $self      = shift;
	my @params    = @_;
	my $time_diff = $params[0];
	foreach my $friction (@{$self->{FFORCES}}) {
		my $magnitude = $friction->{magnitude};
		my $func      = $friction->{implementation};
		if ($func eq 'stokes') {
			# F = 6*PI*eta*r*v
			# magnitude = 6*PI*eta*r
			foreach my $P (@{$self->{p}}) {
				my ($vx, $vy, $vz) = (
					$P->{vx}, $P->{vy}, $P->{vz}
				);

				my $v = sqrt(
						$vx**2 +
						$vy**2 +
						$vz**2
				);

				next if $v == 0;

				my $f = $magnitude * $v;

				$P->{_fx} -= $f * ($vx/$v);
				$P->{_fy} -= $f * ($vy/$v);
				$P->{_fz} -= $f * ($vz/$v);
			}
		}
		elsif ($func eq 'newton') {
			# F = 0.5*c*rho*A*v^2
			# magnitude = c*rho*A
			foreach my $P (@{$self->{p}}) {
				my ($vx, $vy, $vz) = (
					$P->{vx}, $P->{vy}, $P->{vz}
				);

				my $v = sqrt(
						$vx**2 +
						$vy**2 +
						$vz**2
					);

				next if $v == 0;

				my $f = 0.5 * $magnitude * $v**2;

				$P->{_fx} -= $f * ($vx/$v);
				$P->{_fy} -= $f * ($vy/$v);
				$P->{_fz} -= $f * ($vz/$v);
			}
		}
		else {
			foreach my $particle (@{$self->{p}}) {
				my @this_f = $func->({%$particle}, $magnitude);
				$particle->{_fx} += $this_f[0];
				$particle->{_fy} += $this_f[1];
				$particle->{_fz} += $this_f[2];
			}
		}
	}
#	use Data::Dumper;
#	print Dumper $self->{p};
	$self->SUPER::iterate_step(@params);
}

1;
__END__