Physikalisch: Coulombkraft

Herr-Vorragend

Hi,

im Rahmen einer Facharbeit über das Thema Kernfusion muss ich eine Simulation programmieren. In dieser sollen sich zwei Protonen nähern. Nun versuche ich die Abstoßung zu simulieren. Es trat gleich ein Fehler auf, worauf hin ich eine sehr abgespeckte Version programmiert habe, bei der der Fehler immernoch auftrat. Hier der Code (lediglich die Fensterprozedur, der Rest ist 0-8-15 Fenstererstellung):

LRESULT CALLBACK WndProc(HWND hWnd, UINT message, WPARAM wParam, LPARAM lParam)
{
	HDC				hdc;
	PAINTSTRUCT		ps;
	static UINT		offsetX, offsetY;
	static DOUBLE	bulletX, bulletY;
	static UINT		targetX, targetY;
	static DOUBLE	speedX, speedY;
	static DOUBLE	t=0;

	static DOUBLE	alpha;
	static DOUBLE	dx, dy;
	static DOUBLE	r;
	static DOUBLE	F, Fx, Fy;

	switch (message)
	{
	case WM_CREATE:
		offsetX = 50;
		offsetY = 50;

		bulletX = offsetX;
		bulletY = offsetY;

		targetX = 700;
		targetY = 200;

		speedX = 3;
		speedY = 1;

		SetTimer(hWnd, 1, 20, NULL);

		return 0;

	case WM_TIMER:

		dx = targetX - bulletX;
		dy = targetY - bulletY;

		alpha = atan(dy/dx);

		r = sqrt(dx*dx + dy*dy);

		F = 1/(r*r);

		Fx = cos(alpha)*F;
		Fy = sin(alpha)*F;

		bulletX = offsetX + speedX*t	- (0.5 * (Fx/0.001) * (t*t));
		bulletY = offsetY + speedY*t	- (0.5 * (Fy/0.001) * (t*t));
		t++;

		InvalidateRect(hWnd, NULL, TRUE);

		return 0;

	case WM_PAINT:

		hdc = BeginPaint(hWnd, &ps);

		SelectObject(hdc, CreateSolidBrush(RGB(0,0,255)));
		Ellipse(hdc, targetX-5, targetY-5, targetX+5, targetY+5);

		SelectObject(hdc, CreateSolidBrush(RGB(255,0,0)));
		Ellipse(hdc, offsetX + (int)bulletX - 5, offsetY + (int)bulletY - 5, offsetX + (int)bulletX + 5, offsetY + (int)bulletY + 5);

		EndPaint(hWnd, &ps);

		return 0;

	case WM_DESTROY:
		PostQuitMessage(0);
		return 0;
	}

	return DefWindowProc(hWnd, message, wParam, lParam);
}

Das wichtigste steckt in dem WM_TIMER Case-Zweig.
Zuerst wird der horizontale und der vertikale Abstand zwischen den Protonen errechnet (dx und dy). Damit wird der Winkel alpha mittels arcus-tangens und die Hypothenuse r des rechtwinkligen Dreiecks mittels Pythagoras errechnet.
Die Kraft F berechnet sich durch Epsilon*Q1*Q2/r², das hab ich nun mal abgespeckt und durch 1/(r*r) ersetzt. Die vertikale und horizontale Komponente der Kraft errechnet sich nun mit sin und cos des Winkels. Also ergibt sich die neue x und y-koordinate des geschosses durch die Ausgangsposition offsetX/Y plus eine gleichförmige Bewegung durch die Anfangsgeschwindigkeit speedX/speedY minus einer beschleunigten Bewegung, die sich durch 1/2 at² = 1/2 F/m * t² errechnen lässt. Für m hab ich einfach einen passenden Wert eingesetzt.
Ich hoffe soweit ist alles verständlich und korrekt.
Allerdings wird jeder sehen, dass wenn man das ganze testet irgendwas verrückt spielt, sobald das geschoss vollständig abgebremst wurde und umkehren soll. Nun frage ich mich woran das liegt. Irgendwie scheint ein Vorzeichen ständig umgedreht zu werden, nur weiß ich nicht welches. Ich wüsste nicht, wo ich irgend ein Vorzeichen falsch gesetzt haben soll. Liegt es möglicherweise an den Datentypen, dass es vielleicht irgendwo einen Bufferoverflow gibt o.ä.? Oder hab ich einfach einen riesen Denkfehler drin?

Wäre für jede Hilfe dankbar,

MfG, Herr-Vorragend

Walter Z

Division durch 0, weil dx oder r den Wert null erreicht?

Walter Zorn

Herr-Vorragend

Nein daran kanns nicht liegen, dx wird nie 0, sieht man aber auch nur, wenn man sich das ganze kompiliert anschaut. Somit wird r dann auch nie 0. Btw. würde dann das Programm abstürzen und ne Fehlermeldung kommen, oder?

Was hat das mit WinApi zu tun?

Herr-Vorragend

