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 vor 6 Tagen 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 vor 5 Tagen 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 vor 1 Tag 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 vor 1 Tag 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 vor 3 Stunden 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 vor 1 Stunde 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 vor 18 Minuten von Werner-Salomon

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: