Už se to tu kdysi řešilo, ale čas plyne a mě se konečně podařilo Vincentyho metodu přepsat do SQL. S klauzulí WITH už není problém v rámci jednoho dotazu iterovat. Čísla které mi dotaz vrací sedí na kontrolní formulář National Geodetic Survey při NOAA - čímžto si dovolím tvrdit, že mi dotaz kalkuluje čísla správná.
Začal jsem tedy tvořit dotaz, který by mi vrátil chybějící kešky pro doplnění statorového modulu FindsByBearingChallenge - můj výsledek se oproti statorovému rozchází asi o 5 sektorů (337 vs. 342). Můžete říct, že ta přesnost je zběsilá (v azimutech se lišíme v desetinách, ve vzdálenostech jsou to nízké jednotky metrů, případně centimetry), ale matematika je krásná v tom, že by se výsledky kalkulované stejnou metodou lišit neměly a pokud ano, tak velmi málo (na 4 nebo 5 desetinném místě). Toto je pak způsobené zaokrouhlováním a různou interpretací des. čísel v počítačových programech.
Takže "stay tuned", jsem v kontaktu s Gordem a HaLuMa a zjišťuju proč se výpočty liší - první moje podezření je na nedostatečný počet iterací Vincentyho metody v proceduře GG, která ji počítá. Kdysi se pro výpočet VM končilo s iterací kvůli výpočetní složitosti při průběžné diferenci kolem 10^-12, dnes se dá jít až na hranici datového typu při počítačové kalkulaci - mě to v SQLite končí při rozdílu 10^-19 - i tak se počet iterací pohybuje mezi sedmi a osmi. Vycházím ze znalosti inverzní procedury ProjectPoint, která počítá Vincentyho metodou projekci bodu. Tyto procedury autor převzal z již hotového řešení a z pochopitelných důvodů (viz. níže) se nezabýval její přesnou funkčností.
Pokud chcete vidět jak vypadá jeden z použitých vzorečků (vůbec se nedivím, že to autor nekontroloval, to by pozdrželo vývoj GG asi o dva roky - z toho půl roční pobyt v psychiatrické ozdravovně...):
(L + (1 - (f / 16 * (1 - power( (cosU1 * cosU2 * sin(lambdaIter) / (sqrt(power( (cosU2 * sin(lambdaIter) ), 2) + power( (cosU1 * sinU2 - sinU1 * cosU2 * cos(lambdaIter) ), 2) ) ) ), 2) ) * (4 + f * (4 - 3 * (1 - power( (cosU1 * cosU2 * sin(lambdaIter) / (sqrt(power( (cosU2 * sin(lambdaIter) ), 2) + power( (cosU1 * sinU2 - sinU1 * cosU2 * cos(lambdaIter) ), 2) ) ) ), 2) ) ) ) ) ) * f * (cosU1 * cosU2 * sin(lambdaIter) / (sqrt(power( (cosU2 * sin(lambdaIter) ), 2) + power( (cosU1 * sinU2 - sinU1 * cosU2 * cos(lambdaIter) ), 2) ) ) ) * ( (atan( (sqrt(power( (cosU2 * sin(lambdaIter) ), 2) + power( (cosU1 * sinU2 - sinU1 * cosU2 * cos(lambdaIter) ), 2) ) ) / (sinU1 * sinU2 + cosU1 * cosU2 * cos(lambdaIter) ) ) ) + (f / 16 * (1 - power( (cosU1 * cosU2 * sin(lambdaIter) / (sqrt(power( (cosU2 * sin(lambdaIter) ), 2) + power( (cosU1 * sinU2 - sinU1 * cosU2 * cos(lambdaIter) ), 2) ) ) ), 2) ) * (4 + f * (4 - 3 * (1 - power( (cosU1 * cosU2 * sin(lambdaIter) / (sqrt(power( (cosU2 * sin(lambdaIter) ), 2) + power( (cosU1 * sinU2 - sinU1 * cosU2 * cos(lambdaIter) ), 2) ) ) ), 2) ) ) ) ) * (sqrt(power( (cosU2 * sin(lambdaIter) ), 2) + power( (cosU1 * sinU2 - sinU1 * cosU2 * cos(lambdaIter) ), 2) ) ) * ( (CASE/* osetreni pokud se pocita czdalenost na rovniku -> deleni nulou */ WHEN (1 - power( (cosU1 * cosU2 * sin(lambdaIter) / (sqrt(power( (cosU2 * sin(lambdaIter) ), 2) + power( (cosU1 * sinU2 - sinU1 * cosU2 * cos(lambdaIter) ), 2) ) ) ), 2) ) = 0 THEN 0 ELSE ( (sinU1 * sinU2 + cosU1 * cosU2 * cos(lambdaIter) ) - 2 * ( (sinU1 * sinU2) / (1 - power( (cosU1 * cosU2 * sin(lambdaIter) / (sqrt(power( (cosU2 * sin(lambdaIter) ), 2) + power( (cosU1 * sinU2 - sinU1 * cosU2 * cos(lambdaIter) ), 2) ) ) ), 2) ) ) ) END) + (f / 16 * (1 - power( (cosU1 * cosU2 * sin(lambdaIter) / (sqrt(power( (cosU2 * sin(lambdaIter) ), 2) + power( (cosU1 * sinU2 - sinU1 * cosU2 * cos(lambdaIter) ), 2) ) ) ), 2) ) * (4 + f * (4 - 3 * (1 - power( (cosU1 * cosU2 * sin(lambdaIter) / (sqrt(power( (cosU2 * sin(lambdaIter) ), 2) + power( (cosU1 * sinU2 - sinU1 * cosU2 * cos(lambdaIter) ), 2) ) ) ), 2) ) ) ) ) * (sinU1 * sinU2 + cosU1 * cosU2 * cos(lambdaIter) ) * ( -1 + 2 * power( (CASE/* osetreni pokud se pocita czdalenost na rovniku -> deleni nulou */ WHEN (1 - power( (cosU1 * cosU2 * sin(lambdaIter) / (sqrt(power( (cosU2 * sin(lambdaIter) ), 2) + power( (cosU1 * sinU2 - sinU1 * cosU2 * cos(lambdaIter) ), 2) ) ) ), 2) ) = 0 THEN 0 ELSE ( (sinU1 * sinU2 + cosU1 * cosU2 * cos(lambdaIter) ) - 2 * ( (sinU1 * sinU2) / (1 - power( (cosU1 * cosU2 * sin(lambdaIter) / (sqrt(power( (cosU2 * sin(lambdaIter) ), 2) + power( (cosU1 * sinU2 - sinU1 * cosU2 * cos(lambdaIter) ), 2) ) ) ), 2) ) ) ) END), 2) ) ) ) ) - lambda
Samotný kód uveřejním jakmile bude dotestováno a rekonciliováno na Stator, kterýžto pokládám za prubířský kámen a měl by vracet co nejpřesnější výsledky. Pak zase sepíšu nějaký pěkný článek na blog.
)).
GeoGet
URWIGO
CWG