Ich wusste nicht genau wo ich es hinposten soll, hab dann WinAPI genommen, weil das Programm eben mit der WinAPI geschrieben ist und evt. dort irgendwo der Fehler vergraben ist...Wenn jemand anderer Meinung ist bitte verschieben

du machst Haufenweise gdi-leaks:

case WM_PAINT:

        hdc = BeginPaint(hWnd, &ps);

        SelectObject(hdc, CreateSolidBrush(RGB(0,0,255))); // BOOOMMM!!!
        Ellipse(hdc, targetX-5, targetY-5, targetX+5, targetY+5);

        SelectObject(hdc, CreateSolidBrush(RGB(255,0,0)));  // BOOOMMM!!!
        Ellipse(hdc, offsetX + (int)bulletX - 5, offsetY + (int)bulletY - 5, offsetX + (int)bulletX + 5, offsetY + (int)bulletY + 5);

        EndPaint(hWnd, &ps);

        return 0;

du muss den Brush auch wieder entfernen: sowohl außem hdc als auch außem Speicher: DeleteObject(), sag ich nur!
das is zwar nit die Lösung des problems, aber die sollte man trotzdem nit machen

Walter Z

Vermutlich wird die Lösung deiner Gleichung mit der Laufzeit infolge Fehlerfortpflanzung (besser: Ungenauigkeitsfortpflanzung) immer instabiler.

Insbesondere Ausdrücke wie (0.5 * (Fx/0.001) * (t*t)), die dann im nächsten Durchlauf wiederum in Quadrierungen, atan() und sqrt() weiterverarbeitet werden und somit auf sich selbst zurückwirken, lassen das jedenfalls erwarten. Nur eine Vermutung...

Walter Zorn

Herr-Vorragend

Hm, ja das ist ne Möglichkeit....wobei der Fehler ja erst so krass auftritt, wenn das Teilchen vollständig abgebremst wurde. Gibt es irgendwelche "genaueren" Datentypen, bei welchen das Runden vermieden wird?

Das Problem mit den Speicherleaks ist mir eigentlich bekannt, ka wieso ich es vergessen hab. Danke für den Tipp.

MfG

Walter Z

Evtl. die Gleichung so umstellen, dass nirgendwo als Zwischenergebnis betragsmäßig sehr kleine (<< 1.0) oder sehr große Werte auftreten? Falls nötig, sogar mehrere Versionen der Gleichung zur Vefügung stellen, die jeweils auf bestimmte Parametergrößen bzw. Zwischenergebnisse hin optimiert sind und entsprechend aufgerufen werden?

Ich kann mir aber irgendwie nicht vorstellen, dass so etwas in der Facharbeit erwartet wird...

Ups, sorry, inzwischen ziemlich offtopic.

Walter Zorn

Herr-Vorragend

Also in die Facharbeit muss ein praktischer Teil mit rein. Da mein Thema nun mal Kernfusion ist und ich schlecht eine Fusion herstellen kann hab ich beschlossen eine Simulation zu programmieren.
Ich hab nun auch mal den Datentyp long double ausprobiert, hilft aber auch nichts. Ich denk mir wird nichts anderes übrig bleiben also wirklich alles zu erweitern bis keine Kommazahlen mehr rauskommen...

Vielen Dank für die Hilfe, ich meld mich nochmal

Herr-Vorragend

Hi,

Ich hoffe ihr habt den Thread noch nicht vergessen

Ich glaube langsam nicht mehr dran, dass es an den Rundungsfehlern liegt...Ich weiß nimmer was ich machen soll, ich habs in Flash (Actionscript) probiert, ich habs in C probiert, ich habs in Pascal probiert und immer mit den aller größten Datentypen (long double/extended)...irgendwas muss an der Formel nicht stimmen schätz ich mal, ich kann mirs echt anders nimmer erklären...Bitte helft mir Ich muss die Facharbeit in 3 Wochen abgeben und da muss alles funktionieren....

MfG

C++ Forumbot

Dieser Thread wurde von Moderator/in Jochen Kalmbach aus dem Forum WinAPI in das Forum Rund um die Programmierung verschoben.

Im Zweifelsfall bitte auch folgende Hinweise beachten:
C/C++ Forum :: FAQ - Sonstiges :: Wohin mit meiner Frage?

Dieses Posting wurde automatisch erzeugt.

hi,

schmeiß mal den winkel raus: du kannst das ganze viel schöner rein vektoriell machen. es kann sein, dass bei kleinen geschwindigkeiten dy/dx nicht mehr viel mit dem "echten" winkel zu tun hat. obs aber wirklich daran liegt kann ich natürlich nicht sagen

Gregor

Edit: Hier stand Mist.

Chris++ 0

Herr-Vorragend schrieb:

static UINT		offsetX, offsetY;
	static UINT		targetX, targetY;

//...

	bulletX = offsetX + speedX*t	- (0.5 * (Fx/0.001) * (t*t));
	bulletY = offsetY + speedY*t	- (0.5 * (Fy/0.001) * (t*t));

Kann es vieleicht sein, das das Ergebnis aus

