Wie kann ich ermitteln, ob sich ein Punkt (die Sonne) innerhalb eines Polygons auf einer Kugeloberfläche befindet?

Question

Wie kann ich ermitteln, ob sich ein Punkt (die Sonne) innerhalb eines Polygons auf einer Kugeloberfläche befindet?

1 Antwort

Beste Antwort

Hallo,

Willkommen in der Mathelounge!

Wenn ein geschlossenes Polygon sich auf einer endlichen(!) Kugeloberfläche befindet, so teilt es diese endliche Fläche auch nur in zwei endliche Flächen. D.h. man muss zunächst definieren, welche von beiden innen und welche außen sein soll. Weiter beschränke ich mich auf ein Polygon, welches sich nicht selbst schneidet (oder berührt).

Mal angenommen, dass die Fläche, die mathematisch positiv umrundet wird, innen ist. Dann wähle einen beliebiegen Großkreis durch den zu prüfenden Punkt $P$ und wähle eine beliebige Richtung auf diesem Großkreis. Dann berechne vom Punkt ausgehend in der gewählten Richtung den nächsten Schnittpunkt $S$ mit dem Polygon.

Anschließend berechne das Kreuzprodukt aus der lokalen Richtung des Großkreises und der lokalen Richtung des Polygons in $S$. Zeigt das Kreuzprodukt vom Mittelpunkt der Kugel weg, so liegt $P$ innerhalb des Polygons.

Falls es keinen Schnittpunkt gibt oder unendlich viele, so wähle einen neuen Großkreis.

Falls Du dazu noch Fragen hast, so melde Dich bitte.

Gruß Werner

Beantwortet 16 Apr von Werner-Salomon

Hallo Werner,

Zunächst erst einmal ganz vielen lieben Dank, dass du dich mit meinem Problem beschäftigt hast!

Im konkreten Fall geht es um die Programmierung einer Verschattung von Fenstern. Bei zu verschattenden Fenstern sind eigentlich nur die beiden unteren Ecken interessant. Wenn der Schatten dort angelangt ist, muss die Verschattung ausgelöst oder beendet werden. Beide Ecken nehmen wir der Einfachheit halber mit Koordinaten Azimut/Elevation von 0/0 an. Dann nehmen wir einen Winkelmesser und visieren von beiden Fenstereckpunken die Eckpunkte des Polygons an, in welchem die Sonne stehen muss, damit sie auf den jeweiligen Eckpunkt scheint. Die jeweils größeren Werte nehmen wir. Da bekommen wir einen Datensatz von Punkten Azimut/Elevation.

Dieser Datensatz beschreibt jetzt eine endliche Fläche am Himmel, die nicht komplex gestaltet ist (keine Berührung oder Schneidung etc.) auf der Innenseite der Kugel, auf der die Sonne entlang wandert (deren Azimut und Elevation man jetzt auf den betrachteten Teil des Himmels umrechnen muss). Befindet sich die Sonne innerhalb dieser Fläche, dann soll die Verschattung ausgelöst werden.

Im Prinzip ist die Aufgabenstellung vermutlich nicht viel anders als der PiP Test nach Jordan, nur eben auf einer Kugelinnenfläche. Ich vermute, dass dein Lösungsansatz genau das widerspiegelt?

Nun muss ich gestehen, dass meine letzte Mathematikvorlesung im Ingenieurstudium schon über 35 Jahre zurück liegt. Um deinen Ansatz in ein funktionierendes Programm umzusetzen fehlen mir daher noch ein paar Zwischenschritte... Könntest du mir einen groben Ablaufplan mit den entsprechenden mathematischen Formeln skizzieren? Das würde mir unendlich helfen.

Vielen lieben Dank vorab schon einmal!

Kommentiert 16 Apr von Guitardoc

Ich habe mal angefangen die Funktion zu programmieren (in Javascript). tempPolygon beinhaltet dabei schon die aus den Winkeln umgerechneten Werte in kartesischen Koordinaten (x=-cos(h)*cos(A), y=cos(h)*sin(A), z=sin(h))

async function Punkt_in_Polygon(tempPolygon, Sonne_x, Sonne_y, Sonne_z) {
PunktInnerhalb = false;
var i_end = tempPolygon.length;
var i_inc = 1;
if (1 > i_end) {
i_inc = -i_inc;
}
for (i = 1; i_inc >= 0 ? i <= i_end : i >= i_end; i += i_inc) {
// Jeden Polygonpunkt speichern - hier werden die x, y und z Koordinaten nur aus einem String ausgelesen und zugeordnet
tempPolygonpunkt = tempPolygon[(i - 1)].split(',');
tempPolygonpunkt_x = parseFloat((tempPolygonpunkt.slice(0, 1)));
tempPolygonpunkt_y = parseFloat((tempPolygonpunkt.slice(1, 2)));
tempPolygonpunkt_z = parseFloat((tempPolygonpunkt.slice(2, 3)));
// Und den nächsten Polygonpunkt - das gleiche noch einmal für den zu vergleichenden Polygonpunkt
j = (i + 1) % tempPolygon.length + 1;
tempPolygonpunkt_1 = tempPolygon[(j - 1)].split(',');
tempPolygonpunkt_1_x = parseFloat((tempPolygonpunkt_1.slice(0, 1)));
tempPolygonpunkt_1_y = parseFloat((tempPolygonpunkt_1.slice(1, 2)));
tempPolygonpunkt_1_z = parseFloat((tempPolygonpunkt_1.slice(2, 3)));
// Punkt-in-Polygon-Test nach Jordan - hier endet meine Idee, denn der Rest funktioniert nur in der Ebene - im Prinzip müsste man es um die y-Koordinate erweitern, wenn das überhaupt geht
// Hier nur in x-z Ebene!!!
if (((tempPolygonpunkt_z < Sonne_z && tempPolygonpunkt_1_z >= Sonne_z) || (tempPolygonpunkt_1_z < Sonne_z && tempPolygonpunkt_z >= Sonne_z))) {
if ((Sonne_z - tempPolygonpunkt_z) * (tempPolygonpunkt_1_x - tempPolygonpunkt_x) < (Sonne_x - tempPolygonpunkt_x) * (tempPolygonpunkt_1_z - tempPolygonpunkt_z)) {
PunktInnerhalb = !PunktInnerhalb;
}
}
}
return PunktInnerhalb;
}

Kommentiert 21 Apr von Guitardoc

Jetzt geht's aber los ;-)

oben hatte ich es ja schon beschrieben. Ich würde zunächst mal die Ebene bestimmen, die die Sonne von ihrer aktuellen Position überstreicht. Das ist zwear keine Ebene (es sei denn das Fenster befindet sich am Äquator), aber man kann das erstmal annehmen.

Den Richtungsvektor zur Sonne hast Du ja schon. Ich nenne ihn mal $\vec s$. Dann braucht man den Vektor, der in Richtung der Rotationsachse der Erde zeigt (nördliche Orientierung) - ich nenne in mal $\vec{N}$. Daraus berechne ich den Normalenvektor $\vec{N}^*$, der senkrecht auf $\vec s$ steht und nach Norden zeigt.$$\vec{N}^* = \vec N - \frac{\left< \vec{s},\,\vec{N}\right>}{|\vec{s}|^2} \vec{s}$$(siehe auch Orthogonalisierung) Das ist dann auch gleichzeitig der Normalenvektor der Ebene in der sich näherungsweise die Sonne bewegt (der Sonnenebene). Dies ist der Großkreis, den ich in meiner Antwort erwähnt habe.

