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

Question

Aufgabe:

Punkt in sphärischem Polygon

Problem/Ansatz:

Ich habe ein Polygon, welches sich auf einer Kugeloberfläche befindet. Die Koordinaten der Eckpunkte des Polygons sind durch Winkel (Azimut, Elevation) bestimmt. Wie kann ich ermitteln, ob sich ein Punkt (die Sonne) innerhalb dieses Polygons befindet?

Der Ansatz mit dem PiP Test nach Jordan scheitert, da dieser ein Polygon in einer Ebene voraussetzt. Gibt es ein vergleichbares Verfahren auch für Polygone auf einer Kugeloberfläche?

Answer 1

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

Comments

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!

---

Vielen Dank für Dein Feedback. Ich melde mich noch wegen Deiner Fragen ... Du musst Dich aber noch ein wenig gedulden.

---

Das ist kein Problem - ich übe mich gern in Geduld wenn ich geholfen bekomme… :-)

---

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;
}

---

Ich werd‘s mir ansehen und melde mich frühestens Montag Abend.

---

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

---

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;
}

---

Ich denke, hier ist irgendwo ein Fehler. Ich habe den Algorithmus mal mit Testdaten laufen lassen (Versuch macht kluch... ☺):

Zunächst die Koordinaten des Bereichs - da hab ich einen Teil der westlichen Hemisphäre genommen, da die Sonne gerade im Westen steht (die Werte passen vermutlich nicht zum tatsächlich möglichen Sonnenlauf, aber das sollte erstmal keine Rolle spielen vermute ich) :

Azimut, Elevation:

180°, 60°

270°, 45°

300°, 30°

Das ergibt x, y, z:

0.5, 0, 0.86

0, -0.71, 0.71

-0.43, -0.75, 0.5

Die Sonne steht im Moment im Westen:

Azimut, Elevation:

238°, 38°

Das ergibt x, y, z:

0.44, -0.64, 0.63

Es ergeben sich überhaupt keine Winkel, sprich die Liste der Winkel ist leer. Das dürfte eigentlich nicht sein...

---

Ahhh, hab den Fehler entdeckt. Es fehlten natürlich noch die Bereichskoordinaten 180,0 sowie 300,0. Jetzt gibt es auch berechnete Winkel - aber die sind alle positiv...

---

Ich habe mir das alles die Nacht noch mal durch den Kopf gehen lassen. Mein Fehler liegt ja schon bei dem "Vektor, der in Richtung der Rotationsachse der Erde zeigt" (du siehst, meine Mathevorlesungen wegen wirklich weit zurück). Ich gehe davon aus, das damit der Vektor vom Beobachter (also ich an meinem Fenster :-)) weg in Richtung Koordinatenursprung gemeint ist. Das würde bedeuten, dass der Vektor N mit

N.x = cos(h) * cos(A)

N.y = - cos(h) * sin(A)

N.z = - sin (h)

mit

h = Breitengrad

A = 360° - Längengrad

definiert ist. Wenn das so richtig ist, dann gibt es woanders immer noch einen Fehler, denn es ergeben sich nach wie vor nur positive Winkel.

Ich weiß nicht mehr weiter und hoffe noch mal auf deine Hilfe...

---

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... :-)

---

Ich werde mich heute Abend noch melden. Zur Zeit bin ich gerade offline beschäftigt.

---

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?

Nein - nicht richtig. Ich unterstelle, das Szenario spielt nicht direkt am Nordpol ;-)

mal angenommen Du befindest Dich auf dem 50'sten Breitengrad (mitten in Deutschland) bei 50° nördlicher Breite, dann ist \(\vec N\)$$\vec N = \begin{pmatrix} -\cos(50°)\\0 \\ \sin(50°) \end{pmatrix}$$in Deinem System, bei dem \(x\) nach Süden und \(z\) nach oben zeigt.

(geht gleich weiter)

---

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)

---

Die Sonne steht im Moment im Westen:
Azimut, Elevation:
238°, 38°

Das ergibt x, y, z:
0.44, -0.64, 0.63

Bei den Koordinaten des Bereichs bekomme ich die selben Werte wie Du auch. Nur bei der Sonne nicht. Ich habe

bei 238°,38°  → Sonnenvektor: 0.42 -0.67 0.62

nah drann, aber doch daneben !?

---

ich habe ein Ergebnis und auch ein Bild dazu

Die drei schwarzen Vektoren geben den Bereich an (das lila Dreieck). Der gelbe Vektor ist der Vektor zur Sonne und die gelbe Ebene zeigt die momentane Bewegungsebene des Sonnenvektors. Offensichtlich geht diese ganz am lila Dreieck vorbei. Und daher gibt es auch keinen Schnittpunkt.

Du kannst auf das Bild klicken. Dann öffnet sich Geoknecht3D und man kann die Szene mit der Maus rotieren und sich das so besser vorstellen.

ich rechne das morgen (Sonntag) nochmal mit anderen Werten.

Gruß Werner

---

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.

---

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.

---

Hallo Guitardoc,

ich habe den Fehler bei mir gefunden. Ich schrieb oben:

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>$$

Das darf nicht \(\vec N^*\) heißen, sondern \(\vec n_k\)! Es geht ja darum, den Sonnenvektor zu der Ebene zu berechnen, die durch \(p_{k}\) und \(p_{k+1}\) aufgespannt ist.

Aber es ist nicht das einzige Problem. Die Winkel passen noch nicht zur Skizze ... ich suche noch!

---

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.

---

Ich sehe gerade, dass die beiden Winkel falsch sind.

Phi 1 = 133.5668293844902°

Phi 2 = 61.46435809214322°

Da hab ich mich tatsächlich in der Zuordnung der Ausgabe der Werte vertan. Richtig sind die Werte für die Winkel wie in der ausführlichen Liste der Zwischenwerte angegeben.

Ahh - das hat sich mit deiner Antwort überschnitten. OK, dann suche ich erst mal nicht weiter, ob ich einen Fehler im Programm eingebaut habe und warte erst mal ab.

---

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?

---

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.)

---

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... :-)

---

Ich habe hier noch mal den Fall des \(d_{k=4}\) raus gezeichnet (der grüne Vektor)

als ich mir das alles im Kopf vorstellte, lag \(d_{k=4}\) (alias dk40) so rechts neben \(p_4\) - also genau in entgegen gesetzter Orientierung. Dann würden auch die gesuchten Winkel richtig heraus kommen.

Spontan hab ich schon wieder 'ne Menge Ideen wie man das behebt, aber erstmal denke ich das sauber zu Ende, bevor ich wieder ungelegte Eier poste.

---

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?

---

Hallo Guitardoc,

Leider habe ich aktuell den Kopf nicht frei, um mich mit Deinem Problem zu beschäftigen. Der Schlüssel zur Lösung liegt mit Sicherheit in den Vorzeichen/Richtungen der Rotationen.

Ich melde mich erst nach dem 1.Juni wieder - bis dahin bin ich mehr oder weniger offline.

Gruß Werner

---

Hallo Werner,

Solltest du doch noch mal Zeit haben um dich mit dem Problem zu beschäftigen, dann würde mich das sehr freuen. Ich komme einfach nicht weiter…

Viele Grüße!

---

Zeit habe ich nicht viel, aber ich denke noch drüber nach. Ich melde mich spätestens bis Mittwoch Abend.

---

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?

Für Deine Anwendung könnte es reichen. Als allgemeiner Ansatz taugt es nicht. Falls sich die \(p_k\) auf beiden Seiten dieser Projektionsebene befinden, lassen sich Fälle konstruieren, die zu falschen Ergebnissen führen.

(melde mich heute nochmal, ich arbeite noch was aus)

---

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.

---

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}\)?

---

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

---

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.

---

Großes Ooops - diesmal hab ich ein ungelegtes Ei gepostet. Der Sonnenvektor zeigte nach Osten aber das Polygon nach Westen. Mit Sonne.y = -0.9627 stimmt die Berechnung!

Ich teste demnächst mal noch mit anderen Werten, bin jetzt aber auch erst mal wieder offline.

---

Hallo Werner,

Ich habe trotzdem irgendwo einen Fehler in der Berechnung.Manchmal stimmen die Ergebnisse, manchmal aber auch nicht.

Du sagst, du hast den Algorithmus programmiert - könntest du mir diesen zur Verfügung stellen? Dann kann ich genauer verfolgen was wie berechnet wird. Ich finde den Fehler sonst einfach nicht... Vielleicht liegt es auch einfach am Koordinatensystem? Ich hatte angenommen:

Nord: -x

Ost: +y

Süd: +x

West: -y

Zenit: +z

---

... 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

---

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?

---

Hmm, hab schon wieder einen Fall, wo das Ergebnis nicht stimmt. Ich denke, es liegt an der Berechnung von N*. Stutzig macht mich, dass

$${|\vec{s}|^2}$$

bei mir immer 1 ergibt. Irgendwie wird N* anders berechnet als ich es umgesetzt habe...

---

Nächsten Fehler gefunden - ich hab doch tatsächlich $$\lor$$ und $$\land$$ vertauscht... Oh Mann, ich werde alt. Hatte mich ja schon gewundert, weil so macht die Fallunterscheidung ja gar keinen Sinn (siehe meine Bemerkung oben).

Mal sehen, was die Berechnungen jetzt sagen. Werde weiter testen.

---

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.

---

Hallo Werner,

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

Bei mir ergibt sich, dass sich der Punkt nicht innerhalb des Polygons befindet - er sollte aber innerhalb sein. Oder liegt es vielleicht an der einspringenden Ecke?

---

Ich habe noch eine Konstellation gefunden, wo etwas nicht stimmt:

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 - sie sollte es aber sein. Kannst du das bitte auch einmal testen?

---

Die letzte Konstellation ist korrekt, da hab ich mich vertan. Aber bei der vorletzten bin ich ratlos.

---

Ahhh, auch bei der letzten Konstellation stimmt doch was nicht. Die Sonne müsste doch im Bereich liegen, es wird bei mir aber ausgegeben, dass sie außerhalb wäre. Wenn du die Parameter auch mal durch deinen Algorithmus laufen lassen könntest? Das würde mir sehr helfen, den Fehler einzugrenzen. Auch hier gibt es wieder einspringende Bereiche...

---

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

---

Bei der ersten Variante sollten die Punkte eigentlich so aussehen:

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?

---

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

Nein, tut sie nicht! Betrachte dazu den letzten und ersten Punkt im Polygon. Also (170,15) und (10,15). Beide Punkte werden durch einen Großkreis verbunden und nicht durch einen Breitenkreis. D.h. die Projektion dieses Kreissegments auf Deine Ebene (Blick nach Osten) ist keine Gerade, sondern wieder ein Kreis (oder so ähnlich). Und dieser liegt mit der Elevation mittig höher als der Sonnenvektor. Und somit ist letztere draußen.

Schau Dir den kürzesten Flugweg von Rom nach New York an, Er führt nicht(!) über den Breitenkreis von ca. 40° bis 42° sondern an Island vorbei.

Tipp: mehr Zwischenpunkte im Polygon definieren.

---

Stimmt, du hast Recht! Ich war nur auf die oberen Begrenzungen fixiert (welche als Kreise ja den Bereich noch größer machen) und hab die unteren völlig außer acht gelassen! Manchmal sieht man den Wald vor Bäumen nicht - vielen lieben Dank, dass du es dir noch mal angesehen hast! Da werde ich noch mehr Punkte definieren müssen - das hatte ich ohnehin „später“ vor. Aber jetzt bin ich erst mal eine Woche nicht da und kann erst danach zur Tat schreiten…

---

Hallo Werner,

Ich wollte mich noch einmal melden und mich recht herzlich bei dir bedanken. Das System funktioniert, ich habe es ja jetzt mehrere Monate getestet.

Danke noch einmal, dass du dich meinem Problem angenommen und mir mit deinem Wissen geholfen hast. Ich weiß das wirklich sehr zu schätzen!

Viele liebe Grüße!

---

Danke für das Feedback. Es hat mich sehr gefreut, dass Du Dich noch einmal gemeldet hast :-)