しみゅれーしょん☆ぶつり 10

10.pl

use 5.14.0;

my $G   = 1;
my $H   = 3.5;
my $M   = 1;
my $org = 90;
my $dx  = 0.5;
my $dt  = 0.1;
my $nm  = 99;
my $ni  = 50;
my $nk  = 41;

my @m;
$m[0] = [];

while (@{$m[0]} <= 500) {
	my $x = rand(10) - 5;
	my $y = rand(10) - 5;

	push @{$m[0]}, {x => $x+$org, y => $y+$org, vx => $x*$H, vy => $y*$H } if $x*$x + $y*$y <= 5;
}

my @F;
my @ro;
for my $i (0..$nk) {
	$ro[$_->{x}+0.5][$_->{y}+0.5]++ for @{$m[$i]};

	my @phi = map { [0..$nm+1] } 0..$nm+1;
	for (0..$ni) {
		for my $ix (1..$nm) {
			for my $iy (1..$nm) {
				my $p1 = $phi[$ix+1][$iy] + $phi[$ix-1][$iy] + $phi[$ix][$iy+1] + $phi[$ix][$iy-1];
				my $p2 = $G * int $ro[$ix][$iy] * $dx * $dx;
				$phi[$ix][$iy] = ($p1 / 4) - ($p2 / 4);
			}
		}
	}

	for my $ix (1..$nm) {
		for my $iy (1..$nm) {
			$F[$ix][$iy] = {
				x => -1 * ($phi[$ix+1][$iy] - $phi[$ix][$iy]) / $dx,
				y => -1 * ($phi[$ix][$iy+1] - $phi[$ix][$iy]) / $dx,
			};
		}
	}

	my ($Fp, $x, $y, $vx, $vy);
	for my $p (@{$m[$i]}) {
		$Fp = {
			x => $M * $F[$p->{x}+0.5][$p->{y}+0.5]->{x},
			y => $M * $F[$p->{x}+0.5][$p->{y}+0.5]->{y},
		};
		$vx = $p->{vx} + ($Fp->{x}/$M) * $dt;
		$vy = $p->{vy} + ($Fp->{y}/$M) * $dt;
		$x  = $p->{x}  + $vx * $dt;
		$y  = $p->{y}  + $vy * $dt;

		push @{$m[$i+1]}, {x => $x, y => $y, vx => $vx, vy => $vy};
	}
}

for (0..2) {
	open my $fh, '>', "out$_.csv";
	say $fh "$_->{x},$_->{y}" for @{$m[$_*20]};
	close $fh;
}