Diese muss nun mit jedem Ebenenstück, welches sich aus zwei auf einander folgenden Stützpunkte zusammen setzt, zum Schnitt gebracht werden. Da Du wahrscheinlich keine Winkel hast, die 180° überschreiten, reicht es sicher aus, den Vorzeichenwechsel des Skalarprodukts mit $\vec{N}^*$ zu erkennen.

Mal angenommen zwei Punkte des Polygons sind $p_k$ und $p_{k+1}$, dann rechne $$v_k = \operatorname{sgn}\left(\left< p_k,\, \vec{N}^*\right>\right), \quad v_{k+1} = \operatorname{sgn}\left(\left< p_{k+1},\, \vec{N}^*\right>\right)$$Die $\operatorname{sgn}$-Funktion stellt nur das Vorzeichen fest und liefert $1$, $0$ oder $-1$. Ein Vorzeichenwechsel liegt vor, wenn $$v_{k} \cdot v_{k+1} = -1 \quad \lor \quad v_{k}= 0 \land v_{k+1} \ne 0$$dann liegen nämlich die beiden Punkte ober- und unterhalb der Sonnenebene. Nur in diesem Fall ist das Polygonstück relevant!

Die beiden Vektoren spannen dann eine weitere Ebene auf$$E_{k}: \quad \vec{n}_k \vec x = 0 \quad \vec{n}_k = p_k \times p_{k+1}$$Dann gilt es den Winkel zu berechnen, den $\vec{s}$ von dem Schnittgerade der beiden Ebenen entfernt ist. Die Schnittgerade ist$$g_k:\quad \vec{x} = \vec{d}_k\lambda \quad \vec{d}_k = \vec{N}^* \times \vec{n}_{k}$$und der Sinus diese Winkels ist$$\sin\left(\varphi_{k}\right) = \left<\frac{\vec{d}_k \times \vec{s}}{|\vec{d}_k| |\vec{s}|}, \, \frac{\vec{N}^*}{|\vec{N}^*|}\right>$$Du müsstes dann eine gerade Anzahl solcher Winkel erhalten. Ist der Polygonzug konvex, so müssten es genau zwei sein.

Ich glaube, es reicht jetzt aus, die Anzahl der negativen Winkel zu zählen. Ist diese ungerade, scheint die Sonne durch das Polygon auf die Ecke des Fensters. Ist einer der Winkel $=0$ so liegt die Richtung zur Sonne auf der Kante.

Ich hoffe, ich habe mich jetzt nicht mit den Richtungen vertan. Ich habe mir das jetzt alles ohne Zeichnung im Kopf vorgestellt. Versuch es einfach mal nachzuvollziehen

Gruß Werner

Kommentiert 22 Apr von Werner-Salomon

Hallo Werner,

Ganz vielen Dank für deine ausführliche Anleitung! Ich bin (völlig entgegen meiner Erwartung ☺) beim Durcharbeiten nur über einen Punkt gestolpert:

Der Vektor in Richtung der Rotationsachse der Erde müsste $$\vec{N} = \begin{pmatrix} x\\y\\z \end{pmatrix} = \begin{pmatrix} 0\\0\\1 \end{pmatrix}$$ sein, richtig?

Ich habe die Vorgehensweise noch nicht hundertprozentig verstanden, werde es mir heute aber noch mal detailliert zu Gemüte führen. Ich habe den Ablauf schon mal fix programmiert - wenn du magst, kannst du ja mal drüber schauen, dass ich keine Fehler eingebaut habe (bei Skalarpodukten, Kreuzprodukten und Vektorbeträgen kann man schon mal durcheinander kommen...):

Ich hoffe, dass man einigermaßen erkennen kann, was ich programmiert hab... Ansonsten hier noch mal als Javascript:

async function SonneinBereich(tempPolygon, Sonne_x, Sonne_y, Sonne_z) {
SonneInnerhalb = false;
ListePhi = [];
var i_end = tempPolygon.length;
var i_inc = 1;
if (1 > i_end) {
i_inc = -i_inc;
}
for (i = 1; i_inc >= 0 ? i <= i_end : i >= i_end; i += i_inc) {
// Jeden Polygonpunkt speichern
tempPolygonpunkt = tempPolygon[(i - 1)].split(',');
pk_x = parseFloat((tempPolygonpunkt.slice(0, 1)));
pk_y = parseFloat((tempPolygonpunkt.slice(1, 2)));
pk_z = parseFloat((tempPolygonpunkt.slice(2, 3)));
// Und den nächsten Polygonpunkt
j = (i + 1) % tempPolygon.length + 1;
tempPolygonpunkt_1 = tempPolygon[(j - 1)].split(',');
pk_1_x = parseFloat((tempPolygonpunkt_1.slice(0, 1)));
pk_1_y = parseFloat((tempPolygonpunkt_1.slice(1, 2)));
pk_1_z = parseFloat((tempPolygonpunkt_1.slice(2, 3)));
// Vektor N in Richtung der Rotationsachse der Erde (nördliche Ausrichtung)
N_x = 0;
N_y = 0;
N_z = 1;
// Skalarprodukt s,N
s_N = Sonne_x * N_x + Sonne_y * N_y + Sonne_z * N_z;
// Skalarprodukt s,s
s_s = Math.pow(Sonne_x, 2) + Math.pow(Sonne_y, 2) + Math.pow(Sonne_z, 2);
// Vektor Nstern
Nstern_x = N_x - (s_N / s_s) * Sonne_x;
Nstern_y = N_y - (s_N / s_s) * Sonne_y;
Nstern_z = N_z - (s_N / s_s) * Sonne_z;
// vk = Vorzeichen von pk, Nstern -> (-1, 0, +1)
vk = pk_x * Nstern_x + pk_y * Nstern_y + pk_z * Nstern_z;
if (vk < 0) {
vk = -1;
} else if (vk > 0) {
vk = 1;
}
// vk+1 = Vorzeichen von pk+1, Nstern -> (-1, 0, +1)
vk_1 = pk_1_x * Nstern_x + pk_1_y * Nstern_y + pk_1_z * Nstern_z;
if (vk_1 < 0) {
vk_1 = -1;
} else if (vk_1 > 0) {
vk_1 = 1;
}
// Test, ob die Punkte ober- und unterhalb der Sonnenebene liegen
if ((vk * vk_1 == -1 && (vk == 0 || vk_1 != 0))) {
// Kreuzprodukt nk = pk x pk+1
nk_x = pk_y * pk_1_z - pk_z * pk_1_y;
nk_y = pk_z * pk_1_x - pk_x * pk_1_z;
nk_z = pk_x * pk_1_y - pk_y * pk_1_x;
// Kreuzprodukt dk = N* x nk
dk_x = Nstern_y * nk_z - Nstern_z * nk_y;
dk_y = Nstern_z * nk_x - Nstern_x * nk_z;
dk_z = Nstern_x * nk_y - Nstern_y * nk_x;
// Kreuzprodukt dkxs = dk x s
dkxs_x = dk_y * Sonne_z - dk_z * Sonne_y;
dkxs_y = dk_z * Sonne_x - dk_x * Sonne_z;
dkxs_z = dk_x * Sonne_y - dk_y * Sonne_x;
// Skalarprodukte |dk|*|s|, |N*|
_7Cdk_7C__7Cs_7C = Math.sqrt(Math.pow(dk_x, 2) + Math.pow(dk_y, 2) + Math.pow(dk_z, 2)) * Math.sqrt(Math.pow(Sonne_x, 2) + Math.pow(Sonne_y, 2) + Math.pow(Sonne_z, 2));
_7CNstern_7C = Math.sqrt(Math.pow(Nstern_x, 2) + Math.pow(Nstern_y, 2) + Math.pow(Nstern_z, 2));
// sin(Phi) und Phi berechnen
sinPhi = (dkxs_x / _7Cdk_7C__7Cs_7C) * (Nstern_x / _7CNstern_7C) + (dkxs_y / _7Cdk_7C__7Cs_7C) * (Nstern_y / _7CNstern_7C) + (dkxs_z / _7Cdk_7C__7Cs_7C) * (Nstern_z / _7CNstern_7C);
Phi = Math.acos(sinPhi) / Math.PI * 180;
console.debug((['Phi = ',Phi,'°'].join('')));
if (Phi < 0) {
ListePhi.push(Phi);
}
}
}
console.debug(('Liste der negativen Winkel: ' + String(ListePhi)));
SonneInnerhalb = ListePhi.length % 2 === 1;
return SonneInnerhalb;
}

Kommentiert 23 Apr von Guitardoc

Hallo Werner,

Ich komme einfach nicht weiter. Ich denke, es liegt an der Bestimmung des Vektors, der in Richtung der Rotationsachse der Erde zeigt (nördliche Orientierung). Ich vermute, dass meine Umrechnung nicht stimmt, bekomme es aber scheinbar nicht hin... Vielleicht liegt es auch an meiner Annahme des Koordinatensystems?

Mein Richtungsvektor zur Sonne beruht auf einem Horizontsystem mit +x in Richtung Süden (-x in Richtung Norden), +y in Richtung Osten und +z in Richtung Zenit. Das Azimut A verläuft im Uhrzeigersinn mit 0° im Norden und die Elevation h in Richtung Zenit mit 0° am Horizont.

Daraus ergeben sich der Richtungsvektor zur Sonne sowie die Richtungsvektoren zu den Begrenzungspunkten der Fläche:

$$\vec{s} = \begin{pmatrix} x\\y\\z \end{pmatrix} = \begin{pmatrix} -cos(h)cos(A)\\cos(h)sin(A)\\sin(h) \end{pmatrix}$$

Wie errechnet sich denn auf dieser Grundlage der Vektor, der in Richtung der Rotationsachse der Erde zeigt (nördliche Orientierung)? Die Frage ist vermutlich banal, aber ich kann es einfach nicht (mehr)... Könntest du mir bitte noch einmal helfen? Vielen lieben Dank schon mal im Voraus - und ich übe mich selbstverständlich in Geduld, wenn ich weiß, dass ich geholfen bekomme... :-)

Kommentiert 27 Apr von Guitardoc

Ich habe den Ablauf schon mal fix programmiert - wenn du magst, kannst du ja mal drüber schauen, dass ich keine Fehler eingebaut habe (bei Skalarpodukten, Kreuzprodukten und Vektorbeträgen kann man schon mal durcheinander kommen...):

Ja - das 'durcheinander kommen' kann man deutlich reduzieren, wenn man Skalar- und Kreuzprodukt, nebst Vektorbetrag in Unterprogramme verpackt. Oder besser - wenn möglich(!) - objektorientiert in Klassen unterbringt.

Falls das nicht möglich ist, schnitze Dir die Funktionen dot für das Skalarprodukt, cross für das Vektorprodukt, abs für den Betrag und unit für den Einheitsvektor.

Den Code in der obigen Form hier einzustellen, ist nicht sinnvoll. Besser wäre es, wenn Du die Inputparameter angibst (hast Du schon) und dann Zwischenergebnisse postest. Ich kann das dann selber nachrechnen - mit C++ oder Excel, aber nicht mit Javascript.

(ich rechne jetzt mal selber)

Kommentiert 27 Apr von Werner-Salomon

Hallo Werner,

Vielen lieben Dank, dass du dir das noch einmal angesehen hast!

Was ich mir bei

Der Vektor in Richtung der Rotationsachse der Erde müsste $$\vec{N} = \begin{pmatrix} x\\y\\z \end{pmatrix} = \begin{pmatrix} 0\\0\\1 \end{pmatrix}$$sein, richtig?

gedacht habe, weiß ich auch nicht... Können wir das bitte ganz schnell vergessen? :-)

Ja, ich hätte das lieber auch gern objektorientiert programmiert, da sieht man auf den ersten Blick was gehauen und gestochen ist. Die Berechnung muss aber in mein Hausautomationssystem (ioBroker) rein, und das versteht nur Javascript (was ich nicht beherrsche). Es gibt da aber so ein Hilfsmittel, das nennt sich Blockly, da klickert man die Logik mit Bilderchen zusammen (das hatte ich oben mal gepostet) und dann wird das intern in Javascript übersetzt. Unterprogramme lassen sich damit zwar realisieren, aber maximal nur mit einem einzigen Rückgabewert. Und da es in Blockly keine Objekte gibt, muss ich mit diesem Spaghetticode irgendwie klarkommen... Das nächste Problem ist, dass man ohne Aufwand die Daten nicht immer gleichzeitig sieht (man muss dafür erst aufwändige Ausgaben schreiben), daher die Abweichung im Sonnenvektor (da hab ich die Daten erst mit ein paar Minuten Abstand gesehen).

Da ist es wirklich die beste Variante, sich über Zwischenwerte auszutauschen - und das hab ich gerade mal gemacht und die ganzen Ausgaben in den Code reingeschrieben. Damit stimmen sie zeitlich überein. In meinem Test ergeben sich folgende Werte:

Liste der Grenzpunkte (Azimut, Elevation):

180, 0

180, 60

270, 45

300, 30

300, 0

Grenzpunkte umgerechnet in x,y,z Koordinaten:

1, 0, 0

0.5, 0, 0.8660

0, -0.7071, 0.7071

-0.43300, -0.75, 0,5

-0.5, -0.8660, 0

Sonne (Azimut, Elevation):

274, 15.2

Sonne umgerechnet in x,y,z Koordinaten:

-0.0673, -0.9627, 0.2622

Ich bekomme dann zwei Winkel Phi:

Phi 1 = 133.5668293844902°

Phi 2 = 61.46435809214322°

Es müsste jetzt aber einer der beiden Winkel negativ sein, denn die Sonne befindet sich innerhalb der Grenzpunkte? Ich werde mal noch ein paar mehr Ausgaben der Zwischenwerte einprogrammieren, vielleicht erkennt man dann mehr.

PS. Der Geoknecht ist übrigens eine geniale Erfindung. Den hatte ich noch gar nicht entdeckt.

Kommentiert 28 Apr von Guitardoc

So, habe die Ausgaben hinzugefügt. Ich poste die gleich mal im Rohformat (Zeitstempel bitte ignorieren). Es sind die gleichen Ausgangsdaten wie vorhin (gleiche Grenzpunkte, gleicher Sonnenstand):

Verschattung Test Westseite wird geprüft
21:54:51.609

Grenzpunktwinkelliste: 180,0,180,60,270,45,300,30,300,0
21:54:51.610

Grenzpunkte x,y,z: 1,1.2246467991473532e-16,0,0.5000000000000001,6.123233995736767e-17,0.8660254037844386,1.29893408435324e-16,-0.7071067811865476,0.7071067811865475,-0.43301270189221946,-0.75,0.49999999999999994,-0.5000000000000001,-0.8660254037844386,0
21:54:51.610

Sonne x,y,z: -0.0673,0.9627,0.2622
21:54:51.610

Grenzpunkt i: 1, Grenzpunkt j: 3
21:54:51.610

Skalarprodukt s,N: 0.24666583099588252
21:54:51.610

Skalarprodukt s,s: 1.00006942
21:54:51.610

Vektor N* (x,y,z): -0.6089361960364841, -0.23744871181016225, 0.7155242901082678
21:54:51.610

Skalarprodukt vk = pk, N*: -0.6089361960364841
21:54:51.610

Vorzeichen von vk (0, -1, +1): -1
21:54:51.610

Skalarprodukt vk+1 = pk+1, N*: 0.6738536719442225
21:54:51.610

Vorzeichen von vk+1 (0, -1, +1): 1
21:54:51.610

Vorzeichen von vk * vk+1: -1
21:54:51.610

Neuer Winkel Phi 1 wurde gefunden.
21:54:51.610

Kreuzprodukt nk = pk x pk+1: 8.659560562354932e-17, -0.7071067811865475, -0.7071067811865476
21:54:51.610

Kreuzprodukt dk = N* x nk: 0.6738536719442226, -0.43058291352733874, 0.43058291352733874
21:54:51.610

Kreuzprodukt dkxs = dk x s: -0.5274210107796372, -0.20566266286416507, 0.6197406999003132
21:54:51.610

Betrag von dk * Betrag von s: 0.9082617053189251
21:54:51.610

Betrag von N*: 0.9691027764476942
21:54:51.610

sin(Phi 1) = 0.9241542626577512
21:54:51.610

Phi 1 = 22.45883139310086°
21:54:51.610

Grenzpunkt i: 2, Grenzpunkt j: 4
21:54:51.611

Skalarprodukt s,N: 0.24666583099588252
21:54:51.611

Skalarprodukt s,s: 1.00006942
21:54:51.611

Vektor N* (x,y,z): -0.6089361960364841, -0.23744871181016225, 0.7155242901082678
21:54:51.611

Vorzeichen von vk (0, -1, +1): 1
21:54:51.611

Skalarprodukt vk+1 = pk+1, N*: 0.7995257864374837
21:54:51.611

Vorzeichen von vk+1 (0, -1, +1): 1
21:54:51.611

Vorzeichen von vk * vk+1: 1
21:54:51.611

Es wurde kein neuer Winkel Phi gefunden.
21:54:51.611

Skalarprodukt s,N: 0.24666583099588252
21:54:51.611

Skalarprodukt s,s: 1.00006942
21:54:51.611

Vektor N* (x,y,z): -0.6089361960364841, -0.23744871181016225, 0.7155242901082678
21:54:51.611

Vorzeichen von vk (0, -1, +1): 1
21:54:51.611

Vorzeichen von vk+1 (0, -1, +1): 1
21:54:51.611

Es wurde kein neuer Winkel Phi gefunden.
21:54:51.611

Skalarprodukt s,N: 0.24666583099588252
21:54:51.611

Skalarprodukt s,s: 1.00006942
21:54:51.611

Vektor N* (x,y,z): -0.6089361960364841, -0.23744871181016225, 0.7155242901082678
21:54:51.611

Skalarprodukt vk = pk, N*: 0.3151941142403443
21:54:51.611

Vorzeichen von vk (0, -1, +1): 1
21:54:51.611

Skalarprodukt vk+1 = pk+1, N*: 0.7995257864374837
21:54:51.611

Vorzeichen von vk+1 (0, -1, +1): 1
21:54:51.611

Vorzeichen von vk * vk+1: 1
21:54:51.611

Es wurde kein neuer Winkel Phi gefunden.
21:54:51.611

Grenzpunkt i: 3, Grenzpunkt j: 5
21:54:51.611

Skalarprodukt s,N: 0.24666583099588252
21:54:51.611

Skalarprodukt s,s: 1.00006942
21:54:51.611

Vektor N* (x,y,z): -0.6089361960364841, -0.23744871181016225, 0.7155242901082678
21:54:51.611

Skalarprodukt vk = pk, N*: 0.6738536719442225
21:54:51.611

Vorzeichen von vk (0, -1, +1): 1
21:54:51.611

Skalarprodukt vk+1 = pk+1, N*: 0.5101047145417327
21:54:51.611

Vorzeichen von vk+1 (0, -1, +1): 1
21:54:51.611

Vorzeichen von vk * vk+1: 1
21:54:51.611

Es wurde kein neuer Winkel Phi gefunden.
21:54:51.611

Grenzpunkt i: 4, Grenzpunkt j: 1
21:54:51.611

Skalarprodukt s,N: 0.24666583099588252
21:54:51.611

Skalarprodukt s,s: 1.00006942
21:54:51.611

Vektor N* (x,y,z): -0.6089361960364841, -0.23744871181016225, 0.7155242901082678
21:54:51.611

Skalarprodukt vk = pk, N*: 0.7995257864374837
21:54:51.611

Vorzeichen von vk (0, -1, +1): 1
21:54:51.611

Skalarprodukt vk+1 = pk+1, N*: -0.6089361960364841
21:54:51.611

Vorzeichen von vk+1 (0, -1, +1): -1
21:54:51.611

Vorzeichen von vk * vk+1: -1
21:54:51.611

Neuer Winkel Phi 2 wurde gefunden.
21:54:51.611

Kreuzprodukt nk = pk x pk+1: -6.123233995736765e-17, 0.49999999999999994, 0.75
21:54:51.611

Kreuzprodukt dk = N* x nk: -0.5358486789117556, 0.456702147027363, -0.304468098018242
21:54:51.611

Kreuzprodukt dkxs = dk x s: 0.4128587409127361, 0.16099022660729, -0.48512546869340556
21:54:51.611

Betrag von dk * Betrag von s: 0.7671064645971244
21:54:51.611

Betrag von N*: 0.9691027764476942
21:54:51.611

sin(Phi 2) = -0.8565325841395597
21:54:51.611

Phi 2 = 148.92946274571392°
21:54:51.611

Grenzpunkt i: 5, Grenzpunkt j: 2
21:54:51.612

Skalarprodukt s,N: 0.24666583099588252
21:54:51.612

Skalarprodukt s,s: 1.00006942
21:54:51.612

Vektor N* (x,y,z): -0.6089361960364841, -0.23744871181016225, 0.7155242901082678
21:54:51.612

Skalarprodukt vk = pk, N*: 0.5101047145417327
21:54:51.612

Vorzeichen von vk (0, -1, +1): 1
21:54:51.612

Vorzeichen von vk (0, -1, +1): 1
21:54:51.612

Skalarprodukt vk+1 = pk+1, N*: 0.3151941142403443
21:54:51.612

Vorzeichen von vk+1 (0, -1, +1): 1
21:54:51.612

Vorzeichen von vk * vk+1: 1
21:54:51.612

Es wurde kein neuer Winkel Phi gefunden.
21:54:51.612

Liste der negativen Winkel Phi:
21:54:51.612

Test Westseite: Sonne befindet sich nicht im Verschattungsbereich.

Ich hoffe, man kann erkennen, wie der Ablauf ist. Der Punkt i läuft die Liste von Anfang bis Ende durch und der Punkt j ermittelt sich in Anlehnung an das Verfahren von Jordan zu

Länge von tempGrenzpunkte ist dabei die gesamte Anzahl der Punkte.

Kommentiert 28 Apr von Guitardoc

Die Formel oben ist falsch - auch mit $\vec n_k$. Es reicht doch aus, wenn man berechnet auf welcher 'Seite' der durch $p_k$ und $p_{k+1}$ aufgespannten Ebene sich der Sonnenvektor befindet.

Dazu reicht es aus, sich das Vorzeichen des Skalarprodukts $$\operatorname{sgn}\left< \vec s,\, \vec{v}_{k}\right>$$anzusehen (es gibt viele Wege nach Rom!). Eine Berechnung von $\vec{d}_{k}$ ist dann nicht mehr notwendig.

Achtung: das gilt zumindest, wenn der Bereich keine überstumpfen Winkel besitzt.

Sind alle Vorzeichen identisch, dann liegt der Sonnenvektor im Bereich. Das ist hier der Fall. Das Skalarprodukt ist in beiden Fällen positiv. Da die 'Drehrichtung' des Bereichs rechtsdrehend ist, zeigen alle $\vec n_k$ nach innen und die Vorzeichen sind positiv.

Wobei - für eine wirklich schlüssige Lösung komme ich immer wieder auf meine ursprüngliche Idee zurück (s. Antwort). Ich denk' da noch mal drüber nach. Hier noch das Bild zu Deinen Daten:

... musst Du Dir in Geoknecht3D ansehen, sonst sieht man es nicht richtig.

Kommentiert 28 Apr von Werner-Salomon

Ich hab bei mir im Algorithmus noch einen Fehler entdeckt - das habe ich aber erst gesehen, als ich mir die Zeichnung im Geoknecht3D richtig zu Gemüte geführt hab. Ich hatte die Berechnung des Punktes pk+1 aus dem Verfahren von Jordan übernommen, aber das funktioniert nur in der Fläche. Im Raum muss man tatsächlich den nachfolgenden Punkt nehmen. Eigentlich logisch, wenn man es einmal räumlich gesehen hat. Also das konnte nicht funktionieren. Aber jetzt funktioniert es - es gibt genau zwei Schnittpunkte der Ebenen, einmal zwischen P0 und P1 und einmal zwischen P4 und P0. Allerdings ergeben sich immer noch nur positive Winkel. Irgendwo muss sich da noch mal der Fehlerteufel eingeschlichen haben. Soll ich die neuen Zwischenwerte noch mal schicken?

Kommentiert 30 Apr von Guitardoc

ich habe jetzt selber das ganze in C++ programmiert (das war der leichte Teil). Das Problem ist die Orientierung des Vektors $d_{k}$.

Mir ist das ja ein wenig peinlich. Eigentlich müsste ich sowas ohne Probleme hin bekommen. Aber bei den vielen Vektoren kann man schon mal den Überblick verlieren. Und wie ich oben bereits vermutet habe, sind die Vorzeichen das Problem.

Tut mit auch leid, dass mein Antworten immer so zeitverzögert kommen. Irgenwie komme ich z.Zt. zu gar nichts. Und da ich Dir in jedem Fall auch antworten wollte, habe ich dann wohl zu wenig drüber nachgedacht.

ich melde mich in jedem Fall noch ... bis wann benötigst Du denn eine Lösung?

Soll ich die neuen Zwischenwerte noch mal schicken?

Nicht nötig, die habe ich jetzt selber (s.o.)

Kommentiert 1 Mai von Werner-Salomon

Hallo Werner,

Je mehr ich über das Problem und deinen ursprünglichen Lösungsansatz nachdenke, umso besser gefällt er mir. Eigentlich müsste er die richtigen Ergebnisse bringen. Ich bin noch(wieder) ein bisschen langsam, aber ich denke, ich finde raus, warum die Ergebnisse nicht stimmen. Theoretisch müssten sie stimmen. aber praktisch tun sie es nicht - und das bei einem wirklich noch recht einfachen Beispiel. Wahrscheinlich würden uns unsere Studenten auslachen wenn sie wüssten, wer wir hier sind... Das fuchst mich aber total!!!

PS. Und wieder eine Überschneidung! Ich hab keine Deadline für die Lösung. Wenn es warm wird, dann freut sich nur meine Familie, wenn es nicht mehr so warm im Haus ist... aber im Moment ist das noch kein Thema und es lässt sich sicher auch im Juli aushalten. Für eine gute Lösung ist ein bisschen schwitzen ein guter Preis... :-)

Kommentiert 1 Mai von Guitardoc

Hallo Werner,

Ich habe mich jetzt auch noch mal mit dem Thema beschäftigt. Wäre es vielleicht ein Ansatz, direkt in Richtung des Sonnenvektors zu schauen und jetzt die Punkte, die sich aus den Spitzen der einzelnen Grenzpunktvektoren ergeben, auf eine Ebene senkrecht zum Sonnenvektor zu projizieren? Dann könnte man mit dem Verfahren nach Jordan doch berechnen, ob sich die Sonne in diesem jetzt ebenen Polygon befindet? Also quasi so:

PS. Kann den Geoknecht irgendwie nicht dazu bringen, meine Rotation im Link zu speichern - er dreht immer wieder auf deine ursprüngliche Rotation... Ich hoffe aber, dass du erkennen kannst, was ich meine? Im Moment befindet sich die Sonne innerhalb des ebenen Polygons. Wenn man die Sonne weiter wandern lässt, dann wandert sie irgendwann aus dem ebenen Polygon heraus. Die Form des Polygons ändert sich jetzt natürlich ständig - das muss man dann in die Ebene umrechnen - wenn das irgendwie geht?

Kommentiert 13 Mai von Guitardoc

ich denke, ich hab's jetzt. Ich habe ziemlich lange über die Richtung der Vektoren $\vec d_k$ nachgedacht, bis mir 'ne ziemlich simple Lösung eingefallen ist.

Es gilt alles wie oben schon beschrieben bis zur Berechnung der Vorzeichen $v_k$ mit denen die relevanten Paare $\vec p_k,\space \vec p_{k+1}$ gefunden werden. Den Vektor $\vec n_k$ brauchen wir nicht mehr.

Die Schnittvektoren $\vec d_k$ berechnet man dann nach$$d_k' = \left|\left<\vec p_{k+1},\,\vec N^* \right>\right| \vec p_k + \left|\left<\vec p_{k},\,\vec N^*\right>\right|\vec p_{k+1}\\ \vec d_k = \frac{\vec d_k'}{|\vec d_k'|}$$(Bem.: die Normierung ist nicht zwingend nötig) dies ist schlicht der Schnittpunkt der Strecke $\vec p_k \vec p_{k+1}$ mit der aktuellen Sonnenebene.

Der Winkel $\varphi_k$ auf den Vektor $\vec s$ bezogen muss nun in einem Bereich von -180° bis +180° berechnet werden.$$\varphi_k = \operatorname{arctan2}\left(\left<(\vec s \times \vec d_k),\,\vec N^*\right>,\, \left<\vec s,\, \vec d_k\right>\right)$$wobei ich hier davon ausgehe, dass die Parameter arctan2(y,x) lauten. Siehe auch hier.

Von allen Winkeln $\varphi_k$ wird nun derjenige ausgewählt, dessen Absolutwert am kleinsten ist. Dann kann man davon ausgehen, dass zwischen dieser Grenze und $\vec s$ keine weitere Grenze liegt.

Nun muss man nur noch feststellen, ob $\vec s$ 'links' von dieser Grenze liegt. Du erinnerst Dich: Voraussetzung ist ein positiver Drehsinn der $\vec p_k$ um den 'offenen' Bereich herum.

$\vec s$ liegt genau dann 'links', wenn das zugehörige Vorzeichen $v_{k+1}$ identisch mit dem Vorzeichen des Winkels $\varphi_k$ ist. Ist $v_{k+1}=0$, dann vergleiche mit $-v_{k}$. Ist $\varphi_k=0$, dann liegt $\vec s$ auf der Grenze des Polygons.

Ich habe das ganze programmtechnisch überprüft. Das sollte passen ;-)

Achso: und falls kein Paar $\vec p_k,\, \vec p_{k+1}$ gefunden wird, liegt $\vec s$ außerhalb, wenn das Polygon weniger als eine Halbkugel einnimmt.

Kommentiert 5 Jun von Werner-Salomon

Hallo Werner,

Zunächst erst einmal ganz vielen lieben Dank, dass du dir das noch mal angeschaut hast!!!! Ich weiß das wirklich zu schätzen!!!!!

Ich muss gestehen, dass ich erst mal nachschauen musste, was atan2 überhaupt macht. Entweder sind meine grauen Zellen schon so weit degeneriert, dass ich mich an die Mathevorlesung dazu nicht erinnern kann, oder ich bin schon so alt, dass es die Funktion damals noch nicht gab... :-)

Ich hab die Vorgehensweise gerade mit meiner kastrierten Programmiersprache umständlich in das System eingetippt und hoffe, dass ich keine Fehler eingebaut habe. Bisher bin ich mir da noch nicht so ganz sicher... Wenn ich darf, dann würde ich dir einmal meine Ausgangs- und Zwischenwerte schicken?

Eine Frage habe ich aber noch - wie kann es passieren, dass kein Paar $\vec p_k,\, \vec p_{k+1}$ gefunden wird? Wenn man eine Liste von Vektoren hat (und die Liste hat mehrere davon, sonst wird es ja kein Polygon), dann gibt es doch immer ein Paar $\vec p_k,\, \vec p_{k+1}$?

Kommentiert 6 Jun von Guitardoc

Kann den Geoknecht irgendwie nicht dazu bringen, meine Rotation im Link zu speichern

Ja - unten rechts im Bild steht der Befehl der die Kamera positioniert mit der aktuellen Position. Wenn Du eine gute Position gefunden hast, kopiere den einfach mit in das Script oben links hinein.

Ich muss gestehen, dass ich erst mal nachschauen musste, was atan2 überhaupt macht.

Ja - das kennt nicht jeder, daher habe ich auch gleich den Link zum Wiki-Artikel hinzu gefügt.

Ich hab die Vorgehensweise gerade mit meiner kastrierten Programmiersprache umständlich in das System eingetippt

Ok - wäre aber programmtechnisch gar nicht nötig gewesen. Du kannst die Winkel $\varphi_k$ auch nach ihren Cosinus sortieren. Der größte Cosinus-Wert gewinnt das ist der mit dem kleinsten absoluten Winkel und es gilt $$\cos(\varphi_k) = \frac{\left<\vec s,\, \vec d_k\right>}{|\vec s|} \quad |\vec d_k| = 1$$das setzt aber voraus, dass $\vec d_k$ normiert worden ist, sonst kann man die unterschiedlichen Winkel $\varphi_k$ nicht vergleichen. $|\vec s|$ ist ja immer geich, also kann man sich hier die Division auch sparen.

Und im Intervall $(-\pi,\,\pi]$ ist das Vorzeichen des Winkels gleich dem Vorzeichen seines Sinus:$$\operatorname{sgn}(\varphi_k) = \operatorname{sgn}(\sin(\varphi_k)) = \operatorname{sgn}\left(\left< (\vec s\times \vec d_k),\,\vec N^*\right>\right)$$.. hätte ich vielleicht gleich schreiben sollen.

Wenn ich darf, dann würde ich dir einmal meine Ausgangs- und Zwischenwerte schicken?

Äh - schlechte Nachrichten für Dich. Ich bin ab morgen im Laufe des Tages schon wieder für eine Woche so gut wie offline. Ich kann Dir nur empfehlen Dir die Zwischenergebnisse in Geoknecht3D zu veranschaulichen. Ich bin erst ab dem 17./18. Juni wieder erreichbar.

Eine Frage habe ich aber noch - wie kann es passieren, dass kein Paar $\vec p_k,\, \vec p_{k+1}$ gefunden wird?

Stelle Dir ein kleines Loch in Bodennähe vor, welches Du mit einem Polygon beschreibst. Dann kann doch es gut sein, dass die Sonnenebene darüber liegt und es zu keinem relevanten Schnittpunkt mit einem der Polygonsegmente kommt.

Viel Erfolg mit Deinem Programm & Gruß Werner

Kommentiert 6 Jun von Werner-Salomon

Hallo Werner,

Irgendwo habe ich einen Fehler, aber ich finde ihn nicht. Die Sonne befindet sich innerhalb des Polygons, das ergibt der Algorithmus aber nicht... Solltest du doch zwischenzeitlich mal reinschauen, vielleicht siehst du ja auf den ersten Blick wo hier was nicht stimmt:

+x in Richtung Süd, +y in Richtung Ost, +z in Richtung Zenit

Grenzpunktwinkelliste (Azimut, Elevation): 180,0; 180,60; 270,45; 300,30; 300,0

Grenzpunkte (x,y,z):

1, 1.2246467991473532e-16, 0;

0.5000000000000001, 6.123233995736767e-17, 0.8660254037844386;

1.29893408435324e-16, -0.7071067811865476, 0.7071067811865475;

-0.43301270189221946, -0.75, 0.49999999999999994;

-0.5000000000000001, -0.8660254037844386, 0

Vektor Sonne (x,y,z): -0.0673, 0.9627, 0.2622

Vektor N (x,y,z): -0.6255356541281263, 0, 0.7801955815143388

Start der Berechnung der Parameter für alle Grenzpunkte

-------------------- Grenzpunkt i: 1, Grenzpunkt j: 2 --------------------

Skalarprodukt <s,N>: 0.24666583099588252

Skalarprodukt <s,s>: 1.00006942

Vektor Nstern (x,y,z): -0.6089361960364841, -0.23744871181016225, 0.7155242901082678

Skalarprodukt <pk+1,Nstern>: 0.3151941142403443

Skalarprodukt <pk,Nstern>: -0.6089361960364841

Betrag |pk+1,Nstern|: 0.3151941142403443

Betrag |pk,Nstern|: 0.6089361960364841

Vektor d'k (x,y,z): 0.6196622122585864, 7.588673447950457e-17, 0.5273542150514562

Betrag |d'k|: 0.8136852741900537

Vektor dk (x,y,z): 0.7615502355937328, 9.326300584097776e-17, 0.6481058853822655

Kreuzprodukt s*dk (x,y,z): 0.623931535857507, 0.24329599785890318, -0.7331444118060866

Skalarprodukt <s*dk,Nstern>: -0.9622874521168296

Skalarprodukt <s,dk>: 0.11868103229177188

Winkel Phi: -82.96909044086752

Betrag |Phi|: 82.96909044086752

Vorzeichen vk von <pk,Nstern> (0, -1, +1): -1

Vorzeichen vk+1 von <pk+1,Nstern> (0, -1, +1): 1

-------------------- Grenzpunkt i: 2, Grenzpunkt j: 3 --------------------

Skalarprodukt <s,N>: 0.24666583099588252

Skalarprodukt <s,s>: 1.00006942

Vektor Nstern (x,y,z): -0.6089361960364841, -0.23744871181016225, 0.7155242901082678

Skalarprodukt <pk+1,Nstern>: 0.6738536719442225

Skalarprodukt <pk,Nstern>: 0.3151941142403443

Betrag |pk+1,Nstern|: 0.6738536719442225

Betrag |pk,Nstern|: 0.3151941142403443

