Monte-Carlo-Simulation Eierwurf

eierwurf_mcs.pl

#!/usr/bin/perl -w

# ===========================================
#   Laden von Perl-Modulen
# ===========================================

# Trigonometrische Funktionen, Kreiszahl Pi
use Math::Trig;

# Normalverteilte Zufallszahlen
use Math::Random;

# ===========================================
#   Definition der Programm-Parameter
# ===========================================

# Erdbeschleunigung in m s^-2
$g=9.81;

# Dichte Luft in kg m^-3
my $rho_luft=1.2;

# Masse Ei in kg
my $m_ei=0.06;

# Durchmesser des Eis in m
my $d_ei=0.05;

# Aerodynamischer Widerstandsbeiwert des
# Eis
my $c_ei=0.3;

# Breite des Gebaeudes in m
$geb_breite=20;

# Hoehe des Gebaeudes in m
$geb_hoehe=15;

# Schrittweite fuer Runge-Kutta-Integration
# in m s^-1
$h=0.01;

# Max. Iterationen bei Runge-Kutta-Int.
$max_iter_rk=10000;

# Anz. Wuerfe fuer Monte-Carlo-Simulation
my $max_iter_mcs=100000;

# Standardabweichung fuer Abwurfort in m
my $sigma_x0=20;

# Standardabweichung fuer Abwurfgeschw. in
# m s^-1
my $sigma_v0=20;

# Min- und Max-Wert fuer Abwurfwinkel
my $min_beta0=0;
my $max_beta0=pi/2;

# Dateiname fuer Ausgabe der Koordinaten
# der optimalen Bahn
my $outfile="eierwurf_mcs_opt_bahn.dat";

# ===========================================
#   Berechnung von Hilfsgroessen
# ===========================================

# Querschnittsflaeche des Eis in m^2
my $A=0.25*pi*$d_ei**2;

# Hilfsgroesse in m s^-2
$k=$rho_luft*$c_ei*$A/(2*$m_ei);

# Kleinste sinnvolle Abwurfgeschw., damit ein
# vertikal geworfenes Ei ueberhaupt die
# Hoehe des Gebaeudes erreicht
my $v0_min=sqrt(2*$g*$geb_hoehe);

# ===========================================
#   Initialisierungen
# ===========================================

# Initialisierung des Random-Generators mit
# Zeitstring
my $jetzt=localtime time;
random_set_seed_from_phrase($jetzt);

# Array, indem die Parameter der Bahnen
# abgelegt werden, die komplett ueber das
# Gebaeude fuehren
my @bahnen_ueber_geb;

# Laufindex
my $i=0;

# Startzeit
my $start_zeit=localtime time;

# ===========================================
#   Start Hauptprogramm
# ===========================================

# Hauptschleife
while ($i < $max_iter_mcs)
{
	# Parameter fuer Wurfbahn bestimmen
	# Abwurfort x
	my $x0=-abs(random_normal(1,0,$sigma_x0));
	# Abwurfort y (immer Null also vom Boden)
	my $y0=0;
	# Abwurfgeschwindigkeit
	my $v0=$v0_min
		+abs(random_normal(1,0,$sigma_v0));
	# Abwurfwinkel
	my $beta0=random_uniform(1,
		$min_beta0,$max_beta0);

	# Wurfbahn berechnen
	my ($n_rk,
		$ref_x_arr,
		$ref_y_arr,
		$ref_vx_arr,
		$ref_vy_arr)
		=wurfbahn($x0,$y0,$v0,$beta0);

	#Dereferenzieren der uebergebenen Arrays
	my @x=@{$ref_x_arr};
	my @y=@{$ref_y_arr};
	my @vx=@{$ref_vx_arr};
	my @vy=@{$ref_vy_arr};

	# Pruefen, ob die Bahn komplett ueber das
	# Gebaeude fuehrt
	my $drueber=wurfbahn_ueber_geb(\@x,\@y);

	print "$i $x0 $v0 $beta0 $n_rk $drueber\n";
	#ausg_koord_wurfbahn(\@x,\@y);

	# Wenn die Wurfbahn ueber das Gebaeude geht,
	# dann entsprechende Parameter abspeichern
	if ($drueber==1)
	{
		my $par={"x0"=>$x0,
			"v0"=>$v0,
			"beta0"=>$beta0};
		push @bahnen_ueber_geb, $par;
	}

	$i++;
}

# Endzeit
my $end_zeit=localtime time;

# Parameter-Eintraege in Array aufsteigend nach v0
# sortieren. Danach befindet sich die energetisch
# guenstigste Bahn im ersten Element
@bahnen_ueber_geb
	=sort{$a->{'v0'}<=>$b->{'v0'}}
	@bahnen_ueber_geb;

# Parameter zur minimalsten Bahn ausgeben
my $x0_opt=$bahnen_ueber_geb[0]->{'x0'};
my $v0_opt=$bahnen_ueber_geb[0]->{'v0'};
my $beta0_opt=
	$bahnen_ueber_geb[0]->{'beta0'};

print "\nParameter der energetisch guenstigsten Bahn:\n";
print "$x0_opt $v0_opt $beta0_opt\n";

# Ein bisschen Statistik
print "\nIteration MCS: $max_iter_mcs\n";
my $anz_bahnen_ueber_geb=$#bahnen_ueber_geb+1;
print "Anz. Bahnen ueber Geb.: $#bahnen_ueber_geb\n";
print "Startzeit: $start_zeit\n";
print "Endzeit: $end_zeit\n";

# Koordinaten der guenstigsten Bahn berechnen
my ($n_rk,$ref_x_arr,$ref_y_arr,
	$ref_vx_arr,$ref_vy_arr)
	=wurfbahn($x0_opt,0,$v0_opt,
	$beta0_opt);

#Dereferenzieren der uebergebenen Arrays
my @x=@{$ref_x_arr};
my @y=@{$ref_y_arr};
my @vx=@{$ref_vx_arr};
my @vy=@{$ref_vy_arr};

# Koordinaten der guenstigtsen Bahn in Datei
# ausgeben
ausg_koord_wurfbahn(\@x,\@y,\@vx,\@vy,
	$outfile);

#Programmende
print "\n\nProgrammende\n";

# ===========================================
#   Unterroutinen
# ===========================================

sub v_rk4
# Runge-Kutta-Routine zur Integration der
# Bewegungsgleichungen
{
	# Akt. Geschwindigkeit in x
	my $vx_n=shift;
	# Akt. Geschwindigkeit in y
	my $vy_n=shift;

	# 1. Schritt
	my $v_n=sqrt($vx_n**2+$vy_n**2);
	my $abl_vx_n=-$k*$v_n*$vx_n;
	my $vx_a=$vx_n+0.5*$h*$abl_vx_n;
	my $abl_vy_n=-($g+$k*$v_n*$vy_n);
	my $vy_a=$vy_n+0.5*$h*$abl_vy_n;

	# 2. Schritt
	my $v_a=sqrt($vx_a**2+$vy_a**2);
	my $abl_vx_a=-$k*$v_a*$vx_a;
	my $vx_b=$vx_n+0.5*$h*$abl_vx_a;
	my $abl_vy_a=-($g+$k*$v_a*$vy_a);
	my $vy_b=$vy_n+0.5*$h*$abl_vy_a;

	# 3. Schritt
	my $v_b=sqrt($vx_b**2+$vy_b**2);
	my $abl_vx_b=-$k*$v_b*$vx_b;
	my $vx_c=$vx_n+$h*$abl_vx_b;
	my $abl_vy_b=-($g+$k*$v_b*$vy_b);
	my $vy_c=$vy_n+$h*$abl_vy_b;

	# 4. Schritt
	my $v_c=sqrt($vx_c**2+$vy_c**2);
	my $abl_vx_c=-$k*$v_c*$vx_c;
	my $abl_vy_c=-($g+$k*$v_c*$vy_c);
	my $vx_m=$vx_n+$h*1/6*($abl_vx_n
		+2*($abl_vx_a+$abl_vx_b)
		+$abl_vx_c);
	my $vy_m=$vy_n+$h*1/6*($abl_vy_n
		+2*($abl_vy_a+$abl_vy_b)
		+$abl_vy_c);

	# Neue Geschwindigkeitswerte fuer x
	# und y zurueckgeben
	return ($vx_m, $vy_m);
}

sub wurfbahn
# Berechnung der Wurfbahn fuer einen gegebenen
# Satz von Parametern
{
	# Startwert fuer x
	my $x0=shift;
	# Startwert fuer y
	my $y0=shift;
	# Abwurfgeschwindigkeit
	my $v0=shift;
	# Abwurfwinkel
	my $beta0=shift;

	# Initialisierungen
	my (@x_arr, @y_arr, @vx_arr, @vy_arr);
	push @x_arr, $x0;
	push @y_arr, $y0;
	push @vx_arr, $v0*cos($beta0);
	push @vy_arr, $v0*sin($beta0);
	my $i=1;
	my $x=$x0;
	my $y=$y0;
	my $vx=$v0*cos($beta0);
	my $vy=$v0*sin($beta0);

	# Iterative Runge-Kutta-Integration
	while ($y>=0 && $i <= $max_iter_rk)
	{
		($vx, $vy)=v_rk4($vx,$vy,$h,$k,$g);
		$x=$x+$h*$vx;
		$y=$y+$h*$vy;

		# Neue Werte in Arrays ablegen
		push @x_arr, $x;
		push @y_arr, $y;
		push @vx_arr, $vx;
		push @vy_arr, $vy;

		$i++;
	}

	# Arrays zurueckgeben
	return ($i,\@x_arr,\@y_arr,\@vx_arr,\@vy_arr); 
}

sub wurfbahn_ueber_geb
# Ueberpruefen, ob die Wurfbahn auch komplett
# ueber das Gebaeude geht
{
	# x-Koord. Wurfbahn
	my $ref_x_arr=shift;
	my @x_arr=@{$ref_x_arr};
	# y-Koord. Wurfbahn
	my $ref_y_arr=shift;
	my @y_arr=@{$ref_y_arr};

	# Arraygroesse bestimmen
	my $n_arr=$#x_arr;

	my $i=0;
	my $drueber=1;
	# Schleife ueber Bahnpunkte
	while ($i <= $n_arr && $drueber == 1)
	{
		# Fuer jeden Bahnpunkt, dessen Lotpunkt
		# im Gebaeude liegt pruefen, ob er
		# hoeher als das Gebaeude ist. Wenn
		# nicht, setze $drueber zu 1, wodurch
		# die Schleife im naechsten Durchlauf
		# abbricht.
		$drueber=0 if ($x_arr[$i]>=0
			&& $x_arr[$i]<=$geb_breite
			&& $y_arr[$i]<$geb_hoehe);
		$i++;
	}

	# Jetzt muss noch geprueft werden, dass die
	# Wurfbahn ueberhaupt ueber die zweite Dach-
	# kante geht!
	my $x_max=(sort {$b <=> $a} @x_arr)[0];
	$drueber=0 if ($x_max <= $geb_breite);

	return $drueber;
}

sub ausg_koord_wurfbahn
# Ausgabe der Koordinaten der Wurfbahn am
# Bildschirm
{
	# x-Koord. Wurfbahn
	my $ref_x_arr=shift;
	my @x_arr=@{$ref_x_arr};
	# y-Koord. Wurfbahn
	my $ref_y_arr=shift;
	my @y_arr=@{$ref_y_arr};
	# vx Wurfbahn
	my $ref_vx_arr=shift;
	my @vx_arr=@{$ref_vx_arr};
	# vy Wurfbahn
	my $ref_vy_arr=shift;
	my @vy_arr=@{$ref_vy_arr};
	# Ausgabedatei
	my $outfile=shift;

	# Arraygroesse bestimmen
	my $n=$#x_arr;

	# Oeffnen der Ausgabedatei
	open FILE, ">$outfile"
		or die "Konnte $outfile nicht oeffnen: $!\n";

	my $i=0;
	while ($i <= $n)
	{
		# Berechne den momentanen Flugwinkel
		my $beta=asin($vy_arr[$i]
			/sqrt($vx_arr[$i]**2
			+$vy_arr[$i]**2));
		
		# Ausgabe: Zaehler, x, y, vx, vy, beta 
		print FILE "$i $x_arr[$i] $y_arr[$i] $vx_arr[$i] $vy_arr[$i] $beta\n";
		$i++;
	}

	# Schliessen der Ausgabedatei
	close FILE;
}