From 0062067d431922133225aec6f738024b41b4fe42 Mon Sep 17 00:00:00 2001 From: tsteven4 <13596209+tsteven4@users.noreply.github.com> Date: Wed, 10 Nov 2021 19:42:31 -0700 Subject: [PATCH 1/8] Fix kml bounding box, and introduce NVectors. The kml bouning box had problems around the antimeridian, that were unsuccessfully mitigated in 2011. Instead of fixing the obvious mistake in that mitigation we use nvectors to compute the actual center of the bounding box. Nvectors will be useful to avoid problems at discontinuities and singularites that are inherent in the latitude,longitude coordinate system. --- CMakeLists.txt | 4 + GPSBabel.pro | 4 + kml.cc | 22 +- reference/LineStyles.kml | 4 +- reference/bounds-test.kml | 2 +- reference/earth-expertgps-track.kml | 2 +- reference/earth-expertgps.kml | 2 +- reference/earth-gc.kml | 12 +- reference/track/bounds-test-track.kml | 2 +- reference/track/gtrnctr_power-kml.kml | 20 +- reference/track/sim.csv | 49 + reference/track/sim.kml | 1358 +++++++++++++++++++++++++ src/core/nvector.cc | 449 ++++++++ src/core/nvector.h | 94 ++ src/core/vector3d.cc | 171 ++++ src/core/vector3d.h | 66 ++ testo.d/kml.test | 4 + 17 files changed, 2235 insertions(+), 30 deletions(-) create mode 100644 reference/track/sim.csv create mode 100644 reference/track/sim.kml create mode 100644 src/core/nvector.cc create mode 100644 src/core/nvector.h create mode 100644 src/core/vector3d.cc create mode 100644 src/core/vector3d.h diff --git a/CMakeLists.txt b/CMakeLists.txt index 2c1dfbf75..673c7c6fb 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -96,8 +96,10 @@ set(SUPPORT formspec.cc xmltag.cc cet.cc cet_util.cc fatal.cc rgbcolors.cc inifile.cc garmin_fs.cc units.cc gbser.cc gbfile.cc parse.cc session.cc main.cc globals.cc + src/core/nvector.cc src/core/textstream.cc src/core/usasciicodec.cc + src/core/vector3d.cc src/core/xmlstreamwriter.cc ) if(${QT_VERSION_MAJOR} EQUAL "6") @@ -187,8 +189,10 @@ set(HEADERS src/core/datetime.h src/core/file.h src/core/logging.h + src/core/nvector.h src/core/textstream.h src/core/usasciicodec.h + src/core/vector3d.h src/core/xmlstreamwriter.h src/core/xmltag.h ) diff --git a/GPSBabel.pro b/GPSBabel.pro index 0bd2a2696..3ff2771d2 100644 --- a/GPSBabel.pro +++ b/GPSBabel.pro @@ -107,8 +107,10 @@ SUPPORT = route.cc waypt.cc filter_vecs.cc util.cc vecs.cc mkshort.cc \ formspec.cc xmltag.cc cet.cc cet_util.cc fatal.cc rgbcolors.cc \ inifile.cc garmin_fs.cc units.cc gbser.cc \ gbfile.cc parse.cc session.cc main.cc globals.cc \ + src/core/nvector.cc \ src/core/textstream.cc \ src/core/usasciicodec.cc \ + src/core/vector3d.cc \ src/core/xmlstreamwriter.cc versionAtLeast(QT_VERSION, 6.0): SUPPORT += src/core/codecdevice.cc @@ -184,8 +186,10 @@ HEADERS = \ src/core/datetime.h \ src/core/file.h \ src/core/logging.h \ + src/core/nvector.h \ src/core/textstream.h \ src/core/usasciicodec.h \ + src/core/vector3d.h \ src/core/xmlstreamwriter.h \ src/core/xmltag.h diff --git a/kml.cc b/kml.cc index a548e5b85..66652d2f5 100644 --- a/kml.cc +++ b/kml.cc @@ -49,11 +49,13 @@ #include "src/core/datetime.h" // for DateTime #include "src/core/file.h" // for File #include "src/core/logging.h" // for Warning, Fatal +#include "src/core/nvector.h" // for NVector #include "src/core/xmlstreamwriter.h" // for XmlStreamWriter #include "src/core/xmltag.h" // for xml_findfirst, xml_tag, fs_xml, xml_attribute, xml_findnext #include "units.h" // for fmt_setunits, fmt_speed, fmt_altitude, fmt_distance, units_aviation, units_metric, units_nautical, units_statute #include "xmlgeneric.h" // for cb_cdata, cb_end, cb_start, xg_callback, xg_string, xg_cb_type, xml_deinit, xml_ignore_tags, xml_init, xml_read, xg_tag_mapping +using gpsbabel::NVector; // Icons provided and hosted by Google. Used with permission. #define ICON_BASE "https://earth.google.com/images/kml-icons/" @@ -1704,15 +1706,19 @@ void KmlFormat::kml_write_AbstractView() writer->writeEndElement(); // Close gx:TimeSpan tag } -// If our BB spans the antemeridian, flip sign on one. -// This doesn't make our BB optimal, but it at least prevents us from -// zooming to the wrong hemisphere. - if (kml_bounds.min_lon * kml_bounds.max_lon < 0) { - kml_bounds.min_lon = -kml_bounds.max_lon; - } + // Uses NVectors to compute the middle of the bounding box. + // This works even if the BB spans the antimeridian. + // This prevents us from zooming to the wrong hemisphere. + // Gade, K. (2010). A Non-singular Horizontal Position Representation, The Journal of Navigation, Volume 63, Issue 03, pp 395-417, July 2010. + // 5.3.6. Horizontal geographical mean, equation 17. + NVector bb_1(kml_bounds.min_lat, kml_bounds.min_lon); + NVector bb_2(kml_bounds.min_lat, kml_bounds.max_lon); + NVector bb_3(kml_bounds.max_lat, kml_bounds.min_lon); + NVector bb_4(kml_bounds.max_lat, kml_bounds.max_lon); + NVector bb_mid = (bb_1 + bb_2 + bb_3 + bb_4).normalize(); - writer->writeTextElement(QStringLiteral("longitude"), QString::number((kml_bounds.min_lon + kml_bounds.max_lon) / 2, 'f', precision)); - writer->writeTextElement(QStringLiteral("latitude"), QString::number((kml_bounds.min_lat + kml_bounds.max_lat) / 2, 'f', precision)); + writer->writeTextElement(QStringLiteral("longitude"), QString::number(bb_mid.longitude(), 'f', precision)); + writer->writeTextElement(QStringLiteral("latitude"), QString::number(bb_mid.latitude(), 'f', precision)); // It turns out the length of the diagonal of the bounding box gives us a // reasonable guess for setting the camera altitude. diff --git a/reference/LineStyles.kml b/reference/LineStyles.kml index 77c81ca37..8d908bc50 100644 --- a/reference/LineStyles.kml +++ b/reference/LineStyles.kml @@ -4,7 +4,7 @@ GPS device -122.275670 - 37.523278 + 37.523279 5908.813619 @@ -159,7 +159,7 @@ NOTE 1 - 12 test lines - see notes with GPX sample!! four line widths for solid red two line widths for dashed red seven colors for medium width Note that line names have nothing at all to do with line style attributes! These are named for convenient cross-reference, but keep in mind that "red med" is a line label and "Red Medium" is a style label. Waypoints are generated by Topo to exactly match the first track point. Apparently only GPX and KML output from GPSBabel retain track descriptions?? + 12 test lines - see notes with GPX sample!! four line widths for solid red two line widths for dashed red seven colors for medium width Note that line names have nothing at all to do with line style attributes!These are named for convenient cross-reference, but keep in mind that "red med" is a line label and "Red Medium" is a style label. Waypoints are generated by Topo to exactly match the first track point. Apparently only GPX and KML output from GPSBabel retain track descriptions?? #waypoint -122.277670,37.522060 diff --git a/reference/bounds-test.kml b/reference/bounds-test.kml index 65298dd15..8d140e0e4 100644 --- a/reference/bounds-test.kml +++ b/reference/bounds-test.kml @@ -4,7 +4,7 @@ GPS device -117.144015 - 36.438270 + 36.438594 87257.203848 diff --git a/reference/earth-expertgps-track.kml b/reference/earth-expertgps-track.kml index 207b78b3f..ed0f0ad29 100644 --- a/reference/earth-expertgps-track.kml +++ b/reference/earth-expertgps-track.kml @@ -8,7 +8,7 @@ 2002-05-25T19:05:57Z -81.356770 - 36.257086 + 36.698053 2979727.763965 diff --git a/reference/earth-expertgps.kml b/reference/earth-expertgps.kml index fd0b5347a..87ec4654f 100644 --- a/reference/earth-expertgps.kml +++ b/reference/earth-expertgps.kml @@ -8,7 +8,7 @@ 2002-05-25T19:05:57Z -81.356770 - 36.257086 + 36.698053 2979727.763965 diff --git a/reference/earth-gc.kml b/reference/earth-gc.kml index 94d1bfde0..83afb2294 100644 --- a/reference/earth-gc.kml +++ b/reference/earth-gc.kml @@ -8,7 +8,7 @@ 2003-06-29T07:00:00Z -79.930833 - 41.027500 + 41.235699 2109328.437865 @@ -309,7 +309,7 @@ PICTURES of points and of the panels will follow soon. YOU CAN ONLY LOG ONE POIN Good luck!]]> - Found it by Christopher R & Pooh B 2005-07-12
This marker is not in Quebec but it is a Geodesic marker in Clarenville, Newfoundland, Canada! + Found it by Christopher R & Pooh B2005-07-12
This marker is not in Quebec but it is a Geodesic marker in Clarenville, Newfoundland, Canada! Found this one while hunting a traditional cache and thought of this cache right away! @@ -319,7 +319,7 @@ Smiles Pooh Bear Ce marqueur n'est pas au Québec mais c'est un marqueur géodésique dans Clarenville, Terre-Neuve, Canada! -A trouvé celui-ci tandis que chasse une cachette traditionnelle et pensé à cette cachette tout de suite! Elle est située sur la montagne nue dans Clarenville - il y a aactually deux marqueurs à moins de 15 pieds d'un des autres sur la montagne nue... Ours De Pooh De Sourires

Found it by TravelBen 2005-06-26
[:D] 14h22 +A trouvé celui-ci tandis que chasse une cachette traditionnelle et pensé à cette cachette tout de suite! Elle est située sur la montagne nue dans Clarenville - il y a aactually deux marqueurs à moins de 15 pieds d'un des autres sur la montagne nue... Ours De Pooh De Sourires

Found it by TravelBen2005-06-26
[:D] 14h22 Marqueur du Service de la Géodégie (c'est bien un "g" pas un "s") du Québec. @@ -329,7 +329,7 @@ UTM: 18T E 582877 N 5033250 Ce marqueur se trouve dans le ville de Senneville, sur un monument décrivant une page d'histoire du Québec, sur le bas côté avant droit. -Près de la cache: Exo-07 La Jumelle de Loudiver (GCP3VE)

Found it by etasse 2005-06-03
MRN marker 94K4731 in Gatineau, QC. corner of Du Rhone and Gatineau Ave. +Près de la cache: Exo-07 La Jumelle de Loudiver (GCP3VE)

Found it by etasse2005-06-03
MRN marker 94K4731 in Gatineau, QC. corner of Du Rhone and Gatineau Ave. Position Average N 45° 29.5247 W 075° 43.0049 59.49m @@ -341,7 +341,7 @@ UTM 18T 0443996 5037868 This pole has everything: An underground cable warning, a geodesic mark, a bus stop and a garage sale sign. -Judging by the coordinates it looks like the coords should be 45°29'31.5" -75°43'0" I placed the GPS antenna right against the marker, to no avail.

Found it by Katou 2005-06-03
Un bo point géodésique a Lotbinière..en allant faire une nouvelle cache a l'île richelieu ;-)

Found it by Gps_Gulliver&DauphinBleu 2005-05-29
Point Geodesique situe near Port de Plaisance de Longueuil +Judging by the coordinates it looks like the coords should be 45°29'31.5" -75°43'0" I placed the GPS antenna right against the marker, to no avail.

Found it by Katou2005-06-03
Un bo point géodésique a Lotbinière..en allant faire une nouvelle cache a l'île richelieu ;-)

Found it by Gps_Gulliver&DauphinBleu2005-05-29
Point Geodesique situe near Port de Plaisance de Longueuil sur le bord du fleuve st-laurent. Il y a des sentiers et une grande piste cyclable Enjoy !

]]> @@ -426,7 +426,7 @@ Now that it's intuitively obvious to even the most casual observer where the cac Member of Middle Tennessee GeoCachers Club [www.mtgc.org]

]]> - Found it by littlepod 2005-07-03
Enjoyed the puzzle. We seemed to be about 50ft off though. TFTC.

Write note by robertlipe 2005-04-29
TB Drop to show he's hanging out in Nashville until we blast off for Geowoodstock in a few weeks.

Found it by Big Bumblebee 2005-04-18
Found it a while ago. Thanks.

Write note by robertlipe 2005-03-27
I had to renew my permit with the CDC and in doing so, I trolled out here verified that the infectious ooze is fully within specification and industry accepted tolerance. Ooze On!

Found it by Virtual Babe 2004-12-27
This was a great cache, however on this day I considered it a FIFM cache (Fun, Invigorating, Frustrating and Maddening), especially when the cache was not replaced in the proper spot by the previous cacher! Thanks anyway!!

Write note by robertlipe 2004-01-12
I got a complaint from the CDC about oozy rat this weekend. I went out tonight in the dark and verified that the infectious ooze is fully within specification and industry accepted tolerance. (Although I realize now I did misstate the cache container to the reporting officer when confronted. It's, uuuuh, smaller than I said.)

Write note by robertlipe 2003-10-04
In the expectation that this cache will get some traffic in the next 48 hours, Ryan and I checked it earlier today. The Rat is Oozing just as we planned it.

Write note by robertlipe 2003-07-03
It won't earn him a smiley face, but I've confirmed that rickrich would have indeed sunk the battleship! Thanx for playing. You get a copy of the home game and some rice-a-roni...

]]> + Found it by littlepod2005-07-03
Enjoyed the puzzle. We seemed to be about 50ft off though. TFTC.

Write note by robertlipe2005-04-29
TB Drop to show he's hanging out in Nashville until we blast off for Geowoodstock in a few weeks.

Found it by Big Bumblebee2005-04-18
Found it a while ago. Thanks.

Write note by robertlipe2005-03-27
I had to renew my permit with the CDC and in doing so, I trolled out here verified that the infectious ooze is fully within specification and industry accepted tolerance. Ooze On!

Found it by Virtual Babe2004-12-27
This was a great cache, however on this day I considered it a FIFM cache (Fun, Invigorating, Frustrating and Maddening), especially when the cache was not replaced in the proper spot by the previous cacher! Thanks anyway!!

Write note by robertlipe2004-01-12
I got a complaint from the CDC about oozy rat this weekend. I went out tonight in the dark and verified that the infectious ooze is fully within specification and industry accepted tolerance. (Although I realize now I did misstate the cache container to the reporting officer when confronted. It's, uuuuh, smaller than I said.)

Write note by robertlipe2003-10-04
In the expectation that this cache will get some traffic in the next 48 hours, Ryan and I checked it earlier today. The Rat is Oozing just as we planned it.

Write note by robertlipe2003-07-03
It won't earn him a smiley face, but I've confirmed that rickrich would have indeed sunk the battleship! Thanx for playing. You get a copy of the home game and some rice-a-roni...

]]> diff --git a/reference/track/bounds-test-track.kml b/reference/track/bounds-test-track.kml index f0ec3641c..ed0b8cfa1 100644 --- a/reference/track/bounds-test-track.kml +++ b/reference/track/bounds-test-track.kml @@ -4,7 +4,7 @@ GPS device -117.144015 - 36.438270 + 36.438594 87257.203848 diff --git a/reference/track/gtrnctr_power-kml.kml b/reference/track/gtrnctr_power-kml.kml index 10f614624..6667d6337 100644 --- a/reference/track/gtrnctr_power-kml.kml +++ b/reference/track/gtrnctr_power-kml.kml @@ -8,10 +8,10 @@ 2010-05-28T02:41:44Z -122.139608 - 37.382794 + 37.382814 19819.321245 - + - + - + - + - Heart Rate + Heart Rate - Cadence + Cadence - Power + Power diff --git a/reference/track/sim.csv b/reference/track/sim.csv new file mode 100644 index 000000000..8a00c4a00 --- /dev/null +++ b/reference/track/sim.csv @@ -0,0 +1,49 @@ +utc_d,utc_t,lat,lon +2021/11/10,11:00:00 AM,66.000103228632,-179.97904 +2021/11/10,11:00:01 AM,66.0001884601077,-179.97903 +2021/11/10,11:00:02 AM,66.0000230237331,-179.97947 +2021/11/10,11:00:03 AM,66.000045323685,-179.98008 +2021/11/10,11:00:04 AM,66.0001767900345,-179.98048 +2021/11/10,11:00:05 AM,66.0001092928285,-179.98067 +2021/11/10,11:00:06 AM,66.0001668298279,-179.98197 +2021/11/10,11:00:07 AM,66.0000748518391,-179.98320 +2021/11/10,11:00:08 AM,66.0001341555094,-179.98429 +2021/11/10,11:00:09 AM,66.0001213660596,-179.98519 +2021/11/10,11:00:10 AM,66.0001967593428,-179.98558 +2021/11/10,11:00:11 AM,66.0001580625671,-179.98671 +2021/11/10,11:00:12 AM,66.0001378472268,-179.98805 +2021/11/10,11:00:13 AM,66.0001709596168,-179.98906 +2021/11/10,11:00:14 AM,66.0000502908297,-179.98970 +2021/11/10,11:00:15 AM,66.000066383696,-179.99130 +2021/11/10,11:00:16 AM,66.0000462879889,-179.99183 +2021/11/10,11:00:17 AM,66.0001312399228,-179.99340 +2021/11/10,11:00:18 AM,66.0000883525882,-179.99384 +2021/11/10,11:00:19 AM,66.0001647357197,-179.99504 +2021/11/10,11:00:20 AM,66.0001959985441,-179.99596 +2021/11/10,11:00:21 AM,66.0001581122524,-179.99737 +2021/11/10,11:00:22 AM,66.0001028337378,-179.99798 +2021/11/10,11:00:23 AM,66.000105956278,-179.99903 +2021/11/10,11:00:24 AM,66.0000574941283,-179.99956 +2021/11/10,11:00:25 AM,66.0000416369344,179.99915 +2021/11/10,11:00:26 AM,66.0001504925119,179.99831 +2021/11/10,11:00:27 AM,66.0001927654035,179.99675 +2021/11/10,11:00:28 AM,66.0001282065146,179.99621 +2021/11/10,11:00:29 AM,66.0000911889429,179.99493 +2021/11/10,11:00:30 AM,66.0001909788552,179.99362 +2021/11/10,11:00:31 AM,66.0000568057737,179.99293 +2021/11/10,11:00:32 AM,66.0000978383165,179.99211 +2021/11/10,11:00:33 AM,66.0000403599515,179.99138 +2021/11/10,11:00:34 AM,66.0001144849062,179.99045 +2021/11/10,11:00:35 AM,66.0000769781374,179.98854 +2021/11/10,11:00:36 AM,66.0000133158984,179.98783 +2021/11/10,11:00:37 AM,66.0001758676388,179.98699 +2021/11/10,11:00:38 AM,66.0000207460354,179.98606 +2021/11/10,11:00:39 AM,66.0000931231469,179.98538 +2021/11/10,11:00:40 AM,66.0001991309548,179.98429 +2021/11/10,11:00:41 AM,66.0001841979519,179.98311 +2021/11/10,11:00:42 AM,66.0001217081583,179.98251 +2021/11/10,11:00:43 AM,66.0001017396611,179.98182 +2021/11/10,11:00:44 AM,66.0001469324212,179.98183 +2021/11/10,11:00:45 AM,66.0000706021879,179.98175 +2021/11/10,11:00:46 AM,66.0000732617344,179.98178 +2021/11/10,11:00:47 AM,66.0001673249007,179.98162 diff --git a/reference/track/sim.kml b/reference/track/sim.kml new file mode 100644 index 000000000..0b286c82c --- /dev/null +++ b/reference/track/sim.kml @@ -0,0 +1,1358 @@ + + + + GPS device + + + 2021-11-10T11:00:00Z + 2021-11-10T11:00:47Z + + 179.999795 + 66.000106 + 1300.000000 + + + + + + + + normal + #track_n + + + highlight + #track_h + + + + + + + + + normal + #waypoint_n + + + highlight + #waypoint_h + + + + + Tracks + + + + +Distance 1.2 mi +Min Speed 3.1 mph +Max Speed 193.7 mph +Avg Speed 88.4 mph +Start Time2021-11-10T11:00:00Z +End Time2021-11-10T11:00:47Z +]]> + + + 2021-11-10T11:00:00Z + 2021-11-10T11:00:47Z + + + Points + + -0 + + +Longitude: -179.979040 +Latitude: 66.000103 +Heading: 360.0 +Time: 2021-11-10T11:00:00Z + +]]> + + -179.979040 + 66.000103 + 66 + + + 2021-11-10T11:00:00Z + + #track + + -179.979040,66.000103 + + + + -1 + + +Longitude: -179.979030 +Latitude: 66.000188 +Speed: 21.2 mph +Heading: 2.7 +Time: 2021-11-10T11:00:01Z + +]]> + + -179.979030 + 66.000188 + 66 + + + 2021-11-10T11:00:01Z + + #track + + -179.979030,66.000188 + + + + -2 + + +Longitude: -179.979470 +Latitude: 66.000023 +Speed: 60.7 mph +Heading: 227.2 +Time: 2021-11-10T11:00:02Z + +]]> + + -179.979470 + 66.000023 + 66 + + + 2021-11-10T11:00:02Z + + #track + + -179.979470,66.000023 + + + + -3 + + +Longitude: -179.980080 +Latitude: 66.000045 +Speed: 62.0 mph +Heading: 275.1 +Time: 2021-11-10T11:00:03Z + +]]> + + -179.980080 + 66.000045 + 66 + + + 2021-11-10T11:00:03Z + + #track + + -179.980080,66.000045 + + + + -4 + + +Longitude: -179.980480 +Latitude: 66.000177 +Speed: 52.1 mph +Heading: 308.9 +Time: 2021-11-10T11:00:04Z + +]]> + + -179.980480 + 66.000177 + 66 + + + 2021-11-10T11:00:04Z + + #track + + -179.980480,66.000177 + + + + -5 + + +Longitude: -179.980670 +Latitude: 66.000109 +Speed: 25.6 mph +Heading: 228.9 +Time: 2021-11-10T11:00:05Z + +]]> + + -179.980670 + 66.000109 + 66 + + + 2021-11-10T11:00:05Z + + #track + + -179.980670,66.000109 + + + + -6 + + +Longitude: -179.981970 +Latitude: 66.000167 +Speed: 132.4 mph +Heading: 276.2 +Time: 2021-11-10T11:00:06Z + +]]> + + -179.981970 + 66.000167 + 66 + + + 2021-11-10T11:00:06Z + + #track + + -179.981970,66.000167 + + + + -7 + + +Longitude: -179.983200 +Latitude: 66.000075 +Speed: 126.7 mph +Heading: 259.6 +Time: 2021-11-10T11:00:07Z + +]]> + + -179.983200 + 66.000075 + 66 + + + 2021-11-10T11:00:07Z + + #track + + -179.983200,66.000075 + + + + -8 + + +Longitude: -179.984290 +Latitude: 66.000134 +Speed: 111.4 mph +Heading: 277.6 +Time: 2021-11-10T11:00:08Z + +]]> + + -179.984290 + 66.000134 + 66 + + + 2021-11-10T11:00:08Z + + #track + + -179.984290,66.000134 + + + + -9 + + +Longitude: -179.985190 +Latitude: 66.000121 +Speed: 91.2 mph +Heading: 268.0 +Time: 2021-11-10T11:00:09Z + +]]> + + -179.985190 + 66.000121 + 66 + + + 2021-11-10T11:00:09Z + + #track + + -179.985190,66.000121 + + + + -10 + + +Longitude: -179.985580 +Latitude: 66.000197 +Speed: 43.7 mph +Heading: 295.4 +Time: 2021-11-10T11:00:10Z + +]]> + + -179.985580 + 66.000197 + 66 + + + 2021-11-10T11:00:10Z + + #track + + -179.985580,66.000197 + + + + -11 + + +Longitude: -179.986710 +Latitude: 66.000158 +Speed: 114.9 mph +Heading: 265.2 +Time: 2021-11-10T11:00:11Z + +]]> + + -179.986710 + 66.000158 + 66 + + + 2021-11-10T11:00:11Z + + #track + + -179.986710,66.000158 + + + + -12 + + +Longitude: -179.988050 +Latitude: 66.000138 +Speed: 135.8 mph +Heading: 267.9 +Time: 2021-11-10T11:00:12Z + +]]> + + -179.988050 + 66.000138 + 66 + + + 2021-11-10T11:00:12Z + + #track + + -179.988050,66.000138 + + + + -13 + + +Longitude: -179.989060 +Latitude: 66.000171 +Speed: 102.6 mph +Heading: 274.6 +Time: 2021-11-10T11:00:13Z + +]]> + + -179.989060 + 66.000171 + 66 + + + 2021-11-10T11:00:13Z + + #track + + -179.989060,66.000171 + + + + -14 + + +Longitude: -179.989700 +Latitude: 66.000050 +Speed: 71.4 mph +Heading: 245.1 +Time: 2021-11-10T11:00:14Z + +]]> + + -179.989700 + 66.000050 + 66 + + + 2021-11-10T11:00:14Z + + #track + + -179.989700,66.000050 + + + + -15 + + +Longitude: -179.991300 +Latitude: 66.000066 +Speed: 162.1 mph +Heading: 271.4 +Time: 2021-11-10T11:00:15Z + +]]> + + -179.991300 + 66.000066 + 66 + + + 2021-11-10T11:00:15Z + + #track + + -179.991300,66.000066 + + + + -16 + + +Longitude: -179.991830 +Latitude: 66.000046 +Speed: 53.9 mph +Heading: 264.7 +Time: 2021-11-10T11:00:16Z + +]]> + + -179.991830 + 66.000046 + 66 + + + 2021-11-10T11:00:16Z + + #track + + -179.991830,66.000046 + + + + -17 + + +Longitude: -179.993400 +Latitude: 66.000131 +Speed: 160.4 mph +Heading: 277.6 +Time: 2021-11-10T11:00:17Z + +]]> + + -179.993400 + 66.000131 + 66 + + + 2021-11-10T11:00:17Z + + #track + + -179.993400,66.000131 + + + + -18 + + +Longitude: -179.993840 +Latitude: 66.000088 +Speed: 45.8 mph +Heading: 256.5 +Time: 2021-11-10T11:00:18Z + +]]> + + -179.993840 + 66.000088 + 66 + + + 2021-11-10T11:00:18Z + + #track + + -179.993840,66.000088 + + + + -19 + + +Longitude: -179.995040 +Latitude: 66.000165 +Speed: 123.0 mph +Heading: 278.9 +Time: 2021-11-10T11:00:19Z + +]]> + + -179.995040 + 66.000165 + 66 + + + 2021-11-10T11:00:19Z + + #track + + -179.995040,66.000165 + + + + -20 + + +Longitude: -179.995960 +Latitude: 66.000196 +Speed: 93.5 mph +Heading: 274.8 +Time: 2021-11-10T11:00:20Z + +]]> + + -179.995960 + 66.000196 + 66 + + + 2021-11-10T11:00:20Z + + #track + + -179.995960,66.000196 + + + + -21 + + +Longitude: -179.997370 +Latitude: 66.000158 +Speed: 143.1 mph +Heading: 266.2 +Time: 2021-11-10T11:00:21Z + +]]> + + -179.997370 + 66.000158 + 66 + + + 2021-11-10T11:00:21Z + + #track + + -179.997370,66.000158 + + + + -22 + + +Longitude: -179.997980 +Latitude: 66.000103 +Speed: 63.3 mph +Heading: 257.4 +Time: 2021-11-10T11:00:22Z + +]]> + + -179.997980 + 66.000103 + 66 + + + 2021-11-10T11:00:22Z + + #track + + -179.997980,66.000103 + + + + -23 + + +Longitude: -179.999030 +Latitude: 66.000106 +Speed: 106.3 mph +Heading: 270.4 +Time: 2021-11-10T11:00:23Z + +]]> + + -179.999030 + 66.000106 + 66 + + + 2021-11-10T11:00:23Z + + #track + + -179.999030,66.000106 + + + + -24 + + +Longitude: -179.999560 +Latitude: 66.000057 +Speed: 55.0 mph +Heading: 257.3 +Time: 2021-11-10T11:00:24Z + +]]> + + -179.999560 + 66.000057 + 66 + + + 2021-11-10T11:00:24Z + + #track + + -179.999560,66.000057 + + + + -25 + + +Longitude: 179.999150 +Latitude: 66.000042 +Speed: 130.7 mph +Heading: 268.3 +Time: 2021-11-10T11:00:25Z + +]]> + + 179.999150 + 66.000042 + 66 + + + 2021-11-10T11:00:25Z + + #track + + 179.999150,66.000042 + + + + -26 + + +Longitude: 179.998310 +Latitude: 66.000150 +Speed: 89.3 mph +Heading: 287.7 +Time: 2021-11-10T11:00:26Z + +]]> + + 179.998310 + 66.000150 + 66 + + + 2021-11-10T11:00:26Z + + #track + + 179.998310,66.000150 + + + + -27 + + +Longitude: 179.996750 +Latitude: 66.000193 +Speed: 158.4 mph +Heading: 273.8 +Time: 2021-11-10T11:00:27Z + +]]> + + 179.996750 + 66.000193 + 66 + + + 2021-11-10T11:00:27Z + + #track + + 179.996750,66.000193 + + + + -28 + + +Longitude: 179.996210 +Latitude: 66.000128 +Speed: 57.0 mph +Heading: 253.6 +Time: 2021-11-10T11:00:28Z + +]]> + + 179.996210 + 66.000128 + 66 + + + 2021-11-10T11:00:28Z + + #track + + 179.996210,66.000128 + + + + -29 + + +Longitude: 179.994930 +Latitude: 66.000091 +Speed: 130.0 mph +Heading: 265.9 +Time: 2021-11-10T11:00:29Z + +]]> + + 179.994930 + 66.000091 + 66 + + + 2021-11-10T11:00:29Z + + #track + + 179.994930,66.000091 + + + + -30 + + +Longitude: 179.993620 +Latitude: 66.000191 +Speed: 135.0 mph +Heading: 280.6 +Time: 2021-11-10T11:00:30Z + +]]> + + 179.993620 + 66.000191 + 66 + + + 2021-11-10T11:00:30Z + + #track + + 179.993620,66.000191 + + + + -31 + + +Longitude: 179.992930 +Latitude: 66.000057 +Speed: 77.5 mph +Heading: 244.4 +Time: 2021-11-10T11:00:31Z + +]]> + + 179.992930 + 66.000057 + 66 + + + 2021-11-10T11:00:31Z + + #track + + 179.992930,66.000057 + + + + -32 + + +Longitude: 179.992110 +Latitude: 66.000098 +Speed: 83.7 mph +Heading: 277.0 +Time: 2021-11-10T11:00:32Z + +]]> + + 179.992110 + 66.000098 + 66 + + + 2021-11-10T11:00:32Z + + #track + + 179.992110,66.000098 + + + + -33 + + +Longitude: 179.991380 +Latitude: 66.000040 +Speed: 75.3 mph +Heading: 259.0 +Time: 2021-11-10T11:00:33Z + +]]> + + 179.991380 + 66.000040 + 66 + + + 2021-11-10T11:00:33Z + + #track + + 179.991380,66.000040 + + + + -34 + + +Longitude: 179.990450 +Latitude: 66.000114 +Speed: 96.0 mph +Heading: 281.1 +Time: 2021-11-10T11:00:34Z + +]]> + + 179.990450 + 66.000114 + 66 + + + 2021-11-10T11:00:34Z + + #track + + 179.990450,66.000114 + + + + -35 + + +Longitude: 179.988540 +Latitude: 66.000077 +Speed: 193.7 mph +Heading: 267.2 +Time: 2021-11-10T11:00:35Z + +]]> + + 179.988540 + 66.000077 + 66 + + + 2021-11-10T11:00:35Z + + #track + + 179.988540,66.000077 + + + + -36 + + +Longitude: 179.987830 +Latitude: 66.000013 +Speed: 73.6 mph +Heading: 257.6 +Time: 2021-11-10T11:00:36Z + +]]> + + 179.987830 + 66.000013 + 66 + + + 2021-11-10T11:00:36Z + + #track + + 179.987830,66.000013 + + + + -37 + + +Longitude: 179.986990 +Latitude: 66.000176 +Speed: 94.2 mph +Heading: 295.4 +Time: 2021-11-10T11:00:37Z + +]]> + + 179.986990 + 66.000176 + 66 + + + 2021-11-10T11:00:37Z + + #track + + 179.986990,66.000176 + + + + -38 + + +Longitude: 179.986060 +Latitude: 66.000021 +Speed: 101.8 mph +Heading: 247.7 +Time: 2021-11-10T11:00:38Z + +]]> + + 179.986060 + 66.000021 + 66 + + + 2021-11-10T11:00:38Z + + #track + + 179.986060,66.000021 + + + + -39 + + +Longitude: 179.985380 +Latitude: 66.000093 +Speed: 71.2 mph +Heading: 284.7 +Time: 2021-11-10T11:00:39Z + +]]> + + 179.985380 + 66.000093 + 66 + + + 2021-11-10T11:00:39Z + + #track + + 179.985380,66.000093 + + + + -40 + + +Longitude: 179.984290 +Latitude: 66.000199 +Speed: 113.5 mph +Heading: 283.4 +Time: 2021-11-10T11:00:40Z + +]]> + + 179.984290 + 66.000199 + 66 + + + 2021-11-10T11:00:40Z + + #track + + 179.984290,66.000199 + + + + -41 + + +Longitude: 179.983110 +Latitude: 66.000184 +Speed: 119.6 mph +Heading: 268.2 +Time: 2021-11-10T11:00:41Z + +]]> + + 179.983110 + 66.000184 + 66 + + + 2021-11-10T11:00:41Z + + #track + + 179.983110,66.000184 + + + + -42 + + +Longitude: 179.982510 +Latitude: 66.000122 +Speed: 62.7 mph +Heading: 255.6 +Time: 2021-11-10T11:00:42Z + +]]> + + 179.982510 + 66.000122 + 66 + + + 2021-11-10T11:00:42Z + + #track + + 179.982510,66.000122 + + + + -43 + + +Longitude: 179.981820 +Latitude: 66.000102 +Speed: 70.1 mph +Heading: 265.9 +Time: 2021-11-10T11:00:43Z + +]]> + + 179.981820 + 66.000102 + 66 + + + 2021-11-10T11:00:43Z + + #track + + 179.981820,66.000102 + + + + -44 + + +Longitude: 179.981830 +Latitude: 66.000147 +Speed: 11.3 mph +Heading: 5.1 +Time: 2021-11-10T11:00:44Z + +]]> + + 179.981830 + 66.000147 + 66 + + + 2021-11-10T11:00:44Z + + #track + + 179.981830,66.000147 + + + + -45 + + +Longitude: 179.981750 +Latitude: 66.000071 +Speed: 20.7 mph +Heading: 203.1 +Time: 2021-11-10T11:00:45Z + +]]> + + 179.981750 + 66.000071 + 66 + + + 2021-11-10T11:00:45Z + + #track + + 179.981750,66.000071 + + + + -46 + + +Longitude: 179.981780 +Latitude: 66.000073 +Speed: 3.1 mph +Heading: 77.7 +Time: 2021-11-10T11:00:46Z + +]]> + + 179.981780 + 66.000073 + 66 + + + 2021-11-10T11:00:46Z + + #track + + 179.981780,66.000073 + + + + -47 + + +Longitude: 179.981620 +Latitude: 66.000167 +Speed: 28.5 mph +Heading: 325.3 +Time: 2021-11-10T11:00:47Z + +]]> + + 179.981620 + 66.000167 + 66 + + + 2021-11-10T11:00:47Z + + #track + + 179.981620,66.000167 + + + + + Path + #lineStyle + + 1 + +-179.979040,66.000103 +-179.979030,66.000188 +-179.979470,66.000023 +-179.980080,66.000045 +-179.980480,66.000177 +-179.980670,66.000109 +-179.981970,66.000167 +-179.983200,66.000075 +-179.984290,66.000134 +-179.985190,66.000121 +-179.985580,66.000197 +-179.986710,66.000158 +-179.988050,66.000138 +-179.989060,66.000171 +-179.989700,66.000050 +-179.991300,66.000066 +-179.991830,66.000046 +-179.993400,66.000131 +-179.993840,66.000088 +-179.995040,66.000165 +-179.995960,66.000196 +-179.997370,66.000158 +-179.997980,66.000103 +-179.999030,66.000106 +-179.999560,66.000057 +179.999150,66.000042 +179.998310,66.000150 +179.996750,66.000193 +179.996210,66.000128 +179.994930,66.000091 +179.993620,66.000191 +179.992930,66.000057 +179.992110,66.000098 +179.991380,66.000040 +179.990450,66.000114 +179.988540,66.000077 +179.987830,66.000013 +179.986990,66.000176 +179.986060,66.000021 +179.985380,66.000093 +179.984290,66.000199 +179.983110,66.000184 +179.982510,66.000122 +179.981820,66.000102 +179.981830,66.000147 +179.981750,66.000071 +179.981780,66.000073 +179.981620,66.000167 + + + + + + + diff --git a/src/core/nvector.cc b/src/core/nvector.cc new file mode 100644 index 000000000..2a74f89a1 --- /dev/null +++ b/src/core/nvector.cc @@ -0,0 +1,449 @@ +/* + Copyright (C) 2018, 2021 Robert Lipe, gpsbabel.org + + This program is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2 of the License, or + (at your option) any later version. + + This program is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with this program; if not, write to the Free Software + Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. + + */ + +// https://en.wikipedia.org/wiki/N-vector +// http://www.navlab.net/Publications/A_Nonsingular_Horizontal_Position_Representation.pdf + +#include +#include +#include +#include +#include +#include "nvector.h" +#include "vector3d.h" + +namespace gpsbabel +{ + +LatLon::LatLon(double latitude, double longitude) +{ + lat = latitude; + lon = longitude; +} + +NVector::NVector(double latitude_degrees, double longitude_degrees) +{ + // this implements equation 3. + + // The coordinate system choice matches + // A_Nonsingular_Horizontal_Position_Representation.pdf + // The E frame: + // from the earths center + // x points to North pole + // y points towards latitude 0, longitude +90 (east) + // z points to latitude 0, longitude +180, the antimeridian + double latitude_radians = latitude_degrees * M_PI/180.0; + double longitude_radians = longitude_degrees * M_PI/180.0; + x = sin(latitude_radians); + y = sin(longitude_radians)*cos(latitude_radians); + z = -cos(longitude_radians)*cos(latitude_radians); +} + +NVector::NVector(const Vector3D& v) +{ + x = v.getx(); + y = v.gety(); + z = v.getz(); +} + +std::pair PVector::toNVectorAndHeight() const +{ + // this implements equation 23. + constexpr double a = WGS84_SEMI_MAJOR_AXIS_METERS; + constexpr double a2 = a * a; + constexpr double e2 = WGS84_ECCENTRICITY_SQUARED; + constexpr double e4 = e2 * e2; + + double px2 = getx() * getx(); + double py2 = gety() * gety(); + double pz2 = getz() * getz(); + + double q = ((1.0-e2)/(a2)) * px2; + double p = (py2 + pz2) / (a2); + double r = (p+q-e4) / 6.0; + double s = (e4*p*q) / (4.0*r*r*r); + double t = cbrt(1.0 + s + sqrt(s*(2.0+s))); + double u = r * (1.0 + t + (1.0/t)); + double v = sqrt(u*u + e4*q); + double w = e2 * ((u+v-q)/(2.0*v)); + double k = sqrt(u+v+w*w) - w; + double d = k*sqrt(py2 + pz2) / (k+e2); + double sf = 1.0 / (sqrt(d*d + px2)); + double sf2 = k / (k+e2); + + double height = ((k+e2-1.0)/k) * sqrt(d*d+px2); + double nx = sf * getx(); + double ny = sf * sf2 * gety(); + double nz = sf * sf2 * getz(); + + return {Vector3D(nx, ny, nz), height}; +} + +double NVector::latitude() const +{ + // this implements equation 6. + double latitude_radians = atan2(x, sqrt(y*y + z*z)); + double latitude_degrees = latitude_radians * 180.0/M_PI; + return latitude_degrees; +} + +double NVector::longitude() const +{ + // this implements equation 5. + double longitude_radians = atan2(y, -z); + double longitude_degrees = longitude_radians * 180.0/M_PI; + return longitude_degrees; +} + +// great circle distance in radians +double NVector::distance_radians(const NVector& n_EA_E, const NVector& n_EB_E) +{ + // this implements equation 16 using arctan for numerical accuracy. + double result = atan2(crossProduct(n_EA_E, n_EB_E).norm(), + dotProduct(n_EA_E, n_EB_E)); + return result; +} + +// great circle distance in meters +double NVector::distance(const NVector& n_EA_E, const NVector& n_EB_E) +{ + double result = distance_radians(n_EA_E, n_EB_E) * MEAN_EARTH_RADIUS_METERS; + return result; +} + +double NVector::distance(double latitude_a_degrees, double longitude_a_degrees, double latitude_b_degrees, double longitude_b_degrees) +{ + NVector n_EA_E(latitude_a_degrees, longitude_a_degrees); + NVector n_EB_E(latitude_b_degrees, longitude_b_degrees); + double result = distance(n_EA_E, n_EB_E); + return result; +} + +double NVector::azimuth_radians(const NVector& n_EA_E, const NVector& n_EB_E, double height_a, double height_b) +{ + PVector p_EA_E(n_EA_E, height_a); + PVector p_EB_E(n_EB_E, height_b); + Vector3D p_AB_E = p_EB_E - p_EA_E; + + // equation 9 (normalized) + Vector3D kaeast = crossProduct(Vector3D(1.0, 0.0, 0.0), n_EA_E).normalize(); + // equation 10 (normalized, n_EA_E and kaeast are perpendicular unit vectors) + Vector3D kanorth = crossProduct(n_EA_E, kaeast); + + // equation 11 + // Ren = [kanorth(normalized) kaeast(normalized) -n_EA_E]; + // and a rotation from the E frame to the N frame (North, East, Down). + // P_AB_N = R_EN' * p_AB_E + double P_AB_N_x = kanorth*p_AB_E; + double P_AB_N_y = kaeast*p_AB_E; + // double P_AB_N_z = (-n_EA_E)*p_AB_E; + + double azimuth = atan2(P_AB_N_y, P_AB_N_x); + return azimuth; +} + +double NVector::azimuth(const NVector& n_EA_E, const NVector& n_EB_E, double height_a, double height_b) +{ + double azimuth = azimuth_radians(n_EA_E, n_EB_E, height_a, height_b); + double azimuth_degrees = azimuth * 180.0/M_PI; + return azimuth_degrees; +} + +// returns values in the range [0.0,360) +double NVector::azimuth_true_degrees(const NVector& n_EA_E, const NVector& n_EB_E, double height_a, double height_b) +{ + double result = azimuth(n_EA_E, n_EB_E, height_a, height_b); + if (result < 0.0) { + result += 360.0; + } + if (result >= 360.0) { + result = 0.0; + } + return result; +} + +double NVector::azimuth(double latitude_a_degrees, double longitude_a_degrees, + double latitude_b_degrees, double longitude_b_degrees, + double height_a, double height_b) +{ + NVector n_EA_E(latitude_a_degrees, longitude_a_degrees); + NVector n_EB_E(latitude_b_degrees, longitude_b_degrees); + return azimuth(n_EA_E, n_EB_E, height_a, height_b); +} + +// returns values in the range [0.0,360) +double NVector::azimuth_true_degrees(double latitude_a_degrees, double longitude_a_degrees, + double latitude_b_degrees, double longitude_b_degrees, + double height_a, double height_b) +{ + double result = azimuth(latitude_a_degrees, longitude_a_degrees, + latitude_b_degrees, longitude_b_degrees, + height_a, height_b); + if (result < 0.0) { + result += 360.0; + } + if (result >= 360.0) { + result = 0.0; + } + return result; +} + +#if 0 +// This interpolation is nonlinear! +NVector NVector::linepart(const NVector& n_EA_E, const NVector& n_EB_E, double fraction) +{ + Vector3D n_ER_E = (n_EA_E + (n_EB_E - n_EA_E)*fraction).normalize(); + return n_ER_E; +} +#elif 0 +NVector NVector::linepart(const NVector& n_EA_E, const NVector& n_EB_E, double fraction) +{ + // see example 8. + double sab_rad = distance_radians(n_EA_E, n_EB_E); + double sar_rad = fraction*sab_rad; + double az_rad = azimuth_radians(n_EA_E, n_EB_E); + + // equation 9 (normalized) + Vector3D kaeast = crossProduct(Vector3D(1.0, 0.0, 0.0), n_EA_E).normalize(); + // equation 10 (normalized, n_EA_E and kaeast are perpendicular unit vectors) + Vector3D kanorth = crossProduct(n_EA_E, kaeast); + + Vector3D d_EA_E = kanorth*cos(az_rad) + kaeast*sin(az_rad); + Vector3D n_ER_E = n_EA_E*cos(sar_rad) + d_EA_E*sin(sar_rad); + return n_ER_E; +} +#else +NVector NVector::linepart(const NVector& n_EA_E, const NVector& n_EB_E, double fraction) +{ + double dp_a_b = dotProduct(n_EA_E, n_EB_E); + if (dp_a_b >= (1.0-8.0*DBL_EPSILON)) { + // The points are so close we will have trouble constructing a basis. + // Since they are so close the non-linearities in direct n vector + // interpolation are negligble. + Vector3D result = (n_EA_E + (n_EB_E-n_EA_E)*fraction).normalize(); + return result; + } + if (dp_a_b <= (-1.0+8.0*DBL_EPSILON)) { + // The points are so close to be exactly opposite each other that + // there isn't a unique great circle between them. + // Unless fraction is 1.0 or 0.0 there isn't a unique answer. + Vector3D result(nan(""), nan(""), nan("")); + return result; + } + // Form an orthonormal basis with one component perpendicuar to the great + // circle containing A and B, and one component being A. + // Call this the W frame. + // The columns of the rotation matrix from E to W are w1, w2 and w3. + Vector3D w1 = n_EA_E; + Vector3D w2 = crossProduct(n_EB_E, n_EA_E).normalize(); + Vector3D w3 = crossProduct(n_EA_E,w2); + // Rotate A and B to the W frame. + // Vector3D n_EA_W = Vector3D(w1*n_EA_E, w2*n_EA_E, w3*n_EA_E); + Vector3D n_EB_W = Vector3D(w1*n_EB_E, w2*n_EB_E, w3*n_EB_E); + // By construciton n_EA_W.y is (1, 0, 0), + // n_EB_W.y is zero, + // and both n_EA_W and n_EB_W are unit vectors. + // The information is all in the angle between them, + // which is all in n_EB_W.z and n_EB_W.x. + double theta = atan2(n_EB_W.getz(), n_EB_W.getx()); + // We define the interpolated point as the point on the great circle + // between A and B whose great circle distance from A is the given fraction + // of the great circle distance between A and B. For a spheroid the + // distance is proportional to the angle between the vectors. + // The interpolated point is thus: + Vector3D n_EX_W = Vector3D(cos(fraction*theta), 0.0, sin(fraction*theta)); + // Translate the interpolated point back to the E frame. + // We need to invert the matrix composed of the column vectors w1, w2, w3. + // Since this matrix is orthogonal it's inverse equals it's transpose. + Vector3D wt1 = Vector3D(w1.getx(), w2.getx(), w3.getx()); + Vector3D wt2 = Vector3D(w1.gety(), w2.gety(), w3.gety()); + Vector3D wt3 = Vector3D(w1.getz(), w2.getz(), w3.getz()); + Vector3D n_EX_E = Vector3D(wt1*n_EX_W, wt2*n_EX_W, wt3*n_EX_W); + return n_EX_E; +} +#endif + +LatLon NVector::linepart(double latitude_a_degrees, double longitude_a_degrees, + double latitude_b_degrees, double longitude_b_degrees, + double fraction) +{ + NVector n_EA_E(latitude_a_degrees, longitude_a_degrees); + NVector n_EB_E(latitude_b_degrees, longitude_b_degrees); + NVector n_ER_E = linepart(n_EA_E, n_EB_E, fraction); + LatLon result(n_ER_E.latitude(), n_ER_E.longitude()); + return result; +} + +// Find the point of closest approach to Y on the great circle arc AB. +// The great circle arc AB is the shortest of the two possibilities. +#if 0 +// This implementation works and passes regression. +// However, it seems harder to understand than the next implementation. +NVector NVector::crossTrackProjection(const NVector& n_EA_E, const NVector& n_EB_E, const NVector& n_EY_E) +{ + // Compute the normal to the great circle defined by A and B: + Vector3D a_cross_b = crossProduct(n_EA_E, n_EB_E); + // Because a_cross_b is normal to the plane defined by n_EA_E and n_EB_E, + // any great circle defined by the point defined by a_cross_b and + // any other point, e.g. Y, + // will cross the great circle including the points A and B perpendicularly. + // Thus, if we find the intersection of the two great circles we will + // have found the projection of Y onto the great circle defined + // by A and B. + Vector3D x = a_cross_b; + // Compute the normal to the great circle defined by Y and X: + Vector3D x_cross_y = crossProduct(x, n_EY_E); + // A candidate projection from X onto the great circle defined by A and B is + // c_candidate, the other possibity is -c_candidate. + Vector3D c_candidate = crossProduct(a_cross_b, x_cross_y).normalize(); + // pick the candidate closest to Y. + Vector3D c = dotProduct(c_candidate, n_EY_E) >= 0.0 ? c_candidate : -c_candidate; + // find the midpoint between A and B. + Vector3D m = (n_EA_E + n_EB_E).normalize(); + // Now we have a point C on the great circle defined by A and B, + // and the midpoint M of this arc. + // If C is on the arc defined by A and B then the point of closest + // approach to Y is C. But if C is not on the arc, then the point + // of closest approach is either A or B, whichever is close to C. + // Note that A, B, C, M all are all on the great circle defined by A and B, so + // the dot product of any of these two unit vectors monotonically decreases + // the farther apart they are on the great cirle. + Vector3D result; + if (dotProduct(n_EA_E, m) < dotProduct(c, m)) { + result = c; + } else if (dotProduct(n_EA_E, c) < dotProduct(n_EB_E, c)) { + result = n_EB_E; + } else { + result = n_EA_E; + } + return result; +} +#else +NVector NVector::crossTrackProjection(const NVector& n_EA_E, const NVector& n_EB_E, const NVector& n_EY_E) +{ + // Form an orthonormal basis with one component perpendicuar to the great + // circle containing A and B, and one component being A. + // Call this the W frame. + // The columns of the rotation matrix from E to W are w1, w2 and w3. + Vector3D w1 = n_EA_E; + Vector3D w2 = crossProduct(n_EB_E, n_EA_E).normalize(); + Vector3D w3 = crossProduct(n_EA_E,w2); + // Rotate Y to the W frame. + Vector3D n_EY_W = Vector3D(w1*n_EY_E, w2*n_EY_E, w3*n_EY_E); + // By construciton n_EA_W.y is (1, 0, 0), + // n_EB_W.y is zero, + // and both n_EA_W and n_EB_W are unit vectors. + // The projection of Y onto the great cricle defined by A and B + // is just the n_EY_W with the y component set to zero. + Vector3D n_EC_W = Vector3D(n_EY_W.getx(), 0.0, n_EY_W.getz()).normalize(); + // Translate the projected point back to the E frame. + // We need to invert the matrix composed of the column vectors w1, w2, w3. + // Since this matrix is orthogonal it's inverse equals it's transpose. + Vector3D wt1 = Vector3D(w1.getx(), w2.getx(), w3.getx()); + Vector3D wt2 = Vector3D(w1.gety(), w2.gety(), w3.gety()); + Vector3D wt3 = Vector3D(w1.getz(), w2.getz(), w3.getz()); + Vector3D n_EC_E = Vector3D(wt1*n_EC_W, wt2*n_EC_W, wt3*n_EC_W); + // find the midpoint between A and B. + Vector3D n_EM_E = (n_EA_E + n_EB_E).normalize(); + // Now we have a point C on the great circle defined by A and B, + // and the midpoint M of this arc. + // If C is on the arc defined by A and B then the point of closest + // approach to Y is C. But if C is not on the arc, then the point + // of closest approach is either A or B, whichever is close to C. + // Note that A, B, C, M all are all on the great circle defined by A and B, so + // the dot product of any of these two unit vectors monotonically decreases + // the farther apart they are on the great cirle. + Vector3D result; + if (dotProduct(n_EA_E, n_EM_E) < dotProduct(n_EC_E, n_EM_E)) { + result = n_EC_E; + } else if (dotProduct(n_EA_E, n_EC_E) < dotProduct(n_EB_E, n_EC_E)) { + result = n_EB_E; + } else { + result = n_EA_E; + } + return result; +} +#endif + +LatLon NVector::crossTrackProjection(double latitude_a_degrees, double longitude_a_degrees, + double latitude_b_degrees, double longitude_b_degrees, + double latitude_x_degrees, double longitude_x_degrees) +{ + NVector n_EA_E(latitude_a_degrees, longitude_a_degrees); + NVector n_EB_E(latitude_b_degrees, longitude_b_degrees); + NVector n_EX_E(latitude_x_degrees, longitude_x_degrees); + NVector n_EC_E = crossTrackProjection(n_EA_E, n_EB_E, n_EX_E); + LatLon result(NVector(n_EC_E).latitude(), NVector(n_EC_E).longitude()); + return result; +} + +// Find the minimum distance to point X when traveling from A to B. +#if 0 +// This doesn't work if we are more than 90 degrees away +double NVector::crossTrackDistance(const NVector& n_EA_E, const NVector& n_EB_E, const NVector& n_EX_E) +{ + Vector3D c_E = crossProduct(n_EA_E, n_EB_E).normalize(); + double result = fabs((atan2(crossProduct(c_E, n_EX_E).norm(), + dotProduct(c_E, n_EX_E)) - M_PI/2.0)) * MEAN_EARTH_RADIUS_METERS; + return result; +} +#else +double NVector::crossTrackDistance(const NVector& n_EA_E, const NVector& n_EB_E, const NVector& n_EX_E) +{ + NVector n_EP_E = crossTrackProjection(n_EA_E, n_EB_E, n_EX_E); + double result = distance(n_EP_E, n_EX_E); + return result; +} +#endif + +double NVector::crossTrackDistance(double latitude_a_degrees, double longitude_a_degrees, + double latitude_b_degrees, double longitude_b_degrees, + double latitude_x_degrees, double longitude_x_degrees) +{ + NVector n_EA_E(latitude_a_degrees, longitude_a_degrees); + NVector n_EB_E(latitude_b_degrees, longitude_b_degrees); + NVector n_EX_E(latitude_x_degrees, longitude_x_degrees); + double result = crossTrackDistance(n_EA_E, n_EB_E, n_EX_E); + return result; +} + +PVector::PVector(const NVector& n_EA_E, double h=0.0) +{ + // This implements equation 22. + + // a semi-major axis, b semi-minor axis, f flattening + // a/b = 1/(1-f) + // aspect ratio = b/a + constexpr double b = WGS84_SEMI_MINOR_AXIS_METERS; + constexpr double asq_over_bsq = 1.0 / (WGS84_ASPECT_RATIO * WGS84_ASPECT_RATIO); + + double denom = sqrt(n_EA_E.getx()*n_EA_E.getx() + n_EA_E.gety()*n_EA_E.gety()*asq_over_bsq + n_EA_E.getz()*n_EA_E.getz()*asq_over_bsq); + x = (b/denom)*n_EA_E.getx() + h*n_EA_E.getx(); + y = (b/denom)*asq_over_bsq*n_EA_E.gety() + h*n_EA_E.gety(); + z = (b/denom)*asq_over_bsq*n_EA_E.getz() + h*n_EA_E.getz(); +} + +PVector::PVector(const Vector3D& v) +{ + x = v.getx(); + y = v.gety(); + z = v.getz(); +} + +} // namespace gpsbabel diff --git a/src/core/nvector.h b/src/core/nvector.h new file mode 100644 index 000000000..2a175c0f6 --- /dev/null +++ b/src/core/nvector.h @@ -0,0 +1,94 @@ +/* + Copyright (C) 2018, 2021 Robert Lipe, gpsbabel.org + + This program is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2 of the License, or + (at your option) any later version. + + This program is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with this program; if not, write to the Free Software + Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. + + */ + +#ifndef NVECTOR_H +#define NVECTOR_H + +#include "vector3d.h" + +namespace gpsbabel +{ + +#define COMPARE_BEARING_TO_GRTCIRC +#ifdef COMPARE_BEARING_TO_GRTCIRC +constexpr double MEAN_EARTH_RADIUS_METERS = 6378137.0; +#else +constexpr double MEAN_EARTH_RADIUS_METERS = 6371000.0; +#endif +constexpr double WGS84_SEMI_MAJOR_AXIS_METERS = 6378137.0; // a +#ifdef COMPARE_BEARING_TO_GRTCIRC +constexpr double WGS84_FLATTENING = 0.0; +#else +constexpr double WGS84_FLATTENING = 1.0/298.257223563; // (a-b)/a +#endif +constexpr double WGS84_ASPECT_RATIO = 1.0 - WGS84_FLATTENING; // b/a +constexpr double WGS84_SEMI_MINOR_AXIS_METERS = WGS84_SEMI_MAJOR_AXIS_METERS * WGS84_ASPECT_RATIO; // b +constexpr double WGS84_ECCENTRICITY_SQUARED = 1.0 - (WGS84_ASPECT_RATIO * WGS84_ASPECT_RATIO); + +class PVector; + +class LatLon +{ +public: + LatLon(double latitude, double longitude); + + double lat; + double lon; +}; + +class NVector : public Vector3D +{ +// friend class PVector; + +public: + NVector() = default; + NVector(double latitude_degrees, double longitude_degrees); + NVector(const Vector3D& v); + + [[nodiscard]] double latitude() const; + [[nodiscard]] double longitude() const; + + static double distance_radians(const NVector& n_EA_E, const NVector& n_EB_E); + static double distance(const NVector& n_EA_E, const NVector& n_EB_E); + static double distance(double latitude_a_degrees, double longitude_a_degrees, double latitude_b_degrees, double longitude_b_degrees); + static double azimuth_radians(const NVector& n_EA_E, const NVector& n_EB_E, double height_a = 0.0, double height_b = 0.0); + static double azimuth(const NVector& n_EA_E, const NVector& n_EB_E, double height_a = 0.0, double height_b = 0.0); + static double azimuth(double latitude_a_degrees, double longitude_a_degrees, double latitude_b_degrees, double longitude_b_degrees, double height_a = 0.0, double height_b = 0.0); + static double azimuth_true_degrees(const NVector& n_EA_E, const NVector& n_EB_E, double height_a = 0.0, double height_b = 0.0); + static double azimuth_true_degrees(double latitude_a_degrees, double longitude_a_degrees, double latitude_b_degrees, double longitude_b_degrees, double height_a = 0.0, double height_b = 0.0); + static NVector linepart(const NVector& n_EA_E, const NVector& n_EB_E, double fraction); + static LatLon linepart(double latitude_a_degrees, double longitude_a_degrees, double latitude_b_degrees, double longitude_b_degrees, double fraction); + static NVector crossTrackProjection(const NVector& n_EA_E, const NVector& n_EB_E, const NVector& n_EY_E); + static LatLon crossTrackProjection(double latitude_a_degrees, double longitude_a_degrees, double latitude_b_degrees, double longitude_b_degrees, double latitude_x_degrees, double longitude_x_degrees); + static double crossTrackDistance(const NVector& n_EA_E, const NVector& n_EB_E, const NVector& n_EX_E); + static double crossTrackDistance(double latitude_a_degrees, double longitude_a_degrees, double latitude_b_degrees, double longitude_b_degrees, double latitude_x_degrees, double longitude_x_degrees); +}; + +class PVector : public Vector3D +{ +public: + PVector() = default; + PVector(const NVector& n_EA_E, double h); + PVector(const Vector3D& v); + + [[nodiscard]] std::pair toNVectorAndHeight() const; +}; + +} // namespace gpsbabel +#endif // NVECTOR_H diff --git a/src/core/vector3d.cc b/src/core/vector3d.cc new file mode 100644 index 000000000..2328d7600 --- /dev/null +++ b/src/core/vector3d.cc @@ -0,0 +1,171 @@ +/* + Copyright (C) 2018, 2021 Robert Lipe, gpsbabel.org + + This program is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2 of the License, or + (at your option) any later version. + + This program is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with this program; if not, write to the Free Software + Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. + + */ + +// https://en.wikipedia.org/wiki/N-vector +// http://www.navlab.net/Publications/A_Nonsingular_Horizontal_Position_Representation.pdf + +#include +#include +#include "vector3d.h" + +namespace gpsbabel +{ + +Vector3D::Vector3D(double xi, double yi, double zi) +{ + x = xi; + y = yi; + z = zi; +} + +double Vector3D::norm() const +{ + double norm = sqrt(x*x + y*y + z*z); + return norm; +} + +double Vector3D::getx() const +{ + return x; +} + +double Vector3D::gety() const +{ + return y; +} + +double Vector3D::getz() const +{ + return z; +} + +Vector3D& Vector3D::normalize() +{ + *this = *this/this->norm(); + return *this; +} + +Vector3D& Vector3D::operator+=(const Vector3D& rhs) +{ + x += rhs.x; + y += rhs.y; + z += rhs.z; + return *this; +} + +Vector3D& Vector3D::operator-=(const Vector3D& rhs) +{ + x -= rhs.x; + y -= rhs.y; + z -= rhs.z; + return *this; +} + +Vector3D& Vector3D::operator*=(double rhs) +{ + x *= rhs; + y *= rhs; + z *= rhs; + return *this; +} + +Vector3D& Vector3D::operator/=(double rhs) +{ + x /= rhs; + y /= rhs; + z /= rhs; + return *this; +} + +Vector3D Vector3D::operator+ (const Vector3D& rhs) const +{ + Vector3D result = *this; + result += rhs; + return result; +} + +Vector3D Vector3D::operator- (const Vector3D& rhs) const +{ + Vector3D result = *this; + result -= rhs; + return result; +} + +Vector3D Vector3D::operator* (double rhs) const +{ + Vector3D result = *this; + result *= rhs; + return result; +} + +Vector3D Vector3D::operator/ (double rhs) const +{ + Vector3D result = *this; + result /= rhs; + return result; +} + +Vector3D Vector3D::operator- () const +{ + Vector3D result(-x,-y,-z); + return result; +} + +double Vector3D::operator* (const Vector3D& b) const +{ + double result = x*b.x + y*b.y + z*b.z; + return result; +} + +double Vector3D::dotProduct(const Vector3D& a, const Vector3D& b) +{ + double dotproduct = a.x*b.x + a.y*b.y + a.z*b.z; + return dotproduct; +} + +Vector3D Vector3D::crossProduct(const Vector3D& b, const Vector3D& c) +{ + // a = b x c + Vector3D a; + a.x = b.y*c.z - b.z*c.y; + a.y = b.z*c.x - b.x*c.z; + a.z = b.x*c.y - b.y*c.x; + return a; +} + +Vector3D operator* (double lhs, const Vector3D& rhs) +{ + Vector3D result = rhs*lhs; + return result; +} + +std::ostream& operator<< (std::ostream& os, const Vector3D& v) +{ + os << '(' << v.getx() << ", " << v.gety() << ", " << v.getz() << ')'; + return os; +} + +QDebug operator<< (QDebug debug, const Vector3D& v) +{ + QDebugStateSaver saver(debug); + debug.nospace() << '(' << v.getx() << ", " << v.gety() << ", " << v.getz() << ')'; + return debug; +} + +} // namespace gpsbabel diff --git a/src/core/vector3d.h b/src/core/vector3d.h new file mode 100644 index 000000000..aa5eb05f4 --- /dev/null +++ b/src/core/vector3d.h @@ -0,0 +1,66 @@ +/* + Copyright (C) 2018, 2021 Robert Lipe, gpsbabel.org + + This program is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2 of the License, or + (at your option) any later version. + + This program is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with this program; if not, write to the Free Software + Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. + + */ + +#ifndef VECTOR3D_H +#define VECTOR3D_H + +#include + +#include + +namespace gpsbabel +{ + +class Vector3D +{ +public: + Vector3D() = default; + Vector3D(double xi, double yi, double zi); + + [[nodiscard]] double norm() const; + [[nodiscard]] double getx() const; + [[nodiscard]] double gety() const; + [[nodiscard]] double getz() const; + Vector3D& normalize(); + + Vector3D& operator+=(const Vector3D& rhs); + Vector3D& operator-=(const Vector3D& rhs); + Vector3D& operator*=(double rhs); + Vector3D& operator/=(double rhs); + Vector3D operator+ (const Vector3D& rhs) const; + Vector3D operator- (const Vector3D& rhs) const; + Vector3D operator* (double rhs) const; + Vector3D operator/ (double rhs) const; + Vector3D operator- () const; + double operator* (const Vector3D& b) const; + + static Vector3D crossProduct(const Vector3D& b, const Vector3D& c); + static double dotProduct(const Vector3D& a, const Vector3D& b); + +protected: + double x{}; + double y{}; + double z{}; +}; +Vector3D operator* (double lhs, const Vector3D& rhs); +std::ostream& operator<< (std::ostream& os, const Vector3D& v); +QDebug operator<<(QDebug debug, const Vector3D& v); + +} // namespace gpsbabel +#endif // VECTOR3D_H diff --git a/testo.d/kml.test b/testo.d/kml.test index 2863b86bd..2b3664e74 100644 --- a/testo.d/kml.test +++ b/testo.d/kml.test @@ -69,6 +69,10 @@ compare ${REFERENCE}/track/Placemark-Track-1~kml.gpx ${TMPDIR}/Placemark-Track-1 gpsbabel -T -i random,points=20,seed=33,nodelay -f dummy -o kml,track -F ${TMPDIR}/realtime.kml compare ${REFERENCE}/realtime.kml ${TMPDIR}/realtime.kml +# kml bounding box +gpsbabel -t -i unicsv -f ${REFERENCE}/track/sim.csv -o kml -F ${TMPDIR}/sim.kml +compare ${REFERENCE}/track/sim.kml ${TMPDIR}/sim.kml + if [ "${RUNNINGVALGRIND}" != "0" ]; then set -e if which xmllint > /dev/null; From 6a75958a72972c1bd525731e55991bb14b1ca1b2 Mon Sep 17 00:00:00 2001 From: tsteven4 <13596209+tsteven4@users.noreply.github.com> Date: Wed, 10 Nov 2021 20:26:38 -0700 Subject: [PATCH 2/8] restore lost whitespace in kml gc_logs data. --- kml.cc | 2 +- reference/earth-gc.kml | 10 +++++----- 2 files changed, 6 insertions(+), 6 deletions(-) diff --git a/kml.cc b/kml.cc index 66652d2f5..00ad14247 100644 --- a/kml.cc +++ b/kml.cc @@ -1239,7 +1239,7 @@ QString KmlFormat::kml_geocache_get_logs(const Waypoint* wpt) const if (logpart) { gpsbabel::DateTime t = xml_parse_time(logpart->cdata); if (t.isValid()) { - r += t.date().toString(Qt::ISODate); + r = r + " " + t.date().toString(Qt::ISODate); } } diff --git a/reference/earth-gc.kml b/reference/earth-gc.kml index 83afb2294..31a9144b8 100644 --- a/reference/earth-gc.kml +++ b/reference/earth-gc.kml @@ -309,7 +309,7 @@ PICTURES of points and of the panels will follow soon. YOU CAN ONLY LOG ONE POIN Good luck!]]> - Found it by Christopher R & Pooh B2005-07-12
This marker is not in Quebec but it is a Geodesic marker in Clarenville, Newfoundland, Canada! + Found it by Christopher R & Pooh B 2005-07-12
This marker is not in Quebec but it is a Geodesic marker in Clarenville, Newfoundland, Canada! Found this one while hunting a traditional cache and thought of this cache right away! @@ -319,7 +319,7 @@ Smiles Pooh Bear Ce marqueur n'est pas au Québec mais c'est un marqueur géodésique dans Clarenville, Terre-Neuve, Canada! -A trouvé celui-ci tandis que chasse une cachette traditionnelle et pensé à cette cachette tout de suite! Elle est située sur la montagne nue dans Clarenville - il y a aactually deux marqueurs à moins de 15 pieds d'un des autres sur la montagne nue... Ours De Pooh De Sourires

Found it by TravelBen2005-06-26
[:D] 14h22 +A trouvé celui-ci tandis que chasse une cachette traditionnelle et pensé à cette cachette tout de suite! Elle est située sur la montagne nue dans Clarenville - il y a aactually deux marqueurs à moins de 15 pieds d'un des autres sur la montagne nue... Ours De Pooh De Sourires

Found it by TravelBen 2005-06-26
[:D] 14h22 Marqueur du Service de la Géodégie (c'est bien un "g" pas un "s") du Québec. @@ -329,7 +329,7 @@ UTM: 18T E 582877 N 5033250 Ce marqueur se trouve dans le ville de Senneville, sur un monument décrivant une page d'histoire du Québec, sur le bas côté avant droit. -Près de la cache: Exo-07 La Jumelle de Loudiver (GCP3VE)

Found it by etasse2005-06-03
MRN marker 94K4731 in Gatineau, QC. corner of Du Rhone and Gatineau Ave. +Près de la cache: Exo-07 La Jumelle de Loudiver (GCP3VE)

Found it by etasse 2005-06-03
MRN marker 94K4731 in Gatineau, QC. corner of Du Rhone and Gatineau Ave. Position Average N 45° 29.5247 W 075° 43.0049 59.49m @@ -341,7 +341,7 @@ UTM 18T 0443996 5037868 This pole has everything: An underground cable warning, a geodesic mark, a bus stop and a garage sale sign. -Judging by the coordinates it looks like the coords should be 45°29'31.5" -75°43'0" I placed the GPS antenna right against the marker, to no avail.

Found it by Katou2005-06-03
Un bo point géodésique a Lotbinière..en allant faire une nouvelle cache a l'île richelieu ;-)

Found it by Gps_Gulliver&DauphinBleu2005-05-29
Point Geodesique situe near Port de Plaisance de Longueuil +Judging by the coordinates it looks like the coords should be 45°29'31.5" -75°43'0" I placed the GPS antenna right against the marker, to no avail.

Found it by Katou 2005-06-03
Un bo point géodésique a Lotbinière..en allant faire une nouvelle cache a l'île richelieu ;-)

Found it by Gps_Gulliver&DauphinBleu 2005-05-29
Point Geodesique situe near Port de Plaisance de Longueuil sur le bord du fleuve st-laurent. Il y a des sentiers et une grande piste cyclable Enjoy !

]]> @@ -426,7 +426,7 @@ Now that it's intuitively obvious to even the most casual observer where the cac Member of Middle Tennessee GeoCachers Club [www.mtgc.org]

]]> - Found it by littlepod2005-07-03
Enjoyed the puzzle. We seemed to be about 50ft off though. TFTC.

Write note by robertlipe2005-04-29
TB Drop to show he's hanging out in Nashville until we blast off for Geowoodstock in a few weeks.

Found it by Big Bumblebee2005-04-18
Found it a while ago. Thanks.

Write note by robertlipe2005-03-27
I had to renew my permit with the CDC and in doing so, I trolled out here verified that the infectious ooze is fully within specification and industry accepted tolerance. Ooze On!

Found it by Virtual Babe2004-12-27
This was a great cache, however on this day I considered it a FIFM cache (Fun, Invigorating, Frustrating and Maddening), especially when the cache was not replaced in the proper spot by the previous cacher! Thanks anyway!!

Write note by robertlipe2004-01-12
I got a complaint from the CDC about oozy rat this weekend. I went out tonight in the dark and verified that the infectious ooze is fully within specification and industry accepted tolerance. (Although I realize now I did misstate the cache container to the reporting officer when confronted. It's, uuuuh, smaller than I said.)

Write note by robertlipe2003-10-04
In the expectation that this cache will get some traffic in the next 48 hours, Ryan and I checked it earlier today. The Rat is Oozing just as we planned it.

Write note by robertlipe2003-07-03
It won't earn him a smiley face, but I've confirmed that rickrich would have indeed sunk the battleship! Thanx for playing. You get a copy of the home game and some rice-a-roni...

]]> + Found it by littlepod 2005-07-03
Enjoyed the puzzle. We seemed to be about 50ft off though. TFTC.

Write note by robertlipe 2005-04-29
TB Drop to show he's hanging out in Nashville until we blast off for Geowoodstock in a few weeks.

Found it by Big Bumblebee 2005-04-18
Found it a while ago. Thanks.

Write note by robertlipe 2005-03-27
I had to renew my permit with the CDC and in doing so, I trolled out here verified that the infectious ooze is fully within specification and industry accepted tolerance. Ooze On!

Found it by Virtual Babe 2004-12-27
This was a great cache, however on this day I considered it a FIFM cache (Fun, Invigorating, Frustrating and Maddening), especially when the cache was not replaced in the proper spot by the previous cacher! Thanks anyway!!

Write note by robertlipe 2004-01-12
I got a complaint from the CDC about oozy rat this weekend. I went out tonight in the dark and verified that the infectious ooze is fully within specification and industry accepted tolerance. (Although I realize now I did misstate the cache container to the reporting officer when confronted. It's, uuuuh, smaller than I said.)

Write note by robertlipe 2003-10-04
In the expectation that this cache will get some traffic in the next 48 hours, Ryan and I checked it earlier today. The Rat is Oozing just as we planned it.

Write note by robertlipe 2003-07-03
It won't earn him a smiley face, but I've confirmed that rickrich would have indeed sunk the battleship! Thanx for playing. You get a copy of the home game and some rice-a-roni...

]]> From f6e24bae1f3ba8942afc39dc3b8de1d9a9cdf2fd Mon Sep 17 00:00:00 2001 From: tsteven4 <13596209+tsteven4@users.noreply.github.com> Date: Wed, 10 Nov 2021 20:41:22 -0700 Subject: [PATCH 3/8] get nvectors M_PI on MSVC. --- src/core/nvector.cc | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/core/nvector.cc b/src/core/nvector.cc index 2a74f89a1..1bb94cc95 100644 --- a/src/core/nvector.cc +++ b/src/core/nvector.cc @@ -25,6 +25,8 @@ #include #include #include + +#include "defs.h" #include "nvector.h" #include "vector3d.h" From 106b22b74f767e688b09a73627340609b2227039 Mon Sep 17 00:00:00 2001 From: tsteven4 <13596209+tsteven4@users.noreply.github.com> Date: Thu, 11 Nov 2021 07:14:26 -0700 Subject: [PATCH 4/8] improve spelling and capitalization in nvector.cc --- src/core/nvector.cc | 28 ++++++++++++++-------------- 1 file changed, 14 insertions(+), 14 deletions(-) diff --git a/src/core/nvector.cc b/src/core/nvector.cc index 1bb94cc95..2dd425f15 100644 --- a/src/core/nvector.cc +++ b/src/core/nvector.cc @@ -41,7 +41,7 @@ LatLon::LatLon(double latitude, double longitude) NVector::NVector(double latitude_degrees, double longitude_degrees) { - // this implements equation 3. + // This implements equation 3. // The coordinate system choice matches // A_Nonsingular_Horizontal_Position_Representation.pdf @@ -66,7 +66,7 @@ NVector::NVector(const Vector3D& v) std::pair PVector::toNVectorAndHeight() const { - // this implements equation 23. + // This implements equation 23. constexpr double a = WGS84_SEMI_MAJOR_AXIS_METERS; constexpr double a2 = a * a; constexpr double e2 = WGS84_ECCENTRICITY_SQUARED; @@ -99,7 +99,7 @@ std::pair PVector::toNVectorAndHeight() const double NVector::latitude() const { - // this implements equation 6. + // This implements equation 6. double latitude_radians = atan2(x, sqrt(y*y + z*z)); double latitude_degrees = latitude_radians * 180.0/M_PI; return latitude_degrees; @@ -107,7 +107,7 @@ double NVector::latitude() const double NVector::longitude() const { - // this implements equation 5. + // This implements equation 5. double longitude_radians = atan2(y, -z); double longitude_degrees = longitude_radians * 180.0/M_PI; return longitude_degrees; @@ -116,7 +116,7 @@ double NVector::longitude() const // great circle distance in radians double NVector::distance_radians(const NVector& n_EA_E, const NVector& n_EB_E) { - // this implements equation 16 using arctan for numerical accuracy. + // This implements equation 16 using arctan for numerical accuracy. double result = atan2(crossProduct(n_EA_E, n_EB_E).norm(), dotProduct(n_EA_E, n_EB_E)); return result; @@ -237,7 +237,7 @@ NVector NVector::linepart(const NVector& n_EA_E, const NVector& n_EB_E, double f if (dp_a_b >= (1.0-8.0*DBL_EPSILON)) { // The points are so close we will have trouble constructing a basis. // Since they are so close the non-linearities in direct n vector - // interpolation are negligble. + // interpolation are negligible. Vector3D result = (n_EA_E + (n_EB_E-n_EA_E)*fraction).normalize(); return result; } @@ -248,7 +248,7 @@ NVector NVector::linepart(const NVector& n_EA_E, const NVector& n_EB_E, double f Vector3D result(nan(""), nan(""), nan("")); return result; } - // Form an orthonormal basis with one component perpendicuar to the great + // Form an orthonormal basis with one component perpendicular to the great // circle containing A and B, and one component being A. // Call this the W frame. // The columns of the rotation matrix from E to W are w1, w2 and w3. @@ -258,7 +258,7 @@ NVector NVector::linepart(const NVector& n_EA_E, const NVector& n_EB_E, double f // Rotate A and B to the W frame. // Vector3D n_EA_W = Vector3D(w1*n_EA_E, w2*n_EA_E, w3*n_EA_E); Vector3D n_EB_W = Vector3D(w1*n_EB_E, w2*n_EB_E, w3*n_EB_E); - // By construciton n_EA_W.y is (1, 0, 0), + // By construction n_EA_W.y is (1, 0, 0), // n_EB_W.y is zero, // and both n_EA_W and n_EB_W are unit vectors. // The information is all in the angle between them, @@ -312,7 +312,7 @@ NVector NVector::crossTrackProjection(const NVector& n_EA_E, const NVector& n_EB // Compute the normal to the great circle defined by Y and X: Vector3D x_cross_y = crossProduct(x, n_EY_E); // A candidate projection from X onto the great circle defined by A and B is - // c_candidate, the other possibity is -c_candidate. + // c_candidate, the other possibility is -c_candidate. Vector3D c_candidate = crossProduct(a_cross_b, x_cross_y).normalize(); // pick the candidate closest to Y. Vector3D c = dotProduct(c_candidate, n_EY_E) >= 0.0 ? c_candidate : -c_candidate; @@ -325,7 +325,7 @@ NVector NVector::crossTrackProjection(const NVector& n_EA_E, const NVector& n_EB // of closest approach is either A or B, whichever is close to C. // Note that A, B, C, M all are all on the great circle defined by A and B, so // the dot product of any of these two unit vectors monotonically decreases - // the farther apart they are on the great cirle. + // the farther apart they are on the great circle. Vector3D result; if (dotProduct(n_EA_E, m) < dotProduct(c, m)) { result = c; @@ -339,7 +339,7 @@ NVector NVector::crossTrackProjection(const NVector& n_EA_E, const NVector& n_EB #else NVector NVector::crossTrackProjection(const NVector& n_EA_E, const NVector& n_EB_E, const NVector& n_EY_E) { - // Form an orthonormal basis with one component perpendicuar to the great + // Form an orthonormal basis with one component perpendicular to the great // circle containing A and B, and one component being A. // Call this the W frame. // The columns of the rotation matrix from E to W are w1, w2 and w3. @@ -348,10 +348,10 @@ NVector NVector::crossTrackProjection(const NVector& n_EA_E, const NVector& n_EB Vector3D w3 = crossProduct(n_EA_E,w2); // Rotate Y to the W frame. Vector3D n_EY_W = Vector3D(w1*n_EY_E, w2*n_EY_E, w3*n_EY_E); - // By construciton n_EA_W.y is (1, 0, 0), + // By construction n_EA_W.y is (1, 0, 0), // n_EB_W.y is zero, // and both n_EA_W and n_EB_W are unit vectors. - // The projection of Y onto the great cricle defined by A and B + // The projection of Y onto the great circle defined by A and B // is just the n_EY_W with the y component set to zero. Vector3D n_EC_W = Vector3D(n_EY_W.getx(), 0.0, n_EY_W.getz()).normalize(); // Translate the projected point back to the E frame. @@ -370,7 +370,7 @@ NVector NVector::crossTrackProjection(const NVector& n_EA_E, const NVector& n_EB // of closest approach is either A or B, whichever is close to C. // Note that A, B, C, M all are all on the great circle defined by A and B, so // the dot product of any of these two unit vectors monotonically decreases - // the farther apart they are on the great cirle. + // the farther apart they are on the great circle. Vector3D result; if (dotProduct(n_EA_E, n_EM_E) < dotProduct(n_EC_E, n_EM_E)) { result = n_EC_E; From ad0dc3504b41cb0954d40b591f9726341fd9b041 Mon Sep 17 00:00:00 2001 From: tsteven4 <13596209+tsteven4@users.noreply.github.com> Date: Thu, 11 Nov 2021 07:28:57 -0700 Subject: [PATCH 5/8] rename vector3d member variables. --- src/core/nvector.cc | 28 ++++++++++++------------- src/core/vector3d.cc | 50 ++++++++++++++++++++++---------------------- src/core/vector3d.h | 6 +++--- 3 files changed, 42 insertions(+), 42 deletions(-) diff --git a/src/core/nvector.cc b/src/core/nvector.cc index 2dd425f15..1793a49b8 100644 --- a/src/core/nvector.cc +++ b/src/core/nvector.cc @@ -52,16 +52,16 @@ NVector::NVector(double latitude_degrees, double longitude_degrees) // z points to latitude 0, longitude +180, the antimeridian double latitude_radians = latitude_degrees * M_PI/180.0; double longitude_radians = longitude_degrees * M_PI/180.0; - x = sin(latitude_radians); - y = sin(longitude_radians)*cos(latitude_radians); - z = -cos(longitude_radians)*cos(latitude_radians); + x_ = sin(latitude_radians); + y_ = sin(longitude_radians)*cos(latitude_radians); + z_ = -cos(longitude_radians)*cos(latitude_radians); } NVector::NVector(const Vector3D& v) { - x = v.getx(); - y = v.gety(); - z = v.getz(); + x_ = v.getx(); + y_ = v.gety(); + z_ = v.getz(); } std::pair PVector::toNVectorAndHeight() const @@ -100,7 +100,7 @@ std::pair PVector::toNVectorAndHeight() const double NVector::latitude() const { // This implements equation 6. - double latitude_radians = atan2(x, sqrt(y*y + z*z)); + double latitude_radians = atan2(x_, sqrt(y_*y_ + z_*z_)); double latitude_degrees = latitude_radians * 180.0/M_PI; return latitude_degrees; } @@ -108,7 +108,7 @@ double NVector::latitude() const double NVector::longitude() const { // This implements equation 5. - double longitude_radians = atan2(y, -z); + double longitude_radians = atan2(y_, -z_); double longitude_degrees = longitude_radians * 180.0/M_PI; return longitude_degrees; } @@ -436,16 +436,16 @@ PVector::PVector(const NVector& n_EA_E, double h=0.0) constexpr double asq_over_bsq = 1.0 / (WGS84_ASPECT_RATIO * WGS84_ASPECT_RATIO); double denom = sqrt(n_EA_E.getx()*n_EA_E.getx() + n_EA_E.gety()*n_EA_E.gety()*asq_over_bsq + n_EA_E.getz()*n_EA_E.getz()*asq_over_bsq); - x = (b/denom)*n_EA_E.getx() + h*n_EA_E.getx(); - y = (b/denom)*asq_over_bsq*n_EA_E.gety() + h*n_EA_E.gety(); - z = (b/denom)*asq_over_bsq*n_EA_E.getz() + h*n_EA_E.getz(); + x_ = (b/denom)*n_EA_E.getx() + h*n_EA_E.getx(); + y_ = (b/denom)*asq_over_bsq*n_EA_E.gety() + h*n_EA_E.gety(); + z_ = (b/denom)*asq_over_bsq*n_EA_E.getz() + h*n_EA_E.getz(); } PVector::PVector(const Vector3D& v) { - x = v.getx(); - y = v.gety(); - z = v.getz(); + x_ = v.getx(); + y_ = v.gety(); + z_ = v.getz(); } } // namespace gpsbabel diff --git a/src/core/vector3d.cc b/src/core/vector3d.cc index 2328d7600..c2c77bca4 100644 --- a/src/core/vector3d.cc +++ b/src/core/vector3d.cc @@ -29,30 +29,30 @@ namespace gpsbabel Vector3D::Vector3D(double xi, double yi, double zi) { - x = xi; - y = yi; - z = zi; + x_ = xi; + y_ = yi; + z_ = zi; } double Vector3D::norm() const { - double norm = sqrt(x*x + y*y + z*z); + double norm = sqrt(x_*x_ + y_*y_ + z_*z_); return norm; } double Vector3D::getx() const { - return x; + return x_; } double Vector3D::gety() const { - return y; + return y_; } double Vector3D::getz() const { - return z; + return z_; } Vector3D& Vector3D::normalize() @@ -63,33 +63,33 @@ Vector3D& Vector3D::normalize() Vector3D& Vector3D::operator+=(const Vector3D& rhs) { - x += rhs.x; - y += rhs.y; - z += rhs.z; + x_ += rhs.x_; + y_ += rhs.y_; + z_ += rhs.z_; return *this; } Vector3D& Vector3D::operator-=(const Vector3D& rhs) { - x -= rhs.x; - y -= rhs.y; - z -= rhs.z; + x_ -= rhs.x_; + y_ -= rhs.y_; + z_ -= rhs.z_; return *this; } Vector3D& Vector3D::operator*=(double rhs) { - x *= rhs; - y *= rhs; - z *= rhs; + x_ *= rhs; + y_ *= rhs; + z_ *= rhs; return *this; } Vector3D& Vector3D::operator/=(double rhs) { - x /= rhs; - y /= rhs; - z /= rhs; + x_ /= rhs; + y_ /= rhs; + z_ /= rhs; return *this; } @@ -123,19 +123,19 @@ Vector3D Vector3D::operator/ (double rhs) const Vector3D Vector3D::operator- () const { - Vector3D result(-x,-y,-z); + Vector3D result(-x_,-y_,-z_); return result; } double Vector3D::operator* (const Vector3D& b) const { - double result = x*b.x + y*b.y + z*b.z; + double result = x_*b.x_ + y_*b.y_ + z_*b.z_; return result; } double Vector3D::dotProduct(const Vector3D& a, const Vector3D& b) { - double dotproduct = a.x*b.x + a.y*b.y + a.z*b.z; + double dotproduct = a.x_*b.x_ + a.y_*b.y_ + a.z_*b.z_; return dotproduct; } @@ -143,9 +143,9 @@ Vector3D Vector3D::crossProduct(const Vector3D& b, const Vector3D& c) { // a = b x c Vector3D a; - a.x = b.y*c.z - b.z*c.y; - a.y = b.z*c.x - b.x*c.z; - a.z = b.x*c.y - b.y*c.x; + a.x_ = b.y_*c.z_ - b.z_*c.y_; + a.y_ = b.z_*c.x_ - b.x_*c.z_; + a.z_ = b.x_*c.y_ - b.y_*c.x_; return a; } diff --git a/src/core/vector3d.h b/src/core/vector3d.h index aa5eb05f4..f440cc9be 100644 --- a/src/core/vector3d.h +++ b/src/core/vector3d.h @@ -54,9 +54,9 @@ class Vector3D static double dotProduct(const Vector3D& a, const Vector3D& b); protected: - double x{}; - double y{}; - double z{}; + double x_{}; + double y_{}; + double z_{}; }; Vector3D operator* (double lhs, const Vector3D& rhs); std::ostream& operator<< (std::ostream& os, const Vector3D& v); From 68c47b5a7d12c221be4fa74dd2b13472c92b0e76 Mon Sep 17 00:00:00 2001 From: tsteven4 <13596209+tsteven4@users.noreply.github.com> Date: Thu, 11 Nov 2021 08:20:44 -0700 Subject: [PATCH 6/8] kill obsolete comment --- src/core/nvector.h | 1 - 1 file changed, 1 deletion(-) diff --git a/src/core/nvector.h b/src/core/nvector.h index 2a175c0f6..a100b47b5 100644 --- a/src/core/nvector.h +++ b/src/core/nvector.h @@ -54,7 +54,6 @@ class LatLon class NVector : public Vector3D { -// friend class PVector; public: NVector() = default; From 1caaa2ca2cfa4695fc527f082fe5e1ba97448425 Mon Sep 17 00:00:00 2001 From: tsteven4 <13596209+tsteven4@users.noreply.github.com> Date: Fri, 12 Nov 2021 04:50:03 -0700 Subject: [PATCH 7/8] dont use getters within class. --- src/core/nvector.cc | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/src/core/nvector.cc b/src/core/nvector.cc index 1793a49b8..78906a7e9 100644 --- a/src/core/nvector.cc +++ b/src/core/nvector.cc @@ -72,9 +72,9 @@ std::pair PVector::toNVectorAndHeight() const constexpr double e2 = WGS84_ECCENTRICITY_SQUARED; constexpr double e4 = e2 * e2; - double px2 = getx() * getx(); - double py2 = gety() * gety(); - double pz2 = getz() * getz(); + double px2 = x_ * x_; + double py2 = y_ * y_; + double pz2 = z_ * z_; double q = ((1.0-e2)/(a2)) * px2; double p = (py2 + pz2) / (a2); @@ -90,9 +90,9 @@ std::pair PVector::toNVectorAndHeight() const double sf2 = k / (k+e2); double height = ((k+e2-1.0)/k) * sqrt(d*d+px2); - double nx = sf * getx(); - double ny = sf * sf2 * gety(); - double nz = sf * sf2 * getz(); + double nx = sf * x_; + double ny = sf * sf2 * y_; + double nz = sf * sf2 * z_; return {Vector3D(nx, ny, nz), height}; } From 47233e2d8bbed7b9bfb66d1135453a87215be169 Mon Sep 17 00:00:00 2001 From: tsteven4 <13596209+tsteven4@users.noreply.github.com> Date: Sat, 13 Nov 2021 05:55:36 -0700 Subject: [PATCH 8/8] performance tweaks for nvector. --- src/core/nvector.cc | 11 +++++------ src/core/nvector.h | 4 ++++ src/core/vector3d.cc | 15 --------------- src/core/vector3d.h | 6 +++--- 4 files changed, 12 insertions(+), 24 deletions(-) diff --git a/src/core/nvector.cc b/src/core/nvector.cc index 78906a7e9..a35efc691 100644 --- a/src/core/nvector.cc +++ b/src/core/nvector.cc @@ -26,7 +26,6 @@ #include #include -#include "defs.h" #include "nvector.h" #include "vector3d.h" @@ -50,8 +49,8 @@ NVector::NVector(double latitude_degrees, double longitude_degrees) // x points to North pole // y points towards latitude 0, longitude +90 (east) // z points to latitude 0, longitude +180, the antimeridian - double latitude_radians = latitude_degrees * M_PI/180.0; - double longitude_radians = longitude_degrees * M_PI/180.0; + double latitude_radians = latitude_degrees * kRadiansPerDegree; + double longitude_radians = longitude_degrees * kRadiansPerDegree; x_ = sin(latitude_radians); y_ = sin(longitude_radians)*cos(latitude_radians); z_ = -cos(longitude_radians)*cos(latitude_radians); @@ -101,7 +100,7 @@ double NVector::latitude() const { // This implements equation 6. double latitude_radians = atan2(x_, sqrt(y_*y_ + z_*z_)); - double latitude_degrees = latitude_radians * 180.0/M_PI; + double latitude_degrees = latitude_radians * kDegreesPerRadian; return latitude_degrees; } @@ -109,7 +108,7 @@ double NVector::longitude() const { // This implements equation 5. double longitude_radians = atan2(y_, -z_); - double longitude_degrees = longitude_radians * 180.0/M_PI; + double longitude_degrees = longitude_radians * kDegreesPerRadian; return longitude_degrees; } @@ -163,7 +162,7 @@ double NVector::azimuth_radians(const NVector& n_EA_E, const NVector& n_EB_E, do double NVector::azimuth(const NVector& n_EA_E, const NVector& n_EB_E, double height_a, double height_b) { double azimuth = azimuth_radians(n_EA_E, n_EB_E, height_a, height_b); - double azimuth_degrees = azimuth * 180.0/M_PI; + double azimuth_degrees = azimuth * kDegreesPerRadian; return azimuth_degrees; } diff --git a/src/core/nvector.h b/src/core/nvector.h index a100b47b5..04ce9d211 100644 --- a/src/core/nvector.h +++ b/src/core/nvector.h @@ -20,6 +20,7 @@ #ifndef NVECTOR_H #define NVECTOR_H +#include "defs.h" #include "vector3d.h" namespace gpsbabel @@ -41,6 +42,9 @@ constexpr double WGS84_ASPECT_RATIO = 1.0 - WGS84_FLATTENING; // b/a constexpr double WGS84_SEMI_MINOR_AXIS_METERS = WGS84_SEMI_MAJOR_AXIS_METERS * WGS84_ASPECT_RATIO; // b constexpr double WGS84_ECCENTRICITY_SQUARED = 1.0 - (WGS84_ASPECT_RATIO * WGS84_ASPECT_RATIO); +constexpr double kRadiansPerDegree = M_PI/180.0; +constexpr double kDegreesPerRadian = 1.0/kRadiansPerDegree; + class PVector; class LatLon diff --git a/src/core/vector3d.cc b/src/core/vector3d.cc index c2c77bca4..a59b44a72 100644 --- a/src/core/vector3d.cc +++ b/src/core/vector3d.cc @@ -40,21 +40,6 @@ double Vector3D::norm() const return norm; } -double Vector3D::getx() const -{ - return x_; -} - -double Vector3D::gety() const -{ - return y_; -} - -double Vector3D::getz() const -{ - return z_; -} - Vector3D& Vector3D::normalize() { *this = *this/this->norm(); diff --git a/src/core/vector3d.h b/src/core/vector3d.h index f440cc9be..bd9ac5474 100644 --- a/src/core/vector3d.h +++ b/src/core/vector3d.h @@ -34,9 +34,9 @@ class Vector3D Vector3D(double xi, double yi, double zi); [[nodiscard]] double norm() const; - [[nodiscard]] double getx() const; - [[nodiscard]] double gety() const; - [[nodiscard]] double getz() const; + [[nodiscard]] double getx() const {return x_;} + [[nodiscard]] double gety() const {return y_;} + [[nodiscard]] double getz() const {return z_;} Vector3D& normalize(); Vector3D& operator+=(const Vector3D& rhs);