Vektor d'k (x,y,z): 0.33692683597211137, -0.2228758955694348, 0.8064502939065566

Betrag |d'k|: 0.9019730784042433

Vektor dk (x,y,z): 0.37354422658400965, -0.24709816834416315, 0.8940957476616884

Kreuzprodukt s*dk (x,y,z): 0.925535116013747, 0.15811594002795895, -0.3429813202028639

Skalarprodukt <s*dk,Nstern>: -0.8465477247784396

Skalarprodukt <s,dk>: -0.02858902807713498

Winkel Phi: -91.93421859340933

Betrag |Phi|: 91.93421859340933

Vorzeichen vk von <pk,Nstern> (0, -1, +1): 1

Vorzeichen vk+1 von <pk+1,Nstern> (0, -1, +1): 1

-------------------- Grenzpunkt i: 3, Grenzpunkt j: 4 --------------------

Skalarprodukt <s,N>: 0.24666583099588252

Skalarprodukt <s,s>: 1.00006942

Vektor Nstern (x,y,z): -0.6089361960364841, -0.23744871181016225, 0.7155242901082678

Skalarprodukt <pk+1,Nstern>: 0.7995257864374837

Skalarprodukt <pk,Nstern>: 0.6738536719442225

Betrag |pk+1,Nstern|: 0.7995257864374837

Betrag |pk,Nstern|: 0.6738536719442225

Vektor d'k (x,y,z): -0.29178719916856094, -1.0707403592816191, 0.9022769412955634

Betrag |d'k|: 1.4302896089208093

Vektor dk (x,y,z): -0.20400567643690146, -0.7486178691387688, 0.6308351369317119

Kreuzprodukt s*dk (x,y,z): 0.8035925916123443, -0.011035083646251352, 0.24677824729884418

Skalarprodukt <s*dk,Nstern>: -0.31014051929033293

Skalarprodukt <s,dk>: -0.5415598676921943

Winkel Phi: -150.20109474849926

Betrag |Phi|: 150.20109474849926

Vorzeichen vk von <pk,Nstern> (0, -1, +1): 1

Vorzeichen vk+1 von <pk+1,Nstern> (0, -1, +1): 1

-------------------- Grenzpunkt i: 4, Grenzpunkt j: 5 --------------------

Skalarprodukt <s,N>: 0.24666583099588252

Skalarprodukt <s,s>: 1.00006942

Vektor Nstern (x,y,z): -0.6089361960364841, -0.23744871181016225, 0.7155242901082678

Skalarprodukt <pk+1,Nstern>: 0.5101047145417327

Skalarprodukt <pk,Nstern>: 0.7995257864374837

Betrag |pk+1,Nstern|: 0.5101047145417327

Betrag |pk,Nstern|: 0.7995257864374837

Vektor d'k (x,y,z): -0.620644713910417, -1.0749881779418922, 0.2550523572708663

Betrag |d'k|: 1.2672218229533447

Vektor dk (x,y,z): -0.4897680127256357, -0.8483030819628411, 0.20126891176514763

Kreuzprodukt s*dk (x,y,z): 0.4161866494469646, -0.11487177517486723, 0.5285904632670687

Skalarprodukt <s*dk,Nstern>: 0.15206435587038425

Skalarprodukt <s,dk>: -0.7309272810843701

Winkel Phi: 168.24765102810687

Betrag |Phi|: 168.24765102810687

Vorzeichen vk von <pk,Nstern> (0, -1, +1): 1

Vorzeichen vk+1 von <pk+1,Nstern> (0, -1, +1): 1

-------------------- Grenzpunkt i: 5, Grenzpunkt j: 1 --------------------

Skalarprodukt <s,N>: 0.24666583099588252

Skalarprodukt <s,s>: 1.00006942

Vektor Nstern (x,y,z): -0.6089361960364841, -0.23744871181016225, 0.7155242901082678

Skalarprodukt <pk+1,Nstern>: -0.6089361960364841

Skalarprodukt <pk,Nstern>: 0.5101047145417327

Betrag |pk+1,Nstern|: 0.6089361960364841

Betrag |pk,Nstern|: 0.5101047145417327

Vektor d'k (x,y,z): 0.20563661652349058, -0.5273542150514561, 0

Betrag |d'k|: 0.56602905065709

Vektor dk (x,y,z): 0.36329693022782455, -0.9316734087044876, 0

Kreuzprodukt s*dk (x,y,z): 0.24428476776231664, 0.09525645510573559, -0.2870443343245147

Skalarprodukt <s*dk,Nstern>: -0.3767595533344492

Skalarprodukt <s,dk>: -0.9213718739641428

Winkel Phi: -157.7597826894389

Betrag |Phi|: 157.7597826894389

Vorzeichen vk von <pk,Nstern> (0, -1, +1): 1

Vorzeichen vk+1 von <pk+1,Nstern> (0, -1, +1): -1

Ende der Berechnung der Parameter für alle Grenzpunkte

Liste der Winkel Phi: -82.96909044086752,-91.93421859340933,-150.20109474849926,168.24765102810687,-157.7597826894389

Liste der Winkel |Phi|: 82.96909044086752,91.93421859340933,150.20109474849926,168.24765102810687,157.7597826894389

Liste der Vorzeichen vk+1: 1,1,1,1,-1

Liste der Vorzeichen vk: -1,1,1,1,1

Position des betragsmäßig kleinsten Winkels |Phi|: 1

Vorzeichen Phi min des 1. Elements (0, -1, +1): -1

Vorzeichen vk+1 des 1. Elements (0, -1, +1): 1

Vorzeichen vk des 1. Elements (0, -1, +1): -1

Vorzeichen vk+1 ist ungleich 0.

Vorzeichen Phi min ist nicht gleich Vorzeichen vk+1.

Test Westseite: Sonne befindet sich nicht im Verschattungsbereich.

Kommentiert 7 Jun von Guitardoc

... hast Du aber Glück, dass ich nochmal rein geschaut habe ;-)

mein Code (bisher noch unaufgeräumt)

int main()
{
    using namespace std;
    auto breitengrad = 50. * degree;  // ... irgendwo in Deutschland

    initializer_list< tuple< Angle, Angle > > pk = { 
        { 180. * degree, 0. * degree},
        { 180. * degree, 60. * degree},
        { 270. * degree, 45. * degree},
        { 300. * degree, 30. * degree},
        { 300. * degree, 0. * degree}
    };
    auto sonne = make_vector(make_tuple(274. * degree, 15.2 * degree));  // Vektor s
    vector< Vector > pks;
    (cout << fixed).precision(2); // für die Ausgabe in Geoknecht3D auf 2 Nachkommastellen begrenzen
    for (const auto& azimutElevation : pk)
    {
        auto v = make_vector(azimutElevation);
        cout << v << " Geoknecht: " << gk::vektor0(gkFACTOR *v) << endl;
        pks.push_back(v);
    }
    auto N = Vector(-cos(breitengrad), 0, sin(breitengrad));    // dieser Vektor ist parallel zur N-S-Achse der Erde
    cout << "N=" << N << " Geoknecht: " << gk::vektor0(gkFACTOR * N) << endl;
    cout << endl;

    cout << gk::vektor0( 5.*sonne) << "{cc0}" << endl; // Sonnenvektor s in Geoknecht3D

    auto Ns = unit(N - sonne * (sonne * N) / (sonne * sonne));
    cout << gk::vektor0(gkFACTOR * Ns) << "{c0c}" << endl; // N*
    cout << endl;
    double minPhi = 400;
    double phi = 0;
    bool liegtDrin = false;
    for (auto pk = begin(pks); pk != end(pks); ++pk)
    {
        auto pk1 = next(pk, pks); // p_{k+1}
        auto vk = sgn(*pk * Ns);
        auto vk1 = sgn(*pk1 * Ns);
        if (vk * vk1 < 0 || (vk == 0 && vk1 != 0))
        {
            auto dk = unit((abs(*pk1 * Ns) * *pk + abs(*pk * Ns) * *pk1));
            cout << gk::vektor0(gkFACTOR * dk) << "{c00}" << endl;  // 'dk' Geoknecht3D
            cout << gk::dreieck(Vector(), *pk * gkFACTOR, *pk1*gkFACTOR) << "{88f}" << endl;    // Polygonsgmet in Geoknecht3D
            phi = atan2((sonne % dk) * Ns, sonne * dk) * (180 / PI);
            if (abs(phi) < minPhi)
            {
                minPhi = abs(phi);
                auto vd = (vk1 != 0 ? vk1 : -vk);
                liegtDrin = (vd == sgn(phi));
            }
        }
    }
    cout << "Sonne liegt " << (liegtDrin ? "drin": "DRAUSSEN") << endl;

    return 0;
}

