alg/delphi: Koers tussen 2 punten op aarde berekenen

[Megamind]

Ik ben in delphi een applicatie aan het maken, ik kan daarmee al de afstand tussen 2 punten (lat/long) uitrekenen maar nu nog de koers. Het lijkt me heel simpel te doen met wat cosinus formules maar ik wil ook graag grootcircel.

Wie weet wat een goede formule is om te gebruiken?

Op google heb ik al gezocht maar daar kan ik alleen maar vinden van afstand tussen 2 punten...

[Felixje]

Wat bedoel je met grootcircel?

[Megamind]

Nou je kan ook grootcircel vliegen op aarde, lang verhaal, vrij ingewikkeld, maar het gaat erom dat je met een bocht over de aarde vliegt om zo een kortere afstand te krijgen als direct tussen 2 punten (loxodrome koers)

[Felixje]

Dus als ik het goed begrijp, in plaats van een recht lijn over de kaart trekken, er een stuk boven of onder vliegen? Misschien kan je dan de functie van de cos combineren met een kwardratische? Maar hoe dat moet weet ik niet
Zo ver ben ik nog niet met wiskunde

[ Voor 8% gewijzigd door Felixje op 05-02-2005 23:34 ]

[MBV]

Hij bedoelt dat je een touwtje over een bol heen trekt, ipv over een landkaart. Vraag maar aan de KLM hoe ze naar New York vliegen: over groenland, wat je dus echt niet zou verzinnen met alleen een landkaart in je achterhoofd...

[Anoniem: 30225]

Megamind schreef op zaterdag 05 februari 2005 @ 20:47:
Nou je kan ook grootcircel vliegen op aarde, lang verhaal, vrij ingewikkeld, maar het gaat erom dat je met een bocht over de aarde vliegt om zo een kortere afstand te krijgen als direct tussen 2 punten (loxodrome koers)

Kortere afstand?
Ik zou eerder zeggen, langere afstand.

Edit:
Het is maar net hoe je het bekijkt. Als je een routeplanner pakt, en je reist met de auto moet je rekening houden met het wegennet.
Als er een mogelijkheid zou zijn om rechtstreeks van A naar B te rijden dan zou dat korter zijn.
Op grotere hoogtes is de omtrek van de aarde ook groter, dus of de route dan ook korter is is ook nog discutabel. Daarnaast moet je ook nog opstijgen en landen.

[ Voor 35% gewijzigd door Anoniem: 30225 op 06-02-2005 00:35 ]

[Delphi32]

Nou nee hoor, vliegen met een bochtje kan wel degelijk veel korter zijn Ik wilde eerst een heel voorbeeld gaan uittypen, maar toen kwam ik op deze site die het allemaal mooi uitlegt.
Mijn wiskundekennis is nogal stoffig (as in: niet aanwezig), dus op het gebied van koersen blijf ik een antwoord schuldig. Ben wel nieuwsgierig naar het antwoord trouwens.

[Mithrandir]

De grap is, het is eigenljik geen bocht. het is gewoon de korste route op een bol. Op een landkaart (wat eigenlijk een bol geprojecteerd op een plat vlak is) kloppen de hoeken niet. Daardoor LIJKT het alsof een bocht het snelst vliegt, maar in werkelijkheid is het een rechte lijn.

[ Voor 10% gewijzigd door Mithrandir op 06-02-2005 12:45 ]

[.oisyn]

Megamind schreef op zaterdag 05 februari 2005 @ 20:47:
Nou je kan ook grootcircel vliegen op aarde, lang verhaal, vrij ingewikkeld, maar het gaat erom dat je met een bocht over de aarde vliegt om zo een kortere afstand te krijgen als direct tussen 2 punten (loxodrome koers)

Hoe wil je die koers hebben? Heb je voldoende aan het verkrijgen van een rotatie om een bepaalde as?

Het is vrij simpel als je enigszins bekend bent met lineaire algebra. Maak van je begin- en eindpunt vectoren die de (3d) positie op het aardoppervlak aangeven, met het middelpunt van de aarde als de oorsprong. De straal van de aarde maakt hierbij niet echt uit, het is handig als je daar gewoon 1 van maakt, en daar ga ik ook even vanuit in de volgende berekeningen.

Het uitwendig product van die 2 vectoren is de (nog niet genormalizeerde) as waarover de rotatie plaats moet gaan vinden. Bereken de lengte van deze vector. De inverse sinus van deze lengte is de hoek tussen de 2 punten (gezien vanuit de oorsprong). Als je nu het beginpunt roteert over de berekende as met de berekende hoek kom je uit op het eindpunt.

Anoniem: 30225 schreef op zondag 06 februari 2005 @ 00:26:
[...]

Kortere afstand?
Ik zou eerder zeggen, langere afstand.

Ik denk dat de topicstarter bedoelt dat dat korter is dan de rechte lijn die je op landkaarten ziet, omdat de longitudinale en latitudinale (?) lijnen zijn platgeslagen om zo loodrecht op elkaar te liggen, wat in de werkelijkheid natuurlijk niet zo is

[MrBarBarian]

Delphi32 schreef op zondag 06 februari 2005 @ 01:01:
Nou nee hoor, vliegen met een bochtje kan wel degelijk veel korter zijn Ik wilde eerst een heel voorbeeld gaan uittypen, maar toen kwam ik op deze site die het allemaal mooi uitlegt.
Mijn wiskundekennis is nogal stoffig (as in: niet aanwezig), dus op het gebied van koersen blijf ik een antwoord schuldig. Ben wel nieuwsgierig naar het antwoord trouwens.

Nee, das juist fout wat je hier zegt.. Een bol is niet in 2d te projecteren, zonder OF hoekvervorming OF oppervlakte vervorming (of beide natuurlijk).. daarom lijkt de kortste weg in 2d een kromme lijn te zijn, maar is in werkelijkheid kaarsrecht...

Anyway, als je 3 jaar eerder was geweest, had ik je zo die formules gegeven.. nu ben ik ze gelukkig weer vergeten.. Men, wat een drama was dat uitrekenen..

[.oisyn]

MrBarBarian schreef op zondag 06 februari 2005 @ 13:02:
Men, wat een drama was dat uitrekenen..

Drama? Easy as hell!

[Megamind]

Het is wel degelijk een bocht, je moet namelijk om van A naar B te komen bv elk uur een paar graden meer sturen om uit te komen, en toch heb je een flinke bekorting (er zijn wel wat voorwaarden), maar ik heb een aantal formules, alleen ik krijg ze niet in delphi gepropt (er komen andere waarden uit dan mijn navigatie rekenmachine). Ze zijn niet zo eenvoudig helaas, maar volgens mij gaat het erom dat ik mijn data erin zet in graden, minuten en hondersten (31 graden 21 minuten en 10 hondersten) maar delphi rekent het volgens mij met graden, tienden en hondersten (31 graden 89 tienden en 77 hondersten).

Nou wordt het een beetje te moeilijk voor mij

[MSalters]

Helemaal geen bocht, zoals .oisyn uitlegt. Je blijft namelijk in één vlak, het vlak door A,B en het middelpunt van de aarde. Dat is per definitie een grootcirkel, net zoals de evenaar, en je wilt toch niet beweren dat de evenaar slingert?

Delphi rekent ongetwijfeld met radialen.

[MrBarBarian]

Bovendien is het redelijk eenvoudig om van graden naar radialen en datum te converteren..

_{Ik ben niet zo een rekenwonder, en al helemaal niet op een bol, dus het was wel een drama}

[.oisyn]

Megamind schreef op zondag 06 februari 2005 @ 16:57:
Ze zijn niet zo eenvoudig helaas, maar volgens mij gaat het erom dat ik mijn data erin zet in graden, minuten en hondersten (31 graden 21 minuten en 10 hondersten) maar delphi rekent het volgens mij met graden, tienden en hondersten (31 graden 89 tienden en 77 hondersten).

Nou wordt het een beetje te moeilijk voor mij

Heb je mijn uitleg gezien (mijn eerste post in deze draad)? Die is niet zo heel erg moeilijk... Hoe wil je je resultaten eigenlijk hebben?

[Freee!!]

MBV schreef op zaterdag 05 februari 2005 @ 23:35:
Hij bedoelt dat je een touwtje over een bol heen trekt, ipv over een landkaart. Vraag maar aan de KLM hoe ze naar New York vliegen: over groenland, wat je dus echt niet zou verzinnen met alleen een landkaart in je achterhoofd...

Er is nog een bijkomende reden voor die route: op die manier blijven ze altijd in de buurt van land (veiligheid).

[Megamind]

MSalters schreef op zondag 06 februari 2005 @ 17:10:
Helemaal geen bocht, zoals .oisyn uitlegt. Je blijft namelijk in één vlak, het vlak door A,B en het middelpunt van de aarde. Dat is per definitie een grootcirkel, net zoals de evenaar, en je wilt toch niet beweren dat de evenaar slingert?

Delphi rekent ongetwijfeld met radialen.

Sorry je hebt gelijk, je maakt een bocht, op een platte kaart, dus in een merkatorkaart, als je die koers op een gnomomische projectie wilt hebben dan krijg je een juiste rechte lijn.

Ik heb deze functie in delphi geschreven:

code:

var
  LatA,LatB,LonA,LonB,Ver,Diff,Cgc,CgcA,CgcB:Single;
begin
  LatA:=DegToRad(StrToFloat(LatAEdit.text));
  LonA:=DegToRad(StrToFloat(LonAEdit.text));
  LatB:=DegToRad(StrToFloat(LatBEdit.text));
  LonB:=DegToRad(StrToFloat(LonBEdit.text));

  Diff:=DegToRad(StrToFloat(LonBEdit.text)-StrToFloat(LonAEdit.Text));

  Ver:=Sin(LatA)*Sin(LatB)+Cos(LatA)*Cos(LatB)*Cos(diff);
  Ver:=ArcCos(Ver);
  Ver:=RadToDeg(Ver);
  Ver:=Ver*60*0.8684;

  CgcA:=Sin(LatB)-Sin(LatA)*Cos(VeR);
  CgcB:=Cos(LatA)*Sin(Ver);
  Cgc:=CgcA/CgcB;
  Cgc:=ArcCos(Cgc);
  Cgc:=RadToDeg(Cgc);
  Caption:=floattostr(ver)+' '+FloatToStr(Cgc);

Hier komt uit 32,13 nm 62 graden.

Alleen volgens mijn rekenmachine moet het zijn, 25 mijl en 209 graden, met deze coordinaten:
LatA: 52,1759 LonA: 004,4559
LatB: 51,5700 LonB: 004,2700

Als ik dus de RadToDeg weglaat, dan komt ik uit op 1800 mijl

[ Voor 3% gewijzigd door Megamind op 06-02-2005 18:43 ]

[Delphi32]

Als je JavaScript kunt lezen (en vertalen naar Delphi) dan heb je misschien wat aan deze site met veel formules.

[MSalters]

Megamind schreef op zondag 06 februari 2005 @ 18:41:
[...]

Sorry je hebt gelijk, je maakt een bocht, op een platte kaart, dus in een merkatorkaart, als je die koers op een gnomomische projectie wilt hebben dan krijg je een juiste rechte lijn.

Ik heb deze functie in delphi geschreven:

code:

var
  LatA,LatB,LonA,LonB,Ver,Diff,Cgc,CgcA,CgcB:Single;