speedX*t - (0.5 * (Fx/0.001) * (t*t))

nach UINT (OffsetX & OffsetY) konvertiert und somit (vorzeichenlos) gerundet wird?

Probiers mal mit

bulletX = (double)offsetX + speedX*t	- (0.5 * (Fx/0.001) * (t*t));
	bulletY = (double)offsetY + speedY*t	- (0.5 * (Fy/0.001) * (t*t));

oder mit:

bulletX = ( speedX*t	- (0.5 * (Fx/0.001) * (t*t)) ) + (double)offsetX;
	bulletY = ( speedY*t	- (0.5 * (Fy/0.001) * (t*t)) ) + (double)offsetY;

Gregor

Herr-Vorragend schrieb:

case WM_TIMER:

		dx = targetX - bulletX;
		dy = targetY - bulletY;

		alpha = atan(dy/dx);

		r = sqrt(dx*dx + dy*dy);

		F = 1/(r*r);

		Fx = cos(alpha)*F;
		Fy = sin(alpha)*F;

		bulletX = offsetX + speedX*t	- (0.5 * (Fx/0.001) * (t*t));
		bulletY = offsetY + speedY*t	- (0.5 * (Fy/0.001) * (t*t));
		t++;

		InvalidateRect(hWnd, NULL, TRUE);

		return 0;

Ok, mal gucken. Das läuft also schrittweise. In jedem Schritt berechnest Du die momentane relative Position zwischen den beiden Teilchen, dann den momentanen Winkel. dann die momentane Kraft. So, jetzt kommt es Du zählst die Zeit immer weiter hoch. Das heißt, Du lässt die Kraft, die du für einen bestimmten Zeitpunkt berechnet hast, über die gesamte Zeit wirken. Dein Algorithmus soll offensichtlich plötzlich nichtmehr schrittweise funktionieren, sondern stattdessen willst Du einfach die Zeit in eine vorbereitete Formel einsetzen. Da vermischt Du zwei Dinge, die IMHO so nicht zu vermischen sind. Entweder Du musst die Trajektorie komplett vorher berechnen und dann nur noch die Zeit einsetzen oder Du musst komplett schrittweise arbeiten.

Und bitte nicht hau deine Algorithmen nicht in den GUI-Code. Das ist ja schrecklich.

Herr-Vorragend

@Gregor: Ich glaub du hast wirklich Recht, ich hab die Sache mit der Kraft zu sehr vereinfacht...Am einfachsten wär es glaub ich zu der Position ständig etwas hinzuzuaddieren, also so ne Art Momentangeschwindigkeit*Zeitintervall...Nur wie kann ich in diesem Fall die Momentangeschwindigkeit errechnen? Die hängt ja auch von der Kraft ab...KAnn mir da jemand helfen?

Schonmal Tausend Dank für die Tipps!

Herr-Vorragend

Ich hätte da einen Ansatz, aber damit auch ein Problem

Also, ich hab mir gedacht, ich mach einfach eine Endlosschleife und addiere in jedem Durchlauf eine Strecke zur Position des Protons hinzu. Diese ergibt sich durch s = v*t wobei v die Momentangeschwindigkeit ist, die bei jedem Schleifendurchlauf aufgrund des sich ändernden Abstandes variiert. t ist das "Aktualisierungsintervall". Nur wie komm ich an die Momentangeschwindigkeit?

Dazu hab ich mir gedacht, dass das Proton im unendlichen eine kinetische Energie E_kin hat und diese nach und nach in potentielle (elektrische) Energie umgewandelt wird. Also dann

$E_{pot} + E_{kin}^{'} = E_{kin}$

Nach Einsetzen und Auflösen nach v in E_kin' ergibt sich dann:

v = \sqrt{\frac{2(E_{kin} - \frac{Q_{1}\*Q_{2}}{4\*\pi*\epsilon_{0}*r})}{m}}

Ich hoffe da steckt bisher kein Fehler drin und der Ansatz ist richtig...Wenn nicht bitte melden!

Nun ist das Problem, das das ganze zwar schön abgebremst wird, aber das Proton bewegt sich wenn es stehen geblieben ist nicht in die andere Richtung bleibt einfach wo es ist, was ja klar ist, weil E_pot = E_kin und sich das so schnell mit dieser Formel nicht ändert...

Steckt jetzt ein Fehler in meinem Ansatz oder muss ich nur was ergänzen?

MfG, Herr-Vorragend (für jede Hilfe dankbar)

hi,

grob gesagt fehlt die trägheit des teilchens. die 0-8-15-lösung sieht so aus:

F_neu = q₁ q2_q / ( bla* r²)
a_neu = F_neu/m
v_neu = v_alt + Δta_neu
x_neu = x_alt + Δtv_neu

wobei f,a,v,x vektorielle größen sind und jeweils eine x- und eine y-komponente haben. (die formeln gelten für jede komponente einzeln, bis auf die erste: da musst du für r pythagoras benutzen.)

@gregor:
wieso mist? ich dachte, das, was du gepostet hast, benutzt man auch für numerische integration