Bogenlänge einer Ellipse ?
Jürgen HEYN
ich suche eine Formel zur Berrechnung der Bogenlänge einer Ellipse zwischen
zwei beliebigen Winkeln. (Benötigt werden 0° - 90°.)
Im Internet sind einige sehr interessante Abhandlungen über Ellipsen
verfügbar. Aus diesen Veröffentlichungen schließe ich, daß eine solche
Berechnung nur mit Hilfe der elliptischen Integrale zu bewerkstelligen ist.
Leider reicht mein mathematischer Fachverstand nicht aus, um die ganzen Ab-
und Herleitungen nachzuvollziehen. (Ich benötige auch keinen "Beweis", ich
vertraue den Mathematikern ;-))
Gibt es eine "einfache Formel" (die sich in Visual Basic umsetzten läßt), in
die ich nur bekannte Werte einsetzten muß um als Ergebnis die gewünschte
Bogenlänge zu erhalten ?
Folgende Werte sind bekannt oder können errechnet werden:
a Große Halbachse
b Kleine Halbachse
f Abplattung
e numerische Exzentrizität
phi geografische Breite
psi geozentrische Breite
r Radius der Ellipse bei phi/psi
U Umfang der Ellipse
In einem rechtwinkligen Koordinatensystem schneiden sich X-Achse
(waagerecht) und Y-Achse (senkrecht) bei x = 0 und y = 0.
Dieser Nullpunkt ist der Mittelpunkt einer Ellipse mit der großen Halbachse
a = 100 entlang der X-Achse. Die kleine Halbachse b = 65.
Darauf folgt für die Abplattung f = (a - b) / a = 0,35 und für die
numerische Exzentrizität e = Sqr((a^2 - b^2) / a^2) = 0,759934207678533.
An die Ellipse schmiegt sich eine Tangente. Die Senkrechte auf diese
Tangente am Berührungspunkt mit der Ellipse schneidet die X-Achse mit 60°
(geografische Breite phi).
Die geozentrische Breite psi (Winkel vom Mittelpunkt der Ellipse zum
Berührungspunkt der Tangente) errechnet sich psi = Arctan((1-e^2) *
tan(Phi)) = 36,1963473408467°.
Der Radius r für den Winkel psi = ((((1 - e^2) * a) / (Sqr(1 - e^2 *
Sin(phi)^2))) * Sin(phi)) / Sin(psi) = 82,2912982816387.
Der Umfang dieser Ellipse ist der Grenzwert einer Reihe U = 2 * Pi * a * (1-
(1/2)^2 * e^2 - ((1 * 3) / (2 * 4))^2 * e^4 / 3 - ((1 * 3 * 5) / (2* 4 *
6))^2 * e^6 / 5 ...) U = 524,210359956671 (100 Reihenglieder).
Die Bogenlänge einer 1/4 Ellipse (von psi = 0° - 90°) = U / 4 =
131,052589989168.
Vielen Dank für jeden Tipp.
PS: Nachfolgend o.g. Formeln in Visual Basic
Jürgen HEYN
j.h...@gmx.de
Option Explicit
Private Sub Command1_Click()
Dim a As Double 'Große Halbachse
Dim b As Double 'Kleine Halbachse
Dim f As Double 'Abplattung
Dim e As Double 'numerische Exzentrizität
Dim Pi As Double 'Kreiszahl
Dim phi As Double 'geografische Breite
Dim psi As Double 'geozentrische Breite
Dim r As Double 'Radius der Ellipse bei phi/psi
a = 100
b = 65
Pi = 3.14159265358979
phi = 60
f = (a - b) / a
e = Sqr((a ^ 2 - b ^ 2) / a ^ 2)
'alternativ e = Sqr(2 * f - f ^ 2)
'alternativ e = Sqr(1 - (1 - f) ^ 2)
Dim U As Double 'Umfang der Ellipse
Dim n As Integer 'Anzahl der Reihenglieder
Dim tmp(100) As Double 'Variable (Zwischenspeicher für 100
Reihenglieder)
Dim x As Double 'Variable
Dim Odd As Double 'Ungerader Zähler
Dim Evn As Double 'Gerader Nenner
Odd = 1
Evn = 1
x = 1
For n = 1 To 100 Step 1
Odd = Odd * (2 * n - 1)
Evn = Evn * (2 * n)
tmp(n) = ((Odd / Evn) ^ 2) * ((e ^ (2 * n)) / (2 * n - 1))
x = x - tmp(n)
Next n
U = 2 * Pi * a * x
Dim SinPhi As Double 'Sinus phi
Dim SinPsi As Double 'Sinus psi
SinPhi = Sin(Deg2Rad(phi))
SinPsi = Sin(Deg2Rad(CalcPSI(phi, e)))
r = ((((1 - e ^ 2) * a) / (Sqr(1 - e ^ 2 * SinPhi ^ 2))) * SinPhi) /
SinPsi
End Sub
Public Function Deg2Rad(Deg As Double) As Double
Deg2Rad = Deg * (Pi / 180)
End Function
Public Function Rad2Deg(Rad As Double) As Double
Rad2Deg = Rad * (180 / Pi)
End Function
Hermann Kremer
Jürgen HEYN schrieb in Nachricht <9kp1sc$5jlam$1...@ID-94872.news.dfncis.de>...
>Guten Abend,
>
>ich suche eine Formel zur Berrechnung der Bogenlänge einer Ellipse zwischen
>zwei beliebigen Winkeln. (Benötigt werden 0° - 90°.)
Hallo Jürgen,
Du brauchst eine numerische Näherung für das unvollständige elliptische
Integral 2. Art E(k,phi). Dann gilt
Bogenlänge = a * ( E(e, phi2) - E(e, phi1) ) [ mit a, e wie unten ]
Sowas könnte z.B. in "Numerical Recipes in C / Fortran" sein. Das Buch gibt es
wahlweise mit Programmen in C, Fortran-77 und Fortran-90 und kann bei
http://www.ulib.org/webRoot/Books/Numerical_Recipes
kapitelweise als PDF-Dateien heruntergeladen werden.
Falls Du da nichts findest, dann kannst Du auch in den Software-Archiven
http://www.netlib.org
http://gams.nist.gov [ Problem Class C14 ]
nachschauen; das sind größtenteils Fortran- und C-Programme.
Viel Erfolg
Hermann
--
Jan C. Hoffmann
"Jürgen HEYN" <j.h...@gmx.de> schrieb im Newsbeitrag
news:9kp1sc$5jlam$1...@ID-94872.news.dfncis.de...
Sollte man numerisch lösen:
Nach Bronstein: Beispiel für a = 1.5, b = 1, k = 0.74
U = 4 * a * E(0.74, pi/2)
E(0.74, pi/2) =
Integral[0, pi/2] sqrt(1 - 0.74^2 * sin^2(x)) dx
[8.24b]
Eingabesyntax: sqrt(1 - pow(0.74, 2) * pow(sin(x),2))
http://ourworld.compuserve.com/homepages/MTEC/quad.htm
E(0.74, pi/2) = 1.3262755240171515
Ergebnis: U = 7.9576531441029000
Die Rechenungenauigkeit liegt bei < 10^-10
Bronstein gibt auch eine Näherungsformel [3.327c] an:
Ergebnis: U = 7.93
--
Gruss Jan C. Hoffmann (c/o MTEC)
http://ourworld.compuserve.com/homepages/MTEC/
Hermann Kremer
Jan C. Hoffmann <10055...@compuserve.com> schrieb in Nachricht
<9kped7$eei$1...@suaar1ac.prod.compuserve.com>...
>"Jürgen HEYN" <j.h...@gmx.de> schrieb im Newsbeitrag
>news:9kp1sc$5jlam$1...@ID-94872.news.dfncis.de...
>> Guten Abend,
>>
>> ich suche eine Formel zur Berrechnung der Bogenlänge einer Ellipse zwischen
>> zwei beliebigen Winkeln. (Benötigt werden 0° - 90°.)
>>
>
>Nach Bronstein: Beispiel für a = 1.5, b = 1, k = 0.74
>
>U = 4 * a * E(0.74, pi/2)
Vorsicht:
E(k, pi/2) ist das v o l l s t ä n d i g e elliptische Integral 2. Art, Jürgen
benötigt aber eine Näherung für das u n v o l l s t ä n d i g e E(k, phi), da seine
Winkel phi zwischen 0 und pi/2 liegen können ...
E(k,phi) = Integral{0...phi} sqrt(1 - k^2*sin^2(x) ) dx
Die Wurzel kann man wegen k^2 < 1 in eine Taylor-Reihe in sin^2(x) entwickeln und
diese gliedweise integrieren, und für phi = pi/2 ergibt das dann die Formel aus dem
Bronstein.
>
>E(0.74, pi/2) = Integral[0, pi/2] sqrt(1 - 0.74^2 * sin^2(x)) dx
>[8.24b]
>
>Eingabesyntax: sqrt(1 - pow(0.74, 2) * pow(sin(x),2))
>
>http://ourworld.compuserve.com/homepages/MTEC/quad.htm
>
>E(0.74, pi/2) = 1.3262755240171515
>Ergebnis: U = 7.9576531441029000
Ja, das ist der Umfang der gesamten Ellipse ...
>Die Rechenungenauigkeit liegt bei < 10^-10
>
>Bronstein gibt auch eine Näherungsformel [3.327c] an:
>Ergebnis: U = 7.93
Grüße
Hermann
--
Jan C. Hoffmann
"Hermann Kremer" <hermann...@online.de> schrieb im Newsbeitrag
news:9kpn5k$i55$1...@news.online.de...
>
> Jan C. Hoffmann <10055...@compuserve.com> schrieb in Nachricht
> <9kped7$eei$1...@suaar1ac.prod.compuserve.com>...
> >"Jürgen HEYN" <j.h...@gmx.de> schrieb im Newsbeitrag
> >news:9kp1sc$5jlam$1...@ID-94872.news.dfncis.de...
> >> Guten Abend,
> >>
> >> ich suche eine Formel zur Berrechnung der Bogenlänge einer Ellipse zwischen
> >> zwei beliebigen Winkeln. (Benötigt werden 0° - 90°.)
> >>
> >
> >Sollte man numerisch lösen:
> >Nach Bronstein: Beispiel für a = 1.5, b = 1, k = 0.74
> >
> >U = 4 * a * E(0.74, pi/2)
>
> Vorsicht:
> E(k, pi/2) ist das v o l l s t ä n d i g e elliptische Integral 2. Art,
Jürgen
> benötigt aber eine Näherung für das u n v o l l s t ä n d i g e E(k, phi),
da seine
> Winkel phi zwischen 0 und pi/2 liegen können ...
>
> E(k,phi) = Integral{0...phi} sqrt(1 - k^2*sin^2(x) ) dx
>
Hallo Hermann,
a = 1.5
b = 1.0
x1 = 0
x2 = pi/2
L_AB = Integral[x1, x2] sqrt(a^2 - (a^2 - b^2)*sin^2(x)) dx
Eingabesyntax:
sqrt(pow(1.5, 2) - (pow(1.5, 2) - pow(1, 2))*pow(sin(x),2))
http://ourworld.compuserve.com/homepages/MTEC/quad.htm
0 ... pi/2: L = 1.9831799 (U = 4*L = 7.932719795)
0 ... pi/4: L = 1.1157582
pi/4 ... pi/2: L = 0.8674217
[] Bronstein / Semendjajew u.a., Taschenbuch der Mathematik CD, Verlag Harri
Deutsch
Hermann Kremer
>"Hermann Kremer" <hermann...@online.de> schrieb im Newsbeitrag
>news:9kpn5k$i55$1...@news.online.de...
>> Jan C. Hoffmann <10055...@compuserve.com> schrieb in Nachricht
>> <9kped7$eei$1...@suaar1ac.prod.compuserve.com>...
>> >"Jürgen HEYN" <j.h...@gmx.de> schrieb im Newsbeitrag
>> >news:9kp1sc$5jlam$1...@ID-94872.news.dfncis.de...
>> >> Guten Abend,
>> >>
>> >> ich suche eine Formel zur Berrechnung der Bogenlänge einer Ellipse zwischen
>> >> zwei beliebigen Winkeln. (Benötigt werden 0° - 90°.)
>> >>
>> >
>> >Sollte man numerisch lösen:
>> >Nach Bronstein: Beispiel für a = 1.5, b = 1, k = 0.74
>> >
>> >U = 4 * a * E(0.74, pi/2)
>>
>> Vorsicht:
>> E(k, pi/2) ist das v o l l s t ä n d i g e elliptische Integral 2. Art,
>> Jürgen benötigt aber eine Näherung für das u n v o l l s t ä n d i g e
>> E(k, phi), da seine Winkel phi zwischen 0 und pi/2 liegen können ...
>>
>> E(k,phi) = Integral{0...phi} sqrt(1 - k^2*sin^2(x) ) dx
>>
Hallo Jan,
>a = 1.5, b = 1.0, x1 = 0, x2 = pi/2
>
>L_AB = Integral[x1, x2] sqrt(a^2 - (a^2 - b^2)*sin^2(x)) dx
>
>Eingabesyntax:
>sqrt(pow(1.5, 2) - (pow(1.5, 2) - pow(1, 2))*pow(sin(x),2))
>http://ourworld.compuserve.com/homepages/MTEC/quad.htm
>
>0 ... pi/2: L = 1.9831799 (U = 4*L = 7.932719795)
>0 ... pi/4: L = 1.1157582
>pi/4 ... pi/2: L = 0.8674217
>
>[] Bronstein / Semendjajew u.a., Taschenbuch der Mathematik CD, Verlag Harri
>Deutsch
Ja, natürlich - ich wollte Jürgen nur darauf hinweisen, daß es sowohl unvollständige
als auch vollständige elliptische Integrale (1., 2. und 3. Art) gibt ...
Grüße
Hermann
--
Jürgen HEYN
gibt man die Werte in das Rechner-Programm auf Eurer (Deiner?) Homepage ein,
erhält man das lorrekte Ergebnis.
Mich würde interessieren, wie die Rechenroutine aussieht.
Folgende Routine finktioniert, die Ergebnisse weichen aber von den Deinen
ab.
a = 1.5 'große Halbachse
b = 1 'kleine Halbachse
n = 10000 'Anzahl der Schleifendurchläufe
t0 = 0 'Startwinkel
t1 = Pi / 2 'Endwinkel (90° = 1/4 Umfang)
var = 0
For i = 0 To n - 1
var = var + Arc_F(a, b, t0 + (t1 - t0) * (i / n))
Next i
var = var * ((t1 - t0) / n)
Function Arc_F(a, b, x)
Arc_F = Sqr(a^2 * Sin(x)^2 + b^2 * Cos(x)^2)
End Function
'Ergebnis für var = 1,98314067875316
'Dein Ergebnis = 1.9831799
Jeder Vorschlag, der die Genauigkeit erhöht wird gerne akzeptiert.
Grüsse
Jürgen HEYN
j.h...@gmx.de
Hermann Kremer
-----Ursprüngliche Nachricht-----
Von: Jürgen HEYN <j.h...@gmx.de>
An: Hermann Kremer <hermann...@online.de>
Datum: Mittwoch, 8. August 2001 18:00
Betreff: Re: Bogenlänge einer Ellipse ?
Hallo Jürgen,
>Vielen Dank für Ihre Ausführungen zu meiner Frage "Bogenlänge einer Ellipse"
>in der Newsgroup de.sci.mathematik.
>Leider verstehe ich Ihre Antwort nicht.
Ich entnehme Deiner Antwort, daß Du wohl kein Mathematiker bist ;-))
OK, wenn ich das richtig verstanden habe, dann hast Du eine Ellipse mit großer
Halbachse a und kleiner Halbachse b in Normal-Lage im Koordinatensystem,
d.h. mir der großen Ellipsenachse auf der x-Koordinatenachse, der kleinen
Ellipsenachse auf der y-Koordinatenachse und dem Ellipsenmittelpunkt im
Ursprung des Koordinatensystems.
Du hast dann zwei Punkte auf dem Ellipsenumfang, sagen wir A und B, die durch
ihre Öffnungswinkel phiA und phiB (vom Ellipsenmittelpunkt aus gemessene
Winkel zwischen positiver x-Koordinatenachse und jeweiligem Punkt) beschrieben
sind, und Du suchst die Länge des Ellipsenbogens zwischen diesen beiden Punkten.
Stimmt das so ?
Falls ja, dann wird es wohl das einfachste sein, wenn ich Dir ein fertiges Programm
schicke; da ich aber nicht in Basic arbeite, wird das ein C-Programm sein. Wäre
das OK ?
Grüße
Hermann
--
>Ich habe das Internet durchsucht und einige Routinen gefunden. Alle
>errechnen Werte mit Hilfe einer Schleife, die x-mal durchlaufen wird.
>Vermutlich entspricht dies einer Taylor-Reihe, die Sie erwähnen. Wird mit
>Hilfe einer solchen Schleife (Reihe) die Präzision erhöht?
>Wie kann ich für sin^2(x) eine Taylor-Reihe entwickeln.
>Entspricht k dem Wert der numerischen Exzentrizität e ?
>Was ist x für ein Wert und wo befindet er sich in meinem Koordinatensystem ?
>Leider kann ich auch den Ausdruck "Bogenlänge = a * ( E(e, phi2) - E(e,
>phi1) ) [ mit a, e wie unten ]" nicht deuten.
>Was bedeutet "E" und was macht E mit e und phi_ ?
>
>Ich wäre Ihnen sehr verbunden, wenn Sie mir nochmals helfen würden.
>Vielen Dank für Ihre Mühe
>mit freundlichen Grüssen
>Jürgen Heyn
>j.h...@gmx.de
>
Jürgen HEYN
vielen Dank für Deine Antwort.
Du hast es richtig erkannt (ich kann es wohl kaum verheimlichen), ich bin
KEIN Mathematiker. Als Berufshubschrauberpilot habe ich mir als Hobby die
Kartogrphie ausgesucht. (Woher sollte ich denn wissen, daß das alles so
mathematisch ist?)
Die Annahmen meine Ellipse betreffend in Deinen Ausführungen treffen zu.
Es wäre wirklich super toll, wenn Du mir ein C Proramm schicken würdest.
(Jede Zeile (Kommentar) Dokumentation machen es mir leichter.)
Mit freundlichen Grüßen
Jürgen Heyn
j.h...@gmx.de
www.heliheyn.de
PS: Sorry wegen der persönlichen E-Mail, aber ich habe mich wegen der Fragen
etwas geschämt ;-(
Hermann Kremer
Jürgen HEYN schrieb in Nachricht <9kv091$6mhfc$1...@ID-94872.news.dfncis.de>...
>Hallo Hermann,
>
>vielen Dank für Deine Antwort.
>
>Du hast es richtig erkannt (ich kann es wohl kaum verheimlichen), ich bin
>KEIN Mathematiker. Als Berufshubschrauberpilot habe ich mir als Hobby die
... Ist doch äußerst interessant ...
>(Woher sollte ich denn wissen, daß das alles so mathematisch ist?)
Na ja, Mercator, Gauß, Krüger, Lambert, Mollweide, ... waren halt alle Mathematiker
...
>Die Annahmen meine Ellipse betreffend in Deinen Ausführungen treffen zu.
>Es wäre wirklich super toll, wenn Du mir ein C Proramm schicken würdest.
>(Jede Zeile (Kommentar) Dokumentation machen es mir leichter.)
OK, ist notiert. Ich hoffe, daß ich am Wochenende dazu komme - falls nicht ein
anderer
lieber Mensch schon ein fertiges Programm dafür hat - Stichwort:
Bogenlänge einer Ellipse in Normallage als Differenz zweier unvollständiger
elliptischer Integrale 2. Art.
Grüße
Hermann
--
>Mit freundlichen Grüßen
>Jürgen Heyn
>j.h...@gmx.de
>www.heliheyn.de
>
>PS: Sorry wegen der persönlichen E-Mail, aber ich habe mich wegen der Fragen
>etwas geschämt ;-(
Absolut kein Grund zum schämen - ich kann stattdessen keine Hubschrauber fliegen
;-((
Jan C. Hoffmann
> Hallo Jan,
>
> gibt man die Werte in das Rechner-Programm auf Eurer (Deiner?) Homepage ein,
> erhält man das lorrekte Ergebnis.
> Mich würde interessieren, wie die Rechenroutine aussieht.
Die Berechnungsmethode wurde von Gauß entwickelt.
Siehe hierzu
[] H.R. Schwarz, Numerische Mathematik, B.G. Teubner Stuttgart
> Folgende Routine finktioniert, die Ergebnisse weichen aber von den Deinen
> ab.
> a = 1.5 'große Halbachse
> b = 1 'kleine Halbachse
> n = 10000 'Anzahl der Schleifendurchläufe
> t0 = 0 'Startwinkel
> t1 = Pi / 2 'Endwinkel (90° = 1/4 Umfang)
>
> var = 0
> For i = 0 To n - 1
> var = var + Arc_F(a, b, t0 + (t1 - t0) * (i / n))
> Next i
> var = var * ((t1 - t0) / n)
>
> Function Arc_F(a, b, x)
> Arc_F = Sqr(a^2 * Sin(x)^2 + b^2 * Cos(x)^2)
> End Function
>
> 'Ergebnis für var = 1,98314067875316
> 'Dein Ergebnis = 1.9831799
Hallo Jürgen,
mit Deiner Routine
n = 10^6: L = 1,98317955596218
n = 10^7: L = 1,98317990939165 (ca. 1,5 Minuten)
n = 10^8: L = 1,98317994473474 (ca. 15 Minuten)
Gauß Quadratur: L = 1,98317994866132 (ca. 0,1 Sekunden)
n = 10^6 Simpson: L = 1,98317994866136 (ca. 10 Sekunden)
Falls die Gauß Quadratur Methode für Dich ausreicht, kannst Du die Webseite bei
Dir abspeichern und hast damit ein fertiges Programm.
GQ-Test für Kreisumfang
r = a = b = 1
Eingabesyntax:
f(x) = sqrt(1 * pow(sin(x), 2) + 1 * pow(cos(x),2))
L = Integral[0,pi/2] f(x) dx = 1.570796326794899 (GQ)
L = U/4 = d*pi/4 = pi/2 = 1.5707963267949000
Klaus Nagel
Hermann Kremer schrieb:
> Na ja, Mercator, Gauß, Krüger, Lambert, Mollweide, ... waren halt alle Mathematiker
> ...
>
> lieber Mensch schon ein fertiges Programm dafür hat - Stichwort:
> Bogenlänge einer Ellipse in Normallage als Differenz zweier unvollständiger
> elliptischer Integrale 2. Art.
Wenn es sich um kartographische Probleme handelt, dann ist die Exzentrizität gewöhnlich
klein, denn die Erde ist nur schwach abgeplattet. Dann ist es möglich, die Bogenlänge (
vom Äquator aus gezählt) als Funktion der geographischen Breite phi zu entwickeln:
Bogenlänge = a*[a0*phi + a2*sin(2*phi) + a4*sin(4*phi) + ... ]
Dabei hängen die Koeffizienten ai von der Exzentrizität ab. Die allgemeine Formel dazu
weiß ich nicht, doch für die Gauß-Krüger-Transformation habe ich die Werte in einem
alten C-Programm gefunden. Gauß-Krüger benutzt das Besselsche Ellipsoid mit den
Halbachsen
a 6377397.15500
b 6356078.96325
Die Bogenlänge (in Metern) - und das ist in Kilometern der sogenannte Hochwert der
Gauß-Krüger-Koordinaten - errechnet sich zu
hw = 111120.61962 * phi
- 15988.638530 * SIN ( 2. * phi)
+ 16.7299538 * SIN ( 4. * phi)
- 0.02178480 * SIN ( 6. * phi)
+ 0.000030766 * SIN ( 8. * phi)
- 0.0000000452 * SIN (10. * phi)
Dabei sind phi und das Argument von SIN in Grad angegeben. Die Länge des Erdquadranten
ergibt sich daraus zu 111120.61962m * 90 = 10000,856 km.
Das C-Programm zur Umrechnung geographischer Koordinaten in Gauß-Krüger mit einer
Genauigkeit im Zentimeterbereich schicke ich auf Wunsch gern zu.
Gruß,
Klaus Nagel
Hermann Kremer
Klaus Nagel schrieb in Nachricht <3B73AFF8...@T-Online.de>...
>Hermann Kremer schrieb:
>
>> Na ja, Mercator, Gauß, Krüger, Lambert, Mollweide, ... waren halt alle Mathematiker
>> ...
>>
>> lieber Mensch schon ein fertiges Programm dafür hat - Stichwort:
>> Bogenlänge einer Ellipse in Normallage als Differenz zweier unvollständiger
>> elliptischer Integrale 2. Art.
Hallo Klaus,
>Wenn es sich um kartographische Probleme handelt, dann ist die Exzentrizität
gewöhnlich
>klein, denn die Erde ist nur schwach abgeplattet. Dann ist es möglich, die Bogenlänge
(
>vom Äquator aus gezählt) als Funktion der geographischen Breite phi zu entwickeln:
>
>Bogenlänge = a*[a0*phi + a2*sin(2*phi) + a4*sin(4*phi) + ... ]
>
>Dabei hängen die Koeffizienten ai von der Exzentrizität ab. Die allgemeine Formel
dazu
>weiß ich nicht, doch für die Gauß-Krüger-Transformation habe ich die Werte in einem
>alten C-Programm gefunden. Gauß-Krüger benutzt das Besselsche Ellipsoid mit den
>Halbachsen
>
>a 6377397.15500
>b 6356078.96325
>
>Die Bogenlänge (in Metern) - und das ist in Kilometern der sogenannte Hochwert der
>Gauß-Krüger-Koordinaten - errechnet sich zu
>
> hw = 111120.61962 * phi
> - 15988.638530 * SIN ( 2. * phi)
> + 16.7299538 * SIN ( 4. * phi)
> - 0.02178480 * SIN ( 6. * phi)
> + 0.000030766 * SIN ( 8. * phi)
> - 0.0000000452 * SIN (10. * phi)
>
>Dabei sind phi und das Argument von SIN in Grad angegeben.
Auch die Sinus-Argumente in Grad ? Oder ist SIN(x) gleich sin(x) mit Grad->Rad
Umrechnung ?
> Die Länge des Erdquadranten
>ergibt sich daraus zu 111120.61962m * 90 = 10000,856 km.
>Das C-Programm zur Umrechnung geographischer Koordinaten in Gauß-Krüger mit einer
>Genauigkeit im Zentimeterbereich schicke ich auf Wunsch gern zu.
Bei dem von Jürgen in seinem 1. Posting angegebenen Beispiel a = 100, b = 65, e ~=
0.76 handelt es sich aber wohl nicht um ein Erdellipsoid ;-)
Falls die von Dir angegebene Approximation für Jürgens Problem passen sollte, umso
besser. Falls nicht, würde ich ihm ein C-Programm mit den elliptischen Integralen aus
der Cephes-Library basteln.
Grüße
Hermann
--
>
>Gruß,
>
>Klaus Nagel
>
Klaus Nagel
Hermann Kremer schrieb:
>
> >Die Bogenlänge (in Metern) - und das ist in Kilometern der sogenannte Hochwert der
> >Gauß-Krüger-Koordinaten - errechnet sich zu
> >
> > hw = 111120.61962 * phi
> > - 15988.638530 * SIN ( 2. * phi)
> > + 16.7299538 * SIN ( 4. * phi)
> > - 0.02178480 * SIN ( 6. * phi)
> > + 0.000030766 * SIN ( 8. * phi)
> > - 0.0000000452 * SIN (10. * phi)
> >
> >Dabei sind phi und das Argument von SIN in Grad angegeben.
>
> Auch die Sinus-Argumente in Grad ? Oder ist SIN(x) gleich sin(x) mit Grad->Rad
> Umrechnung ?
>
Für phi = 0° und phi = 90° müssen alle Korrekturterme verschwinden. Damit müßte es klar
sein.
SIN(x) ist bei mir definiert als
#define SIN(x) sin(0.01745329252 * x)
> > Die Länge des Erdquadranten
> >ergibt sich daraus zu 111120.61962m * 90 = 10000,856 km.
> Bei dem von Jürgen in seinem 1. Posting angegebenen Beispiel a = 100, b = 65, e ~=
> 0.76 handelt es sich aber wohl nicht um ein Erdellipsoid ;-)
Das hat mich auch stutzig gemacht, diese Exzentrizität und das Wort Abwicklung eines
Kegelstumpfes erinnert mehr an Klempnerarbeiten. Aber welcher Klempner spricht dabei schon
von geographischer und geozentrischer Breite? Ich nehme an, daß Jürgen zum Zeichnen solch
große Exzentrizitäten benutzt hat, weil man bei kleinen Werten kaum einen Unterschied zum
Kreis erkennen kann.
>
> Falls die von Dir angegebene Approximation für Jürgens Problem passen sollte, umso
> besser. Falls nicht, würde ich ihm ein C-Programm mit den elliptischen Integralen aus
> der Cephes-Library basteln.
Andere Exzentrizitäten erfordern andere Koeffizienten, die allgemeine Formel für a_i[e]
habe ich (noch) nicht abgeleitet, aber für e ~ 0.75 wird die Konvergenz schlecht sein.
Gruß,
Klaus Nagel
Jürgen HEYN
z.Zt versuche ich hinter eine sogenannte "Polyconic Projection" (Mehrkegel
Projektion) zu steigen.
Auf jeden geografischen Breitengrad phi (Punkt P) wird ein Kegel aufgesetzt
dessen Seitenlänge der Strecke 1/Tan(psi) entspricht. (psi ist die
geozentrische Breite Winkel zwischen x = 0 und y = 0 nach P.
Wickelt man den Kegel ab, so erhält man einen Kreisausschnitt (Sektor).
Ich lege jeden geografischen Punkt auf die Line die der Kegelboden
beschreibt.
Korrekt, ich versuche erste Ansätze in CorelDRAW zu zeichen. Eine Grafik
(ausschließlich Breitengrade) für eine Kugel habe ich bereits erstellt.
Der nächste Schritt wäre die Umrechnung der geografischen Koordinaten, um
alle Ländergrenzen einzeichen zu können.
Nur wegen des korrekten Rechenweges hätte ich gerne die Formel für die
Errechnung der Bogenlänge, da o.g. Kreisbogen jeweils den korrekten Abstand
haben soll.
Nochmals vielen Dank für Alle, die bereit sind mir zu helfen.
Mit freundlichen Grüssen
Jürgen Heyn
j.h...@gmx.de
www.heliheyn.de
Jürgen HEYN
vielen Dank für den Tipp.
Ich hätte erwartet, daß der Rechner ein Java Applet ist.
Erst jetzt habe ich bemerkt, daß es ein Script ist, welches in HTML
integriert wurde.
Mal sehen, ob ich Java entschlüsselt bekomme ...
Ansonsten gibt es sicherlich eine kompetente Newsgroup.
Nochmals Danke
Grüsse
Jürgen Heyn
j.h...@gmx.de
www.heliheyn.de
Hermann Kremer
Jürgen HEYN schrieb in Nachricht <9l14mn$6vft0$1...@ID-94872.news.dfncis.de>...
>Hallo Leute (Kollegen ;-),
Hallo Jürgen,
>z.Zt versuche ich hinter eine sogenannte "Polyconic Projection" (Mehrkegel
>Projektion) zu steigen.
>
>Auf jeden geografischen Breitengrad phi (Punkt P) wird ein Kegel aufgesetzt
>dessen Seitenlänge der Strecke 1/Tan(psi) entspricht. (psi ist die
>geozentrische Breite Winkel zwischen x = 0 und y = 0 nach P.
>Wickelt man den Kegel ab, so erhält man einen Kreisausschnitt (Sektor).
>Ich lege jeden geografischen Punkt auf die Line die der Kegelboden
>beschreibt.
>
>Korrekt, ich versuche erste Ansätze in CorelDRAW zu zeichen. Eine Grafik
>(ausschließlich Breitengrade) für eine Kugel habe ich bereits erstellt.
>Der nächste Schritt wäre die Umrechnung der geografischen Koordinaten, um
>alle Ländergrenzen einzeichen zu können.
>
>Nur wegen des korrekten Rechenweges hätte ich gerne die Formel für die
>Errechnung der Bogenlänge, da o.g. Kreisbogen jeweils den korrekten Abstand
>haben soll.
OK, ganz kurz dazu: Die Bogenlänge einer beliebigen Ellipse läßt sich leider nicht
mit den elementaren Funktionen Sinus, Cosinus usw. beschreiben, sondern man
benötigt dazu eine spezielle Funktion E(k, phi) mit den zwei Argumenten k und
phi, die den schönen Namen "unvollständiges elliptisches Integral 2. Art" trägt
(es gibt auch ein solches "... 1. Art", das wird z.B. bei der Berechnung von Pendeln
mit großer Auslenkung gebraucht, und eins "... 3. Art", das wird auch irgendwo
gebraucht :-). Neben diesen "unvollständigen ..." gibt es dann auch noch die
"vollständigen ... 1., 2. und 3. Art", die werden irgendwo auch gebraucht.
Wenn Du jetzt auf dem Umfang einer Ellipse mit Halbachsen a, b, a >= b und damit
der numerischen Exzentrizität e = sqrt(a^2 - b^2) / a zwei Punkte mit den vom
Äquator (große Achse) aus gemessenen geozentrischen Winkeln psi_1 und psi_2 hast,
dann hat das Ellipsenstück zwischen diesen beiden Punkten die Bogenlänge:
B = a*[ E(e, pi/2 - arctan( (a/b)*tan(psi_1) ) -
- E(e, pi/2 - arctan( (a/b)*tan(psi_2) ) ]
also (Große Halbachse) mal (Differenz zweier elliptischer Integrale). Wenn allerdings
die Exzentrizität e sehr klein ist, wie beim Erdellipsoid, dann kann man diese
elliptischen Integrale durch Reihenentwicklungen approximieren und erhält damit die
von Klaus angegebe Näherungsformel für die Bogenlänge. Bei großer Exzentrizität (e >
0.1) braucht man jedoch sehr viel mehr Terme in der Näherungsformel, und dann ist es
meistens genauer, die elliptischen Integrale mittels spezieller numerischer
Algorithmen direkt auszurechnen. Solche Algorithmen habe ich, und ich würde Dir ein
C-Programm damit basteln, falls Du genaue Werte für stark exzentrische Ellipsen
wirklich brauchst - dieses Programm wird dann aber leider ziemlich kompliziert.
Falls Du aber nur mit dem Erdellipsoid rechnen willst, dann reichen die Formeln von
Klaus m.E. völlig aus.
PS: Deine Homepage ist übrigens hervorragend gemacht (Die kleine Bildergeschichte ist
schwärzester Humor, und die Peters-Projektion ist wirklich spannend), nur lädt sie
halt fürchterlich langsam wegen der vielen schönen Bilder ...
Grüße
Hermann
--
Jürgen HEYN
vielen Dank für Eure Angebote, mir entsprechende C Routinen zur Verfügung zu
stellen.
Ich denke, daß ich zuerst versuche das Java Script von der Internetseite von
Jan C. Hoffmann umzusetzten.
Die Berechnung der Bogenlänge einer beliebigen Ellipse mit Hilfe der Gauß
Quadratur erscheint mir derzeit als (relativ) einfache Lösung mit
ausreichender Genauigkeit. (siehe letztes Posting von Jan)
Falls ich doch auf Euer Angebot zurückkommen möchte, würde ich mich in
dieser Newsgroup melden.
Nochmals vielen Dank und ein schönes Wochenende
Jan C. Hoffmann
Hallo Hermann,
Zusammenfassung:
x = a*cos(t)
y = b*sin(t)
L = Integral dl
dl^2 = dx^2 + dy^2
dx = -a*sin(t)
dy = b*cos(t)
L = integral[x1, x2] sqrt((a*sin(t))^2 + (b*cos(t))^2) dt
a = 6377397.15500 m
b = 6356078.96325 m
x1 = 0
x2 = pi/2
Eingabesyntax:
sqrt(pow(6377397.15500*sin(x), 2) + pow(6356078.96325*cos(x), 2))
http://ourworld.compuserve.com/homepages/MTEC/quad.htm
L = 10000855.764771388 m (Klaus Nagel: 10000855,7658 m)
Übereinstimmung auf einen cm !!!
Beispiel 1° bei pi/4
x1 = pi/4 = 0.785398163397448
x2 = pi/4 + 1*pi/180 = 0.802851455917392
L = 111123.94402206751 m (Erdumfang = 40004619.8479441 m)
Genauigkeit ca. 12 Stellen (µm-Bereich!)
L = 111123.944022 m
Hermann Kremer
>"Hermann Kremer" <hermann...@online.de> schrieb im Newsbeitrag
>news:9l1bk1$2hu$1...@news.online.de...
>> Jürgen HEYN schrieb in Nachricht <9l14mn$6vft0$1...@ID-94872.news.dfncis.de>...
Hallo Jan,
>Zusammenfassung:
>
>x = a*cos(t)
>y = b*sin(t)
>
>L = Integral dl
>
>dl^2 = dx^2 + dy^2
>
>dx = -a*sin(t)
>dy = b*cos(t)
>
>L = integral[x1, x2] sqrt((a*sin(t))^2 + (b*cos(t))^2) dt
Ja, OK, das stimmt schon - nur ist t nicht der geozentrische Ellipsenwinkel,
sondern ein Hilfswinkel, der bei der bekannten Konstruktion mittels affiner
Kreisabbildung auftritt:
Konzentrische Kreise mit Radien a und b,
Schnittpunkte S_b, S_a der beiden Kreise mit einer Geraden durch den
Ursprung und Steigungswinkel t,
Parallele zur x-Achse durch S_b, zur y-Achse durch S_a. Der Schnittpunkt
der beiden Parallelen ist der Ellipsenpunkt mit dem geozentrischen
Winkel phi,
--> tan(phi) = (b/a)*tan(t) --> t = arctan(a * tan(phi) / b) .
Für phi = 0 und phi = pi/2 ist natürlich t = phi, und bei a/b ~= 1 ist
ebenfalls
t ~= phi, aber nicht bei großer Exzentrizität. Dafür mußt Du die Integralgrenzen
noch
umrechnen:
L(phi1, phi2) = Integral{arctan(a*tan(phi1)/b) ... arctan(a*tan(phi2)/b)} ... dt
>a = 6377397.15500 m
>b = 6356078.96325 m
>
>x1 = 0
>x2 = pi/2
>
>Eingabesyntax:
>sqrt(pow(6377397.15500*sin(x), 2) + pow(6356078.96325*cos(x), 2))
>http://ourworld.compuserve.com/homepages/MTEC/quad.htm
>
>L = 10000855.764771388 m (Klaus Nagel: 10000855,7658 m)
>
>Übereinstimmung auf einen cm !!!
>
>Beispiel 1° bei pi/4
>
>x1 = pi/4 = 0.785398163397448
>x2 = pi/4 + 1*pi/180 = 0.802851455917392
>
>L = 111123.94402206751 m (Erdumfang = 40004619.8479441 m)
>
>Genauigkeit ca. 12 Stellen (µm-Bereich!)
>
>L = 111123.944022 m
Ja, klar, beim Erdellipsoid ist a ~= b ...
Grüße
Hermann
--
Jan C. Hoffmann
Hallo Hermann,
Jürgen würde sein Ziel bei 111116 m um ~ 7 m Meter verfehlen.
Neuer Ansatz für eine Punktlandung [1, Abb. 8.13 b]:
dl^2 = rho^2 * d_phi^2 + d_rho^2
L = Integral[phi1, phi2] sqr(rho^2 + (d_rho/d_phi)^2) * d_phi
x^2/a^2 + y^2/b^2 = 1
rho^2 = x^2 + y^2
tan^2(phi) = y^2/x^2
Numerisch differenziert (d_rho/d_phi) und integriert über GQ:
phi1 = pi/4 = 0.785398163397448
phi2 = pi/4 + 1*pi/180 = 0.802851455917392
L = 111116,828151609 m für delta_phi = 1°
Erdumfang über die Pole = 40003310,9337449 m
Fehlerabschätzung: Auf ~ 12 Stellen genau
L = 111116,828151 m
[1] Bronstein / Semendjajew u.a., Taschenbuch der Mathematik CD, Verlag Harri
Deutsch
--
Gruss Jan C. Hoffmann (c/o MTEC)
http://ourworld.compuserve.com/homepages/MTEC/
Die zu integrierende Funktion für Jürgen:
Function Funktion(ByVal phi As Double) As Double
Dim a, b, xQ, f1, f2, fStr, rho, d_phi As Double
a = 6377397.155
b = 6356078.96325
xQ = 1 / (1 / a ^ 2 + (Tan(phi)) ^ 2 / b ^ 2) 'xQ = x^2
rho = (xQ + b ^ 2 * (1 - xQ / a ^ 2)) ^ (1 / 2)
d_phi = 10 ^ -30
f1 = rho - d_phi / 2
f2 = rho + d_phi / 2
fStr = (f2 - f1) / d_phi
Funktion = Sqr(rho ^ 2 + fStr ^ 2)
End Function
Klaus Nagel
http://home.wtal.de/aw/excel/source/excel_koordinaten.txt
steht die Reihenentwicklung der Bogenlänge nach der geographischen Breite phi
(Bogenmaß!) in Abhängigkeit von den Halbachsen a und b angegeben:
N = (a - b) / (a + b)
alpha = (a + b) / 2 * (1 + 1 / 4 * N ^ 2 + 1 / 64 * N ^ 4)
beta = -3 / 2 * N + 9 / 16 * N ^ 3 - 3 / 32 * N ^ 5
gamma = 15 / 16 * N ^ 2 - 15 / 32 * N ^ 4
delta = -35 / 48 * N ^ 3 + 105 / 256 * N ^ 5
epsilon = 315 / 512 * N ^ 4
mb = alpha * (phi + beta * sin(2 * phi) + gamma * sin(4 * phi) + delta *
sin(6 * phi) + epsilon * sin(8 * phi))
Für das Besselsche Ellipsoid stimmen die Ergebnisse mit der von mir angegebenen
Formel bis auf eine maximale Differenz von 1 mm bei phi=90° überein. Die Formel
stammte, wenn ich mich recht erinnere, aus einem Buch von Krüger.
Erdquadrant : 10000855,765 m
Erdumfang über die Pole: 40003423,059 m
Meridianlänge zwischen 45° und 46° : 111129,192 m
Die Werte stimmen gut überein mit denen meines GPS-Gerätes.
"Jan C. Hoffmann" schrieb:
> L = 111116,828151609 m für delta_phi = 1°
> Erdumfang über die Pole = 40003310,9337449 m
>
> Fehlerabschätzung: Auf ~ 12 Stellen genau
>
> L = 111116,828151 m
Woher kommen die Abweichungen, beim Erdumfang immerhin 112 m?
Gruß,
Klaus Nagel
Jan C. Hoffmann
"Klaus Nagel" <nagel...@T-Online.de> schrieb im Newsbeitrag
news:3B764E4A...@T-Online.de...
Hallo Klaus,
Auch wenn man i.A. unter Numerik eine Approximation vesteht, ist es doch so, daß
die Genauigkeit beliebig hoch getrieben werden kann.
> > Fehlerabschätzung: Auf ~ 12 Stellen genau
> >
> > L = 111116,828151 m
>
> Woher kommen die Abweichungen, beim Erdumfang immerhin 112 m?
In bezug auf meine Aussage, daß auf 12 Stellen gerechnet wurde, eine gute Frage.
Kreis gewählt für Algorithmen-Test:
a = b = r = 10000000 m
Sub IntegrationsTest()
Dim pi As Double
pi = Application.WorksheetFunction.pi
Cells(5, 5).Value = 4 * GaussIntegration(0, pi / 2)
Cells(6, 5).Value = 4 * pi / 2 * 10000000 'Referenz
Cells(7, 5).Value = 360 * GaussIntegration(pi / 4, pi / 4 + pi / 180)
End Sub
Ergebnisse:
Cells(5, 5).Value = 62831853,0717958
Cells(6, 5).Value = 62831853,0717959
Cells(7, 5).Value = 62831853,0717959
Mathematische Geauigkeit im µm-Bereich !!!
Die Berechnungs-Methode:
'Funktion für GQ
Function Funktion(ByVal phi As Double) As Double
Dim f1, f2, fStr, d_phi As Double
d_phi = 10 ^ -30
f1 = rho(phi - d_phi / 2)
f2 = rho(phi + d_phi / 2)
fStr = (f2 - f1) / d_phi
Funktion = Sqr(rho(phi) ^ 2 + fStr ^ 2)
End Function
Function rho(phi As Double) As Double
Dim a, b, xQ As Double
a = 10000000 '6377397.15508
b = 10000000 '6356078.9629
xQ = 1 / (1 / a ^ 2 + (Tan(phi)) ^ 2 / b ^ 2) 'xQ = x^2
rho = (xQ + b ^ 2 * (1 - xQ / a ^ 2)) ^ (1 / 2)
End Function
Wenn obige 2 Funktionen (2. Methode) ok sind, müssen meine Ergebnisse stimmen.
Jan C. Hoffmann
"Jan C. Hoffmann" <10055...@compuserve.com> schrieb im Newsbeitrag
news:9l5vgr$fha$1...@suaar1aa.prod.compuserve.com...
Ein weiterer Test:
Ellipse
Sub IntegrationsTest()
Dim pi As Double
Cells.Select
Selection.ClearContents
Range("A1").Select
pi = Application.WorksheetFunction.pi
Cells(6, 5).Value = GaussIntegration(0, pi / 4)
Cells(7, 5).Value = GaussIntegration(pi / 4, pi / 2)
Cells(8, 5).Value = GaussIntegration(0, pi / 2)
End Sub
Function Funktion(ByVal phi As Double) As Double
Dim f1, f2, fStr, d_phi As Double
d_phi = 10 ^ -60
f1 = rho(phi - d_phi / 2)
f2 = rho(phi + d_phi / 2)
fStr = (f2 - f1) / d_phi
Funktion = Sqr(rho(phi) ^ 2 + fStr ^ 2)
End Function
Function rho(phi As Double) As Double
Dim a, b, xQ As Double
a = 15000000 '6377397.15508
b = 10000000 '6356078.9629
xQ = 1 / (1 / a ^ 2 + (Tan(phi)) ^ 2 / b ^ 2) 'xQ = x^2
rho = (xQ + b ^ 2 * (1 - xQ / a ^ 2)) ^ (1 / 2)
End Function
Ergebnis:
10729678,7679426 = L(0...pi/4)
8312735,4015074 = L(p/4...pi/2)
19042414,1694500 = L(0...p/2)
Klaus Nagel
"Jan C. Hoffmann" schrieb:
> > Die Berechnungs-Methode:
> >
> > 'Funktion für GQ
> >
> > Function Funktion(ByVal phi As Double) As Double
> > Dim f1, f2, fStr, d_phi As Double
> > d_phi = 10 ^ -30
> > f1 = rho(phi - d_phi / 2)
> > f2 = rho(phi + d_phi / 2)
> > fStr = (f2 - f1) / d_phi
Die numerische Differentiation (f(x+delta/2)-f(x-delta/2))/ delta mit delta = 10
^-30 erscheint mir höchst fragwürdig.
> > Funktion = Sqr(rho(phi) ^ 2 + fStr ^ 2)
> > End Function
> >
> > Function rho(phi As Double) As Double
> > Dim a, b, xQ As Double
> > a = 10000000 '6377397.15508
> > b = 10000000 '6356078.9629
> > xQ = 1 / (1 / a ^ 2 + (Tan(phi)) ^ 2 / b ^ 2) 'xQ = x^2
> > rho = (xQ + b ^ 2 * (1 - xQ / a ^ 2)) ^ (1 / 2)
Das ist richtig, aber sehr kompliziert gemacht, rho ergibt sich nach Pythagoras aus
x=a*cos(phi) und y = b*sin(phi). Der Winkel phi ist nicht die geographische Breite,
sondern der Winkel zwischen x-Achse und dem um den Faktor a/b in y-Richtung
gestreckten Ellipsenpunkt vom Mittelpunkt aus gesehen.
> Ergebnis:
>
> 10729678,7679426 = L(0...pi/4)
> 8312735,4015074 = L(p/4...pi/2)
> 19042414,1694500 = L(0...p/2)
Die Ergebnisse können nicht stimmen, denn der Meridianbogen vom Äquator bis zur
geographischen Breite 45° ist wegen der Abplattung kürzer als der Rest bis zum Pol.
Das liegt wohl daran, daß Du Hermanns Hinweis noch nicht beachtet hast.
Aber auch für den Quadranten (L(0...pi/2) weicht das Ergebnis erheblich von meinem
ab, ich erhalte 19831799,4866.
Ich bin inzwischen sicher, daß meine Ergebnisse stimmen, denn ich erhalte auf fünf
verschiedene Weisen Werte, die um höchstens 1 mm abweichen, sowohl beim
Erdquadranten, als auch für den Meridianbogen zwischen den Breiten 45° und 46°.
Die Methoden die ich angewandt habe sind:
1. Direkte Berechnung der Funktion a*(EllipticE(1,e)-EllipticE(cos(t*rad),e)) mit
Maple.
2. Gauß-Integration mit n=6 Stützstellen
3. Numerische Integration nach der Keplerschen Faßregel (Simpson)
4. Reihenentwicklung für das Besselsche Ellipsoid ( siehe Beitrag vom 10.8. 17:48)
5. Reihenentwicklung für beliebige Exzentrizität (siehe Beitrag von heute 11:37)
Die Länge des Quadranten für a=1.5 und b = 1.0 für die verschiedenen Verfahren:
1. 1.9831799487
2. 1.9831799842
3. 1.9832604775 ( 4 Intervalle)
1.9831799487 (12 Intervalle)
1.9831799487 (96 Intervalle)
4. Nicht anzuwenden
5. 1.9831794500
Es hat mich überrascht, daß auch die Reihenentwicklung bei der hohen Exzentrizität
von 0.745 noch solch gute Werte liefert.
Gruß,
Klaus Nagel
Jan C. Hoffmann
"Klaus Nagel" <nagel...@T-Online.de> schrieb im Newsbeitrag
Hallo Klaus,
genau diese Ergebnisse hatte ich auch schon am 10.08.2001.
<Zitat vom 10.08.2001>
n = 10^6: L = 1,98317955596218
n = 10^7: L = 1,98317990939165 (ca. 1,5 Minuten)
n = 10^8: L = 1,98317994473474 (ca. 15 Minuten)
Gauß Quadratur: L = 1,98317994866132 (ca. 0,1 Sekunden)
n = 10^6 Simpson: L = 1,98317994866136 (ca. 10 Sekunden)
</Zitat vom 10.ß8.2001>
Basis dafür war:
L = Integral[0,pi/2] Sqr(a^2 * Sin(x)^2 + b^2 * Cos(x)^2) dx [1]
Ist das richtig?
Die Integrationen sind nicht das Problem, sondern der Ansatz. Da wählte ich im
2. Anlauf den Bronstein, Seite 443, 8.60c, Abb 8.13b. Das weicht von [1] ab.
Jetzt brauchen wir einen kompetenten Schiedsrichter, der uns/mir erklärt, daß
der Ansatz im Bronstein die falsche Methode und [1] besser ist.
TIA
Klaus Nagel
Hallo Jan,
"Jan C. Hoffmann" schrieb:
> Basis dafür war:
> L = Integral[0,pi/2] Sqr(a^2 * Sin(x)^2 + b^2 * Cos(x)^2) dx [1]
>
> Ist das richtig?
Ja
> Die Integrationen sind nicht das Problem, sondern der Ansatz. Da wählte ich im
> 2. Anlauf den Bronstein, Seite 443, 8.60c, Abb 8.13b. Das weicht von [1] ab.
>
> Jetzt brauchen wir einen kompetenten Schiedsrichter, der uns/mir erklärt, daß
> der Ansatz im Bronstein die falsche Methode und [1] besser ist.
Den Bronstein besitze ich nicht, habe aber die Integration von Wurzel(rho^2 +
delta_rho^2) ausgeführt und komme auf die gleichen Werte für die Länge des
Quadranten (1.9831799). Man muß dabei wieder darauf achten, daß der Hilfswinkel
nicht die geographische Breite ist (Hermanns Hinweis). Wenn phi der Winkel zwischen
x-Achse und Ellipsenpunkt ist, vom Mittelpunkt aus gesehen, dann gilt rho =
b/Wurzel(1-(e*cos(phi))^2), auch die von Dir benutzte Formel läßt sich so
vereinfachen. Die Ableitung nach phi ergibt einen recht komplizierten Ausdruck,
deswegen hast Du vermutlich die etwas fragwürdige numerische Differentiation
gewählt. Damit ließen sich die Abweichungen bei der Integration erklären.
Etwas einfacher wird es, wenn man rho vom Brennpunkt aus mißt. Dann ergibt sich rho
= p/(1+e*cos(phi)), mit p = a*(1-e^2); auch so komme ich zum gleichen Ergebnis.
Gruß,
Klaus Nagel
Jürgen HEYN
einstweilen vielen Dank für die Ausführungen.
In der Funktion für die GaussIntegration werden 32 Array Variablen x und w
definiert:
x[ 1] = 0.024350292663424
x[ 2] = 0.072993121787799
x[ 3] = 0.12146281929612
x[ 4] = 0.16964442042399
x[ 5] = 0.21742364374001
...
x[32] = 0.99930504173577
w[ 1] = 0.04869095700914
w[ 2] = 0.048575467441503
w[ 3] = 0.048344762234803
w[ 4] = 0.047999388596458
w[ 5] = 0.04754016571483
...
w[32] = 1.7832807216964E-03
Ich nehme mal an, daß diese Variablen irgendeine Abhängigkeit von der
mathematischen Konstanten e = 2.7182818284590452354 sind.
Es wäre sehr interessant zu wissen, wie diese Variablen berrechnet werden.
Danke
Jürgen Heyn
j.h...@gmx.de
www.heliheyn.de
Klaus Nagel
"Jürgen HEYN" schrieb:
Bei numerischer Integration approximiert man die zu integrierende Funktion f(x)
durch ein Polynom und berechnet das Integral dieses Polynoms. Der berechnete
Wert hängt nur ab von den Werten von f(x) an den Stützstellen. Der einfachste
Fall ist die Annäherung durch eine Konstante, nämlich dem Wert von f(x) in der
Mitte des Intervalls.
Beispiel: f(x)=x^2 im Intervall [0,1] würde den ziemlich schlechten Wert f(1/2)
= 1/4 liefern.
Die Keplersche Faßregel benutzt 3 Punkte, nämlich die Intervallenden (x1 und
x3) und die Intervallmitte (x2). Dadurch kann man eine Parabel legen und wenn
man die integriert erhält man (f(x0) + 4*f(x1) + f(x2))/6. Für das Beispiel
erhält man (0 + 4*1/4 + 1)/6 = 1/3, also den richtigen Wert. Durch mehr
Stützstellen und Polynome höherer Ordnung erhält man bessere
Integrationsformeln, die auch für kompliziertere Funktionen gute Werte liefern.
Bei der Gaußschen Integration wählt man die Abszissen der Stützstellen nicht
äquidistant sondern so, daß das Integral für Polynome bis zu einer möglichst
hohen Ordnung verschwindet. Ich erkläre es am Beispiel mit 6 Stützstellen (,
das ich auch für den Ellipsenbogen benutzt habe). Man betrachtet das
Integrationsintervall von [-1,1], die Lage der Stützstellen in einem Intervall
[a,b] ergibt sich daraus durch eine lineare Transformation. 6 Stützstellen
bestimmen eindeutig ein Polynom 5ten Grades; dieses würde also richtig
berechnet, ganz gleich, wo die Stützstellen liegen. Die Integrale x^n für
ungerades n sind Null. Wende ich die Methode auf ein Polynom P6 6ten Grades an,
so erhalte ich eine Näherungspolynom P5 vom Grad 5. Die Differenz P6-P5 ist ein
Polynom Q6 mit Nullstellen an den Stützstellen, es gilt also Q6 = a*x * R5(x)
mit R5 = (x-x1)*(x-x2)*..*(x-x6). Wählt man die Stützstellen so, daß das
Integral von x*R5 verschwindet, dann liefert die Integralformel auch für
Polynome 6ten Grades korrekte Werte. Die Stützstellen liegen sicher
symmetrisch, man kann daher drei Werte frei wählen. Damit lassen sich dann auch
die Integrale x^3*R5 und x^5*R5 zu Null machen. Weil x^11 sowieso verschwindet,
erhält man also mit 6 Stützpunkten eine Formel, die für alle Polynome bis zum
11ten Grad den exakten Wert liefert.
Dein Beispiel oben gilt für 64 Stützstellen, wegen der Symmetrie ist nur die
Hälfte angegeben.
Gruß,
Klaus Nagel
Jan C. Hoffmann
"Klaus Nagel" <nagel...@T-Online.de> schrieb im Newsbeitrag
Hallo Klaus,
Private Sub IntegrationsTest()
Dim pi As Double
Cells.Select
Selection.ClearContents
pi = Application.WorksheetFunction.pi
Cells(6, 5).Value = GaussIntegration(0, pi / 4)
Cells(7, 5).Value = GaussIntegration(pi / 4, pi / 2)
Cells(8, 5).Value = GaussIntegration(0, pi / 2)
Cells(9, 5).Value = 1.98317994866132
Range("A1").Select
End Sub
Ergebnisse:
1,12825922320136
0,85492073925379
1,98318000977101
1,98317994866132
Basis:
Private Function Funktion(ByVal phi As Double) As Double
Dim f1, f2, fStr, d_phi As Double
d_phi = 10 ^ -10
f1 = r(phi - d_phi / 2)
f2 = r(phi + d_phi / 2)
fStr = (f2 - f1) / d_phi
Funktion = Sqr(r(phi) ^ 2 + fStr ^ 2)
End Function
Private Function r(phi As Double) As Double
Dim a, b, eps As Double
a = 1.5 '6377397.15508
b = 1 '6356078.9629
eps = Sqr(a ^ 2 - b ^ 2) / a
r = Sqr(b ^ 2 / (1 - (eps * Cos(phi)) ^ 2))
End Function
Erdumfang = 40003423,0529789 m
Dein Erg. = 40003423,059 m
Das reicht für eine Punktlandung. Ich hab mich mit der Numerik wieder versöhnt.
Danke für die Hinweise.
Klaus Nagel
Bei der Beschreibung habe ich einen Fehler gemacht. 6 Funktionswerte legen zwar ein
Polynom 5ten Grades fest, aber 6 Nullstellen ein Polynom 6ten Grades bis auf einen
konstanten Faktor. Der spielt aber hier keine Rolle, da das Integral zu Null
gemacht werden soll. Ich wiederhole den jetzt hoffentlich korrekten Absatz:
> Bei der Gaußschen Integration wählt man die Abszissen der Stützstellen nicht
> äquidistant sondern so, daß das Integral für Polynome bis zu einer möglichst
> hohen Ordnung verschwindet. Ich erkläre es am Beispiel mit 6 Stützstellen (,
> das ich auch für den Ellipsenbogen benutzt habe). Man betrachtet das
> Integrationsintervall von [-1,1], die Lage der Stützstellen in einem Intervall
> [a,b] ergibt sich daraus durch eine lineare Transformation. 6 Stützstellen
> bestimmen eindeutig ein Polynom 5ten Grades; dieses würde also richtig
> berechnet, ganz gleich, wo die Stützstellen liegen. Die Integrale x^n für
> ungerades n sind Null. Wende ich die Methode auf ein Polynom P6 6ten Grades an,
> so erhalte ich eine Näherungspolynom P5 vom Grad 5. Die Differenz P6-P5 ist ein
> Polynom Q6 mit Nullstellen an den Stützstellen, es gilt also Q6 = a* R6(x)
> mit R6 = (x-x1)*(x-x2)*..*(x-x6). Wählt man die Stützstellen so, daß das
> Integral von R6 verschwindet, dann liefert die Integralformel auch für
> Polynome 6ten Grades korrekte Werte. Die Stützstellen liegen sicher
> symmetrisch, man kann daher drei Werte frei wählen. Damit lassen sich dann auch
> die Integrale x^2*R6 und x^4*R6 zu Null machen. Weil x^11 sowieso verschwindet,
> erhält man also mit 6 Stützpunkten eine Formel, die für alle Polynome bis zum
> 11ten Grad den exakten Wert liefert.
>
Nennt man die Stützstellen +u, -u, +v, -v, +w und -w, so ergeben sich deren Werte
aus den Bedingungen
Integral(x=-1 bis x=1 von (x^2-u^2)*(x^2-v^2)*(x^2-w^2)) = 0
Integral(x=-1 bis x=1 von (x^2-u^2)*(x^2-v^2)*(x^2-w^2) * x^2) = 0
Integral(x=-1 bis x=1 von (x^2-u^2)*(x^2-v^2)*(x^2-w^2) * x^4) = 0
Man erhält:
u = .238619
v = .661209
w= .932470
Das sind Nullstellen von Polynomen, sie haben nichts mit der Zahl e zu tun.
Die Gewichte erhält man aus der Forderung, daß die Formel für x^6, x^8 und x^10 die
bekannten Integralwerte liefern muß.
Gruß,
Klaus Nagel
Hermann Kremer
<9l6nkq$1vu$1...@suaar1aa.prod.compuserve.com>...
>"Klaus Nagel" <nagel...@T-Online.de> schrieb im Newsbeitrag
>news:3B76CB99...@T-Online.de...
Hallo Jan,
>Die Integrationen sind nicht das Problem, sondern der Ansatz. Da wählte ich im
>2. Anlauf den Bronstein, Seite 443, 8.60c, Abb 8.13b. Das weicht von [1] ab.
>
>Jetzt brauchen wir einen kompetenten Schiedsrichter, der uns/mir erklärt, daß
>der Ansatz im Bronstein die falsche Methode und [1] besser ist.
Für die Ellipse gibt es mehrere richtige Formeln, die nur unterschiedlich definiert
sind:
1) Parameterform mit Hilfswinkel t
------------------------------------------------
x = a*cos(t)
y = b*sin (t)
2) Parameterform mit geozentrischem Winkel phi
-------------------------------------------------------------------
x = a*b*cos(phi)/sqrt( a^2*sin^2(phi) + b^2*cos^2(phi) )
y = a*b*sin (phi)/sqrt( a^2*sin^2(phi) + b^2*cos^2(phi) )
3) Polarform bezogen auf Brennpunkt, Brennpunktwinkel chi
---------------------------------------------------------------------------------
rho = p/( 1 + e*cos(chi) )
4) Polarform bezogen auf Ellipsenzentrum, geozentr. Winkel phi
-------------------------------------------------------------------------------------
r = a*b/sqrt( a^2*sin^2(phi) + b^2*cos^2(phi) )
Klaus hat in seinen Formeln wohl überall den geozentrischen Winkel phi.
Noch was: Wenn Du die Bogenlänge für eine __Viertelellipse__ ausrechnest,
dann liefern die Integrale mit den Def. nach 1) und 2) wegen
t = 0 <--> phi = 0 , t = pi/2 <--> phi = pi/2
natürlich den gleichen Wert. Die größten Abweichungen erhältst Du für einen
kurzen Ellipsenbogen bei phi = pi/4 +- w...
Grüße
Hermann
--
>
Klaus Nagel
Hermann Kremer schrieb:
>
>
> Klaus hat in seinen Formeln wohl überall den geozentrischen Winkel phi.
>
Die Reihenentwicklungen arbeiten direkt mit der geographischen Breite (Schnittwinkel der
Normalen mit der x-Achse), dadurch sind sie bestens für die praktische Berechnung
geeignet. Auch die Integrationen beziehen sich auf die geographische Breite, zur
Integration der Meridianbogenlängen habe ich die Grenzen für die Hilfswinkel entsprechend
angepaßt.
Gruß,
Klaus Nagel
Jan C. Hoffmann
Hallo Hermann,
letztlich ist es eine hohe Anforderung an die Numerik:
> r = a*b/sqrt( a^2*sin^2(phi) + b^2*cos^2(phi) )
Umfang: 40003423,0626111 m
>Die größten Abweichungen erhältst Du für einen
>kurzen Ellipsenbogen bei phi = pi/4 +- w...
phi = 45° - 1°: L = -111123,944115787 m
phi = 45° + 1°: L = 111117,450813601 m
Genauigkeit dürfte bei >= 9 Stellen liegen:
z.B. 111117,450 m
--
Gruss Jan C. Hoffmann (c/o MTEC)
http://ourworld.compuserve.com/homepages/MTEC/
Basis:
Private Function Funktion(ByVal phi As Double) As Double
Dim f1, f2, fStr, d_phi As Double
d_phi = 10 ^ -10
f1 = r(phi - d_phi / 2)
f2 = r(phi + d_phi / 2)
fStr = (f2 - f1) / d_phi
Funktion = Sqr(r(phi) ^ 2 + fStr ^ 2)
End Function
Private Function r(phi As Double) As Double
Dim a, b As Double
a = 6377397.15508
b = 6356078.9629
r = a * b / Sqr(a ^ 2 * Sin(phi) ^ 2 + b ^ 2 * Cos(phi) ^ 2)
End Function
Jan C. Hoffmann
> Hallo Jan,
>
> einstweilen vielen Dank für die Ausführungen.
>
> In der Funktion für die GaussIntegration werden 32 Array Variablen x und w
> definiert:
>
> x[ 1] = 0.024350292663424
> x[ 2] = 0.072993121787799
> x[ 3] = 0.12146281929612
> x[ 4] = 0.16964442042399
> x[ 5] = 0.21742364374001
> ...
> x[32] = 0.99930504173577
>
> w[ 1] = 0.04869095700914
> w[ 2] = 0.048575467441503
> w[ 3] = 0.048344762234803
> w[ 4] = 0.047999388596458
> w[ 5] = 0.04754016571483
> ...
> w[32] = 1.7832807216964E-03
>
Hallo Joerg,
da empfehle ich folgende Lteratur:
[1] H. Selder, Einführung in die Numerische Mathematik für Ingenieure, Hanser -
[2] H.R. Schwarz, Numerische Mathematik, B.G. Teubner Stuttgart
[1] ist einfacher zu verstehen.
Klaus Nagel
"Jan C. Hoffmann" schrieb:
> letztlich ist es eine hohe Anforderung an die Numerik:
Hallo Jan,
das finde ich nicht, der Integrand ist sehr harmlos, für kleine
Exzentrizitäten fast konstant. Das sieht man auch daran, daß das
Simpsonverfahren mit nur 12 Intervallen den gesamten Erdquadranten auf
Mikrometer genau liefert und Gauß-Integration mit 64 Stützstellen heißt
mit Kanonen auf Spatzen schießen. Deine Gradangaben beziehen sich auf
den Hilfswinkel, nicht auf die von Jürgen gewünschte geographische
Breite. Ich füge meine Ergebnisse bei, für Hilfswinkel und geographische
Breite.
> >Die größten Abweichungen erhältst Du für einen
> >kurzen Ellipsenbogen bei phi = pi/4 +- w...
>
> phi = 45° - 1°: L = -111123,944115787 m
> phi = 45° + 1°: L = 111117,450813601 m
>
Meridianbogen zwischen 44.00° 45.00°
(Bezogen auf Hilfswinkel)
Direkt 111117.4507979489642393
Gauss (n=6) 111117.4507979489642393
Simpson (n=12) 111117.4507979489634432
Meridianbogen zwischen 45.00° 46.00°
(Bezogen auf Hilfswinkel)
Direkt 111123.9440220625195308
Gauss (n=6) 111123.9440220625195308
Simpson (n=12) 111123.9440220625213173
Meridianbogen zwischen 0.00° 90.00°
(Bezogen auf Hilfswinkel)
Direkt 10000855.7647713774049398
Gauss (n=6) 10000855.7647672021632249
Simpson (n=12) 10000855.7647713774049398
Meridianbogen zwischen 44.00° 45.00°
(bezogen auf geographische Breite)
Direkt 111109.7128632500099358
Gauss (n=6) 111109.7128632500099358
Simpson (n=12) 111109.7128632500088922
Reihe (festes e) 111109.7128744436156210
Reide (bel. e) 111109.7128632507207942
Meridianbogen zwischen 45.00° 46.00°
(bezogen auf geographische Breite)
Direkt 111129.1923173713741983
Gauss (n=6) 111129.1923173713741983
Simpson (n=12) 111129.1923173713757377
Reihe (festes e) 111129.1923289540033686
Reide (bel. e) 111129.1923173707075088
Gruß,
Klaus Nagel
Jan C. Hoffmann
"Klaus Nagel" <nagel...@T-Online.de> schrieb im Newsbeitrag
>
>
> "Jan C. Hoffmann" schrieb:
>
> > letztlich ist es eine hohe Anforderung an die Numerik:
>
> Hallo Jan,
> das finde ich nicht, der Integrand ist sehr harmlos, für kleine
> Exzentrizitäten fast konstant. Das sieht man auch daran, daß das
> Simpsonverfahren mit nur 12 Intervallen den gesamten Erdquadranten auf
> Mikrometer genau liefert und Gauß-Integration mit 64 Stützstellen heißt
> mit Kanonen auf Spatzen schießen. Deine Gradangaben beziehen sich auf
> den Hilfswinkel, nicht auf die von Jürgen gewünschte geographische
> Breite. Ich füge meine Ergebnisse bei, für Hilfswinkel und geographische
> Breite.
Hallo Klaus,
ich habe absichtlich nur die Numerik (Differentiation und Integration) einer
Ellipse behandelt, nachdem Du folgende Webseite hier bekanntgegeben hast.
http://home.wtal.de/aw/excel/excel_koordinaten.htm
http://home.wtal.de/aw/excel/source/excel_koordinaten.txt
Das hab ich mal durchgespielt:
Sub Meridianbogen0()
Dim L1, L2, pi As Double
pi = Application.WorksheetFunction.pi
L2 = Meridianbogen(pi / 4, 6377397.15508, 6356078.9629)
L1 = Meridianbogen(pi / 4 - 3.14159 / 180, 6377397.15508, 6356078.9629)
Cells(5, 5).Value = L2 - L1
L1 = Meridianbogen(pi / 4, 6377397.15508, 6356078.9629)
L2 = Meridianbogen(pi / 4 + 3.14159 / 180, 6377397.15508, 6356078.9629)
Cells(6, 5).Value = L2 - L1
End Sub
Function Meridianbogen(phi, a, b) ' in m
Dim N, alpha, beta, gamma, delta, epsilon
N = (a - b) / (a + b)
alpha = (a + b) / 2 * (1 + 1 / 4 * N ^ 2 + 1 / 64 * N ^ 4)
beta = -3 / 2 * N + 9 / 16 * N ^ 3 - 3 / 32 * N ^ 5
gamma = 15 / 16 * N ^ 2 - 15 / 32 * N ^ 4
delta = -35 / 48 * N ^ 3 + 105 / 256 * N ^ 5
epsilon = 315 / 512 * N ^ 4
Meridianbogen = alpha * (phi + beta * Sin(2 * phi) + gamma * Sin(4 * phi) +
delta * Sin(6 * phi) + epsilon * Sin(8 * phi))
End Function
111109,6190185180
111129,0984401270
ok.
Dafür braucht man keine Integration mehr und spätestens dann war die Frage
diesbezüglich beantwortet. Auf dieser Webseite findet man aber nicht, wie man
die Bogenlänge einer Ellipse berechnet.
Meine Berechnungsmethode für Ellipse:
Pol im Mittelpunkt, große Achse ist Polarachse.
Klaus Nagel
"Jan C. Hoffmann" schrieb:
> Auf dieser Webseite findet man aber nicht, wie man
> die Bogenlänge einer Ellipse berechnet.
>
Hallo Jan,
diese Aussage verstehe ich nicht. Die Funktion Meridianbogen(phi, a, b)
liefert doch genau die Bogenlänge der Ellipse mit Halbachsen a und b
zwischen der geographischen Breite 0° (Äquator, große Achse) und der
geographischen Breite phi. Jeden beliebigen Ellipsenbogen kann man damit
als Differenz zweier Werte ermitteln.
Oder bemängelst Du, daß nicht angegeben ist, wie man zu dieser Formel
gelangt?
Die numerischen Integrationen habe ich nur gemacht, um zu prüfen, wie
gut die Reihenentwicklung ist. Sie liefert für alle Breiten exakte Werte
mit der Stellenzahl, mit der a und b gegeben sind, und mehr ist nicht
sinnvoll.
Gruß,
Klaus Nagel
Jan C. Hoffmann
"Klaus Nagel" <nagel...@T-Online.de> schrieb im Newsbeitrag
>
>
> "Jan C. Hoffmann" schrieb:
>
> > Auf dieser Webseite findet man aber nicht, wie man
> > die Bogenlänge einer Ellipse berechnet.
> >
> Hallo Jan,
> diese Aussage verstehe ich nicht. Die Funktion Meridianbogen(phi, a, b)
> liefert doch genau die Bogenlänge der Ellipse mit Halbachsen a und b
> zwischen der geographischen Breite 0° (Äquator, große Achse) und der
> geographischen Breite phi. Jeden beliebigen Ellipsenbogen kann man damit
> als Differenz zweier Werte ermitteln.
Hallo Klaus,
ist das wirklich noch eine 'mathematische' Ellipse oder ein 'getrimmtes'
Gebilde?
> Oder bemängelst Du, daß nicht angegeben ist, wie man zu dieser Formel
> gelangt?
Sieh Dir u.a. auch noch einmal das Posting vom 9.8. von Jürgen an. Interessant
ist vielleicht auch folgende Webseite:
http://ourworld.compuserve.com/homepages/MTEC/simpson-trapezoidal.htm
- aus einer früheren Diskussion -
> Die numerischen Integrationen habe ich nur gemacht, um zu prüfen, wie
> gut die Reihenentwicklung ist. Sie liefert für alle Breiten exakte Werte
> mit der Stellenzahl, mit der a und b gegeben sind, und mehr ist nicht
> sinnvoll.
Das ist richtig.
Hermann Kremer
>"Klaus Nagel" <nagel...@T-Online.de> schrieb im Newsbeitrag
>news:3B796664...@T-Online.de...
>> "Jan C. Hoffmann" schrieb:
>>
>> > Auf dieser Webseite findet man aber nicht, wie man
>> > die Bogenlänge einer Ellipse berechnet.
>> >
>> Hallo Jan,
>> diese Aussage verstehe ich nicht. Die Funktion Meridianbogen(phi, a, b)
>> liefert doch genau die Bogenlänge der Ellipse mit Halbachsen a und b
>> zwischen der geographischen Breite 0° (Äquator, große Achse) und der
>> geographischen Breite phi. Jeden beliebigen Ellipsenbogen kann man damit
>> als Differenz zweier Werte ermitteln.
>
>Hallo Klaus,
>
>ist das wirklich noch eine 'mathematische' Ellipse oder ein 'getrimmtes'
>Gebilde?
Hallo Jan,
es ist eine mathematische Ellipse. Stell es Dir so vor: Man pickt einen
Ellipsenpunkt heraus, berechnet dafür die Normale und den Punkt, in dem
die Normale die x-Achse schneidet, bestimmt den Winkel des Vektors
von Schnittpunkt nach Ellipsenpunkt und nennt diesen Winkel "geographische
Breite". Da wird doch nichts dabei "getrimmt". Es gibt halt viele Möglichkeiten
zur Definition der gleichen Ellipse ...
>> Oder bemängelst Du, daß nicht angegeben ist, wie man zu dieser Formel
>> gelangt?
Die exakte Formel für den Breitenwinkel als Funktion z.B. des geozentrischen
Winkels ist nicht allzu schwierig herzuleiten, und dann kann man diese für sehr
kleine Exzentrizitäten geeignet approximieren.
>Sieh Dir u.a. auch noch einmal das Posting vom 9.8. von Jürgen an.
Ja, er benutzt geographische Breiten.
>Interessant ist vielleicht auch folgende Webseite:
>http://ourworld.compuserve.com/homepages/MTEC/simpson-trapezoidal.htm
Bei großen Exzentrizitäten brauchst Du sicherlich irgend eine numerische
Quadratur für die elliptischen Integrale - da gibt es aber speziell dafür
zugeschnittene Algorithmen, z.B. in
http://www.netlib.org
Grüße
Hermann
--
Klaus Nagel
Hermann Kremer schrieb:
> es ist eine mathematische Ellipse. Stell es Dir so vor: Man pickt einen
> Ellipsenpunkt heraus, berechnet dafür die Normale und den Punkt, in dem
> die Normale die x-Achse schneidet, bestimmt den Winkel des Vektors
> von Schnittpunkt nach Ellipsenpunkt und nennt diesen Winkel "geographische
> Breite".
Die Bezeichnung ist sinnvoll, denn sie entspricht der Polhöhe, das ist
der Winkel zwischen der Horizontalen und dem Himmelspol. Ganz grob kann
man - auf der nördlichen Halbkugel - seine geographische Breite an der
Höhe des Polarsterns erkennen; genauer wird es, wenn man die
Winkelhalbierende zwischen der oberen und unteren Kulmination des
Polarsterns ( oder eines Zirkumpolalarsterns ) ermittelt.
> >> Oder bemängelst Du, daß nicht angegeben ist, wie man zu dieser Formel
> >> gelangt?
>
> Die exakte Formel für den Breitenwinkel als Funktion z.B. des geozentrischen
> Winkels ist nicht allzu schwierig herzuleiten, und dann kann man diese für sehr
> kleine Exzentrizitäten geeignet approximieren.
>
Ich hatte eigentlich an die Koeffizienten der Reihenentwicklung gedacht
und da ist mir noch nicht klar, wie man auf die kommt. Für die
Hilfswinkel weiß ich es, nicht aber für die geographische Breite.
Gruß,
Klaus
Jan C. Hoffmann
Hallo Hermann, Hallo Klaus
Ellipse: (x/a)^2 + (y/b)^2 = 1
Ich habe die Meridianfunktion L [1] bei 20 und 70° differenziert. Dann a und b
berechnet. Das gleiche bei 0 und 45° und a und b berechnet.
Wäre die Meridianfunktion eine Ellipse, müßte a und b konstant bleiben. Nach
Bronstein hat das wirklich nach einem getrimmten Gebilde ausgesehen.
phi = 20° und 70°:
a = 6390550,145616
0,000000 = L_Str(20*PI()/180;a;b) - c_1
c_1 = 6388127,803802 = L'(20°)_Meridian
b = 6369930,880496
0,000000 = L_Str(70*PI()/180;a;b) - c_2
c_2 = 6372332,572937 = L'(70°)_Meridian
phi = 0° und 45°:
a = 6334832,032602
0,000000 = L_Str(0*PI()/180;a;b) - c_1
c_1 = 6334832,032602 = L'(0°)_Meridian
b = 6427532,386276
0,000000 = L_Str(45*PI()/180;a;b) - c_2
c_2 = 6380677,223206 = L'(45°)_Meridian
Abweichungen:
delta-a = 6390550,145616 - 6334832,032602 = 55718,113014
delta-b = 6369930,880496 - 6427532,386276 = -57601,50578
delta-a und delta-b sollten Null sein.
Basis:
Berechnet wurden a und b über 2 nichtlineare Gleichungen mit den Funktionen, wie
zuvor angewandt:
Function L_Str(ByVal phi, ByVal a, ByVal b As Double) As Double
Dim f1, f2, fStr, d_phi As Double
d_phi = 10 ^ -30
f1 = r(phi - d_phi / 2, a, b)
f2 = r(phi + d_phi / 2, a, b)
fStr = (f2 - f1) / d_phi
L_Str = Sqr(r(phi, a, b) ^ 2 + fStr ^ 2)
End Function
Function r(ByVal phi, ByVal a, ByVal b As Double) As Double
r = a * b / Sqr(a ^ 2 * Sin(phi) ^ 2 + b ^ 2 * Cos(phi) ^ 2)
End Function
Function Meridianbogen(phi) As Double ' phi in ° Erg in m
Dim N, alpha, beta, gamma, delta, epsilon, a, b As Double
a = 6377397.15508
b = 6356078.9629
N = (a - b) / (a + b)
alpha = (a + b) / 2 * (1 + 1 / 4 * N ^ 2 + 1 / 64 * N ^ 4)
beta = -3 / 2 * N + 9 / 16 * N ^ 3 - 3 / 32 * N ^ 5
gamma = 15 / 16 * N ^ 2 - 15 / 32 * N ^ 4
delta = -35 / 48 * N ^ 3 + 105 / 256 * N ^ 5
epsilon = 315 / 512 * N ^ 4
Meridianbogen = alpha * (phi + beta * Sin(2 * phi) + gamma * Sin(4 * phi) +
delta * Sin(6 * phi) + epsilon * Sin(8 * phi))
End Function
Klaus Nagel
"Jan C. Hoffmann" schrieb:
> Ellipse: (x/a)^2 + (y/b)^2 = 1
>
> Ich habe die Meridianfunktion L [1] bei 20 und 70° differenziert. Dann a und b
> berechnet. Das gleiche bei 0 und 45° und a und b berechnet.
>
> Wäre die Meridianfunktion eine Ellipse, müßte a und b konstant bleiben. Nach
> Bronstein hat das wirklich nach einem getrimmten Gebilde ausgesehen.
>
>
>
> Function L_Str(ByVal phi, ByVal a, ByVal b As Double) As Double
> Dim f1, f2, fStr, d_phi As Double
> d_phi = 10 ^ -30
> f1 = r(phi - d_phi / 2, a, b)
> f2 = r(phi + d_phi / 2, a, b)
> fStr = (f2 - f1) / d_phi
> L_Str = Sqr(r(phi, a, b) ^ 2 + fStr ^ 2)
> End Function
>
> Function r(ByVal phi, ByVal a, ByVal b As Double) As Double
> r = a * b / Sqr(a ^ 2 * Sin(phi) ^ 2 + b ^ 2 * Cos(phi) ^ 2)
> End Function
>
> Function Meridianbogen(phi) As Double ' phi in ° Erg in m
> Dim N, alpha, beta, gamma, delta, epsilon, a, b As Double
> a = 6377397.15508
> b = 6356078.9629
> N = (a - b) / (a + b)
> alpha = (a + b) / 2 * (1 + 1 / 4 * N ^ 2 + 1 / 64 * N ^ 4)
> beta = -3 / 2 * N + 9 / 16 * N ^ 3 - 3 / 32 * N ^ 5
> gamma = 15 / 16 * N ^ 2 - 15 / 32 * N ^ 4
> delta = -35 / 48 * N ^ 3 + 105 / 256 * N ^ 5
> epsilon = 315 / 512 * N ^ 4
> Meridianbogen = alpha * (phi + beta * Sin(2 * phi) + gamma * Sin(4 * phi) +
> delta * Sin(6 * phi) + epsilon * Sin(8 * phi))
> End Function
Hallo Jan,
es kommt zwar in allen Funktionen der Winkel phi vor, er bedeutet aber
nicht immer das gleiche. In L_Str und r bezeichnet phi den
geozentrischen Winkel, bei der Funktion Meridianbogen ist er die
geographische Breite (in Bogenmaß, nicht in Grad wie der Kommentar
sagt). Das erklärt die Abweichungen. Die Zusammenhänge zwischen den
Winkeln macht man sich am besten klar, indem man eine Ellipse zeichnet
mit einer Tangente in einem Ellipsenpunkt P im ersten Quadranten. Die
Tangente schneidet die x-Achse im Punkt Q, irgendwo rechts der Ellipse.
Der Schnittwinkel ist komplementär zur geographischen Breite. Streckt
man die Zeichnung um den Faktor a/b in y-Richtung, dann geht die Ellipse
in einen Kreis mit Radius a über. Das Bild der Tangente tangiert den
Kreis in P', der Schnittpunkt Q bleibt unverändert. Aus dieser Zeichnung
kann mit einfachen geometrischen Überlegungen die Beziehungen zwischen
den Winkeln ermitteln.
Die geographische Breite ist ein sinnvoller Parameter, um Orte auf der
Erde (Ellipsoid) zu bezeichen, denn sie kann aus der Messung von
Sternhöhen direkt gemessen werden. Die geographische Breite wurde schon
bestimmt, ehe die wahre Gestalt der Erde bekannt war, ob sie eine Kugel
ist, plattgedrückt (Oblatum) oder an den Polen in die Länge gezogen
(Oblongum). Im Gegenteil, um 1740 schickte die Pariser Akademie zwei
Expeditionen nach Lappland und Südamerika, um durch genaue Messungen der
Bogenlängen für weit entfernte geographische Breiten die Erdform zu
berechnen.
Gruß,
Klaus Nagel
Hermann Kremer
>"Jan C. Hoffmann" schrieb:
Hallo Jan,
>> Ellipse: (x/a)^2 + (y/b)^2 = 1
>>
>> [ SNIP ]
Als kleine Ergänzung zum Beitrag von Klaus:
>Die Zusammenhänge zwischen den
>Winkeln macht man sich am besten klar, indem man eine Ellipse zeichnet
>mit einer Tangente in einem Ellipsenpunkt P im ersten Quadranten. Die
>Tangente schneidet die x-Achse im Punkt Q, irgendwo rechts der Ellipse.
>Der Schnittwinkel ist komplementär zur geographischen Breite. Streckt
>man die Zeichnung um den Faktor a/b in y-Richtung, dann geht die Ellipse
>in einen Kreis mit Radius a über. Das Bild der Tangente tangiert den
>Kreis in P', der Schnittpunkt Q bleibt unverändert. Aus dieser Zeichnung
>kann mit einfachen geometrischen Überlegungen die Beziehungen zwischen
>den Winkeln ermitteln.
>Die geographische Breite ist ein sinnvoller Parameter, um Orte auf der
>Erde (Ellipsoid) zu bezeichen, denn sie kann aus der Messung von
>Sternhöhen direkt gemessen werden. Die geographische Breite wurde schon
>bestimmt, ehe die wahre Gestalt der Erde bekannt war, ob sie eine Kugel
>ist, plattgedrückt (Oblatum) oder an den Polen in die Länge gezogen
>(Oblongum). Im Gegenteil, um 1740 schickte die Pariser Akademie zwei
>Expeditionen nach Lappland und Südamerika, um durch genaue Messungen der
>Bogenlängen für weit entfernte geographische Breiten die Erdform zu
>berechnen.
Näheres dazu z.B. in
http://www.geodesy.tu-berlin.de/~kerstin/sites/erd3/vortraege/Gradmess.htm
Für die Ellipse (x/a)^2 + (y/b)^2 = 1 lautet in einem Ellipsenpunkt P(xP, yP)
die Tangentengleichung
(x * xP)/a^2 + (y * yP)/b^2 = 1
und die Normalengleichung
(y/yP - 1)*b^2 = (x/xP - 1)*a^2 .
Für den Schnittpunkt N der Normalen mit der x-Achse findet man daraus
xN = e^2 * xP ,
wobei e^2 = (a^2 - b^2)/a^2 das Quadrat der numer. Exzentrizität ist.
Für den g e o z e n t r i s c h e n Winkel psi (Winkel zwischen
Ellipsenmittelpunkt
im Koordinatenursprung und Punkt P) ergibt sich damit der g e o g r a p h i s c h
e
Winkel phi (Winkel zwischen Normalenschnittpunkt N und P), d.h. die geographische
Breite, zu
tan(phi) = yP/(xP - xN) = (yP/xP)/(1 - e^2) = tan(psi)/(1 - e^2) ,
d.h.
phi = arctan( tan(psi) / (1 - e^2) ) = arctan( (a/b)^2 * tan(psi) ) .
Für sehr kleine Exzentrizitäten e sind die beiden Winkel fast gleich, nicht aber
für größere Exzentrizitäten. Wenn man also Ellipsenstücke ausrechnet, muß man
immer peinlichst darauf achten, von welchem Punkt auf der x-Achse man den
Polarwinkel zählt:
Koordinatenursprung (geozentrisch),
Schnittpunkt N (geographisch),
Brennpunkt ( gibt es dafür eine eigene Bezeichnung ? )
Dazu kommt dann noch der Hilfswinkel t in der Parameterform
x = a*cos(t), y = b*sin(t) der Ellipse, das wäre dann der vierte Winkel ;-)
Grüße
Hermann
--
>Gruß,
>
>Klaus Nagel
Jan C. Hoffmann
"Klaus Nagel" <nagel...@T-Online.de> schrieb im Newsbeitrag
Hallo Klaus, Hallo Hermann,
leider habe ich über dieses spezielle Thema keine Literatur und habe angenommen,
daß das 'Gebilde' den Pol im Mittelpunkt hat.
Vielleicht kannst Du meine Funktionen korrigieren, falls Du die passende
Literatur oder Kenntnis dazu hast? Ein 'spezieller' Hinweis auf eine
Internetseite wäre willkommen.
Für mich zu interessant, um die Diskussion abzubrechen.
Klaus Nagel
"Jan C. Hoffmann" schrieb:
> leider habe ich über dieses spezielle Thema keine Literatur und habe angenommen,
> daß das 'Gebilde' den Pol im Mittelpunkt hat.
>
> Vielleicht kannst Du meine Funktionen korrigieren, falls Du die passende
> Literatur oder Kenntnis dazu hast? Ein 'spezieller' Hinweis auf eine
> Internetseite wäre willkommen.
Hallo Jan,
hier ist ein Verfahren, mit numerischer Integration die Bogenlänge zu
bestimmen. Es benutzt nicht Deine Methode mit Polarkoordinaten, der
Hilfswinkel psi ist daher etwas anderes als Dein Hilfswinkel phi.
Ellipse mit den Halbachsen a und b. Die Bogenlänge zwischen den
geographischen Breiten phi0 und phi1 berechnet sich wie folgt:
e^2 = 1 - (b/a)^2
psi0 = arctan(tan(phi0)*a/b)
psi1 = arctan(tan(phi1)*a/b)
f(psi)=a*sqrt(1-e^2*cos(psi)^2)
bogen = Integral(psi = psi0 .. psi1, f(psi))
Für die stark exzentrische Ellipse mit a=10 und b=1 lieferten
verschiedene numerische Integrationen folgende Werte für die Länge eines
Quadranten:
Gauß (n= 6) 10.1596114684
Simpson (n=12) 10.1620449110
Simpson (n=20) 10.1601372702
Simpson (n=40) 10.1599367110
Simpson (n=60) 10.1599354625
Simpson (n=80) 10.1599354504
Direkt 10.1599354502
Direkt ist der von Maple ermittelte Wert des elliptischen Integrals, ich
weiß nicht nach welchem Verfahren. Gauß mit n=12 liefert vermutlich ein
gutes Ergebnis, aber ich hatte keine Lust die Koeffizienten einzutippen
oder zu berechnen.
Für kleine Exzentrizitäten ist die Reihenentwicklung vorzuziehen.
Gruß,
Klaus Nagel
Paul Ebermann
> Erst jetzt habe ich bemerkt, daß es ein Script ist, welches in HTML
> integriert wurde.
>
> Mal sehen, ob ich Java entschlüsselt bekomme ...
> Ansonsten gibt es sicherlich eine kompetente Newsgroup.
Bevor du an die falsche Newsgroup schreibst:
- Die in HTML eingebauten Scripte sind in JavaScript
programmiert.
- Applets sind in Java programmiert.
Diese beiden Programmiersprachen haben außer den ersten
4 Buchstaben und Ähnlichkeiten in der Syntax nichts gemeinsam.
Es gibt auch unterschiedliche Newsgroups:
Java: de.comp.lang.java (Homepage: http://www.dclj.de)
JavaScript: de.comp.lang.javascript
Nähere Infos dazu auch in den regelmäßigen Postings von
Hubert Partl in beiden Newsgroups ("Java != JavaScript")
Paul