begin
  LatA:=DegToRad(StrToFloat(LatAEdit.text));
  LonA:=DegToRad(StrToFloat(LonAEdit.text));
  LatB:=DegToRad(StrToFloat(LatBEdit.text));
  LonB:=DegToRad(StrToFloat(LonBEdit.text));

  Diff:=DegToRad(StrToFloat(LonBEdit.text)-StrToFloat(LonAEdit.Text));

  Ver:=Sin(LatA)*Sin(LatB)+Cos(LatA)*Cos(LatB)*Cos(diff);
  Ver:=ArcCos(Ver);
  Ver:=RadToDeg(Ver);
  Ver:=Ver*60*0.8684;

  CgcA:=Sin(LatB)-Sin(LatA)*Cos(VeR);
  CgcB:=Cos(LatA)*Sin(Ver);
  Cgc:=CgcA/CgcB;
  Cgc:=ArcCos(Cgc);
  Cgc:=RadToDeg(Cgc);
  Caption:=floattostr(ver)+' '+FloatToStr(Cgc);

Hier komt uit 32,13 nm 62 graden.

toon volledige bericht

Zonder dat ik Delphi ken weet ik al dat dit niet goed kan zijn: Ik mis vrij belangrijke variabelen als de straal van de aarde, De eerste assignment aan ver lijkt me incorrect, mis je daar haakjes of zo? (als LatA=LatB=arcsin(1) dan is Ver=1, zelfs als ook LonA=LonB )

[.oisyn]

Voor de derde keer, wat wil je nou precies uitrekenen?

MSalters: er klopt wel meer niet, zo doet hij op regel 13 namelijk een RadToDeg, terwijl hij het resultaat daarvan (na wat vermenigvuldigingen) weer gebruikt als invoer voor Sin op regel 17. Aangezien de Sin op regel 17 zijn invoer nog altijd in radialen verwacht klopt die RadToDeg gewoon niet.

[ Voor 72% gewijzigd door .oisyn op 06-02-2005 23:26 ]

[Megamind]

.oisyn schreef op zondag 06 februari 2005 @ 23:23:
Voor de derde keer, wat wil je nou precies uitrekenen?

MSalters: er klopt wel meer niet, zo doet hij op regel 13 namelijk een RadToDeg, terwijl hij het resultaat daarvan (na wat vermenigvuldigingen) weer gebruikt als invoer voor Sin op regel 17. Aangezien de Sin op regel 17 zijn invoer nog altijd in radialen verwacht klopt die RadToDeg gewoon niet.

Ik wil uitrekenen als ik in punt A sta wat de koerst is tov het noorden is om bij punt B uit te komen. Er zijn inderdaad een paar dingen niet goed gegaan in die formule, ik heb hem nu verbeterd met:

code:

function Distance(ALat1, ALon1, ALat2, ALon2: Extended):Extended;
var
 LTheta : Extended;
 LDist  : Extended;
begin
 LTheta := ALon1 - ALon2;
 LDist := Sin(DegToRad(ALat1)) * Sin(DegToRad(ALat2)) + Cos(DegToRad(ALat1)) *
          Cos(DegToRad(ALat2)) * Cos(DegToRad(LTheta));
 LDist := ACos(LDist);
 LDist := Rad2Deg(LDist);
 Result := LDist * 60;
end;

Label1.Caption := FloatToStr(RoundTo(Course(052.308056, 004.764167, 051.956944, 004.437222),-4));

Deze waarden komen wel aardig in de buurt, het zijn nu ook andere coordinaten, niet meer met minuten. Alleen dit is dus een afstand in nm. Als ik hiermee een koers formule gebruik, dan komt er een getal iets van 1,3... uit..

Ik snap dus dat Delphi met radialen werkt, en dit is de koers formule:

code:

  Crs := (Sin(Deg2Rad(ALon2)) - Sin(Deg2Rad(ALon1)) * Cos (LDist) )/
         (Cos(Deg2Rad(ALon1)) * Sin(LDist));

  Crs := ArcCos(Crs);

Ik moet dus hier de graden niet naar radialen omtoveren of juist radialen naar graden..

* Megamind was niet zo goed in algebra

[ Voor 13% gewijzigd door Megamind op 07-02-2005 07:49 ]

[kasper_vk]

MBV schreef op zaterdag 05 februari 2005 @ 23:35:
Hij bedoelt dat je een touwtje over een bol heen trekt, ipv over een landkaart. Vraag maar aan de KLM hoe ze naar New York vliegen: over groenland, wat je dus echt niet zou verzinnen met alleen een landkaart in je achterhoofd...

Om het vraagstukje nog ietsje moeilijker te maken: de aarde is zelfs geen exacte bol, maar een ietsje platgeslagen bol. (vandaar de zwaartekracht op de polen minder sterk is als op de evenaar.)
Ook hierdoor zal een route die ietje meer naar de polen buigt korter zijn.
(Uiteindelijk werkt de methode die .oisyn noemde dan natuurlijk nog steeds: je hebt hebt de 'snijlijn' nodig van de aarde en het vlak door het middelpunt, beginpunt en eindpunt.)

[curry684]

Kun je niet gewoon Mercatorprojectie gebruiken? Op Wikipedia staan de omrekenformules daarvoor

[.oisyn]

De hoek met het noorden kun je wel berekenen door het noordelijkste punt (0, 1, 0) en het eindpunt te projecteren op een vlak wat parallel loopt aan het vlak dat de aarde raakt in het beginpunt. De normaal van dit vlak is gelijk aan de positie van het beginpunt op de eenheidsbol (aarde met straal 1 ) Je kunt de hoek weer berekenen door het uitwendige product te nemen van die 2 geprojecteerde punten, en deze evt. negatief te maken als dat uitwendig product de aarde in wijst ipv uit.

In pseudocode:

C++:

vector A = polarToVector (latA, lonA);    // A = beginpunt
vector B = polarToVector (latB, lonB);    // B = eindpunt
scalar distance = acos (dot (A, B)) * AARDE_STRAAL;   // afstand AB

// nu beide punten projecteren op een vlak dat A als normaal heeft.
vector N = (0, 1, 0);   // het noorden
vector NprojA = N - N * dot (A, N);   // N geprojecteerd op A
vector BprojA = B - B * dot (A, B);   // B geprojecteerd op A
vector axis = cross (NprojA, BprojA);  // bereken het uitwendig product
scalar angle = asin (length (axis));   // de hoek met het noorden
if (dot (A, axis) < 0)    // kijk of het uitwendig product de aarde in wijst
    angle = -angle;  // zo ja, dan is de hoek negatief

Let wel dat hier nog geen enkele controle is of de opgegeven coordinaten wel kloppen. Als A en B aan weerskanten van de aarde liggen dan maakt het natuurlijk niet uit in welke richting je reist. Eveneens als A op de zuidpool ligt dan is elke kant richting het noorden.

[ Voor 10% gewijzigd door .oisyn op 07-02-2005 11:44 ]

[SchizoDuckie]

.oisyn schreef op zondag 06 februari 2005 @ 12:51:
[...]
Het uitwendig product van die 2 vectoren is de (nog niet genormalizeerde) as waarover de rotatie plaats moet gaan vinden. Bereken de lengte van deze vector. De inverse sinus van deze lengte is de hoek tussen de 2 punten (gezien vanuit de oorsprong). Als je nu het beginpunt roteert over de berekende as met de berekende hoek kom je uit op het eindpunt.

Holy crap

Nou had ik niet meer dan een 5 voor wiskunde meestal, dus ik ben je bij die vectoren al kwijtgeraakt... Heb je dit ook ergens in jip-en-janneke taal voor de minder wiskundig begaafden onder ons?

[.oisyn]

De functies in mijn vorige stukje code zijn overigens als volgt:

C++:

vector polarToVector (scalar latitude, scalar longitude)
{
    vector r;
    r.x = -cos (latitude) * sin (longitude);
    r.y = sin (latitude);
    r.z = cos (latitude) * cos (longitude);
    return r;
}

scalar dot (vector a, vector b)
{
    return a.x * b.x + a.y * b.y + a.z * b.z;
}

vector cross (vector a, vector b)
{
    vector r;
    r.x = a.y * b.z - a.z * b.y;
    r.y = a.z * b.x - a.x * b.z;
    r.z = a.x * b.y - a.y * b.x;
    return r;
}

scalar length (vector v)
{
    return sqrt (v.x * v.x + v.y * v.y + v.z * v.z);
}