Das ist in C++. Die Unterprogramme und die Klasse CVevtor und CAngle sind nicht dabei, aber ich hoffe es ist einigermaßen selbst erklärend.

Hier der Link zum Geoknecht3D-Bild. Die $\vec d_k$ sind die roten Vektoren in den lila Polygonsegmenten. Die Vektoren wurden alle mit $\operatorname{gFAKTOR = 6}$ multipliziert, damit man es besser sehen kann.

Gruß Werner

Kommentiert 7 Jun von Werner-Salomon

Hallo Werner,

Danke, dass du noch mal reingeschaut und den Quelltext gepostet hast - der hat mir sehr geholfen und ich habe meinen Fehler möglicherweise gefunden - dein Quelltext hat dazu die Lösung gebracht! Ich hatte die Prüfung, ob das Polygonstück überhaupt relevant ist, unter den Tisch fallen lassen. Irgendwie muss ich das mit gelöscht haben... Daher stimmten die Ergebnisse manchmal, manchmal aber auch nicht.

Ich werde mal noch weiter testen, aber im Moment stimmen die Ergebnisse.

Ich habe aber trotzdem noch zwei Fragen. Du hast geschrieben

$$\vec{N}^* = \vec N - \frac{\left< \vec{s},\,\vec{N}\right>}{|\vec{s}|^2} \vec{s}$$

In deinem Quelltext steht dafür

auto Ns = unit(N - sonne * (sonne * N) / (sonne * sonne))

Ich gehe davon aus, dass die Prozedur dann intern die Berechnung entsprechend auseinandernimmt, also das was in Klammern steht als Skalarprodukte berechnet und dann den Rest, der nicht in Klammern steht, als Vektoren betrachtet? Nicht dass ich hier schon einen Fehler drin habe.

Dann noch eine Frage hierzu:

$$v_{k} \cdot v_{k+1} = -1 \quad \lor \quad v_{k}= 0 \land v_{k+1} \ne 0$$

Die Bedingung muss erfüllt sein - aber wenn

$$v_{k}= 0$$ oder $$v_{k+1}= 0$$

dann kann das doch nie erfüllt sein, insofern müsste doch

$$v_{k} \cdot v_{k+1} = -1$$

als Test ausreichen? Oder hab ich schon wieder einen Denkfehler?

Kommentiert 8 Jun von Guitardoc

Ich vermute, ich habe das Problem entdeckt. Kann es sein, dass die Richtungen der Vektoren, welche das Polygon definieren, nicht berücksichtigt werden?

Mein Haus ist exakt nach den vier Himmelsrichtungen ausgerichtet. Bei bestimmten Sonnenständen wurde nicht nur angezeigt, dass die Sonne im Westen ins Fenster rein scheint, sondern auch in manche Fenster im Osten. Da die Sonne im Moment tagsüber recht hoch steht, betraf das nur manchmal die Fenster auf der Ostseite (wenn die Sonne abends tiefer stand).

Es dürften eigentlich nur Polygone berücksichtig werden, wo der Sonnenvektor entgegen der Vektoren, welche den Verschattungsbereich definieren, gerichtet ist.

Für meinen Anwendungsfall ist es aber vermutlich am einfachsten, zu bestimmen, ob sich die Sonne in einem Bereich von 180° Azimut um das bestreffende Fenster herum befindet. Und nur dann wird der Verschattungstest überhaupt gemacht. Dann sollte dieser Fehler ausgeschlossen sein.

Kommentiert 9 Jun von Guitardoc

Hallo, da bin ich wieder ...

Ich habe scheinbar immer noch irgendwo einen Fehler drin. Könntest du bitte einmal dein Programm mit folgenden Daten laufen lassen:

Definition des Polygons (Azimut, Elevation): 10,15;10,40;90,65;150,35;145,30;170,15
Sonne (Azimut, Elevation): 109,44

Der Polygonzug sieht ein wenig wild aus. Bist Du sicher, dass das Polygon so korrekt ist? Die roten Vektoren sind die $\vec d_k$. Die Öffnng ist dan der Stelle seht schmal. Guckst Du hier. Die grüne Ebene ist die Sonnenebene

Bei mir ergibt sich, dass sich der Punkt nicht innerhalb des Polygons befindet

Mein Programm sagt dasselbe und man sieht es auch ganz klar in der Szene!

Definition des Polygons (Azimut, Elevation): 195,25;270,55;305,40;305,20;252,20;252,40;245,40;245,25

Sonne (Azimut, Elevation): 281,21.5
Auch hier ergibt der Algorithmus, dass die Sonne nicht im Verschattungsbereich wäre

Die Sonne liegt auch hier draußen, wenn auch nur knapp. Auch hier sieht das Polygon etwas 'wild' aus; siehe hier. In dieser Szene gibt es vier Vektoren $\vec d_k$ (rot)

Bestätige doch bitte zunächst mal, ob Deine Polygone so sind, wie sie in den beiden Geoknecht3D-Szenen aussehen.

Gruß Werner

Kommentiert 17 Jun von Werner-Salomon

Bei der ersten Variante sollten die Punkte eigentlich so aussehen:

Screenshot at Jun 18 09-36-58.png
Ich schaue dabei direkt nach Osten aus dem Fenster. Die jeweils oberen Werte sind die aus dem Fenster heraus gemessenen Winkel (horizontal, vertikal) und in Klammern darunter umgerechnet in Azimut, Elevation.

Die Sonne sollte sich mit 109,44 doch aber innerhalb der (hier mit Geraden verbundenen) umrandeten Fläche befinden?

Bei der zweiten Variante sieht das Polygon so aus:

Hier schaue ich jetzt direkt nach Westen. Oben sind wieder die aus dem Fenster gemessenen Winkel (horizontal, vertikal) und darunter in Klammern die umgerechneten Punkte in Azimut, Elevation.

Auch hier sollte sich die Sonne mit 281, 21.5 doch gerade noch so in der umrandeten Fläche befinden?

Oder hab ich da schon wieder einen Denkfehler drin?

Kommentiert 18 Jun von Guitardoc

Ein anderes Problem?

Stell deine Frage

Wie kann ich ermitteln, ob sich ein Punkt (die Sonne) innerhalb eines Polygons auf einer Kugeloberfläche befindet?

1 Antwort

Ähnliche Fragen

Eingabetools:

Beliebte Fragen:

Heiße Lounge-Fragen: