• Willkommen im Geoclub - dem größten deutschsprachigen Geocaching-Forum. Registriere dich kostenlos, um alle Inhalte zu sehen und neue Beiträge zu erstellen.

Schnittpunkt zweier Orthodrome?

jmsanta

jmsanta

Geoguru
muza schrieb:
Code: 
P - Lon:-171.0 Lat:49.00432103313932
Q - Lon:8.999999999999993 Lat:-49.00432103313933
Sieht doch ganz gut aus.
Zum Vergrößern anklicken....
bis du dir gaaaanz sicher, daß die 15. Stelle hinter dem Komma bei Q - Lon wirklich eine Drei ist? :???:
 
OP
M

muza

Geocacher
Ja, habs direkt von der Ausgabekonsole kopiert.
Denke mal, dass das ein Rundungsfehler in Java bei der Berechnung ist.

Hier mal ein Beispiel mit Schnitt am Pol
Code: 
A - Lon:8.0 Lat:49.0
B - Lon:-172.0 Lat:49.0
C - Lon:30.0 Lat:5.0
D - Lon:-150.0 Lat:5.0

P - Lon:14.036243467926477 Lat:90.0
Q - Lon:-165.96375653207352 Lat:-90.0
Oder hier mal eins
Code: 
A - Lon:-25.0 Lat:-33.25
B - Lon:155.0 Lat:80.0
C - Lon:20.0 Lat:50.0
D - Lon:-160.0 Lat:5.0

P - Lon:26.56505117707799 Lat:90.0
Q - Lon:-153.43494882292202 Lat:-90.0


Wenn x bei atan2( y, x ) gegen 0 geht, dann muss man wahrscheinlich anhand des Vorzeichens von y bestimmen, ob die Länge 90° oder -90° beträgt, oder?

Denke mal, wegen der Ungenauigkeit wird hier nichts abgefangen
Code: 
A - Lon:90.0 Lat:-30.0
B - Lon:90.0 Lat:30.0
C - Lon:80.0 Lat:0.0
D - Lon:100.0 Lat:10.0

P - Lon:-90.00000000000001 Lat:5.11568442019437
Q - Lon:89.99999999999999 Lat:-5.115684420194341
 
t31

t31

Geowizard
ich habe das mal in Excel versucht nachzubasteln - nicht erschrecken, Excel-VBA ist halt anders :eek:ps:
Code: 
Public Const r = 6371000.8
Public Const pi = 3.14159265358979
Public Const rad = 57.2957795130823
Public Function WGS84toXYZ(i, Lon, Lat)
    Lon = Lon + 180
    Lat = Lat + 90
    Select Case i
        Case 0: WGS84toXYZ = r * Sin(Lat / rad) * Cos(Lon / rad)
        Case 1: WGS84toXYZ = r * Sin(Lat / rad) * Sin(Lon / rad)
        Case 2: WGS84toXYZ = r * Cos(Lat / rad)
    End Select
End Function
Public Function XYZtoWGS84(i, x, y, z)
    Select Case i
        Case 0: XYZtoWGS84 = rad * IIf(x, Atn(y / x) - (x > 0) * pi, pi / 2 + (y > 0) * pi) - 180
        Case 1: XYZtoWGS84 = rad * WorksheetFunction.Acos(z / r) - 90
    End Select
End Function
Sub test()
    LatA = 49
    LonA = 8
    LatB = 49
    LonB = 10
    LatC = 50
    LonC = 9
    LatD = 48
    LonD = 9
    Dim a(2), b(2), C(2), D(2)
    For i = 0 To 2
        a(i) = WGS84toXYZ(i, LonA, LatA)
        b(i) = WGS84toXYZ(i, LonB, LatB)
        C(i) = WGS84toXYZ(i, LonC, LatC)
        D(i) = WGS84toXYZ(i, LonD, LatD)
    Next i
    e1 = a(1) * b(2) - a(2) * b(1)
    e2 = a(2) * b(0) - a(0) * b(2)
    e3 = a(0) * b(1) - a(1) * b(0)
    f1 = e1 * C(0) + e2 * C(1) + e3 * C(2)
    f2 = e1 * D(0) + e2 * D(1) + e3 * D(2)
    n = -(f2 / f1)
    g1 = n * C(0) + D(0)
    g2 = n * C(1) + D(1)
    g3 = n * C(2) + D(2)
    n = Sqr(g1 * g1 + g2 * g2 + g3 * g3)
    P1 = r * g1 / n
    P2 = r * g2 / n
    P3 = r * g3 / n
    Q1 = -P1
    Q2 = -P2
    Q3 = -P3
    Dim P(1), Q(1)
    For i = 0 To 1
        P(i) = XYZtoWGS84(i, P1, P2, P3)
        Q(i) = XYZtoWGS84(i, Q1, Q2, Q3)
    Next i
    Debug.Print "P - Lon: " & P(0) & " Lat: " & P(1)
    Debug.Print "Q - Lon: " & Q(0) & " Lat: " & Q(1)
End Sub
Als Ergebnis erhalte ich jedoch:
P - Lon: -169,672375413926 Lat: -48,8924602352394
Q - Lon: 10,3276245860744 Lat: 48,8924602352394
Zum Vergrößern anklicken....
Keine Ahnung, ob hier die Rechengenauigkeit (Rundungsfehler) oder meine Schussligkeit bei der Umsetzung in VBA zugeschlagen hat.
 
t31

t31

Geowizard
Dank http://www.office-loesung.de klappt es nun auch unter Excel.
P - Lon: -171 Lat: 49,0043210331393
Q - Lon: 9 Lat: -49,0043210331393
Zum Vergrößern anklicken....
Code: 
Public Const r = 6371000.8
Public Const pi = 3.14159265358979
Public Const rad = 57.2957795130823
Public Function WGS84toXYZ(Lon, Lat) As Variant
Dim arr(2)
    Lon = Lon + 180
    Lat = Lat + 90
    arr(0) = r * Sin(Lat / rad) * Cos(Lon / rad)
    arr(1) = r * Sin(Lat / rad) * Sin(Lon / rad)
    arr(2) = r * Cos(Lat / rad)
    WGS84toXYZ = arr
End Function
Public Function XYZtoWGS84(x, y, z) As Variant
    Dim arr(1)
    arr(0) = rad * IIf(x, Atn(y / x) - (x > 0) * pi, pi / 2 + (y > 0) * pi) - 180
    arr(1) = rad * WorksheetFunction.Acos(z / r) - 90
    XYZtoWGS84 = arr
End Function
Sub test()
    LatA = 49
    LonA = 8
    LatB = 49
    LonB = 10
    LatC = 50
    LonC = 9
    LatD = 48
    LonD = 9
    a = WGS84toXYZ(LonA, LatA)
    b = WGS84toXYZ(LonB, LatB)
    C = WGS84toXYZ(LonC, LatC)
    D = WGS84toXYZ(LonD, LatD)
    e1 = a(1) * b(2) - a(2) * b(1)
    e2 = a(2) * b(0) - a(0) * b(2)
    e3 = a(0) * b(1) - a(1) * b(0)
    f1 = e1 * C(0) + e2 * C(1) + e3 * C(2)
    f2 = e1 * D(0) + e2 * D(1) + e3 * D(2)
    n = -(f2 / f1)
    g1 = n * C(0) + D(0)
    g2 = n * C(1) + D(1)
    g3 = n * C(2) + D(2)
    n = Sqr(g1 * g1 + g2 * g2 + g3 * g3)
    P1 = r * g1 / n
    P2 = r * g2 / n
    P3 = r * g3 / n
    P = XYZtoWGS84(P1, P2, P3)
    Q = XYZtoWGS84(-P1,- P2,- P3)
    Debug.Print "P - Lon: " & P(0) & " Lat: " & P(1)
    Debug.Print "Q - Lon: " & Q(0) & " Lat: " & Q(1)
End Sub
 
KoenigDickBauch

KoenigDickBauch

Geoguru
So und hier kommt die Antwort auf die zweite Frage:

Code: 
program Orthodrome;
var
  LatA,LonA: real;
  LatB,LonB: real;
  LatC,LonC: real;
  LatD,LonD: real;

  A1,A2,A3: real;
  B1,B2,B3: real;
  C1,C2,C3: real;
  D1,D2,D3: real;

  e1,e2,e3: real;
  f1,f2: real;
  g1,g2,g3: real;

  v1,v2,v3: real;
  M11,M12,M13: real;
  M21,M22,M23: real;
  M31,M32,M33: real;

  n: real;
  r: real;
  w: real;
  CW,SW: real;

  P1,P2,P3: real;
  Q1,Q2,Q3: real;

  A1P,A2P,A3P,AW: real;
  B1P,B2P,B3P,BW: real;
  P1P,P2P,P3P,PW: real;
  Q1P,Q2P,Q3P,QW: real;

  A1PP,A2PP,A3PP: real;
  B1PP,B2PP,B3PP: real;
  P1PP,P2PP,P3PP: real;
  Q1PP,Q2PP,Q3PP: real;

  E1T,E2T,E3T: real;
  E1TT,E2TT,E3TT: real;

  LatP,LonP,HP: real;
  LatQ,LonQ,HQ: real;

begin
   r := 6371000.8;

   LatA:=49;
   LonA:=8;
   LatB:=49;
   LonB:=10;
   LatC:=50;
   LonC:=9;
   LatD:=48;
   LonD:=9;

   Writeln('A: ',WGS84(LatA,LonA));
   Writeln('B: ',WGS84(LatB,LonB));
   Writeln('C: ',WGS84(LatC,LonC));
   Writeln('D: ',WGS84(LatD,LonD));

   WGS84toXYZ(LatA,LonA,0,A1,A2,A3);
   WGS84toXYZ(LatB,LonB,0,B1,B2,B3);
   WGS84toXYZ(LatC,LonC,0,C1,C2,C3);
   WGS84toXYZ(LatD,LonD,0,D1,D2,D3);

   e1:=a2*b3-a3*b2;
   e2:=a3*b1-a1*b3;
   e3:=a1*b2-a2*b1;

   f1:=e1*c1+e2*c2+e3*c3;
   f2:=e1*d1+e2*d2+e3*d3;

   n:=-(f2/f1);

   g1:=n*c1+d1;
   g2:=n*c2+d2;
   g3:=n*c3+d3;

   n := sqrt(g1*g1+g2*g2+g3*g3);
   g1:=g1/n;
   g2:=g2/n;
   g3:=g3/n;
   
   P1:=r*g1;
   P2:=r*g2;
   P3:=r*g3;

   Q1:=-P1;
   Q2:=-P2;
   Q3:=-P3;

(* auf XY abbilden *)
   n := sqrt(e1*e1+e2*e2+e3*e3);
   e1:=e1/n;
   e2:=e2/n;
   e3:=e3/n;

(* um z drehen *)
   w:=-arctan2(E2,E1);
   cw:=cos(w);
   sw:=sin(w);

   M11:=cw; M12:=-sw; M13:=0;
   M21:=sw; M22:=cw;  M23:=0;
   M31:=0;  M32:=0;   M33:=1;

   E1T:= M11*E1 + M12*E2 + M13* E3;
   E2T:= M21*E1 + M22*E2 + M23* E3;
   E3T:= M31*E1 + M32*E2 + M33* E3;

   A1P:= M11*A1 + M12*A2 + M13* A3;
   A2P:= M21*A1 + M22*A2 + M23* A3;
   A3P:= M31*A1 + M32*A2 + M33* A3;

   B1P:= M11*B1 + M12*B2 + M13* B3;
   B2P:= M21*B1 + M22*B2 + M23* B3;
   B3P:= M31*B1 + M32*B2 + M33* B3;

   P1P:= M11*P1 + M12*P2 + M13* P3;
   P2P:= M21*P1 + M22*P2 + M23* P3;
   P3P:= M31*P1 + M32*P2 + M33* P3;

   Q1P:= M11*Q1 + M12*Q2 + M13* Q3;
   Q2P:= M21*Q1 + M22*Q2 + M23* Q3;
   Q3P:= M31*Q1 + M32*Q2 + M33* Q3;

(* drehen um y  *)
   w:=arctan2(E3T,E1T)-pi/2;
   cw:=cos(w);
   sw:=sin(w);

   M11:=cw;  M12:=0;  M13:=sw;
   M21:=0;   M22:=1;  M23:=0;
   M31:=-sw; M32:=0;  M33:=cw;

   E1TT:= M11*E1T + M12*E2T + M13* E3T;
   E2TT:= M21*E1T + M22*E2T + M23* E3T;
   E3TT:= M31*E1T + M32*E2T + M33* E3T;

   A1PP:= M11*A1P + M12*A2P + M13* A3P;
   A2PP:= M21*A1P + M22*A2P + M23* A3P;
   A3PP:= M31*A1P + M32*A2P + M33* A3P;

   B1PP:= M11*B1P + M12*B2P + M13* B3P;
   B2PP:= M21*B1P + M22*B2P + M23* B3P;
   B3PP:= M31*B1P + M32*B2P + M33* B3P;

   P1PP:= M11*P1P + M12*P2P + M13* P3P;
   P2PP:= M21*P1P + M22*P2P + M23* P3P;
   P3PP:= M31*P1P + M32*P2P + M33* P3P;

   Q1PP:= M11*Q1P + M12*Q2P + M13* Q3P;
   Q2PP:= M21*Q1P + M22*Q2P + M23* Q3P;
   Q3PP:= M31*Q1P + M32*Q2P + M33* Q3P;

   writeln;
   writeln('E: ',E1 :12:6,E2 :12:6,E3 :12:6);
   writeln;
   writeln('A: ',A1:12:2,A2:12:2,A3:12:2);
   writeln('E: ',B1:12:2,B2:12:2,B3:12:2);
   writeln('P: ',P1:12:2,P2:12:2,P3:12:2);
   writeln('Q: ',Q1:12:2,Q2:12:2,Q3:12:2);
   writeln;
   writeln('E: ',E1T:12:6,E2T:12:6,E3T:12:6);
   writeln;
   writeln('A: ',A1P:12:2,A2P:12:2,A3P:12:2);
   writeln('B: ',B1P:12:2,B2P:12:2,B3P:12:2);
   writeln('P: ',P1P:12:2,P2P:12:2,P3P:12:2);
   writeln('Q: ',Q1P:12:2,Q2P:12:2,Q3P:12:2);

   writeln;
   writeln('E: ',E1TT:12:6,E2TT:12:6,E3TT:12:6);
   writeln;
   writeln('A: ',A1PP:12:2,A2PP:12:2,A3PP:12:2);
   writeln('B: ',B1PP:12:2,B2PP:12:2,B3PP:12:2);
   writeln('P: ',P1PP:12:2,P2PP:12:2,P3PP:12:2);
   writeln('Q: ',Q1PP:12:2,Q2PP:12:2,Q3PP:12:2);

   AW:=arctan2(A2PP,A1PP);
   BW:=arctan2(B2PP,B1PP);
   PW:=arctan2(P2PP,P1PP);
   QW:=arctan2(Q2PP,Q1PP);

   if AW>BW then begin
      BW:=BW+2*PI
   end;
   if AW>PW then begin
      PW:=PW+2*PI
   end;
   if AW>QW then begin
      QW:=QW+2*PI
   end;

   writeln;
   writeln('A: ',AW/Pi*180:7:2);
   writeln('B: ',BW/Pi*180:7:2);
   writeln('P: ',PW/Pi*180:7:2);
   writeln('Q: ',QW/Pi*180:7:2);
   writeln;

   if (AW<=PW) and (PW<=BW) then begin
      writeln('P liegt zwischen A und B')
   end;
   if (AW<=QW) and (QW<=BW) then begin
      writeln('Q liegt zwischen A und B')
   end;

   XYZtoWGS84(P1,P2,P3,LatP,LonP,HP);
   XYZtoWGS84(Q1,Q2,Q3,LatQ,LonQ,HQ);

   writeln;
   Writeln('P: ',WGS84(LatP,LonP));
   Writeln('Q: ',WGS84(LatQ,LonQ));
   writeln;
   writeln('HP: ',HP:7:0,' m');
   writeln('HQ: ',HQ:7:0,' m');

end.
Kann man natürlich viel schöner schreiben.

und das Ergebnis:
Code: 
A: N 49° 0.0000 E 8° 0.0000
B: N 49° 0.0000 E 10° 0.0000
C: N 50° 0.0000 E 9° 0.0000
D: N 48° 0.0000 E 9° 0.0000

E:    -0.743307   -0.117728    0.658510

A:   4151634.45   583474.17  4790558.75
E:   4128742.43   728008.69  4790558.75
P:   4143714.06   656299.83  4794638.20
Q:  -4143714.06  -656299.83 -4794638.20

E:     0.752572    0.000000    0.658510

A:  -4191796.41    73168.08  4790558.75
B:  -4191796.41   -73168.08  4790558.75
P:  -4195365.98       -0.00  4794638.20
Q:   4195365.98        0.00 -4794638.20

E:     0.000000    0.000000    1.000000

A:  -6365580.12    73168.08       -0.00
B:  -6365580.12   -73168.08       -0.00
P:  -6371000.80       -0.00       -0.00
Q:   6371000.80        0.00        0.00

A:  179.34
B:  180.66
P:  180.00
Q:  360.00

P liegt zwischen A und B

P: N 49° 0.2503 E 9° 0.0000
Q: S 49° 0.2503 W 171° 0.0000

HP:    5002 m
HQ:    5002 m
 
KoenigDickBauch

KoenigDickBauch

Geoguru
Hallo,

und nun ist es Bestandteil der Mopsos Sriptsprache:

Code: 
program T16;
var
  r: real;

  LatA,LonA: real;
  LatB,LonB: real;
  LatC,LonC: real;
  LatD,LonD: real;

  LatP,LonP,HP: real;
  LatQ,LonQ,HQ: real;

  Inside: integer;

begin
   r := 6371000.8;

   LatA:=49;
   LonA:=8;
   LatB:=49;
   LonB:=10;
   LatC:=50;
   LonC:=9;
   LatD:=48;
   LonD:=9;

   Orthodrome(LatA,LonA,LatB,LonB,LatC,LonC,LatD,LonD,
     LatP,LonP,LatQ,LonQ, Inside);

   Writeln('P: ',WGS84(LatP,LonP));
   Writeln('Q: ',WGS84(LatQ,LonQ));

   if (Inside = 1) then begin
      writeln('P liegt zwischen A und B')
   end;
   if (Inside = 2) then begin
      writeln('Q liegt zwischen A und B')
   end;

end.

Antwort ist:
Code: 
P: N 49° 0.2503 E 9° 0.0000
Q: S 49° 0.2503 W 171° 0.0000
P liegt zwischen A und B
 
Wallraff

Wallraff

Geocacher
Hallo,

programmtechnisch kann man Thomas nichts mehr am Zeug flicken...
Ich glaube, bei Dir ist unter WGStoXYZ sogar WGS84 drin, da eine Höhe von 5 km rausgerechnet wird.
Bei muza wohl nur eine Kugel ... ?

Ich kann nur noch mit der Höheren Geodäsie zuschlagen. Auf dem Ellipsoid sind die Großkreise Geodätische. Und die schlängeln sich je nach Startpunkt, Länge, Azimut unterschiedlich s-förmig übers Ellipsoid. Damit liegen sie nicht in einer aus Anfangs- End- und Mittelpunkt aufgespannten Ebene.
Wie groß mag der Unterschied sein ? Na, das löse ich auch nicht.

Als Buße für meine hiesigen Beiträge könnte ich höchstens versuchen, dass von mir genannte Verfahren in ein Basic-Programm zu fassen. Basic. Aber all die Sonderfälle, die abgefangen werden wollen ...


Grüße Wallraff
 
KoenigDickBauch

KoenigDickBauch

Geoguru
Wallraff schrieb:
programmtechnisch kann man Thomas nichts mehr am Zeug flicken...
Zum Vergrößern anklicken....
Doch sicher, das geht.
Wallraff schrieb:
Ich glaube, bei Dir ist unter WGStoXYZ sogar WGS84 drin, da eine Höhe von 5 km rausgerechnet wird.
Zum Vergrößern anklicken....
Ist der der erste Schritt für die Umrechnung von WGS84toUTM

Wallraff schrieb:
Bei muza wohl nur eine Kugel ... ?
Zum Vergrößern anklicken....
Interessanter Stoff hierzu:

http://www.ottmarlabonde.de/L1/Rotationsellipsoide.html

Wallraff schrieb:
Ich kann nur noch mit der Höheren Geodäsie zuschlagen. Auf dem Ellipsoid sind die Großkreise Geodätische. Und die schlängeln sich je nach Startpunkt, Länge, Azimut unterschiedlich s-förmig übers Ellipsoid. Damit liegen sie nicht in einer aus Anfangs- End- und Mittelpunkt aufgespannten Ebene.
Zum Vergrößern anklicken....
Ja, das ist sicher eine spannende Aufgabe. Hmm....

Wallraff schrieb:
Wie groß mag der Unterschied sein ? Na, das löse ich auch nicht.
Zum Vergrößern anklicken....
Das wäre die Frage zu einem D5 Mystery. :D

Und wenn du dich schon mit Ellipsen beschäftigst. Dann mach mal diesen Mystery. Wird Zeit, das er mal wieder gefunden wird.

Wallraff schrieb:
Als Buße für meine hiesigen Beiträge könnte ich höchstens versuchen, dass von mir genannte Verfahren in ein Basic-Programm zu fassen. Basic. Aber all die Sonderfälle, die abgefangen werden wollen ...
Zum Vergrößern anklicken....
Mir würde adhoc nur einer einfallen alle Punkte liegen auf einem Großkreis. Das wird man nur in den Winkel der beiden Ebenennormalen herausbekommen. Nur was heißt gleich?
 
Wallraff

Wallraff

Geocacher
Hallo,

wenn der polnächste Punkt gesucht wird, wird's bei Orthodromen am Äquator kritisch.
Auch der Anblick des Nenners in der zitierten Formel weckt üble Erinnerungen,
(IF Divide by Zero THEN Ergebnis ziemlich groß)

Ich hab jetzt zwar frei und fast alle Weihnachtsgeschenke, aber vor dem Fest herumprogrammieren...

Gute Lösungen sind ja da.

Gute Nacht
Wallraff
 
Du musst dich einloggen oder registrieren, um hier zu antworten.
Oben