Die Lösung für ein Überflüssiges Problem: Welcher Computer kann einen Schuss in Echtzeit verfolgen?

Heute mal wieder ein Blog dessen Thema wohl wenige interessieren würde, es ist etwas auf das ich aus Neugier gestoßen bin. Ausgangsbasis war dieses Video der Computer-Chronicles das über den 386-Prozessor geht. Der Moderator steht in der Eingangszene an einem Schießstand und doziert darüber, dass in 1/20 Sekunde die Kugel das Ziel trifft und in dieser Zeit der Prozessor 200.000 Operationen durchführt.

Ich dachte mir: „Der 386 kann wahrscheinlich schneller berechnen, wo genau die Kugel auftrifft, als diese für den Flug braucht“ und das inspirierte mich zu meinem heutigen Blog. Ich dachte nach: Der 386 ist schon recht schnell, ich habe mit einem Computer mit TMS 9900 angefangen, dann mit einem Z80A und dann hatte ich einen 286-er. Ich bin mir relativ sicher, dass ein 286-er diese Aufgabe lösen kann, bestimmt nicht der TMS 9900, denn der Rechner war ziemlich lahm. Aber ein Z80? Vielleicht ein 8086, aber ein Z80 ist mit nur 8 Bit pro Operation, ohne Multiplikations- und Divisionsbefehle doch sehr lahm. Ich wollte es ausprobieren.

Randbedingungen

Zuerst muss man mal definieren, was man berechnen und vor allem wie genau man es berechnen will. Ich habe es mal in Alltagssprache so formuliert: Der Computer soll in der Zeit die angegeben ist (1/20 Sekunde) die Koordinaten und Geschwindigkeit berechnen, wo eine Kugel auftrifft. Diese Koordinaten sollen von den realen Koordinaten (berechnet mit einer höheren Zeitauflösung) maximal so weit abweichen, wie ein guter Schütze seine Schüsse platziert.

Das letzte Kriterium ist dehnbar. Denn natürlich trifft auch der Schütze nicht immer in den Mittelkreis. Ich habe mich mal informiert und Wettkampfschießschieben haben einen Durchmesser von 26 cm und 5 Kreise und die Mitte. Daraus kann man als erforderliche Genauigkeit etwa 2-3 cm ableiten. Als Entfernung habe ich 25 m angenommen, das passt am besten zu der Zeitangabe von 1/20 Sekunde.

Die Physik

Die Physik hinter dem Problem ist relativ einfach. Ich gehe von einem Schießstand aus, das heißt es gibt keinen Seitenwind. Dann wirken nur zwei Kräfte auf die Kugel ein, nachdem sie den Lauf verlassen hat. Die einfachere Kraft ist die Erdbeschleunigung. Sie zieht die Kugel Richtung Erdmittelpunkt. Bei einer so kurzen Strecke und einer geringen Änderung der Entfernung vom Erdmittelpunkt ist dies relativ einfach zu berechnen:

Die Geschwindigkeit Richtung Erdmittelpunkt, die ich im Folgenden als Y-Achse definiere, ist berechenbar nach:

v_y =v_y – g * t

g ist die Erdbeschleunigung von im Mittel 9,80665 m/s². Allerdings ist sie örtlich leicht schwankend, weil die Erde nicht rund ist und zudem noch lokale Vertiefungen und Buckel ist. Am nächsten dem Erdmiteilpunkt ist man am Pol und am weitesten auf einem Berg in Äquatornähe (z.b. dem Kilimandscharo) entfernt. Ich gehe aber von der mittleren Erdbeschleunigung aus. Es wäre kein Problem das Programm so abzuändern das man den geografisch genauen Wert eingibt, wenn man ihn den kennt.

Die Kugel wird durch die Luft abgebremst. Dies ist komplexer. Die Abbremsung hängt von der Geschwindigkeit, der Luftdichte, der Fläche, Masse und einem Luftwiderstandsbeiwert cw ab. Für die Kraft gilt folgende Formel:

F_L = ½ v² * A * ρ * cw

F_L ist die Kraft die der Luftwiderstand ausübt.

v die momentane Geschwindigkeit

A die Fläche der Kugel auf die die Kraft einwirkt

ρ : Die Luftdichte. Auf Meereshöhe sind dies 1,203 kg/m³

cw: Der Luftwiderstandsbeiwert. Ich habe mal geschaut was typische Werte für Kugeln sind und kam auf Werte von 0,25 bis 0,3, ich habe mit 0,275 gerechnet.

Die Abbremsung ist aber auf die Masse bezogen:

v = v- F_l/m

Logisch, ein massiver Körper wird weniger stark abgebremst als ein hohler Körper aus demselben Material und eine Kugel aus Metall stärker als eine Paintballkugel

Die Abbremsung durch den Luftwiderstand ist richtungsunabhängig, zumindest solange es keinen Wind gibt. Die Abbremsung durch die Gravitation wirkt sich aber nur auf die Geschwindigkeit Richtung Erdmittelpunkt aus, das heißt die Geschwindigkeit in Y-Richtung sinkt, bzw. wenn eine Kugel anfangs keine Geschwindigkeit in Y-Richtung hatte (horizontaler Schuss) so bekommt sie nun eine negative Geschwindigkeit, was sich darin äußert, das sie fällt und unterhalb des Zielpunktes aufschlägt.

Das heißt wir müssen die Geschwindigkeit in zwei Komponenten aufteilen, eine horizontale (Index x) und eine vertikale (Index y). Nach den Gesetzen der Mathematik gilt:

v_x = v * cos(winkel)

v_y = v * cos(winkel)

Der Winkel ist der Schusswinkel relativ zum Horizont also einer Linie parallel zur Erdoberfläche. Bei einem Horizontalen Schuss ist er 0 Grad. Für eine längere Strecke wird der Schütze die Erdanziehung kompensieren und mit einem Winkel größer Null leicht „nach oben“ schießen.

Aus der Physik wissen wir auch:

v = Wurzel(v_x² + v_y²)

Was letztlich auf dem Satz des Pythagoras beruht.

Die Routine selbst berechnet nun einfach in einer Schleife ausgehend vom Weg 0 und den Geschwindigkeiten in x und y-Richtung (Vx und Vy) die Änderungen (Abbremsung durch den Luftwiderstand und Beschleunigung durch die Gravitation in einem kleinen Zeitintervall. Je kleiner es ist, desto genauer ist das Ergebnis. Das Herausfinden der Schrittgröße dt für die geforderte Genauigkeit von < 3 cm auf 50 m ist auch eine Aufgabe, die man durch Probieren herausfindet. Durch Addition der Veränderungen bekommt man eine neue Geschwindigkeit für den nächsten Durchlauf und analog addiert man die Geschwindigkeiten um den zurückgelegten Weg zu erhalten.

Integrieren

Nun wird die Kugel durch den Luftwiderstand immer langsamer, gleichzeitig nimmt die Y-Geschwindigkeit durch die Anziehung der Erde zu. Wäre die Abbremsung nicht, so könnte man den Weg relativ einfach nach der bekannten Formel S = ½ v-a*t² berechnen. So ist der Weg aber ein Integral einer komplexen Formel. Numerisch gibt es einige Verfahren dies zu lösen. Ich habe das einfachste genommen, das der eine oder andere vielleicht noch vom Mathematikunterricht kennt. Wenn man da die Fläche unter einer Kurve – die Kurve ist in unserem falle die Flugbahn und die Fläche ist der Weg – bestimmen will, so wendet man einfach die obigen Formeln an und teilt die Berechnung in sehr kleine Intervalle auf, deren Summe dann das Integral ist. In der Physik ist das Integral der Geschwindigkeit der Weg.

Mathematisch gilt:

Integral von F(x) = F(x)+F(x+dx) / 2 *dx

In unserem Fall:

Weg = V_x(t) + V_x(t+dt) / 2 * dt

Real läuft dies so ab: Wir rechnen in einer Schleife die aktuelle Geschwindigkeit aus – wobei wir bei jedem Durchgang den Luftwiderstand und die Anziehung der Erde neu berechnen, basierend auf den alten Werten. Die Summe der Geschwindigkeiten in x Richtung ist dann der zurückgelegte weg. Wir stoppen, wenn wir die Zieldistanz erreicht haben.

Optimieren

Da wir nun von einem viel langsameren Computer als heute reden, will ich einige Worte über das Optimieren verlieren. Das naheliegendste ist es Operationen zusammenzufassen. In der Berechnung des Luftwiderstandes gibt es nur eine variable Größe: Die Geschwindigkeit. Alle anderen Größen können einmal berechnet werden und als eine Konstante betrachtet werden inklusive des Zeitschlitzes dt.

Bei der Gravitationskonstante war ich verleitet sie gar nicht erst in der Schleife zu berechnen. Im Video war von einer Flugzeit von ¹/₂₀ Sekunde die Rede. Praktisch dürfte nach den Gesetzen der Physik die Kugel so maximal 1,22 cm tief fallen, was unterhalb der 3 cm Grenze ist. Das hätte einen sehr großen Vorteil, nämlich das man die Geschwindigkeit nicht in x und y Komponente aufsplitten müsste, was den Rechenaufwand auf ein Drittel reduziert.

Aber es soll ja physikalisch korrekt sein und auch noch korrekte Ergebnisse geben, wenn man einen Präzisionsschuss über eine sehr viel längere Strecke berechnet.

Eine andere Optimierung liegt im Weg und der neuen Geschwindigkeit. Wir errechnen ja die Geschwindigkeitsänderung in x und y Richtung. Es bietet sich an diese zuerst mit 0,5 zu multiplizieren, denn dann reduziert sich die Formel für den Weg auf:

S_x = S_x +(v_x + 0.5 dx)

V_x neu = V_x + dx

Wir müssen, wenn wir den Faktor 2 in die Berechnung einbeziehen dann nur die Änderung dx einmal beim Weg und zweimal bei der Berechnung von Vx neu einbeziehen. Der Faktor 2 steht oben unterm Bruchstrich. Aber Divisionen brauchen bei Rechnern ohne Fließkommaprozessor, also alles vor dem 486 deutlich länger als Multiplikation und die doppelte Addition geht noch schneller als die Multiplikation.

Mein Programm hat beide Schleifentypen integriert. Also die ausformulierte und die optimierte Form. Daneben läuft als Vergleich noch eine Berechnung mit hoher Zeitauflösung von 1 Mikrosekunde (Lazarus) bzw. 100 Mikrosekunden (DOS/CP/M). Das ist nötig, weil die Plattformen langsamer sind und auch keinen Coprozessor haben. Man muss zudem bedenken das je kleiner ich den Zeitschritt mache um so mehr spielt die Genauigkeit eine Rolle, weil die Änderungen dx,dy immer kleiner werden nach wie vor aber in den Termen die normalen Geschwindigkeiten und der Weg auftauchen. Der alte Datentyp Real hat 12 Stellen Genauigkeit, dagegen kann man bei Lazarus mit Extended als Datentyp mit bis zu 20 Stellen Genauigkeit rechnen.

Ich nehme normalerweise sprechende Variablen. Das wir hier Kürzel wie „sy“ finden, hat mit der Physik zu tun. Die Größen für den Weg sind „S“ und für die Geschwindigkeit“v“. Der zweite Index gibt dann die Komponente (x,y) an.

Tests

Die Frage ist natürlich, wie man das Programm nun testet. Ich habe zwei Z80 Computer – genauer gesagt, einen Rechner mit einem Z180 und einem mit einem eZ80. Zumindest beim Z180 kann man relativ einfach die Geschwindigkeit auf einen Z80 extrapolieren, beim eZ80 ist das deutlich schwerer, weil er intern stark optimiert ist. Aber beide haben nur Echtzeituhren mit nur 1 Sekunde Zeitauflösung. Das ist also nicht zielführend.

DOX Box ist als Ziel für einen 8086 oder 80286 sinnvoll. Aber die Dokumentation von DOXBox weist, aus das man einstellt, wie viele Instruktionen pro Sekunde abgearbeitet werden. Nun brauchen die Instruktionen aber unterschiedlich lange zur Ausführung und gerade die hier benötigten Multiplikationen und Divisionen brauchen etwa beim 8086 10-mal mehr Zeit als eine Addition/Subtraktion. Es ist aber die einzige Emulatorplattform die ich habe.

Ich habe das Programm daher auf einem Emulator für den Amstrad CPC unter CP/M laufen lassen. Der Emulator muss akkurat sein, denn sonst wären Spiele nicht spielbar. Das liefert einem einen Anhaltspunkt für die minimale Performance und man kann aus Benchmarkvergleichen wie dem ct’ Highlevel Benchmark dann auf andere Prozessoren schließen. Der Vorteil ist natürlich auch das ich es in Pascal entwickeln kann und nur zum Compilieren auf die Zielplattform gehen muss. Der Z80A im Armstrad CPC läuft nominell mit 4 MHz. Allerdings muss er den Speicherzugriff mit einem ASIC Teilen, das die 16.000 Bytes des Bildschirmspeichers 50 bzw. 60-mal pro Sekunde an den Monitor überträgt. Man nutzte aus, das der Speicherzugriff bei einem normalen Z80 immer in den ersten zwei von vier Takten einen Maschinenzyklus erfolgte. Allerdings mussten dann alle Befehle immer Vielfache von 4 Zyklen haben, auch wenn dies in der Realität nicht so war. So lief er effektiv nur mit 3,2 MHz. Dafür wurde ein Waitimpuls an die Z80 angelegt, wenn das Gate Array auf den Speicher zugriff.

Die Ergebnisse

Aus Bequemlichkeitsgründen habe ich die Masse einer Patrone und ihren Durchmesser schon im Programm verankert. Wer will, kann diese leicht eingeben lassen. Sie entsprechen den Wikipediaangaben für die Standardpatrone 9 x 19 mm Parabellum.

Also moderne PC’s sind beim test außen vor. Ich habe bei einem CPC Emulator es mal getestet und dann unter DOS-BOX mit der Einstellung eines 8088 mit 4,77 MHz also des Prozessors, der im IBM PC verbaut war. Tests unter Lazarus ergaben das bei einem Zeitschlitz von 5 ms (0.005 als Eingabe) die Abweichung in der Y-Richtung kleiner als 3 cm ist.

Prozessor	MHz	Laufzeit 0.005 s	Laufzeit 0.0001 s
Z80	3,2	0,26 s / 0.67	11,3 s / 30,65
8086	4,77 MHz	0,17 s / 0,39 s	6,26 s /17,25 s

 

Ein Tipp: Damit ihr nicht dauernd alles neu eingeben müsst, könnt ihr die Parameter auf der Kommandozeile eingeben:

Schuss 365 0.005 25 0

Wäre die Kombination, die ich für Tests nutze (365 m/s, 0.005 s Zeitauflösung 25 m Distanz Winkel 0 Grad). Würde man nun ein Gewehr simulieren, das eine höhere Mündungsgeschwindigkeit hat, so wäre die Abweichung kleiner aber auch die Flugzeit geringer. Da die 25 m in 0,071 s zurückgelegt werden – in etwa passend zur 1/20 Sekunde in dem Video – sind beide Berechnungen langsamer als die Realität. Die Optimierung bringt aber etwas, sie halbiert in etwa die Laufzeit. Ich denke ein Z80 mit 12 MHz – den gab es in den Achtzigern nicht oder ein 8086 mit 8 MHz würden Realzeit errechnen. Alle Prozessoren ab dem 80286 wären dann schnell genug um den Einschlagspunkt schneller zu berechnen als die Kugel fliegt.

Im Zip Archiv findet ihr die Files:

„Schuss.pas“ für die Dos-BOX (Compiler: Turbo Pascal 3)

„Schuss2.pas“ für einen CPC-Emulator (Compiler: Turbo Pascal 3)

„timecpc.inl“ enthält den Maschinencode für die Zeitabfrage für den CPC als InlineCode. Wichtig: es gibt zwei Routinen für den CPC 464 und 6128, in Schuss2.pas wird derzeit die für den 6128 (CP/M 3) verwendet.

„Schuss.lpi, .lpr, .lps: Projektdateien für Lazarus und als .exe

Wie viel Overhead dazukommt, merkt man auch daran, dass die compilierte .exe unter Lazarus 415 KByte groß ist, unter TP 3 noch 11 KByte.